BOSS 7.1.1
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 const int RUN_BEGIN_DATA_23 = 74031;
241 const int RUN_END_DATA_23 = 78536;
242 const int RUN_BEGIN_MC_23 = -78536;
243 const int RUN_END_MC_23 = -74031;
244
245
246 if (run_no >= RUN_BEGIN_DATA_10 && run_no <= RUN_END_DATA_10)
247 return 0;
248 else if (run_no >= RUN_BEGIN_DATA_11 && run_no <= RUN_END_DATA_11)
249 return 1;
250 else if (run_no >= RUN_BEGIN_MC_10 && run_no <= RUN_END_MC_10)
251 return 2;
252 else if (run_no >= RUN_BEGIN_MC_11 && run_no <= RUN_END_MC_11)
253 return 3;
254
255 else if (run_no >= RUN_BEGIN_DATA_22 && run_no <= RUN_END_DATA_22)
256 return 1;
257 else if (run_no >= RUN_BEGIN_MC_22 && run_no <= RUN_END_MC_22)
258 return 3;
259 else if (run_no >= RUN_BEGIN_DATA_23 && run_no <= RUN_END_DATA_23)
260 return 1;
261 else if (run_no >= RUN_BEGIN_MC_23 && run_no <= RUN_END_MC_23)
262 return 3;
263 else
264 return -1;
265
266}
267
268void SimplePIDSvc::loadTOFInfo(EvtRecTrack *track)
269{
270 //Initialization
271 for (unsigned int i = 0; i < 8; i++)
272 {
273 for (unsigned int j = 0; j < 5; j++)
274 m_tof_dt[i][j] = -99.;
275 m_tof_ph[i] = -99.;
276 }
277 for (unsigned int i = 0; i < 2; i++)
278 {
279 m_tof_zr[i] = -9999.;
280 m_tof_counter[i] = -1;
281 }
282 for (unsigned int i = 0; i < 5; i++)
283 {
284 m_tof_chi[i] = -99.;
285 }
286 m_tof_barrel = -1;
287
288 if (!track->isExtTrackValid() || !track->isTofTrackValid()) return;
289
290 SmartRefVector<RecTofTrack> tofTrk = track->tofTrack();
291 SmartRefVector<RecTofTrack>::iterator it;
292 RecExtTrack* extTrk = track->extTrack();
293 double zrhit[2];
294 zrhit[0] = extTrk->tof1Position().z();
295 zrhit[1] = extTrk->tof2Position().z();
296
297 TofHitStatus *hitst = new TofHitStatus;
298
299 for (it = tofTrk.begin(); it != tofTrk.end(); it++)
300 {
301 unsigned int st = (*it)->status();
302 hitst->setStatus(st);
303 if (hitst->is_raw()) continue; //empty TOF hit
304 bool barrel = hitst->is_barrel();
305 bool readout = hitst->is_readout();
306 bool counter = hitst->is_counter();
307 bool cluster = hitst->is_cluster();
308 int layer = hitst->layer();
309 double tof = (*it)->tof();
310 double ph = (*it)->ph();
311 m_tof_counter[layer-1] = (*it)->tofID();
312
313 if (barrel)
314 {
315 m_tof_barrel = 1;
316 }
317 else
318 {
319 m_tof_barrel = 0;
320 zrhit[0] = extTrk->tof1Position().rho();
321 }
322 m_tof_zr[0] = zrhit[0];
323 m_tof_zr[1] = zrhit[1];
324
325 int idx = -1;
326 if (readout)
327 {
328 if (barrel)
329 idx = ((st & 0xC0) >> 5) + (((st ^ 0x20) & 0x20) >> 5) - 2;
330 else
331 idx = 7;
332 }
333 else if (counter)
334 {
335 idx = layer + 3;
336 }
337 else if (cluster)
338 {
339 idx = 6;
340 }
341 if (idx == -1) continue;
342 m_tof_ph[idx] = ph;
343 for (unsigned int i = 0; i < 5; i++)
344 {
345 double offset = (*it)->toffset(i);
346 double texp = (*it)->texp(i);
347 if (texp < 0.0) continue;
348 double dt = tof - offset - texp;
349 m_tof_dt[idx][i] = dt;
350 }
351 }
352 delete hitst;
353}
354
355void SimplePIDSvc::tofBarrelCorrection()
356{
357 const double EPS = 1e-4;
358 const double BG_LOW = 0.20;
359 const double BG_UP = 7.40;
360 const double COST_LOW = -0.81;
361 const double COST_UP = 0.81;
362 const double Q_LOW = 0.;
363 const double Q_UP = 9000.;
364 const double P_LOW = 0.2;
365 const double P_UP = 1.3;
366 const int BIN_BG = 15;
367 const int BIN_COST = 15;
368 const int BIN_P = 15;
369 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};
370 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};
371 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};
372 int idx = getRunIdx(m_run);
373
374 if (idx != -1)
375 {
376 for (unsigned int i = 0; i < 4; i++)
377 {// only correct e, pi, K
378 double bg;
379 int bin_bg, bin_cost, bin_wgt;
380 int pid;
381 if (i == 0)
382 {
383 bg = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
384 bin_bg = findBin(P, BIN_P, bg);
385 pid = 0;
386 }
387 else if (i == 2 || i == 3)
388 {
389 bg = max(BG_LOW+EPS, min(BG_UP-EPS, m_betagamma[i]));
390 bin_bg = findBin(BG, BIN_BG, bg);
391 pid = 1;
392 }
393 else
394 {
395 continue;
396 }
397 double cost = m_cost[i];
398 int charge = m_charge[i];
399 double t[5], q;
400 double offset, sigma;
401 double offset_q, offset_bgcost;
402 int flag[4] = {0, 0, 0, 0, };
403 cost = max(COST_LOW+EPS, min(COST_UP-EPS, cost));
404 bin_cost = findBin(COST, BIN_COST, cost);
405 if (bin_bg == -1 || bin_cost == -1) continue;
406
407 //corrections
408 for (unsigned int j = 0; j < 4; j++)
409 {
410 t[j] = m_tof_dt[j][i];
411 if (fabs(t[j] + 99.) < EPS)//no readout
412 flag[j] = 0;
413 else
414 flag[j] = 1;
415 q = m_tof_ph[j];
416 q = max(Q_LOW+EPS, min(Q_UP-EPS, q));
417 if (charge == 1)
418 {
419 offset_q = h_tof_p_q_offset[pid][idx][j]->Interpolate( q );
420 offset_bgcost = h_tof_p_bgcost_offset[pid][idx][j]->Interpolate( bg, cost );
421 t[j] = t[j] - offset_q - offset_bgcost;
422 }
423 else
424 {
425 offset_q = h_tof_m_q_offset[pid][idx][j]->Interpolate( q );
426 offset_bgcost = h_tof_m_bgcost_offset[pid][idx][j]->Interpolate( bg, cost );
427 t[j] = t[j] - offset_q - offset_bgcost;
428 }
429 }
430 bin_wgt = flag[0]*8 + flag[1]*4 + flag[2]*2 + flag[3] - 1;
431 if (bin_wgt == -1) continue;
432 t[4] = 0;
433 for (unsigned int j = 0; j < 4; j++)
434 {
435 if (charge == 1)
436 t[4] += t[j] * h_tof_p_wgt[pid][idx][bin_wgt][j]->GetBinContent( bin_bg+1, bin_cost+1 );
437 else
438 t[4] += t[j] * h_tof_m_wgt[pid][idx][bin_wgt][j]->GetBinContent( bin_bg+1, bin_cost+1 );
439 }
440 if (charge == 1)
441 {
442 t[4] /= h_tof_p_wgt[pid][idx][bin_wgt][4]->GetBinContent( bin_bg+1, bin_cost+1 );
443 offset = h_tof_p_final_offset[pid][idx][bin_wgt]->Interpolate( bg, cost );
444 sigma = h_tof_p_final_sigma[pid][idx][bin_wgt]-> Interpolate( bg, cost );
445 m_tof_chi[i] = (t[4] - offset) / sigma;
446 }
447 else
448 {
449 t[4] /= h_tof_m_wgt[pid][idx][bin_wgt][4]->GetBinContent( bin_bg+1, bin_cost+1 );
450 offset = h_tof_m_final_offset[pid][idx][bin_wgt]->Interpolate( bg, cost );
451 sigma = h_tof_m_final_sigma[pid][idx][bin_wgt]-> Interpolate( bg, cost );
452 m_tof_chi[i] = (t[4] - offset) / sigma;
453 }
454 }
455 }
456}
457
458void SimplePIDSvc::tofEndcapCorrection()
459{
460 const double EPS = 1e-4;
461 const double BG_LOW = 0.30;
462 const double BG_UP = 7.40;
463 const double Q_LOW = 0.;
464 const double Q_UP = 6000.;
465 const double COST_EAST_LOW = 0.720;
466 const double COST_EAST_UP = 0.930;
467 const double COST_WEST_LOW = -0.930;
468 const double COST_WEST_UP = -0.720;
469 const double P_LOW = 0.2;
470 const double P_UP = 1.3;
471
472 int idx = getRunIdx(m_run);
473
474 if (idx != -1)
475 {
476 for (unsigned int i = 0; i < 4; i++)
477 {// only correct e, pi, K
478 int pid;
479 double bg;
480 if (i == 0)
481 {
482 bg = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
483 pid = 0;
484 }
485 else if (i == 2 || i == 3)
486 {
487 bg = max(BG_LOW+EPS, min(BG_UP-EPS, m_betagamma[i]));
488 pid = 1;
489 }
490 else
491 {
492 continue;
493 }
494
495 int flag; //0:east, 1:west
496 double cost = m_cost[i];
497 int charge = m_charge[i];
498 double t = m_tof_dt[7][i];
499 double q = m_tof_ph[7];
500 double off_q, off_bg, off_cost;
501 double sg_q, sg_bg, sg_cost;
502 if (cost > 0)
503 {
504 flag = 0;
505 cost = max(COST_EAST_LOW+EPS, min(COST_EAST_UP-EPS, cost));
506 }
507 else
508 {
509 flag = 1;
510 cost = max(COST_WEST_LOW+EPS, min(COST_WEST_UP-EPS, cost));
511 }
512 q = max(Q_LOW+EPS, min(Q_UP-EPS, q));
513
514 //corrections
515 if (charge == 1)
516 {
517 off_q = h_tofec_p_q_offset[pid][idx][flag] ->Interpolate( q );
518 sg_q = h_tofec_p_q_sigma[pid][idx][flag] ->Interpolate( q );
519 off_bg = h_tofec_p_bg_offset[pid][idx][flag] ->Interpolate( bg );
520 sg_bg = h_tofec_p_bg_sigma[pid][idx][flag] ->Interpolate( bg );
521 off_cost = h_tofec_p_cost_offset[pid][idx][flag]->Interpolate( cost );
522 sg_cost = h_tofec_p_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 else
526 {
527 off_q = h_tofec_m_q_offset[pid][idx][flag] ->Interpolate( q );
528 sg_q = h_tofec_m_q_sigma[pid][idx][flag] ->Interpolate( q );
529 off_bg = h_tofec_m_bg_offset[pid][idx][flag] ->Interpolate( bg );
530 sg_bg = h_tofec_m_bg_sigma[pid][idx][flag] ->Interpolate( bg );
531 off_cost = h_tofec_m_cost_offset[pid][idx][flag]->Interpolate( cost );
532 sg_cost = h_tofec_m_cost_sigma[pid][idx][flag] ->Interpolate( cost );
533 m_tof_chi[i] = (((t - off_q) / sg_q - off_bg) / sg_bg - off_cost) / sg_cost;
534 }
535 }
536 }
537}
538
539void SimplePIDSvc::loadDedxInfo(EvtRecTrack *track)
540{
541 if (track->isMdcDedxValid())
542 {
543 RecMdcDedx* dedx_trk = track->mdcDedx();
544 for (unsigned int i = 0; i < 5; i++)
545 m_dedx_chi[i] = dedx_trk->chi(i);
546 }
547 else
548 {
549 for (unsigned int i = 0; i < 5; i++)
550 m_dedx_chi[i] = -99;
551 }
552}
553
554void SimplePIDSvc::dedxCorrection()
555{
556 int idx = getRunIdx(m_run);
557 const double EPS = 1e-4;
558 const double BG_LOW = 0.20;
559 const double BG_UP = 7.40;
560 const double COST_LOW = -0.93;
561 const double COST_UP = 0.93;
562 const double P_LOW = 0.2;
563 const double P_UP = 1.3;
564 if (idx != -1)
565 {
566 double offset, sigma;
567 for (unsigned int i = 0; i < 4; i++)
568 {// only correct e, pi, K
569 double bg;
570 int pid;
571 if (i == 0)
572 {
573 bg = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
574 pid = 0;
575 }
576 else if (i == 2 || i == 3)
577 {
578 bg = max(BG_LOW+EPS, min(BG_UP-EPS, m_betagamma[i]));
579 pid = 1;
580 }
581 else
582 {
583 continue;
584 }
585 double cost = m_cost[i];
586 double charge = m_charge[i];
587 cost = max(COST_LOW+EPS, min(COST_UP-EPS, cost));
588 if (charge == 1)
589 {
590 offset = h_dedx_p_offset[pid][idx]->Interpolate( bg, cost );
591 sigma = h_dedx_p_sigma[pid][idx] ->Interpolate( bg, cost );
592 m_dedx_chi[i] = (m_dedx_chi[i] - offset) / sigma;
593 }
594 else
595 {
596 offset = h_dedx_m_offset[pid][idx]->Interpolate( bg, cost );
597 sigma = h_dedx_m_sigma[pid][idx] ->Interpolate( bg, cost );
598 m_dedx_chi[i] = (m_dedx_chi[i] - offset) / sigma;
599 }
600 }
601 }
602}
603
604void SimplePIDSvc::tofBarrelSecondCorrection()
605{
606
607 int idx = getRunIdx(m_run);
608 const double EPS = 1e-4;
609 const double P_LOW = 0.0;
610 const double P_UP = 1.3;
611 const int BIN_P = 10;
612
613 const int RUN_BEGIN_DATA_10 = 11414;
614 const int RUN_END_DATA_10 = 14604;
615 const int RUN_BEGIN_MC_10 = -14604;
616 const int RUN_END_MC_10 = -11414;
617 const int RUN_BEGIN_DATA_11 = 20448;
618 const int RUN_END_DATA_11 = 23454;
619 const int RUN_BEGIN_MC_11 = -23454;
620 const int RUN_END_MC_11 = -20448;
621
622 const int RUN_BEGIN_DATA_22 = 70521;
623 const int RUN_END_DATA_22 = 73930;
624 const int RUN_BEGIN_MC_22 = -73930;
625 const int RUN_END_MC_22 = -70521;
626
627 const int RUN_BEGIN_DATA_23 = 74031;
628 const int RUN_END_DATA_23 = 78536;
629 const int RUN_BEGIN_MC_23 = -78536;
630 const int RUN_END_MC_23 = -74031;
631
632 double 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};
633 if (idx != -1)
634 {
635
636 for (unsigned int i = 2; i < 4; i++){// second correct for tof of pi k
637
638 int aa=99,bb=99;
639
640 int bin_p;
641 double ptk = m_p[i];
642 ptk = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
643 bin_p = findBin(P, BIN_P, ptk);
644
645 if(i==2){// pi
646
647 if((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))
648 {aa=2; bb=2;}// round0304 data
649
650 if((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))
651 {aa=2; bb=3;}// round0304 inc
652
653 if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
654 {aa=3; bb=2;}// round15 data
655
656 if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
657 {aa=3; bb=3;}// round15 inc
658
659 if(m_run>=RUN_BEGIN_DATA_23 && m_run<=RUN_END_DATA_23 && idx == 1)
660 {aa=7; bb=2;}// round16 data added by Yijia Zeng
661
662 if(m_run>=RUN_BEGIN_MC_23 && m_run<=RUN_END_MC_23 && idx == 3)
663 {aa=7; bb=3;}// round16 inc added by Yijia Zeng
664
665 }
666
667 if(i==3){// k
668
669 if((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))
670 {aa=2; bb=0;}
671
672 if((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))
673 {aa=2; bb=1;}
674
675 if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
676 {aa=3; bb=0;}
677
678 if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
679 {aa=3; bb=1;}
680
681 if(m_run>=RUN_BEGIN_DATA_23 && m_run<=RUN_END_DATA_23 && idx == 1)
682 {aa=7; bb=0;}// round16 data added by Yijia Zeng
683
684 if(m_run>=RUN_BEGIN_MC_23 && m_run<=RUN_END_MC_23 && idx == 3)
685 {aa=7; bb=1;}// round16 inc added by Yijia Zeng
686
687 }
688
689 if(m_tof_chi[i]!=-99 && m_run>0){
690
691 //NO Correction on the higher SIX momentum region from Kaikai He
692 if(!(m_run>=RUN_BEGIN_DATA_23 && m_run<=RUN_END_DATA_23)){
693 if(bin_p<4){
694 m_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]);
695 }
696 if(bin_p>=4){
697 m_tof_chi[i] = (m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p]);
698 }
699 }
700 //NO Correction on the higher FIVE momentum region from Yijia for round 16 data sample
701 //The Corrected Distribution is somehow worse than the Original one
702 else{
703 if(bin_p<5){
704 m_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]);
705 }
706 if(bin_p>=5){
707 m_tof_chi[i] = (m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p]);
708 }
709 }
710
711 }
712 if(m_tof_chi[i]!=-99 && m_run<0){
713 m_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]);
714 }
715
716 }
717
718 }//end idx!= -1
719
720}
721
722void SimplePIDSvc::tofEndcapSecondCorrection()
723{
724
725 int idx = getRunIdx(m_run);
726 const double EPS = 1e-4;
727 const double P_LOW = 0.0;
728 const double P_UP = 1.3;
729 const int BIN_P = 10;
730
731 const int RUN_BEGIN_DATA_10 = 11414;
732 const int RUN_END_DATA_10 = 14604;
733 const int RUN_BEGIN_MC_10 = -14604;
734 const int RUN_END_MC_10 = -11414;
735 const int RUN_BEGIN_DATA_11 = 20448;
736 const int RUN_END_DATA_11 = 23454;
737 const int RUN_BEGIN_MC_11 = -23454;
738 const int RUN_END_MC_11 = -20448;
739
740 const int RUN_BEGIN_DATA_22 = 70521;
741 const int RUN_END_DATA_22 = 73930;
742 const int RUN_BEGIN_MC_22 = -73930;
743 const int RUN_END_MC_22 = -70521;
744
745 const int RUN_BEGIN_DATA_23 = 74031;
746 const int RUN_END_DATA_23 = 78536;
747 const int RUN_BEGIN_MC_23 = -78536;
748 const int RUN_END_MC_23 = -74031;
749
750 double 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};
751 if (idx != -1)
752 {
753
754 for (unsigned int i = 2; i < 4; i++){// second correct for tof of pi k
755
756 int aa=99,bb=99;
757
758 int bin_p;
759 double ptk = m_p[i];
760 ptk = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
761 bin_p = findBin(P, BIN_P, ptk);
762
763 if(i==2){// pi
764
765 if((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))
766 {aa=4; bb=2;}// round0304 data endcap
767
768 if((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))
769 {aa=4; bb=3;}// round0304 inc endcap
770
771 if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
772 {aa=5; bb=2;}// round15 data endcap
773
774 if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
775 {aa=5; bb=3;}// round15 inc endcap
776
777 if(m_run>=RUN_BEGIN_DATA_23 && m_run<=RUN_END_DATA_23 && idx == 1)
778 {aa=8; bb=2;}// round16 data endcap added by Yijia Zeng
779
780 if(m_run>=RUN_BEGIN_MC_23 && m_run<=RUN_END_MC_23 && idx == 3)
781 {aa=8; bb=3;}// round16 inc endcap added by Yijia Zeng
782
783 }
784
785 if(i==3){// k
786
787 if((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))
788 {aa=4; bb=0;}
789
790 if((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))
791 {aa=4; bb=1;}
792
793 if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
794 {aa=5; bb=0;}
795
796 if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
797 {aa=5; bb=1;}
798
799 if(m_run>=RUN_BEGIN_DATA_23 && m_run<=RUN_END_DATA_23 && idx == 1)
800 {aa=8; bb=0;}// round16 data endcap added by Yijia Zeng
801
802 if(m_run>=RUN_BEGIN_MC_23 && m_run<=RUN_END_MC_23 && idx == 3)
803 {aa=8; bb=1;}// round16 inc endcap added by Yijia Zeng
804
805 }
806
807 if(m_tof_chi[i]!=-99 && ( aa==5 || aa==8 )){// for round15 and round16 data and inc
808
809 m_tof_chi[i] = (m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p]) / (60./110.);
810
811 }
812
813 if(m_tof_chi[i]!=-99 && aa==4){// for round0304 data and inc
814
815 m_tof_chi[i] = (m_tof_chi[i] - m_gaussion_mean[aa][bb][bin_p]);
816
817 }
818
819
820 }
821
822 }//end idx!= -1
823
824
825}
826
827void SimplePIDSvc::dedxSecondCorrection()
828
829{
830
831 int idx = getRunIdx(m_run);
832 const double EPS = 1e-4;
833 const double P_LOW = 0.0;
834 const double P_UP = 1.3;
835 const int BIN_P = 10;
836
837 const int RUN_BEGIN_DATA_10 = 11414;
838 const int RUN_END_DATA_10 = 14604;
839 const int RUN_BEGIN_MC_10 = -14604;
840 const int RUN_END_MC_10 = -11414;
841 const int RUN_BEGIN_DATA_11 = 20448;
842 const int RUN_END_DATA_11 = 23454;
843 const int RUN_BEGIN_MC_11 = -23454;
844 const int RUN_END_MC_11 = -20448;
845
846 const int RUN_BEGIN_DATA_22 = 70521;
847 const int RUN_END_DATA_22 = 73930;
848 const int RUN_BEGIN_MC_22 = -73930;
849 const int RUN_END_MC_22 = -70521;
850
851 const int RUN_BEGIN_DATA_23 = 74031;
852 const int RUN_END_DATA_23 = 78536;
853 const int RUN_BEGIN_MC_23 = -78536;
854 const int RUN_END_MC_23 = -74031;
855
856 double 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};
857 if (idx != -1)
858 {
859
860 for (unsigned int i = 2; i < 4; i++){// second correct for dedx of pi k
861
862 int aa=99,bb=99;
863
864 int bin_p;
865 double ptk = m_p[i];
866 ptk = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
867 bin_p = findBin(P, BIN_P, ptk);
868
869 if(i==2){// pi
870
871 if((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))
872 {aa=0; bb=2;}// round0304 data
873
874 if((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))
875 {aa=0; bb=3;}// round0304 inc
876
877 if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
878 {aa=1; bb=2;}// round15 data
879
880 if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
881 {aa=1; bb=3;}// round15 inc
882
883 if(m_run>=RUN_BEGIN_DATA_23 && m_run<=RUN_END_DATA_23 && idx == 1)
884 {aa=6; bb=2;}// round16 data added by Yijia Zeng
885
886 if(m_run>=RUN_BEGIN_MC_23 && m_run<=RUN_END_MC_23 && idx == 3)
887 {aa=6; bb=3;}// round16 inc added by Yijia Zeng
888
889 }
890
891 if(i==3){// k
892
893 if((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))
894 {aa=0; bb=0;}
895
896 if((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))
897 {aa=0; bb=1;}
898
899 if(m_run>=RUN_BEGIN_DATA_22 && m_run<=RUN_END_DATA_22 && idx == 1)
900 {aa=1; bb=0;}
901
902 if(m_run>=RUN_BEGIN_MC_22 && m_run<=RUN_END_MC_22 && idx == 3)
903 {aa=1; bb=1;}
904
905 if(m_run>=RUN_BEGIN_DATA_23 && m_run<=RUN_END_DATA_23 && idx == 1)
906 {aa=6; bb=0;}// round16 data added by Yijia Zeng
907
908 if(m_run>=RUN_BEGIN_MC_23 && m_run<=RUN_END_MC_23 && idx == 3)
909 {aa=6; bb=1;}// round16 inc added by Yijia Zeng
910
911 }
912
913 if(m_dedx_chi[i]!=-99 && m_run>0){
914 m_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]);
915 }
916 if(m_dedx_chi[i]!=-99 && m_run<0){
917 m_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]);
918 }
919
920 }
921
922 }//end idx!= -1
923
924}
925
926void SimplePIDSvc::loadSecondPar()
927{
928 string path = getenv("SIMPLEPIDSVCROOT");
929
930 for(int i=0; i<9; i++){
931
932 const char *dir;
933
934 if(i == 0) dir = "round0304_dedx";
935 else if(i == 1) dir = "round15_dedx";
936 else if(i == 2) dir = "round0304_tof";
937 else if(i == 3) dir = "round15_tof";
938 else if(i == 4) dir = "round0304_tof_endcap";
939 else if(i == 5) dir = "round15_tof_endcap";
940 //Added by Yijia Zeng for round 16 data sample
941 else if(i == 6) dir = "round16_dedx";
942 else if(i == 7) dir = "round16_tof";
943 else if(i == 8) dir = "round16_tof_endcap";
944
945 else{
946 cout << "Boundary Error! " << endl;
947 exit(1);
948 }
949
950 for(int j=0; j<4; j++){
951
952 const char *name;
953
954 if(j == 0) name = "data_K";
955 else if(j == 1) name = "inc_K";
956 else if(j == 2) name = "data_pi";
957 else if(j == 3) name = "inc_pi";
958 else{
959 cout << "Boundary Error! " << endl;
960 exit(1);
961 }
962
963 ifstream second_cor( path + Form("/share/second_correct/%s/%s_%s.dat",dir,dir,name));
964
965 for(int m=0; m<10; m++){
966
967 second_cor>>m_gaussion_mean[i][j][m]>>m_gaussion_sigma[i][j][m]>>m_gaussion_sigmab[i][j][m];
968
969 }
970
971 second_cor.close();
972
973 }
974 }
975
976
977}
978
979void SimplePIDSvc::loadHistogram()
980{
981 string path = getenv("SIMPLEPIDSVCROOT");
982 vector<string> filename;
983 for (unsigned int idx = 0; idx < 2; idx++)
984 {
985 const char *dir;
986 if (idx == 0)
987 dir = "electron";
988 else if (idx == 1)
989 dir = "kpi";
990 else
991 {
992 cout << "Boundary Error! " << endl;
993 exit(1);
994 }
995
996 //dedx
997 filename.clear();
998 filename.push_back( path + Form("/share/%s/dedx/dedx_d10.root", dir) );
999 filename.push_back( path + Form("/share/%s/dedx/dedx_d11.root", dir) );
1000 filename.push_back( path + Form("/share/%s/dedx/dedx_m10.root", dir) );
1001 filename.push_back( path + Form("/share/%s/dedx/dedx_m11.root", dir) );
1002 for (unsigned int i = 0; i < filename.size(); i++)
1003 {
1004 f_dedx[idx][i] = new TFile(filename[i].c_str(), "READ");
1005 const char *name;
1006 if (i == 0)
1007 name = "d10";
1008 else if (i == 1)
1009 name = "d11";
1010 else if (i == 2)
1011 name = "m10";
1012 else if (i == 3)
1013 name = "m11";
1014 else
1015 {
1016 cout << "Boundary Error! " << endl;
1017 exit(1);
1018 }
1019 h_dedx_p_offset[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_p_offset_%s", name) );
1020 h_dedx_p_sigma[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_p_sigma_%s" , name) );
1021 h_dedx_m_offset[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_m_offset_%s", name) );
1022 h_dedx_m_sigma[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_m_sigma_%s" , name) );
1023 }
1024 //tof_barrel q
1025 filename.clear();
1026 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_d10.root", dir) );
1027 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_d11.root", dir) );
1028 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_m10.root", dir) );
1029 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_m11.root", dir) );
1030 for (unsigned int i = 0; i < filename.size(); i++)
1031 {
1032 f_tof_q[idx][i] = new TFile(filename[i].c_str(), "READ");
1033 const char *name;
1034 if (i == 0)
1035 name = "d10";
1036 else if (i == 1)
1037 name = "d11";
1038 else if (i == 2)
1039 name = "m10";
1040 else if (i == 3)
1041 name = "m11";
1042 else
1043 {
1044 cout << "Boundary Error! " << endl;
1045 exit(1);
1046 }
1047 for (unsigned int j = 0; j < 4; j++)
1048 {
1049 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) );
1050 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) );
1051 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) );
1052 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) );
1053 }
1054 }
1055 //tof_barrel bg&cost
1056 filename.clear();
1057 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_d10.root", dir) );
1058 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_d11.root", dir) );
1059 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_m10.root", dir) );
1060 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_m11.root", dir) );
1061 for (unsigned int i = 0; i < filename.size(); i++)
1062 {
1063 f_tof_bgcost[idx][i] = new TFile(filename[i].c_str(), "READ");
1064 const char *name;
1065 if (i == 0)
1066 name = "d10";
1067 else if (i == 1)
1068 name = "d11";
1069 else if (i == 2)
1070 name = "m10";
1071 else if (i == 3)
1072 name = "m11";
1073 else
1074 {
1075 cout << "Boundary Error! " << endl;
1076 exit(1);
1077 }
1078 for (unsigned int j = 0; j < 4; j++)
1079 {
1080 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) );
1081 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) );
1082 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) );
1083 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) );
1084 }
1085 }
1086 //tof_barrel wgt
1087 filename.clear();
1088 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_d10.root", dir) );
1089 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_d11.root", dir) );
1090 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_m10.root", dir) );
1091 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_m11.root", dir) );
1092 for (unsigned int i = 0; i < filename.size(); i++)
1093 {
1094 f_tof_wgt[idx][i] = new TFile(filename[i].c_str(), "READ");
1095 const char *name;
1096 if (i == 0)
1097 name = "d10";
1098 else if (i == 1)
1099 name = "d11";
1100 else if (i == 2)
1101 name = "m10";
1102 else if (i == 3)
1103 name = "m11";
1104 else
1105 {
1106 cout << "Boundary Error! " << endl;
1107 exit(1);
1108 }
1109 for (unsigned int j = 0; j < 15; j++)
1110 {
1111 for (unsigned int k = 0; k < 5; k++)
1112 {
1113 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) );
1114 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) );
1115 }
1116 }
1117 }
1118 //tof_barrel corr
1119 filename.clear();
1120 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_d10.root", dir) );
1121 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_d11.root", dir) );
1122 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_m10.root", dir) );
1123 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_m11.root", dir) );
1124 for (unsigned int i = 0; i < filename.size(); i++)
1125 {
1126 f_tof_final[idx][i] = new TFile(filename[i].c_str(), "READ");
1127 const char *name;
1128 if (i == 0)
1129 name = "d10";
1130 else if (i == 1)
1131 name = "d11";
1132 else if (i == 2)
1133 name = "m10";
1134 else if (i == 3)
1135 name = "m11";
1136 else
1137 {
1138 cout << "Boundary Error! " << endl;
1139 exit(1);
1140 }
1141 for (unsigned int j = 0; j < 15; j++)
1142 {
1143 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) );
1144 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) );
1145 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) );
1146 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) );
1147 }
1148 }
1149 //tof_endcap q
1150 filename.clear();
1151 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_d10.root", dir) );
1152 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_d11.root", dir) );
1153 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_m10.root", dir) );
1154 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_m11.root", dir) );
1155 for (unsigned int i = 0; i < filename.size(); i++)
1156 {
1157 f_tofec_q[idx][i] = new TFile(filename[i].c_str(), "READ");
1158 const char *name;
1159 if (i == 0)
1160 name = "d10";
1161 else if (i == 1)
1162 name = "d11";
1163 else if (i == 2)
1164 name = "m10";
1165 else if (i == 3)
1166 name = "m11";
1167 else
1168 {
1169 cout << "Boundary Error! " << endl;
1170 exit(1);
1171 }
1172 for (unsigned int j = 0; j < 2; j++)
1173 {
1174 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) );
1175 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) );
1176 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) );
1177 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) );
1178 }
1179 }
1180 //tof_endcap bg
1181 filename.clear();
1182 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_d10.root", dir) );
1183 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_d11.root", dir) );
1184 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_m10.root", dir) );
1185 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_m11.root", dir) );
1186 for (unsigned int i = 0; i < filename.size(); i++)
1187 {
1188 f_tofec_bg[idx][i] = new TFile(filename[i].c_str(), "READ");
1189 const char *name;
1190 if (i == 0)
1191 name = "d10";
1192 else if (i == 1)
1193 name = "d11";
1194 else if (i == 2)
1195 name = "m10";
1196 else if (i == 3)
1197 name = "m11";
1198 else
1199 {
1200 cout << "Boundary Error! " << endl;
1201 exit(1);
1202 }
1203 for (unsigned int j = 0; j < 2; j++)
1204 {
1205 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) );
1206 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) );
1207 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) );
1208 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) );
1209 }
1210 }
1211 //tof_endcap cost
1212 filename.clear();
1213 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_d10.root", dir) );
1214 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_d11.root", dir) );
1215 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_m10.root", dir) );
1216 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_m11.root", dir) );
1217 for (unsigned int i = 0; i < filename.size(); i++)
1218 {
1219 f_tofec_cost[idx][i] = new TFile(filename[i].c_str(), "READ");
1220 const char *name;
1221 if (i == 0)
1222 name = "d10";
1223 else if (i == 1)
1224 name = "d11";
1225 else if (i == 2)
1226 name = "m10";
1227 else if (i == 3)
1228 name = "m11";
1229 else
1230 {
1231 cout << "Boundary Error! " << endl;
1232 exit(1);
1233 }
1234 for (unsigned int j = 0; j < 2; j++)
1235 {
1236 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) );
1237 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) );
1238 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) );
1239 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) );
1240 }
1241 }
1242 }
1243 for (unsigned int idx = 0; idx < 3; idx++)
1244 {
1245 const char *dir;
1246 if (idx == 0)
1247 dir = "electron/emc";
1248 else if (idx == 1)
1249 dir = "kpi/emc_pion";
1250 else if (idx == 2)
1251 dir = "kpi/emc_kaon";
1252 else
1253 {
1254 cout << "Boundary Error! " << endl;
1255 exit(1);
1256 }
1257 //emc
1258 filename.clear();
1259 filename.push_back( path + Form("/share/%s/emc_d10.root", dir) );
1260 filename.push_back( path + Form("/share/%s/emc_d11.root", dir) );
1261 filename.push_back( path + Form("/share/%s/emc_m10.root", dir) );
1262 filename.push_back( path + Form("/share/%s/emc_m11.root", dir) );
1263 for (unsigned int i = 0; i < filename.size(); i++)
1264 {
1265 f_emc[idx][i] = new TFile(filename[i].c_str(), "READ");
1266 const char *name;
1267 if (i == 0)
1268 name = "d10";
1269 else if (i == 1)
1270 name = "d11";
1271 else if (i == 2)
1272 name = "m10";
1273 else if (i == 3)
1274 name = "m11";
1275 else
1276 {
1277 cout << "Boundary Error! " << endl;
1278 exit(1);
1279 }
1280 for (unsigned int j = 0; j < 15; j++)
1281 {
1282 for (unsigned int k = 0; k < 25; k++)
1283 {
1284 h_emc_ep[idx][i][j][k] = (TH1D*)f_emc[idx][i]->Get( Form("h_ep_%s_%d_%d", name, j, k) );
1285 h_emc_e35[idx][i][j][k] = (TH1D*)f_emc[idx][i]->Get( Form("h_e35_%s_%d_%d", name, j, k) );
1286 }
1287 }
1288 }
1289 }
1290 cout << "Successfully Return from Loading Initializations by package SimplePIDSvc ... " << endl;
1291}
1292
1293int SimplePIDSvc::findBin(double *a, int length, double value)
1294{
1295 for (int i = 0; i < length; i++)
1296 {
1297 if (value > a[i] && value <= a[i+1])
1298 {
1299 return i;
1300 }
1301 }
1302 if (value < a[0])
1303 {
1304 return 0;
1305 }
1306 else
1307 {
1308 return length;
1309 }
1310}
1311
1312double SimplePIDSvc::getChi2(int i)
1313{
1314 return pow(m_dedx_chi[i], 2) + pow(m_tof_chi[i], 2);
1315}
1316
1317void SimplePIDSvc::loadEMCInfo(EvtRecTrack *track)
1318{
1319 //Initialization
1320 for (unsigned int i = 0; i < 5; i++)
1321 {
1322 m_emc_eop[i] = -99.;
1323 m_emc_likelihood[i] = -99.;
1324 }
1325 m_emc_e = -99.;
1326 m_emc_e13 = -99.;
1327 m_emc_e35 = -99.;
1328 m_emc_sec = -99.;
1329 m_emc_lat = -99.;
1330 m_lh_electron = -99.;
1331
1332 if (!track->isEmcShowerValid()) return;
1333
1334 RecEmcShower* emc_trk = track->emcShower();
1335 m_emc_e = emc_trk->energy();
1336 for (unsigned int i = 0; i < 5; i++)
1337 {
1338 m_emc_eop[i] = m_emc_e / fabs(m_p[i]);
1339 }
1340 double eseed = emc_trk->eSeed();
1341 double e3 = emc_trk->e3x3();
1342 double e5 = emc_trk->e5x5();
1343 m_emc_sec = emc_trk->secondMoment() / 1000.;
1344 m_emc_lat = emc_trk->latMoment();
1345 if (e3 != 0)
1346 {
1347 m_emc_e13 = eseed / e3;
1348 }
1349 if (e5 != 0)
1350 {
1351 m_emc_e35 = e3 / e5;
1352 }
1353}
1354
1355bool SimplePIDSvc::calEMCLikelihood()
1356{
1357 if (m_emc_eop[0] < 0)
1358 return false;
1359
1360 int idx = getRunIdx(m_run);
1361 const Int_t BIN_P = 15;
1362 const Int_t BIN_COST = 25;
1363 const Int_t BIN_PID = 3;
1364 const double EPS = 1e-4;
1365 //electron
1366 double P[BIN_PID][BIN_P + 1] = {
1367 {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},
1368 {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},
1369 {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}, };
1370 double COST[BIN_PID][BIN_COST + 1] = {
1371 {-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},
1372 {-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},
1373 {-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}, };
1374
1375 double vp, vcost;
1376 int pid;
1377 int bin_p, bin_cost;
1378 for (unsigned int i = 0; i < 4; i++)
1379 {
1380 //only e, pi ,K
1381 if (i == 0)
1382 pid = 0;
1383 else if (i == 2)
1384 pid = 1;
1385 else if (i == 3)
1386 pid = 2;
1387 else
1388 continue;
1389
1390 vp = max(P[pid][0]+EPS, min(P[pid][BIN_P]-EPS, m_p[i]));
1391 vcost = max(COST[pid][0]+EPS, min(COST[pid][BIN_COST]-EPS, m_cost[i]));
1392 bin_p = findBin(P[pid], BIN_P, vp);
1393 bin_cost = findBin(COST[pid], BIN_COST, vcost);
1394
1395 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);
1396 }
1397 double a = m_prob[0] > 0 ? m_prob[0] : 0;
1398 double b = m_prob[2] > 0 ? m_prob[2] : 0;
1399 double c = m_prob[3] > 0 ? m_prob[3] : 0;
1400 double sum = a * m_emc_likelihood[0] + b * m_emc_likelihood[2] + c * m_emc_likelihood[3];
1401
1402 if (sum > 0 && m_prob[0] > 0)
1403 {
1404 m_lh_electron = m_prob[0] * m_emc_likelihood[0] / sum;
1405 return true;
1406 }
1407 else
1408 {
1409 return false;
1410 }
1411}
1412
1413bool SimplePIDSvc::ispion()
1414{
1415 if (m_prob[2] > 0.00 && m_prob[2] > m_prob[3])
1416 return true;
1417 else
1418 return false;
1419}
1420
1421bool SimplePIDSvc::iskaon()
1422{
1423 if (m_prob[3] > 0.00 && m_prob[3] > m_prob[2])
1424 return true;
1425 else
1426 return false;
1427}
1428
1429//bool SimplePIDSvc::iselectron_(bool emc)
1430//{
1431// if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
1432// {
1433// if (!emc)
1434// return true;
1435// else if (fabs(m_cost[0]) < 0.7 && m_emc_eop[0] > 0 && m_emc_eop[0] < 0.8)
1436// return false;
1437// 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)
1438// return false;
1439// else if (fabs(m_cost[0]) > 0.85 && m_emc_eop[0] > 0 && m_emc_eop[0] < 0.6)
1440// return false;
1441// else
1442// return true;
1443// }
1444// else
1445// return false;
1446//}
1447
1448bool SimplePIDSvc::iselectron(bool emc)
1449{
1450 if (!emc)
1451 {
1452 if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
1453 return true;
1454 else
1455 return false;
1456 }
1457 else
1458 {
1459 if (calEMCLikelihood())
1460 {
1461 if (m_lh_electron > m_eid_ratio)
1462 return true;
1463 else
1464 return false;
1465 }
1466 else
1467 {
1468 if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
1469 return true;
1470 else
1471 return false;
1472 }
1473 }
1474}
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:1312
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:1448
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
charge
float charge
Definition prepare_data.cxx:28
bg
float bg
Definition prepare_data.cxx:28
b
const double b
Definition slope.cxx:9