BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
SimplePIDSvc.cxx
Go to the documentation of this file.
1#include "GaudiKernel/IDataProviderSvc.h"
2#include "GaudiKernel/MsgStream.h"
3#include "GaudiKernel/SmartDataPtr.h"
4#include "EventModel/EventModel.h"
5#include "EvtRecEvent/EvtRecEvent.h"
6#include "EventModel/EventHeader.h"
7#include "SimplePIDSvc/SimplePIDSvc.h"
8
9#include "DstEvent/TofHitStatus.h"
10#include "TMath.h"
11#include "TFile.h"
12#include "TMatrixD.h"
13#include "TArray.h"
14#include <fstream>
15#include <iostream>
16#include <cstdlib>
17#include <cmath>
18using namespace std;
19DECLARE_COMPONENT(SimplePIDSvc)
20
21SimplePIDSvc::SimplePIDSvc(const std::string& name, ISvcLocator* svcLoc) : base_class(name, svcLoc)
22{
23 declareProperty("DedxChiCut", m_dedx_chi_cut = 4);
24 declareProperty("TofChiCut", m_tof_chi_cut = 4);
25 declareProperty("IsTofCorr", m_tof_corr = true);
26 declareProperty("IsDedxCorr", m_dedx_corr = true);
27 declareProperty("EidRatio", m_eid_ratio = 0.80);
28}
29
30SimplePIDSvc::~SimplePIDSvc() {;}
31
32StatusCode SimplePIDSvc::initialize()
33{
34 MsgStream log(messageService(), name());
35 log << MSG::INFO << "in SimplePIDSvc initialize()" << endreq;
36
37 StatusCode sc = Service::initialize();
38
39 sc = serviceLocator()->service("EventDataSvc", eventSvc_, true);
40
41 loadHistogram();
42 loadSecondPar();//the second correct factor
43
44 return sc;
45}
46
47StatusCode SimplePIDSvc::finalize()
48{
49 MsgStream log(messageService(), name());
50 log << MSG::INFO << "in SimplePIDSvc finalize()" << endreq;
51
52 StatusCode sc = Service::finalize();
53
54 for (unsigned int i = 0; i < 2; i++)
55 {
56 for (unsigned int j = 0; j < 4; j++)
57 {
58 f_dedx[i][j]->Close();
59 f_tof_q[i][j]->Close();
60 f_tof_bgcost[i][j]->Close();
61 f_tof_wgt[i][j]->Close();
62 f_tof_final[i][j]->Close();
63 f_tofec_q[i][j]->Close();
64 f_tofec_bg[i][j]->Close();
65 f_tofec_cost[i][j]->Close();
66 }
67 }
68 for (unsigned int i = 0; i < 3; i++)
69 {
70 for (unsigned int j = 0; j < 4; j++)
71 {
72 f_emc[i][j]->Close();
73 }
74 }
75
76 return sc;
77}
78
79/*StatusCode SimplePIDSvc::queryInterface(const InterfaceID& riid, void** ppvIF)
80{
81 if (ISimplePIDSvc::interfaceID().versionMatch(riid))
82 {
83 *ppvIF = dynamic_cast<ISimplePIDSvc*>(this);
84 }
85 else
86 {
87 return Service::queryInterface(riid, ppvIF);
88 }
89 addRef();
90 return StatusCode::SUCCESS;
91}
92*/
93
94
95void SimplePIDSvc::preparePID(EvtRecTrack* track)
96{
97
98 SmartDataPtr<Event::EventHeader> eventHeaderpid(eventSvc_, "/Event/EventHeader");
99 m_run = eventHeaderpid->runNumber();
100
101 if (track->isMdcKalTrackValid())
102 {
103 RecMdcKalTrack *mdckalTrk = track->mdcKalTrack();
104 RecMdcKalTrack::PidType trk_type[5] = {
105 RecMdcKalTrack::electron,
106 RecMdcKalTrack::muon,
107 RecMdcKalTrack::pion,
108 RecMdcKalTrack::kaon,
109 RecMdcKalTrack::proton,
110 };
111 double mass[5] = {
112 0.000511,
113 0.105658,
114 0.13957,
115 0.493677,
116 0.938272,
117 };
118 for(unsigned int pid = 0; pid < 5; pid++)
119 {
120 mdckalTrk->setPidType(trk_type[pid]);
121 m_p[pid] = mdckalTrk->p();
122 m_betagamma[pid] = m_p[pid] / mass[pid];
123 m_charge[pid] = mdckalTrk->charge();
124 m_cost[pid] = cos(mdckalTrk->theta());
125 }
126 }
127 else
128 {
129 for(unsigned int i = 0; i < 5; i++)
130 {
131 m_p[i] = -99;
132 m_betagamma[i] = -99;
133 m_cost[i] = -99;
134 m_charge[i] = 0;
135 }
136
137 }
138
139 //dE/dx PID
140 loadDedxInfo(track);
141 if (m_dedx_corr)
142 {
143 dedxCorrection();
144 dedxSecondCorrection();
145 }
146 //TOF PID
147 loadTOFInfo(track);
148 if (m_tof_corr)
149 {
150 if (m_tof_barrel == 1)
151 {
152 tofBarrelCorrection();
153 tofBarrelSecondCorrection();
154 }
155 else if (m_tof_barrel == 0)
156 {
157 tofEndcapCorrection();
158 tofEndcapSecondCorrection();
159 }
160 }
161 //EMC
162 loadEMCInfo(track);
163
164 calprob();
165}
166
167void SimplePIDSvc::calprob()
168{
169 bool usededx = false;
170 bool usetof = false;
171
172 for (unsigned int i = 0; i < 5 ;i++)
173 {
174 if (!usededx && fabs(m_dedx_chi[i]) < m_dedx_chi_cut)
175 usededx = true;
176 if (!usetof && fabs(m_tof_chi[i]) < m_tof_chi_cut)
177 usetof = true;
178
179 m_dedx_only[i] = false;
180 }
181 if (!usededx)
182 {
183 for(unsigned int i = 0; i < 5; i++)
184 m_dedx_chi[i] = -99;
185 }
186 if (!usetof)
187 {
188 for(unsigned int i = 0; i < 5; i++)
189 m_tof_chi[i] = -99;
190 }
191
192 for (unsigned int i = 0; i < 5; i++)
193 {
194 m_prob[i] = -99;
195 double chi2 = 0;
196 int ndf = 0;
197
198 if (usededx && usetof)
199 {
200 chi2 = pow(m_dedx_chi[i], 2) + pow(m_tof_chi[i], 2);
201 ndf = 2;
202 }
203 else if (usededx && !usetof)
204 {
205 chi2 = pow(m_dedx_chi[i], 2);
206 ndf = 1;
207 m_dedx_only[i] = true;
208 }
209 else if (!usededx && usetof)
210 {
211 chi2 = pow(m_tof_chi[i],2);
212 ndf = 1;
213 }
214 if (ndf > 0 && chi2 > 0)
215 m_prob[i] = TMath::Prob(chi2, ndf);
216 }
217}
218
219int SimplePIDSvc::getRunIdx(int run_no)
220{
221 // -1: beyond correction region
222 // 0: 2010 psi(3770) data
223 // 1: 2011 psi(3770) data
224 // 2: 2010 psi(3770) mc
225 // 3: 2011 psi(3770) mc
226 const int RUN_BEGIN_DATA_10 = 11414;
227 const int RUN_END_DATA_10 = 14604;
228 const int RUN_BEGIN_MC_10 = -14604;
229 const int RUN_END_MC_10 = -11414;
230 const int RUN_BEGIN_DATA_11 = 20448;
231 const int RUN_END_DATA_11 = 23454;
232 const int RUN_BEGIN_MC_11 = -23454;
233 const int RUN_END_MC_11 = -20448;
234
235 const int RUN_BEGIN_DATA_22 = 70521;
236 const int RUN_END_DATA_22 = 73930;
237 const int RUN_BEGIN_MC_22 = -73930;
238 const int RUN_END_MC_22 = -70521;
239
240 if (run_no >= RUN_BEGIN_DATA_10 && run_no <= RUN_END_DATA_10)
241 return 0;
242 else if (run_no >= RUN_BEGIN_DATA_11 && run_no <= RUN_END_DATA_11)
243 return 1;
244 else if (run_no >= RUN_BEGIN_MC_10 && run_no <= RUN_END_MC_10)
245 return 2;
246 else if (run_no >= RUN_BEGIN_MC_11 && run_no <= RUN_END_MC_11)
247 return 3;
248
249 else if (run_no >= RUN_BEGIN_DATA_22 && run_no <= RUN_END_DATA_22)
250 return 1;
251 else if (run_no >= RUN_BEGIN_MC_22 && run_no <= RUN_END_MC_22)
252 return 3;
253 else
254 return -1;
255
256}
257
258void SimplePIDSvc::loadTOFInfo(EvtRecTrack *track)
259{
260 //Initialization
261 for (unsigned int i = 0; i < 8; i++)
262 {
263 for (unsigned int j = 0; j < 5; j++)
264 m_tof_dt[i][j] = -99.;
265 m_tof_ph[i] = -99.;
266 }
267 for (unsigned int i = 0; i < 2; i++)
268 {
269 m_tof_zr[i] = -9999.;
270 m_tof_counter[i] = -1;
271 }
272 for (unsigned int i = 0; i < 5; i++)
273 {
274 m_tof_chi[i] = -99.;
275 }
276 m_tof_barrel = -1;
277
278 if (!track->isExtTrackValid() || !track->isTofTrackValid()) return;
279
280 SmartRefVector<RecTofTrack> tofTrk = track->tofTrack();
281 SmartRefVector<RecTofTrack>::iterator it;
282 RecExtTrack* extTrk = track->extTrack();
283 double zrhit[2];
284 zrhit[0] = extTrk->tof1Position().z();
285 zrhit[1] = extTrk->tof2Position().z();
286
287 TofHitStatus *hitst = new TofHitStatus;
288
289 for (it = tofTrk.begin(); it != tofTrk.end(); it++)
290 {
291 unsigned int st = (*it)->status();
292 hitst->setStatus(st);
293 if (hitst->is_raw()) continue; //empty TOF hit
294 bool barrel = hitst->is_barrel();
295 bool readout = hitst->is_readout();
296 bool counter = hitst->is_counter();
297 bool cluster = hitst->is_cluster();
298 int layer = hitst->layer();
299 double tof = (*it)->tof();
300 double ph = (*it)->ph();
301 m_tof_counter[layer-1] = (*it)->tofID();
302
303 if (barrel)
304 {
305 m_tof_barrel = 1;
306 }
307 else
308 {
309 m_tof_barrel = 0;
310 zrhit[0] = extTrk->tof1Position().rho();
311 }
312 m_tof_zr[0] = zrhit[0];
313 m_tof_zr[1] = zrhit[1];
314
315 int idx = -1;
316 if (readout)
317 {
318 if (barrel)
319 idx = ((st & 0xC0) >> 5) + (((st ^ 0x20) & 0x20) >> 5) - 2;
320 else
321 idx = 7;
322 }
323 else if (counter)
324 {
325 idx = layer + 3;
326 }
327 else if (cluster)
328 {
329 idx = 6;
330 }
331 if (idx == -1) continue;
332 m_tof_ph[idx] = ph;
333 for (unsigned int i = 0; i < 5; i++)
334 {
335 double offset = (*it)->toffset(i);
336 double texp = (*it)->texp(i);
337 if (texp < 0.0) continue;
338 double dt = tof - offset - texp;
339 m_tof_dt[idx][i] = dt;
340 }
341 }
342 delete hitst;
343}
344
345void SimplePIDSvc::tofBarrelCorrection()
346{
347 const double EPS = 1e-4;
348 const double BG_LOW = 0.20;
349 const double BG_UP = 7.40;
350 const double COST_LOW = -0.81;
351 const double COST_UP = 0.81;
352 const double Q_LOW = 0.;
353 const double Q_UP = 9000.;
354 const double P_LOW = 0.2;
355 const double P_UP = 1.3;
356 const int BIN_BG = 15;
357 const int BIN_COST = 15;
358 const int BIN_P = 15;
359 double BG[BIN_BG + 1] = {0.20, 0.87, 1.11, 1.35, 1.55, 1.72, 1.91, 2.17, 2.63, 3.05, 3.47, 3.93, 4.50, 5.27, 6.00, 7.40};
360 double COST[BIN_COST + 1] = {-0.81, -0.64, -0.53, -0.43, -0.33, -0.23, -0.13, -0.04, 0.05, 0.14, 0.24, 0.34, 0.44, 0.54, 0.65, 0.81};
361 double P[BIN_P + 1] = {0.20, 0.47, 0.56, 0.65, 0.72, 0.79, 0.86, 0.92, 0.98, 1.03, 1.08, 1.13, 1.17, 1.22, 1.26, 1.30};
362 int idx = getRunIdx(m_run);
363
364 if (idx != -1)
365 {
366 for (unsigned int i = 0; i < 4; i++)
367 {// only correct e, pi, K
368 double bg;
369 int bin_bg, bin_cost, bin_wgt;
370 int pid;
371 if (i == 0)
372 {
373 bg = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
374 bin_bg = findBin(P, BIN_P, bg);
375 pid = 0;
376 }
377 else if (i == 2 || i == 3)
378 {
379 bg = max(BG_LOW+EPS, min(BG_UP-EPS, m_betagamma[i]));
380 bin_bg = findBin(BG, BIN_BG, bg);
381 pid = 1;
382 }
383 else
384 {
385 continue;
386 }
387 double cost = m_cost[i];
388 int charge = m_charge[i];
389 double t[5], q;
390 double offset, sigma;
391 double offset_q, offset_bgcost;
392 int flag[4] = {0, 0, 0, 0, };
393 cost = max(COST_LOW+EPS, min(COST_UP-EPS, cost));
394 bin_cost = findBin(COST, BIN_COST, cost);
395 if (bin_bg == -1 || bin_cost == -1) continue;
396
397 //corrections
398 for (unsigned int j = 0; j < 4; j++)
399 {
400 t[j] = m_tof_dt[j][i];
401 if (fabs(t[j] + 99.) < EPS)//no readout
402 flag[j] = 0;
403 else
404 flag[j] = 1;
405 q = m_tof_ph[j];
406 q = max(Q_LOW+EPS, min(Q_UP-EPS, q));
407 if (charge == 1)
408 {
409 offset_q = h_tof_p_q_offset[pid][idx][j]->Interpolate( q );
410 offset_bgcost = h_tof_p_bgcost_offset[pid][idx][j]->Interpolate( bg, cost );
411 t[j] = t[j] - offset_q - offset_bgcost;
412 }
413 else
414 {
415 offset_q = h_tof_m_q_offset[pid][idx][j]->Interpolate( q );
416 offset_bgcost = h_tof_m_bgcost_offset[pid][idx][j]->Interpolate( bg, cost );
417 t[j] = t[j] - offset_q - offset_bgcost;
418 }
419 }
420 bin_wgt = flag[0]*8 + flag[1]*4 + flag[2]*2 + flag[3] - 1;
421 if (bin_wgt == -1) continue;
422 t[4] = 0;
423 for (unsigned int j = 0; j < 4; j++)
424 {
425 if (charge == 1)
426 t[4] += t[j] * h_tof_p_wgt[pid][idx][bin_wgt][j]->GetBinContent( bin_bg+1, bin_cost+1 );
427 else
428 t[4] += t[j] * h_tof_m_wgt[pid][idx][bin_wgt][j]->GetBinContent( bin_bg+1, bin_cost+1 );
429 }
430 if (charge == 1)
431 {
432 t[4] /= h_tof_p_wgt[pid][idx][bin_wgt][4]->GetBinContent( bin_bg+1, bin_cost+1 );
433 offset = h_tof_p_final_offset[pid][idx][bin_wgt]->Interpolate( bg, cost );
434 sigma = h_tof_p_final_sigma[pid][idx][bin_wgt]-> Interpolate( bg, cost );
435 m_tof_chi[i] = (t[4] - offset) / sigma;
436 }
437 else
438 {
439 t[4] /= h_tof_m_wgt[pid][idx][bin_wgt][4]->GetBinContent( bin_bg+1, bin_cost+1 );
440 offset = h_tof_m_final_offset[pid][idx][bin_wgt]->Interpolate( bg, cost );
441 sigma = h_tof_m_final_sigma[pid][idx][bin_wgt]-> Interpolate( bg, cost );
442 m_tof_chi[i] = (t[4] - offset) / sigma;
443 }
444 }
445 }
446}
447
448void SimplePIDSvc::tofEndcapCorrection()
449{
450 const double EPS = 1e-4;
451 const double BG_LOW = 0.30;
452 const double BG_UP = 7.40;
453 const double Q_LOW = 0.;
454 const double Q_UP = 6000.;
455 const double COST_EAST_LOW = 0.720;
456 const double COST_EAST_UP = 0.930;
457 const double COST_WEST_LOW = -0.930;
458 const double COST_WEST_UP = -0.720;
459 const double P_LOW = 0.2;
460 const double P_UP = 1.3;
461
462 int idx = getRunIdx(m_run);
463
464 if (idx != -1)
465 {
466 for (unsigned int i = 0; i < 4; i++)
467 {// only correct e, pi, K
468 int pid;
469 double bg;
470 if (i == 0)
471 {
472 bg = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
473 pid = 0;
474 }
475 else if (i == 2 || i == 3)
476 {
477 bg = max(BG_LOW+EPS, min(BG_UP-EPS, m_betagamma[i]));
478 pid = 1;
479 }
480 else
481 {
482 continue;
483 }
484
485 int flag; //0:east, 1:west
486 double cost = m_cost[i];
487 int charge = m_charge[i];
488 double t = m_tof_dt[7][i];
489 double q = m_tof_ph[7];
490 double off_q, off_bg, off_cost;
491 double sg_q, sg_bg, sg_cost;
492 if (cost > 0)
493 {
494 flag = 0;
495 cost = max(COST_EAST_LOW+EPS, min(COST_EAST_UP-EPS, cost));
496 }
497 else
498 {
499 flag = 1;
500 cost = max(COST_WEST_LOW+EPS, min(COST_WEST_UP-EPS, cost));
501 }
502 q = max(Q_LOW+EPS, min(Q_UP-EPS, q));
503
504 //corrections
505 if (charge == 1)
506 {
507 off_q = h_tofec_p_q_offset[pid][idx][flag] ->Interpolate( q );
508 sg_q = h_tofec_p_q_sigma[pid][idx][flag] ->Interpolate( q );
509 off_bg = h_tofec_p_bg_offset[pid][idx][flag] ->Interpolate( bg );
510 sg_bg = h_tofec_p_bg_sigma[pid][idx][flag] ->Interpolate( bg );
511 off_cost = h_tofec_p_cost_offset[pid][idx][flag]->Interpolate( cost );
512 sg_cost = h_tofec_p_cost_sigma[pid][idx][flag] ->Interpolate( cost );
513 m_tof_chi[i] = (((t - off_q) / sg_q - off_bg) / sg_bg - off_cost) / sg_cost;
514 }
515 else
516 {
517 off_q = h_tofec_m_q_offset[pid][idx][flag] ->Interpolate( q );
518 sg_q = h_tofec_m_q_sigma[pid][idx][flag] ->Interpolate( q );
519 off_bg = h_tofec_m_bg_offset[pid][idx][flag] ->Interpolate( bg );
520 sg_bg = h_tofec_m_bg_sigma[pid][idx][flag] ->Interpolate( bg );
521 off_cost = h_tofec_m_cost_offset[pid][idx][flag]->Interpolate( cost );
522 sg_cost = h_tofec_m_cost_sigma[pid][idx][flag] ->Interpolate( cost );
523 m_tof_chi[i] = (((t - off_q) / sg_q - off_bg) / sg_bg - off_cost) / sg_cost;
524 }
525 }
526 }
527}
528
529void SimplePIDSvc::loadDedxInfo(EvtRecTrack *track)
530{
531 if (track->isMdcDedxValid())
532 {
533 RecMdcDedx* dedx_trk = track->mdcDedx();
534 for (unsigned int i = 0; i < 5; i++)
535 m_dedx_chi[i] = dedx_trk->chi(i);
536 }
537 else
538 {
539 for (unsigned int i = 0; i < 5; i++)
540 m_dedx_chi[i] = -99;
541 }
542}
543
544void SimplePIDSvc::dedxCorrection()
545{
546 int idx = getRunIdx(m_run);
547 const double EPS = 1e-4;
548 const double BG_LOW = 0.20;
549 const double BG_UP = 7.40;
550 const double COST_LOW = -0.93;
551 const double COST_UP = 0.93;
552 const double P_LOW = 0.2;
553 const double P_UP = 1.3;
554 if (idx != -1)
555 {
556 double offset, sigma;
557 for (unsigned int i = 0; i < 4; i++)
558 {// only correct e, pi, K
559 double bg;
560 int pid;
561 if (i == 0)
562 {
563 bg = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
564 pid = 0;
565 }
566 else if (i == 2 || i == 3)
567 {
568 bg = max(BG_LOW+EPS, min(BG_UP-EPS, m_betagamma[i]));
569 pid = 1;
570 }
571 else
572 {
573 continue;
574 }
575 double cost = m_cost[i];
576 double charge = m_charge[i];
577 cost = max(COST_LOW+EPS, min(COST_UP-EPS, cost));
578 if (charge == 1)
579 {
580 offset = h_dedx_p_offset[pid][idx]->Interpolate( bg, cost );
581 sigma = h_dedx_p_sigma[pid][idx] ->Interpolate( bg, cost );
582 m_dedx_chi[i] = (m_dedx_chi[i] - offset) / sigma;
583 }
584 else
585 {
586 offset = h_dedx_m_offset[pid][idx]->Interpolate( bg, cost );
587 sigma = h_dedx_m_sigma[pid][idx] ->Interpolate( bg, cost );
588 m_dedx_chi[i] = (m_dedx_chi[i] - offset) / sigma;
589 }
590 }
591 }
592}
593
594void SimplePIDSvc::tofBarrelSecondCorrection()
595{
596
597int idx = getRunIdx(m_run);
598const double EPS = 1e-4;
599const double P_LOW = 0.0;
600const double P_UP = 1.3;
601const int BIN_P = 10;
602
603const int RUN_BEGIN_DATA_10 = 11414;
604const int RUN_END_DATA_10 = 14604;
605const int RUN_BEGIN_MC_10 = -14604;
606const int RUN_END_MC_10 = -11414;
607const int RUN_BEGIN_DATA_11 = 20448;
608const int RUN_END_DATA_11 = 23454;
609const int RUN_BEGIN_MC_11 = -23454;
610const int RUN_END_MC_11 = -20448;
611
612const int RUN_BEGIN_DATA_22 = 70521;
613const int RUN_END_DATA_22 = 73930;
614const int RUN_BEGIN_MC_22 = -73930;
615const int RUN_END_MC_22 = -70521;
616
617double P[BIN_P + 1] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.3};
618if (idx != -1)
619{
620
621for (unsigned int i = 2; i < 4; i++){// second correct for tof of pi k
622
623int aa=99,bb=99;
624
625int bin_p;
626double ptk = m_p[i];
627ptk = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
628bin_p = findBin(P, BIN_P, ptk);
629
630if(i==2){// pi
631
632if((m_run>=RUN_BEGIN_DATA_10 && m_run<=RUN_END_DATA_10 && idx == 0)||(m_run>=RUN_BEGIN_DATA_11 && m_run<=RUN_END_DATA_11 && idx == 1))
633{aa=2; bb=2;}// round0304 data
634
635if((m_run>=RUN_BEGIN_MC_10 && m_run<=RUN_END_MC_10 && idx == 2)||(m_run>=RUN_BEGIN_MC_11 && m_run<=RUN_END_MC_11 && idx == 3))
636{aa=2; bb=3;}// round0304 inc
637
638if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
639{aa=3; bb=2;}// round15 data
640
641if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
642{aa=3; bb=3;}// round15 inc
643
644}
645
646if(i==3){// k
647
648if((m_run>=RUN_BEGIN_DATA_10 && m_run<=RUN_END_DATA_10 && idx == 0)||(m_run>=RUN_BEGIN_DATA_11 && m_run<=RUN_END_DATA_11 && idx == 1))
649{aa=2; bb=0;}
650
651if((m_run>=RUN_BEGIN_MC_10 && m_run<=RUN_END_MC_10 && idx == 2)||(m_run>=RUN_BEGIN_MC_11 && m_run<=RUN_END_MC_11 && idx == 3))
652{aa=2; bb=1;}
653
654if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
655{aa=3; bb=0;}
656
657if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
658{aa=3; bb=1;}
659
660}
661
662if(m_tof_chi[i]!=-99 && m_run>0){
663
664if(bin_p<4){
665m_tof_chi[i] = (m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p]) / (m_gaussion_sigmab[aa][bb][bin_p] / m_gaussion_sigmab[aa][bb+1][bin_p]);
666}
667if(bin_p>=4){
668m_tof_chi[i] = (m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p]);
669}
670
671}
672if(m_tof_chi[i]!=-99 && m_run<0){
673m_tof_chi[i] = (m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p]) / (m_gaussion_sigmab[aa][bb][bin_p] / m_gaussion_sigmab[aa][bb][bin_p]);
674}
675
676}
677
678}//end idx!= -1
679
680}
681
682void SimplePIDSvc::tofEndcapSecondCorrection()
683{
684
685int idx = getRunIdx(m_run);
686const double EPS = 1e-4;
687const double P_LOW = 0.0;
688const double P_UP = 1.3;
689const int BIN_P = 10;
690
691const int RUN_BEGIN_DATA_10 = 11414;
692const int RUN_END_DATA_10 = 14604;
693const int RUN_BEGIN_MC_10 = -14604;
694const int RUN_END_MC_10 = -11414;
695const int RUN_BEGIN_DATA_11 = 20448;
696const int RUN_END_DATA_11 = 23454;
697const int RUN_BEGIN_MC_11 = -23454;
698const int RUN_END_MC_11 = -20448;
699
700const int RUN_BEGIN_DATA_22 = 70521;
701const int RUN_END_DATA_22 = 73930;
702const int RUN_BEGIN_MC_22 = -73930;
703const int RUN_END_MC_22 = -70521;
704
705double P[BIN_P + 1] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.3};
706if (idx != -1)
707{
708
709for (unsigned int i = 2; i < 4; i++){// second correct for tof of pi k
710
711int aa=99,bb=99;
712
713int bin_p;
714double ptk = m_p[i];
715ptk = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
716bin_p = findBin(P, BIN_P, ptk);
717
718if(i==2){// pi
719
720if((m_run>=RUN_BEGIN_DATA_10 && m_run<=RUN_END_DATA_10 && idx == 0)||(m_run>=RUN_BEGIN_DATA_11 && m_run<=RUN_END_DATA_11 && idx == 1))
721{aa=4; bb=2;}// round0304 data endcap
722
723if((m_run>=RUN_BEGIN_MC_10 && m_run<=RUN_END_MC_10 && idx == 2)||(m_run>=RUN_BEGIN_MC_11 && m_run<=RUN_END_MC_11 && idx == 3))
724{aa=4; bb=3;}// round0304 inc endcap
725
726if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
727{aa=5; bb=2;}// round15 data endcap
728
729if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
730{aa=5; bb=3;}// round15 inc endcap
731
732}
733
734if(i==3){// k
735
736if((m_run>=RUN_BEGIN_DATA_10 && m_run<=RUN_END_DATA_10 && idx == 0)||(m_run>=RUN_BEGIN_DATA_11 && m_run<=RUN_END_DATA_11 && idx == 1))
737{aa=4; bb=0;}
738
739if((m_run>=RUN_BEGIN_MC_10 && m_run<=RUN_END_MC_10 && idx == 2)||(m_run>=RUN_BEGIN_MC_11 && m_run<=RUN_END_MC_11 && idx == 3))
740{aa=4; bb=1;}
741
742if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
743{aa=5; bb=0;}
744
745if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
746{aa=5; bb=1;}
747
748}
749
750if(m_tof_chi[i]!=-99 && aa==5){// for round15 data and inc
751
752m_tof_chi[i] = (m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p]) / (60./110.);
753
754}
755
756if(m_tof_chi[i]!=-99 && aa==4){// for round0304 data and inc
757
758m_tof_chi[i] = (m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p]);
759
760}
761
762
763}
764
765}//end idx!= -1
766
767
768}
769
770void SimplePIDSvc::dedxSecondCorrection()
771
772{
773
774int idx = getRunIdx(m_run);
775const double EPS = 1e-4;
776const double P_LOW = 0.0;
777const double P_UP = 1.3;
778const int BIN_P = 10;
779
780const int RUN_BEGIN_DATA_10 = 11414;
781const int RUN_END_DATA_10 = 14604;
782const int RUN_BEGIN_MC_10 = -14604;
783const int RUN_END_MC_10 = -11414;
784const int RUN_BEGIN_DATA_11 = 20448;
785const int RUN_END_DATA_11 = 23454;
786const int RUN_BEGIN_MC_11 = -23454;
787const int RUN_END_MC_11 = -20448;
788
789const int RUN_BEGIN_DATA_22 = 70521;
790const int RUN_END_DATA_22 = 73930;
791const int RUN_BEGIN_MC_22 = -73930;
792const int RUN_END_MC_22 = -70521;
793
794double P[BIN_P + 1] = {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.3};
795if (idx != -1)
796{
797
798for (unsigned int i = 2; i < 4; i++){// second correct for dedx of pi k
799
800int aa=99,bb=99;
801
802int bin_p;
803double ptk = m_p[i];
804ptk = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
805bin_p = findBin(P, BIN_P, ptk);
806
807if(i==2){// pi
808
809if((m_run>=RUN_BEGIN_DATA_10 && m_run<=RUN_END_DATA_10 && idx == 0)||(m_run>=RUN_BEGIN_DATA_11 && m_run<=RUN_END_DATA_11 && idx == 1))
810{aa=0; bb=2;}// round0304 data
811
812if((m_run>=RUN_BEGIN_MC_10 && m_run<=RUN_END_MC_10 && idx == 2)||(m_run>=RUN_BEGIN_MC_11 && m_run<=RUN_END_MC_11 && idx == 3))
813{aa=0; bb=3;}// round0304 inc
814
815if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
816{aa=1; bb=2;}// round15 data
817
818if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
819{aa=1; bb=3;}// round15 inc
820
821}
822
823if(i==3){// k
824
825if((m_run>=RUN_BEGIN_DATA_10 && m_run<=RUN_END_DATA_10 && idx == 0)||(m_run>=RUN_BEGIN_DATA_11 && m_run<=RUN_END_DATA_11 && idx == 1))
826{aa=0; bb=0;}
827
828if((m_run>=RUN_BEGIN_MC_10 && m_run<=RUN_END_MC_10 && idx == 2)||(m_run>=RUN_BEGIN_MC_11 && m_run<=RUN_END_MC_11 && idx == 3))
829{aa=0; bb=1;}
830
831if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
832{aa=1; bb=0;}
833
834if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
835{aa=1; bb=1;}
836
837}
838
839if(m_dedx_chi[i]!=-99 && m_run>0){
840m_dedx_chi[i] = (m_dedx_chi[i] - m_gaussion_mean[aa][bb][bin_p]) / (m_gaussion_sigmab[aa][bb][bin_p] / m_gaussion_sigmab[aa][bb+1][bin_p]);
841}
842if(m_dedx_chi[i]!=-99 && m_run<0){
843m_dedx_chi[i] = (m_dedx_chi[i] - m_gaussion_mean[aa][bb][bin_p]) / (m_gaussion_sigmab[aa][bb][bin_p] / m_gaussion_sigmab[aa][bb][bin_p]);
844}
845
846}
847
848}//end idx!= -1
849
850}
851
852void SimplePIDSvc::loadSecondPar()
853{
854 string path = getenv("SIMPLEPIDSVCROOT");
855
856for(int i=0; i<6; i++){
857
858const char *dir;
859
860if(i == 0) dir = "round0304_dedx";
861else if(i == 1) dir = "round15_dedx";
862else if(i == 2) dir = "round0304_tof";
863else if(i == 3) dir = "round15_tof";
864else if(i == 4) dir = "round0304_tof_endcap";
865else if(i == 5) dir = "round15_tof_endcap";
866else{
867cout << "Boundary Error! " << endl;
868exit(1);
869}
870
871for(int j=0; j<4; j++){
872
873const char *name;
874
875if(j == 0) name = "data_K";
876else if(j == 1) name = "inc_K";
877else if(j == 2) name = "data_pi";
878else if(j == 3) name = "inc_pi";
879else{
880cout << "Boundary Error! " << endl;
881exit(1);
882}
883
884ifstream second_cor( path + Form("/share/second_correct/%s/%s_%s.dat",dir,dir,name));
885
886for(int m=0; m<10; m++){
887
888second_cor>>m_gaussion_mean[i][j][m]>>m_gaussion_sigma[i][j][m]>>m_gaussion_sigmab[i][j][m];
889
890}
891
892second_cor.close();
893
894}
895}
896
897
898}
899
900void SimplePIDSvc::loadHistogram()
901{
902 string path = getenv("SIMPLEPIDSVCROOT");
903 vector<string> filename;
904 for (unsigned int idx = 0; idx < 2; idx++)
905 {
906 const char *dir;
907 if (idx == 0)
908 dir = "electron";
909 else if (idx == 1)
910 dir = "kpi";
911 else
912 {
913 cout << "Boundary Error! " << endl;
914 exit(1);
915 }
916
917 //dedx
918 filename.clear();
919 filename.push_back( path + Form("/share/%s/dedx/dedx_d10.root", dir) );
920 filename.push_back( path + Form("/share/%s/dedx/dedx_d11.root", dir) );
921 filename.push_back( path + Form("/share/%s/dedx/dedx_m10.root", dir) );
922 filename.push_back( path + Form("/share/%s/dedx/dedx_m11.root", dir) );
923 for (unsigned int i = 0; i < filename.size(); i++)
924 {
925 f_dedx[idx][i] = new TFile(filename[i].c_str(), "READ");
926 const char *name;
927 if (i == 0)
928 name = "d10";
929 else if (i == 1)
930 name = "d11";
931 else if (i == 2)
932 name = "m10";
933 else if (i == 3)
934 name = "m11";
935 else
936 {
937 cout << "Boundary Error! " << endl;
938 exit(1);
939 }
940 h_dedx_p_offset[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_p_offset_%s", name) );
941 h_dedx_p_sigma[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_p_sigma_%s" , name) );
942 h_dedx_m_offset[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_m_offset_%s", name) );
943 h_dedx_m_sigma[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_m_sigma_%s" , name) );
944 }
945 //tof_barrel q
946 filename.clear();
947 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_d10.root", dir) );
948 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_d11.root", dir) );
949 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_m10.root", dir) );
950 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_m11.root", dir) );
951 for (unsigned int i = 0; i < filename.size(); i++)
952 {
953 f_tof_q[idx][i] = new TFile(filename[i].c_str(), "READ");
954 const char *name;
955 if (i == 0)
956 name = "d10";
957 else if (i == 1)
958 name = "d11";
959 else if (i == 2)
960 name = "m10";
961 else if (i == 3)
962 name = "m11";
963 else
964 {
965 cout << "Boundary Error! " << endl;
966 exit(1);
967 }
968 for (unsigned int j = 0; j < 4; j++)
969 {
970 h_tof_p_q_offset[idx][i][j] = (TH1D*)f_tof_q[idx][i]->Get( Form("h_tof_p_q_offset_%s_%d", name, j) );
971 h_tof_m_q_offset[idx][i][j] = (TH1D*)f_tof_q[idx][i]->Get( Form("h_tof_m_q_offset_%s_%d", name, j) );
972 h_tof_p_q_sigma[idx][i][j] = (TH1D*)f_tof_q[idx][i]->Get( Form("h_tof_p_q_sigma_%s_%d" , name, j) );
973 h_tof_m_q_sigma[idx][i][j] = (TH1D*)f_tof_q[idx][i]->Get( Form("h_tof_m_q_sigma_%s_%d" , name, j) );
974 }
975 }
976 //tof_barrel bg&cost
977 filename.clear();
978 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_d10.root", dir) );
979 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_d11.root", dir) );
980 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_m10.root", dir) );
981 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_m11.root", dir) );
982 for (unsigned int i = 0; i < filename.size(); i++)
983 {
984 f_tof_bgcost[idx][i] = new TFile(filename[i].c_str(), "READ");
985 const char *name;
986 if (i == 0)
987 name = "d10";
988 else if (i == 1)
989 name = "d11";
990 else if (i == 2)
991 name = "m10";
992 else if (i == 3)
993 name = "m11";
994 else
995 {
996 cout << "Boundary Error! " << endl;
997 exit(1);
998 }
999 for (unsigned int j = 0; j < 4; j++)
1000 {
1001 h_tof_p_bgcost_offset[idx][i][j] = (TH2D*)f_tof_bgcost[idx][i]->Get( Form("h_tof_p_bgcost_offset_%s_%d", name, j) );
1002 h_tof_m_bgcost_offset[idx][i][j] = (TH2D*)f_tof_bgcost[idx][i]->Get( Form("h_tof_m_bgcost_offset_%s_%d", name, j) );
1003 h_tof_p_bgcost_sigma[idx][i][j] = (TH2D*)f_tof_bgcost[idx][i]->Get( Form("h_tof_p_bgcost_sigma_%s_%d" , name, j) );
1004 h_tof_m_bgcost_sigma[idx][i][j] = (TH2D*)f_tof_bgcost[idx][i]->Get( Form("h_tof_m_bgcost_sigma_%s_%d" , name, j) );
1005 }
1006 }
1007 //tof_barrel wgt
1008 filename.clear();
1009 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_d10.root", dir) );
1010 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_d11.root", dir) );
1011 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_m10.root", dir) );
1012 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_m11.root", dir) );
1013 for (unsigned int i = 0; i < filename.size(); i++)
1014 {
1015 f_tof_wgt[idx][i] = new TFile(filename[i].c_str(), "READ");
1016 const char *name;
1017 if (i == 0)
1018 name = "d10";
1019 else if (i == 1)
1020 name = "d11";
1021 else if (i == 2)
1022 name = "m10";
1023 else if (i == 3)
1024 name = "m11";
1025 else
1026 {
1027 cout << "Boundary Error! " << endl;
1028 exit(1);
1029 }
1030 for (unsigned int j = 0; j < 15; j++)
1031 {
1032 for (unsigned int k = 0; k < 5; k++)
1033 {
1034 h_tof_p_wgt[idx][i][j][k] = (TH2D*)f_tof_wgt[idx][i]->Get( Form("h_tof_p_wgt_%s_%d_%d", name, j, k) );
1035 h_tof_m_wgt[idx][i][j][k] = (TH2D*)f_tof_wgt[idx][i]->Get( Form("h_tof_m_wgt_%s_%d_%d", name, j, k) );
1036 }
1037 }
1038 }
1039 //tof_barrel corr
1040 filename.clear();
1041 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_d10.root", dir) );
1042 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_d11.root", dir) );
1043 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_m10.root", dir) );
1044 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_m11.root", dir) );
1045 for (unsigned int i = 0; i < filename.size(); i++)
1046 {
1047 f_tof_final[idx][i] = new TFile(filename[i].c_str(), "READ");
1048 const char *name;
1049 if (i == 0)
1050 name = "d10";
1051 else if (i == 1)
1052 name = "d11";
1053 else if (i == 2)
1054 name = "m10";
1055 else if (i == 3)
1056 name = "m11";
1057 else
1058 {
1059 cout << "Boundary Error! " << endl;
1060 exit(1);
1061 }
1062 for (unsigned int j = 0; j < 15; j++)
1063 {
1064 h_tof_p_final_offset[idx][i][j] = (TH2D*)f_tof_final[idx][i]->Get( Form("h_tof_p_final_offset_%s_%d", name, j) );
1065 h_tof_m_final_offset[idx][i][j] = (TH2D*)f_tof_final[idx][i]->Get( Form("h_tof_m_final_offset_%s_%d", name, j) );
1066 h_tof_p_final_sigma[idx][i][j] = (TH2D*)f_tof_final[idx][i]->Get( Form("h_tof_p_final_sigma_%s_%d" , name, j) );
1067 h_tof_m_final_sigma[idx][i][j] = (TH2D*)f_tof_final[idx][i]->Get( Form("h_tof_m_final_sigma_%s_%d" , name, j) );
1068 }
1069 }
1070 //tof_endcap q
1071 filename.clear();
1072 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_d10.root", dir) );
1073 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_d11.root", dir) );
1074 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_m10.root", dir) );
1075 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_m11.root", dir) );
1076 for (unsigned int i = 0; i < filename.size(); i++)
1077 {
1078 f_tofec_q[idx][i] = new TFile(filename[i].c_str(), "READ");
1079 const char *name;
1080 if (i == 0)
1081 name = "d10";
1082 else if (i == 1)
1083 name = "d11";
1084 else if (i == 2)
1085 name = "m10";
1086 else if (i == 3)
1087 name = "m11";
1088 else
1089 {
1090 cout << "Boundary Error! " << endl;
1091 exit(1);
1092 }
1093 for (unsigned int j = 0; j < 2; j++)
1094 {
1095 h_tofec_p_q_offset[idx][i][j] = (TH1D*)f_tofec_q[idx][i]->Get( Form("h_tofec_p_q_offset_%s_%d", name, j) );
1096 h_tofec_m_q_offset[idx][i][j] = (TH1D*)f_tofec_q[idx][i]->Get( Form("h_tofec_m_q_offset_%s_%d", name, j) );
1097 h_tofec_p_q_sigma[idx][i][j] = (TH1D*)f_tofec_q[idx][i]->Get( Form("h_tofec_p_q_sigma_%s_%d" , name, j) );
1098 h_tofec_m_q_sigma[idx][i][j] = (TH1D*)f_tofec_q[idx][i]->Get( Form("h_tofec_m_q_sigma_%s_%d" , name, j) );
1099 }
1100 }
1101 //tof_endcap bg
1102 filename.clear();
1103 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_d10.root", dir) );
1104 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_d11.root", dir) );
1105 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_m10.root", dir) );
1106 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_m11.root", dir) );
1107 for (unsigned int i = 0; i < filename.size(); i++)
1108 {
1109 f_tofec_bg[idx][i] = new TFile(filename[i].c_str(), "READ");
1110 const char *name;
1111 if (i == 0)
1112 name = "d10";
1113 else if (i == 1)
1114 name = "d11";
1115 else if (i == 2)
1116 name = "m10";
1117 else if (i == 3)
1118 name = "m11";
1119 else
1120 {
1121 cout << "Boundary Error! " << endl;
1122 exit(1);
1123 }
1124 for (unsigned int j = 0; j < 2; j++)
1125 {
1126 h_tofec_p_bg_offset[idx][i][j] = (TH1D*)f_tofec_bg[idx][i]->Get( Form("h_tofec_p_bg_offset_%s_%d", name, j) );
1127 h_tofec_m_bg_offset[idx][i][j] = (TH1D*)f_tofec_bg[idx][i]->Get( Form("h_tofec_m_bg_offset_%s_%d", name, j) );
1128 h_tofec_p_bg_sigma[idx][i][j] = (TH1D*)f_tofec_bg[idx][i]->Get( Form("h_tofec_p_bg_sigma_%s_%d" , name, j) );
1129 h_tofec_m_bg_sigma[idx][i][j] = (TH1D*)f_tofec_bg[idx][i]->Get( Form("h_tofec_m_bg_sigma_%s_%d" , name, j) );
1130 }
1131 }
1132 //tof_endcap cost
1133 filename.clear();
1134 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_d10.root", dir) );
1135 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_d11.root", dir) );
1136 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_m10.root", dir) );
1137 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_m11.root", dir) );
1138 for (unsigned int i = 0; i < filename.size(); i++)
1139 {
1140 f_tofec_cost[idx][i] = new TFile(filename[i].c_str(), "READ");
1141 const char *name;
1142 if (i == 0)
1143 name = "d10";
1144 else if (i == 1)
1145 name = "d11";
1146 else if (i == 2)
1147 name = "m10";
1148 else if (i == 3)
1149 name = "m11";
1150 else
1151 {
1152 cout << "Boundary Error! " << endl;
1153 exit(1);
1154 }
1155 for (unsigned int j = 0; j < 2; j++)
1156 {
1157 h_tofec_p_cost_offset[idx][i][j] = (TH1D*)f_tofec_cost[idx][i]->Get( Form("h_tofec_p_cost_offset_%s_%d", name, j) );
1158 h_tofec_m_cost_offset[idx][i][j] = (TH1D*)f_tofec_cost[idx][i]->Get( Form("h_tofec_m_cost_offset_%s_%d", name, j) );
1159 h_tofec_p_cost_sigma[idx][i][j] = (TH1D*)f_tofec_cost[idx][i]->Get( Form("h_tofec_p_cost_sigma_%s_%d" , name, j) );
1160 h_tofec_m_cost_sigma[idx][i][j] = (TH1D*)f_tofec_cost[idx][i]->Get( Form("h_tofec_m_cost_sigma_%s_%d" , name, j) );
1161 }
1162 }
1163 }
1164 for (unsigned int idx = 0; idx < 3; idx++)
1165 {
1166 const char *dir;
1167 if (idx == 0)
1168 dir = "electron/emc";
1169 else if (idx == 1)
1170 dir = "kpi/emc_pion";
1171 else if (idx == 2)
1172 dir = "kpi/emc_kaon";
1173 else
1174 {
1175 cout << "Boundary Error! " << endl;
1176 exit(1);
1177 }
1178 //emc
1179 filename.clear();
1180 filename.push_back( path + Form("/share/%s/emc_d10.root", dir) );
1181 filename.push_back( path + Form("/share/%s/emc_d11.root", dir) );
1182 filename.push_back( path + Form("/share/%s/emc_m10.root", dir) );
1183 filename.push_back( path + Form("/share/%s/emc_m11.root", dir) );
1184 for (unsigned int i = 0; i < filename.size(); i++)
1185 {
1186 f_emc[idx][i] = new TFile(filename[i].c_str(), "READ");
1187 const char *name;
1188 if (i == 0)
1189 name = "d10";
1190 else if (i == 1)
1191 name = "d11";
1192 else if (i == 2)
1193 name = "m10";
1194 else if (i == 3)
1195 name = "m11";
1196 else
1197 {
1198 cout << "Boundary Error! " << endl;
1199 exit(1);
1200 }
1201 for (unsigned int j = 0; j < 15; j++)
1202 {
1203 for (unsigned int k = 0; k < 25; k++)
1204 {
1205 h_emc_ep[idx][i][j][k] = (TH1D*)f_emc[idx][i]->Get( Form("h_ep_%s_%d_%d", name, j, k) );
1206 h_emc_e35[idx][i][j][k] = (TH1D*)f_emc[idx][i]->Get( Form("h_e35_%s_%d_%d", name, j, k) );
1207 }
1208 }
1209 }
1210 }
1211 cout << "Successfully Return from Loading Initializations by package SimplePIDSvc ... " << endl;
1212}
1213
1214int SimplePIDSvc::findBin(double *a, int length, double value)
1215{
1216 for (int i = 0; i < length; i++)
1217 {
1218 if (value > a[i] && value <= a[i+1])
1219 {
1220 return i;
1221 }
1222 }
1223 if (value < a[0])
1224 {
1225 return 0;
1226 }
1227 else
1228 {
1229 return length;
1230 }
1231}
1232
1233double SimplePIDSvc::getChi2(int i)
1234{
1235 return pow(m_dedx_chi[i], 2) + pow(m_tof_chi[i], 2);
1236}
1237
1238void SimplePIDSvc::loadEMCInfo(EvtRecTrack *track)
1239{
1240 //Initialization
1241 for (unsigned int i = 0; i < 5; i++)
1242 {
1243 m_emc_eop[i] = -99.;
1244 m_emc_likelihood[i] = -99.;
1245 }
1246 m_emc_e = -99.;
1247 m_emc_e13 = -99.;
1248 m_emc_e35 = -99.;
1249 m_emc_sec = -99.;
1250 m_emc_lat = -99.;
1251 m_lh_electron = -99.;
1252
1253 if (!track->isEmcShowerValid()) return;
1254
1255 RecEmcShower* emc_trk = track->emcShower();
1256 m_emc_e = emc_trk->energy();
1257 for (unsigned int i = 0; i < 5; i++)
1258 {
1259 m_emc_eop[i] = m_emc_e / fabs(m_p[i]);
1260 }
1261 double eseed = emc_trk->eSeed();
1262 double e3 = emc_trk->e3x3();
1263 double e5 = emc_trk->e5x5();
1264 m_emc_sec = emc_trk->secondMoment() / 1000.;
1265 m_emc_lat = emc_trk->latMoment();
1266 if (e3 != 0)
1267 {
1268 m_emc_e13 = eseed / e3;
1269 }
1270 if (e5 != 0)
1271 {
1272 m_emc_e35 = e3 / e5;
1273 }
1274}
1275
1276bool SimplePIDSvc::calEMCLikelihood()
1277{
1278 if (m_emc_eop[0] < 0)
1279 return false;
1280
1281 int idx = getRunIdx(m_run);
1282 const Int_t BIN_P = 15;
1283 const Int_t BIN_COST = 25;
1284 const Int_t BIN_PID = 3;
1285 const double EPS = 1e-4;
1286 //electron
1287 double P[BIN_PID][BIN_P + 1] = {
1288 {0.20, 0.47, 0.56, 0.65, 0.72, 0.79, 0.86, 0.92, 0.98, 1.03, 1.08, 1.13, 1.17, 1.22, 1.26, 1.30},
1289 {0.20, 0.26, 0.31, 0.35, 0.39, 0.42, 0.46, 0.49, 0.53, 0.57, 0.62, 0.67, 0.73, 0.80, 0.88, 1.05},
1290 {0.20, 0.33, 0.39, 0.43, 0.48, 0.52, 0.56, 0.61, 0.67, 0.73, 0.76, 0.81, 0.85, 0.90, 0.96, 1.05}, };
1291 double COST[BIN_PID][BIN_COST + 1] = {
1292 {-0.930, -0.910, -0.905, -0.897, -0.890, -0.881, -0.871, -0.858, -0.775, -0.732, -0.669, -0.561, -0.330, 0.199, 0.515, 0.645, 0.718, 0.766, 0.804, 0.870, 0.882, 0.891, 0.898, 0.906, 0.913, 0.930},
1293 {-0.930, -0.810, -0.728, -0.648, -0.574, -0.501, -0.431, -0.364, -0.295, -0.228, -0.161, -0.096, -0.031, 0.035, 0.100, 0.167, 0.234, 0.301, 0.370, 0.439, 0.510, 0.580, 0.655, 0.733, 0.813, 0.930},
1294 {-0.930, -0.804, -0.721, -0.643, -0.568, -0.497, -0.429, -0.362, -0.293, -0.228, -0.161, -0.096, -0.029, 0.035, 0.100, 0.166, 0.233, 0.298, 0.365, 0.432, 0.500, 0.571, 0.644, 0.722, 0.805, 0.930}, };
1295
1296 double vp, vcost;
1297 int pid;
1298 int bin_p, bin_cost;
1299 for (unsigned int i = 0; i < 4; i++)
1300 {
1301 //only e, pi ,K
1302 if (i == 0)
1303 pid = 0;
1304 else if (i == 2)
1305 pid = 1;
1306 else if (i == 3)
1307 pid = 2;
1308 else
1309 continue;
1310
1311 vp = max(P[pid][0]+EPS, min(P[pid][BIN_P]-EPS, m_p[i]));
1312 vcost = max(COST[pid][0]+EPS, min(COST[pid][BIN_COST]-EPS, m_cost[i]));
1313 bin_p = findBin(P[pid], BIN_P, vp);
1314 bin_cost = findBin(COST[pid], BIN_COST, vcost);
1315
1316 m_emc_likelihood[i] = h_emc_ep[pid][idx][bin_p][bin_cost]->Interpolate(m_emc_eop[i]) * h_emc_e35[pid][idx][bin_p][bin_cost]->Interpolate(m_emc_e35);
1317 }
1318 double a = m_prob[0] > 0 ? m_prob[0] : 0;
1319 double b = m_prob[2] > 0 ? m_prob[2] : 0;
1320 double c = m_prob[3] > 0 ? m_prob[3] : 0;
1321 double sum = a * m_emc_likelihood[0] + b * m_emc_likelihood[2] + c * m_emc_likelihood[3];
1322
1323 if (sum > 0 && m_prob[0] > 0)
1324 {
1325 m_lh_electron = m_prob[0] * m_emc_likelihood[0] / sum;
1326 return true;
1327 }
1328 else
1329 {
1330 return false;
1331 }
1332}
1333
1334bool SimplePIDSvc::ispion()
1335{
1336 if (m_prob[2] > 0.00 && m_prob[2] > m_prob[3])
1337 return true;
1338 else
1339 return false;
1340}
1341
1342bool SimplePIDSvc::iskaon()
1343{
1344 if (m_prob[3] > 0.00 && m_prob[3] > m_prob[2])
1345 return true;
1346 else
1347 return false;
1348}
1349
1350//bool SimplePIDSvc::iselectron_(bool emc)
1351//{
1352// if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
1353// {
1354// if (!emc)
1355// return true;
1356// else if (fabs(m_cost[0]) < 0.7 && m_emc_eop[0] > 0 && m_emc_eop[0] < 0.8)
1357// return false;
1358// else if (fabs(m_cost[0]) >= 0.7 && fabs(m_cost[0])<0.8 && m_emc_eop[0] > 0 && m_emc_eop[0] < -7.5*fabs(m_cost[0])+6.05)
1359// return false;
1360// else if (fabs(m_cost[0]) > 0.85 && m_emc_eop[0] > 0 && m_emc_eop[0] < 0.6)
1361// return false;
1362// else
1363// return true;
1364// }
1365// else
1366// return false;
1367//}
1368
1369bool SimplePIDSvc::iselectron(bool emc)
1370{
1371 if (!emc)
1372 {
1373 if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
1374 return true;
1375 else
1376 return false;
1377 }
1378 else
1379 {
1380 if (calEMCLikelihood())
1381 {
1382 if (m_lh_electron > m_eid_ratio)
1383 return true;
1384 else
1385 return false;
1386 }
1387 else
1388 {
1389 if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
1390 return true;
1391 else
1392 return false;
1393 }
1394 }
1395}
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
P
double P(RecMdcKalTrack *trk)
Definition: CalibEventSelect.cxx:87
mass
double mass
Definition: CosmicGenerator.cxx:138
sigma
TTree * sigma
Definition: CurConst_dEdXCal.cxx:6
q
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
Definition: KKsem.h:33
dt
TGraph2DErrors * dt
Definition: McCor.cxx:45
EPS
const double EPS
Definition: TRunge.cxx:43
t
TTree * t
Definition: binning.cxx:23
DstEmcShower::latMoment
double latMoment() const
Definition: DstEmcShower.h:52
DstEmcShower::eSeed
double eSeed() const
Definition: DstEmcShower.h:47
DstEmcShower::e3x3
double e3x3() const
Definition: DstEmcShower.h:48
DstEmcShower::secondMoment
double secondMoment() const
Definition: DstEmcShower.h:51
DstEmcShower::e5x5
double e5x5() const
Definition: DstEmcShower.h:49
DstEmcShower::energy
double energy() const
Definition: DstEmcShower.h:45
DstExtTrack::tof1Position
const Hep3Vector tof1Position() const
Definition: DstExtTrack.h:58
DstExtTrack::tof2Position
const Hep3Vector tof2Position() const
Definition: DstExtTrack.h:94
DstMdcDedx::chi
double chi(int i) const
Definition: DstMdcDedx.h:58
DstMdcKalTrack::theta
const double theta() const
Definition: DstMdcKalTrack.h:70
DstMdcKalTrack::setPidType
static void setPidType(PidType pidType)
Definition: DstMdcKalTrack.h:59
DstMdcKalTrack::p
const double p() const
Definition: DstMdcKalTrack.h:82
DstMdcKalTrack::charge
const int charge() const
Definition: DstMdcKalTrack.h:65
EvtRecTrack::isMdcDedxValid
bool isMdcDedxValid()
Definition: EvtRecTrack.h:45
EvtRecTrack::mdcDedx
RecMdcDedx * mdcDedx()
Definition: EvtRecTrack.h:55
EvtRecTrack::isExtTrackValid
bool isExtTrackValid()
Definition: EvtRecTrack.h:49
EvtRecTrack::extTrack
RecExtTrack * extTrack()
Definition: EvtRecTrack.h:56
EvtRecTrack::isMdcKalTrackValid
bool isMdcKalTrackValid()
Definition: EvtRecTrack.h:44
EvtRecTrack::tofTrack
SmartRefVector< RecTofTrack > tofTrack()
Definition: EvtRecTrack.h:57
EvtRecTrack::isTofTrackValid
bool isTofTrackValid()
Definition: EvtRecTrack.h:46
EvtRecTrack::emcShower
RecEmcShower * emcShower()
Definition: EvtRecTrack.h:58
EvtRecTrack::isEmcShowerValid
bool isEmcShowerValid()
Definition: EvtRecTrack.h:47
EvtRecTrack::mdcKalTrack
RecMdcKalTrack * mdcKalTrack()
Definition: EvtRecTrack.h:54
SimplePIDSvc::~SimplePIDSvc
virtual ~SimplePIDSvc()
Definition: SimplePIDSvc.cxx:30
SimplePIDSvc::getChi2
double getChi2(int i)
Definition: SimplePIDSvc.cxx:1233
SimplePIDSvc::preparePID
void preparePID(EvtRecTrack *track)
Definition: SimplePIDSvc.cxx:95
SimplePIDSvc::initialize
virtual StatusCode initialize()
Definition: SimplePIDSvc.cxx:32
SimplePIDSvc::iselectron
bool iselectron(bool emc=true)
Definition: SimplePIDSvc.cxx:1369
SimplePIDSvc::finalize
virtual StatusCode finalize()
Definition: SimplePIDSvc.cxx:47
TofHitStatus::is_barrel
bool is_barrel() const
Definition: TofHitStatus.h:26
TofHitStatus::layer
unsigned int layer() const
Definition: TofHitStatus.h:28
TofHitStatus::is_cluster
bool is_cluster() const
Definition: TofHitStatus.h:25
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition: TofHitStatus.cxx:134
TofHitStatus::is_counter
bool is_counter() const
Definition: TofHitStatus.h:24
TofHitStatus::is_readout
bool is_readout() const
Definition: TofHitStatus.h:23
TofHitStatus::is_raw
bool is_raw() const
Definition: TofHitStatus.h:22
EvtCyclic3::c_str
char * c_str(Index i)
Definition: EvtCyclic3.cc:252
bpkt::ifstream
std::ifstream ifstream
Definition: bpkt_streams.h:44
charge
float charge
Definition: prepare_data.cxx:28
bg
float bg
Definition: prepare_data.cxx:28
b
const double b
Definition: slope.cxx:9