BOSS 7.0.5
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;
19
20SimplePIDSvc::SimplePIDSvc(const std::string& name, ISvcLocator* svcLoc) : Service(name, svcLoc)
21{
22 declareProperty("DedxChiCut", m_dedx_chi_cut = 4);
23 declareProperty("TofChiCut", m_tof_chi_cut = 4);
24 declareProperty("IsTofCorr", m_tof_corr = true);
25 declareProperty("IsDedxCorr", m_dedx_corr = true);
26 declareProperty("EidRatio", m_eid_ratio = 0.80);
27}
28
29SimplePIDSvc::~SimplePIDSvc() {;}
30
31StatusCode SimplePIDSvc::initialize()
32{
33 MsgStream log(messageService(), name());
34 log << MSG::INFO << "in SimplePIDSvc initialize()" << endreq;
35
36 StatusCode sc = Service::initialize();
37
38 sc = serviceLocator()->service("EventDataSvc", eventSvc_, true);
39
40 loadHistogram();
41
42 return sc;
43}
44
45StatusCode SimplePIDSvc::finalize()
46{
47 MsgStream log(messageService(), name());
48 log << MSG::INFO << "in SimplePIDSvc finalize()" << endreq;
49
50 StatusCode sc = Service::finalize();
51
52 for (unsigned int i = 0; i < 2; i++)
53 {
54 for (unsigned int j = 0; j < 4; j++)
55 {
56 f_dedx[i][j]->Close();
57 f_tof_q[i][j]->Close();
58 f_tof_bgcost[i][j]->Close();
59 f_tof_wgt[i][j]->Close();
60 f_tof_final[i][j]->Close();
61 f_tofec_q[i][j]->Close();
62 f_tofec_bg[i][j]->Close();
63 f_tofec_cost[i][j]->Close();
64 }
65 }
66 for (unsigned int i = 0; i < 3; i++)
67 {
68 for (unsigned int j = 0; j < 4; j++)
69 {
70 f_emc[i][j]->Close();
71 }
72 }
73
74 return sc;
75}
76
77StatusCode SimplePIDSvc::queryInterface(const InterfaceID& riid, void** ppvIF)
78{
79 if (ISimplePIDSvc::interfaceID().versionMatch(riid))
80 {
81 *ppvIF = dynamic_cast<ISimplePIDSvc*>(this);
82 }
83 else
84 {
85 return Service::queryInterface(riid, ppvIF);
86 }
87 addRef();
88 return StatusCode::SUCCESS;
89}
90
91void SimplePIDSvc::preparePID(EvtRecTrack* track)
92{
93
94 SmartDataPtr<Event::EventHeader> eventHeaderpid(eventSvc_, "/Event/EventHeader");
95 m_run = eventHeaderpid->runNumber();
96
97 if (track->isMdcKalTrackValid())
98 {
99 RecMdcKalTrack *mdckalTrk = track->mdcKalTrack();
100 RecMdcKalTrack::PidType trk_type[5] = {
101 RecMdcKalTrack::electron,
102 RecMdcKalTrack::muon,
103 RecMdcKalTrack::pion,
104 RecMdcKalTrack::kaon,
105 RecMdcKalTrack::proton,
106 };
107 double mass[5] = {
108 0.000511,
109 0.105658,
110 0.13957,
111 0.493677,
112 0.938272,
113 };
114 for(unsigned int pid = 0; pid < 5; pid++)
115 {
116 mdckalTrk->setPidType(trk_type[pid]);
117 m_p[pid] = mdckalTrk->p();
118 m_betagamma[pid] = m_p[pid] / mass[pid];
119 m_charge[pid] = mdckalTrk->charge();
120 m_cost[pid] = cos(mdckalTrk->theta());
121 }
122 }
123 else
124 {
125 for(unsigned int i = 0; i < 5; i++)
126 {
127 m_p[i] = -99;
128 m_betagamma[i] = -99;
129 m_cost[i] = -99;
130 m_charge[i] = 0;
131 }
132
133 }
134
135 //dE/dx PID
136 loadDedxInfo(track);
137 if (m_dedx_corr)
138 {
139 dedxCorrection();
140 }
141 //TOF PID
142 loadTOFInfo(track);
143 if (m_tof_corr)
144 {
145 if (m_tof_barrel == 1)
146 {
147 tofBarrelCorrection();
148 }
149 else if (m_tof_barrel == 0)
150 tofEndcapCorrection();
151 }
152 //EMC
153 loadEMCInfo(track);
154
155 calprob();
156}
157
158void SimplePIDSvc::calprob()
159{
160 bool usededx = false;
161 bool usetof = false;
162
163 for (unsigned int i = 0; i < 5 ;i++)
164 {
165 if (!usededx && fabs(m_dedx_chi[i]) < m_dedx_chi_cut)
166 usededx = true;
167 if (!usetof && fabs(m_tof_chi[i]) < m_tof_chi_cut)
168 usetof = true;
169
170 m_dedx_only[i] = false;
171 }
172 if (!usededx)
173 {
174 for(unsigned int i = 0; i < 5; i++)
175 m_dedx_chi[i] = -99;
176 }
177 if (!usetof)
178 {
179 for(unsigned int i = 0; i < 5; i++)
180 m_tof_chi[i] = -99;
181 }
182
183 for (unsigned int i = 0; i < 5; i++)
184 {
185 m_prob[i] = -99;
186 double chi2 = 0;
187 int ndf = 0;
188
189 if (usededx && usetof)
190 {
191 chi2 = pow(m_dedx_chi[i], 2) + pow(m_tof_chi[i], 2);
192 ndf = 2;
193 }
194 else if (usededx && !usetof)
195 {
196 chi2 = pow(m_dedx_chi[i], 2);
197 ndf = 1;
198 m_dedx_only[i] = true;
199 }
200 else if (!usededx && usetof)
201 {
202 chi2 = pow(m_tof_chi[i],2);
203 ndf = 1;
204 }
205 if (ndf > 0 && chi2 > 0)
206 m_prob[i] = TMath::Prob(chi2, ndf);
207 }
208}
209
210
211
212int SimplePIDSvc::getRunIdx(int run_no)
213{
214 // -1: beyond correction region
215 // 0: 2010 psi(3770) data
216 // 1: 2011 psi(3770) data
217 // 2: 2010 psi(3770) mc
218 // 3: 2011 psi(3770) mc
219 const int RUN_BEGIN_DATA_10 = 11414;
220 const int RUN_END_DATA_10 = 14604;
221 const int RUN_BEGIN_MC_10 = -14604;
222 const int RUN_END_MC_10 = -11414;
223 const int RUN_BEGIN_DATA_11 = 20448;
224 const int RUN_END_DATA_11 = 23454;
225 const int RUN_BEGIN_MC_11 = -23454;
226 const int RUN_END_MC_11 = -20448;
227 if (run_no >= RUN_BEGIN_DATA_10 && run_no <= RUN_END_DATA_10)
228 return 0;
229 else if (run_no >= RUN_BEGIN_DATA_11 && run_no <= RUN_END_DATA_11)
230 return 1;
231 else if (run_no >= RUN_BEGIN_MC_10 && run_no <= RUN_END_MC_10)
232 return 2;
233 else if (run_no >= RUN_BEGIN_MC_11 && run_no <= RUN_END_MC_11)
234 return 3;
235 else
236 return -1;
237}
238
239void SimplePIDSvc::loadTOFInfo(EvtRecTrack *track)
240{
241 //Initialization
242 for (unsigned int i = 0; i < 8; i++)
243 {
244 for (unsigned int j = 0; j < 5; j++)
245 m_tof_dt[i][j] = -99.;
246 m_tof_ph[i] = -99.;
247 }
248 for (unsigned int i = 0; i < 2; i++)
249 {
250 m_tof_zr[i] = -9999.;
251 m_tof_counter[i] = -1;
252 }
253 for (unsigned int i = 0; i < 5; i++)
254 {
255 m_tof_chi[i] = -99.;
256 }
257 m_tof_barrel = -1;
258
259 if (!track->isExtTrackValid() || !track->isTofTrackValid()) return;
260
261 SmartRefVector<RecTofTrack> tofTrk = track->tofTrack();
262 SmartRefVector<RecTofTrack>::iterator it;
263 RecExtTrack* extTrk = track->extTrack();
264 double zrhit[2];
265 zrhit[0] = extTrk->tof1Position().z();
266 zrhit[1] = extTrk->tof2Position().z();
267
268 TofHitStatus *hitst = new TofHitStatus;
269
270 for (it = tofTrk.begin(); it != tofTrk.end(); it++)
271 {
272 unsigned int st = (*it)->status();
273 hitst->setStatus(st);
274 if (hitst->is_raw()) continue; //empty TOF hit
275 bool barrel = hitst->is_barrel();
276 bool readout = hitst->is_readout();
277 bool counter = hitst->is_counter();
278 bool cluster = hitst->is_cluster();
279 int layer = hitst->layer();
280 double tof = (*it)->tof();
281 double ph = (*it)->ph();
282 m_tof_counter[layer-1] = (*it)->tofID();
283
284 if (barrel)
285 {
286 m_tof_barrel = 1;
287 }
288 else
289 {
290 m_tof_barrel = 0;
291 zrhit[0] = extTrk->tof1Position().rho();
292 }
293 m_tof_zr[0] = zrhit[0];
294 m_tof_zr[1] = zrhit[1];
295
296 int idx = -1;
297 if (readout)
298 {
299 if (barrel)
300 idx = ((st & 0xC0) >> 5) + (((st ^ 0x20) & 0x20) >> 5) - 2;
301 else
302 idx = 7;
303 }
304 else if (counter)
305 {
306 idx = layer + 3;
307 }
308 else if (cluster)
309 {
310 idx = 6;
311 }
312 if (idx == -1) continue;
313 m_tof_ph[idx] = ph;
314 for (unsigned int i = 0; i < 5; i++)
315 {
316 double offset = (*it)->toffset(i);
317 double texp = (*it)->texp(i);
318 if (texp < 0.0) continue;
319 double dt = tof - offset - texp;
320 m_tof_dt[idx][i] = dt;
321 }
322 //cout << "barrel: " << barrel << "\t" << m_tof_barrel << "\t";
323 //cout << "idx: " << idx << "\t";
324 //cout << "m_dt:" << m_tof_dt[idx][2] << endl;
325 }
326 delete hitst;
327}
328
329void SimplePIDSvc::tofBarrelCorrection()
330{
331 const double EPS = 1e-4;
332 const double BG_LOW = 0.20;
333 const double BG_UP = 7.40;
334 const double COST_LOW = -0.81;
335 const double COST_UP = 0.81;
336 const double Q_LOW = 0.;
337 const double Q_UP = 9000.;
338 const double P_LOW = 0.2;
339 const double P_UP = 1.3;
340 const int BIN_BG = 15;
341 const int BIN_COST = 15;
342 const int BIN_P = 15;
343 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};
344 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};
345 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};
346 int idx = getRunIdx(m_run);
347
348 if (idx != -1)
349 {
350 for (unsigned int i = 0; i < 4; i++)
351 {// only correct e, pi, K
352 double bg;
353 int bin_bg, bin_cost, bin_wgt;
354 int pid;
355 if (i == 0)
356 {
357 bg = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
358 bin_bg = findBin(P, BIN_P, bg);
359 pid = 0;
360 }
361 else if (i == 2 || i == 3)
362 {
363 bg = max(BG_LOW+EPS, min(BG_UP-EPS, m_betagamma[i]));
364 bin_bg = findBin(BG, BIN_BG, bg);
365 pid = 1;
366 }
367 else
368 {
369 continue;
370 }
371 double cost = m_cost[i];
372 int charge = m_charge[i];
373 double t[5], q;
374 double offset, sigma;
375 double offset_q, offset_bgcost;
376 int flag[4] = {0, 0, 0, 0, };
377 cost = max(COST_LOW+EPS, min(COST_UP-EPS, cost));
378 bin_cost = findBin(COST, BIN_COST, cost);
379 if (bin_bg == -1 || bin_cost == -1) continue;
380
381 //corrections
382 for (unsigned int j = 0; j < 4; j++)
383 {
384 t[j] = m_tof_dt[j][i];
385 if (fabs(t[j] + 99.) < EPS)//no readout
386 flag[j] = 0;
387 else
388 flag[j] = 1;
389 q = m_tof_ph[j];
390 q = max(Q_LOW+EPS, min(Q_UP-EPS, q));
391 if (charge == 1)
392 {
393 offset_q = h_tof_p_q_offset[pid][idx][j]->Interpolate( q );
394 offset_bgcost = h_tof_p_bgcost_offset[pid][idx][j]->Interpolate( bg, cost );
395 t[j] = t[j] - offset_q - offset_bgcost;
396 }
397 else
398 {
399 offset_q = h_tof_m_q_offset[pid][idx][j]->Interpolate( q );
400 offset_bgcost = h_tof_m_bgcost_offset[pid][idx][j]->Interpolate( bg, cost );
401 t[j] = t[j] - offset_q - offset_bgcost;
402 }
403 }
404 bin_wgt = flag[0]*8 + flag[1]*4 + flag[2]*2 + flag[3] - 1;
405 if (bin_wgt == -1) continue;
406 t[4] = 0;
407 for (unsigned int j = 0; j < 4; j++)
408 {
409 if (charge == 1)
410 t[4] += t[j] * h_tof_p_wgt[pid][idx][bin_wgt][j]->GetBinContent( bin_bg+1, bin_cost+1 );
411 else
412 t[4] += t[j] * h_tof_m_wgt[pid][idx][bin_wgt][j]->GetBinContent( bin_bg+1, bin_cost+1 );
413 }
414 if (charge == 1)
415 {
416 t[4] /= h_tof_p_wgt[pid][idx][bin_wgt][4]->GetBinContent( bin_bg+1, bin_cost+1 );
417 offset = h_tof_p_final_offset[pid][idx][bin_wgt]->Interpolate( bg, cost );
418 sigma = h_tof_p_final_sigma[pid][idx][bin_wgt]-> Interpolate( bg, cost );
419 m_tof_chi[i] = (t[4] - offset) / sigma;
420 }
421 else
422 {
423 t[4] /= h_tof_m_wgt[pid][idx][bin_wgt][4]->GetBinContent( bin_bg+1, bin_cost+1 );
424 offset = h_tof_m_final_offset[pid][idx][bin_wgt]->Interpolate( bg, cost );
425 sigma = h_tof_m_final_sigma[pid][idx][bin_wgt]-> Interpolate( bg, cost );
426 m_tof_chi[i] = (t[4] - offset) / sigma;
427 }
428 }
429 }
430}
431
432void SimplePIDSvc::tofEndcapCorrection()
433{
434 const double EPS = 1e-4;
435 const double BG_LOW = 0.30;
436 const double BG_UP = 7.40;
437 const double Q_LOW = 0.;
438 const double Q_UP = 6000.;
439 const double COST_EAST_LOW = 0.720;
440 const double COST_EAST_UP = 0.930;
441 const double COST_WEST_LOW = -0.930;
442 const double COST_WEST_UP = -0.720;
443 const double P_LOW = 0.2;
444 const double P_UP = 1.3;
445
446 int idx = getRunIdx(m_run);
447
448 if (idx != -1)
449 {
450 for (unsigned int i = 0; i < 4; i++)
451 {// only correct e, pi, K
452 int pid;
453 double bg;
454 if (i == 0)
455 {
456 bg = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
457 pid = 0;
458 }
459 else if (i == 2 || i == 3)
460 {
461 bg = max(BG_LOW+EPS, min(BG_UP-EPS, m_betagamma[i]));
462 pid = 1;
463 }
464 else
465 {
466 continue;
467 }
468
469 int flag; //0:east, 1:west
470 double cost = m_cost[i];
471 int charge = m_charge[i];
472 double t = m_tof_dt[7][i];
473 double q = m_tof_ph[7];
474 double off_q, off_bg, off_cost;
475 double sg_q, sg_bg, sg_cost;
476 if (cost > 0)
477 {
478 flag = 0;
479 cost = max(COST_EAST_LOW+EPS, min(COST_EAST_UP-EPS, cost));
480 }
481 else
482 {
483 flag = 1;
484 cost = max(COST_WEST_LOW+EPS, min(COST_WEST_UP-EPS, cost));
485 }
486 q = max(Q_LOW+EPS, min(Q_UP-EPS, q));
487
488 //corrections
489 if (charge == 1)
490 {
491 off_q = h_tofec_p_q_offset[pid][idx][flag] ->Interpolate( q );
492 sg_q = h_tofec_p_q_sigma[pid][idx][flag] ->Interpolate( q );
493 off_bg = h_tofec_p_bg_offset[pid][idx][flag] ->Interpolate( bg );
494 sg_bg = h_tofec_p_bg_sigma[pid][idx][flag] ->Interpolate( bg );
495 off_cost = h_tofec_p_cost_offset[pid][idx][flag]->Interpolate( cost );
496 sg_cost = h_tofec_p_cost_sigma[pid][idx][flag] ->Interpolate( cost );
497 m_tof_chi[i] = (((t - off_q) / sg_q - off_bg) / sg_bg - off_cost) / sg_cost;
498 }
499 else
500 {
501 off_q = h_tofec_m_q_offset[pid][idx][flag] ->Interpolate( q );
502 sg_q = h_tofec_m_q_sigma[pid][idx][flag] ->Interpolate( q );
503 off_bg = h_tofec_m_bg_offset[pid][idx][flag] ->Interpolate( bg );
504 sg_bg = h_tofec_m_bg_sigma[pid][idx][flag] ->Interpolate( bg );
505 off_cost = h_tofec_m_cost_offset[pid][idx][flag]->Interpolate( cost );
506 sg_cost = h_tofec_m_cost_sigma[pid][idx][flag] ->Interpolate( cost );
507 m_tof_chi[i] = (((t - off_q) / sg_q - off_bg) / sg_bg - off_cost) / sg_cost;
508 }
509 }
510 }
511}
512
513void SimplePIDSvc::loadDedxInfo(EvtRecTrack *track)
514{
515 if (track->isMdcDedxValid())
516 {
517 RecMdcDedx* dedx_trk = track->mdcDedx();
518 for (unsigned int i = 0; i < 5; i++)
519 m_dedx_chi[i] = dedx_trk->chi(i);
520 }
521 else
522 {
523 for (unsigned int i = 0; i < 5; i++)
524 m_dedx_chi[i] = -99;
525 }
526}
527
528void SimplePIDSvc::dedxCorrection()
529{
530 int idx = getRunIdx(m_run);
531 const double EPS = 1e-4;
532 const double BG_LOW = 0.20;
533 const double BG_UP = 7.40;
534 const double COST_LOW = -0.93;
535 const double COST_UP = 0.93;
536 const double P_LOW = 0.2;
537 const double P_UP = 1.3;
538 if (idx != -1)
539 {
540 double offset, sigma;
541 for (unsigned int i = 0; i < 4; i++)
542 {// only correct e, pi, K
543 double bg;
544 int pid;
545 if (i == 0)
546 {
547 bg = max(P_LOW+EPS, min(P_UP-EPS, m_p[i]));
548 pid = 0;
549 }
550 else if (i == 2 || i == 3)
551 {
552 bg = max(BG_LOW+EPS, min(BG_UP-EPS, m_betagamma[i]));
553 pid = 1;
554 }
555 else
556 {
557 continue;
558 }
559 double cost = m_cost[i];
560 double charge = m_charge[i];
561 cost = max(COST_LOW+EPS, min(COST_UP-EPS, cost));
562 if (charge == 1)
563 {
564 offset = h_dedx_p_offset[pid][idx]->Interpolate( bg, cost );
565 sigma = h_dedx_p_sigma[pid][idx] ->Interpolate( bg, cost );
566 m_dedx_chi[i] = (m_dedx_chi[i] - offset) / sigma;
567 }
568 else
569 {
570 offset = h_dedx_m_offset[pid][idx]->Interpolate( bg, cost );
571 sigma = h_dedx_m_sigma[pid][idx] ->Interpolate( bg, cost );
572 m_dedx_chi[i] = (m_dedx_chi[i] - offset) / sigma;
573 }
574 }
575 }
576}
577
578void SimplePIDSvc::loadHistogram()
579{
580 string path = getenv("SIMPLEPIDSVCROOT");
581 vector<string> filename;
582 for (unsigned int idx = 0; idx < 2; idx++)
583 {
584 const char *dir;
585 if (idx == 0)
586 dir = "electron";
587 else if (idx == 1)
588 dir = "kpi";
589 else
590 {
591 cout << "Boundary Error! " << endl;
592 exit(1);
593 }
594
595 //dedx
596 filename.clear();
597 filename.push_back( path + Form("/share/%s/dedx/dedx_d10.root", dir) );
598 filename.push_back( path + Form("/share/%s/dedx/dedx_d11.root", dir) );
599 filename.push_back( path + Form("/share/%s/dedx/dedx_m10.root", dir) );
600 filename.push_back( path + Form("/share/%s/dedx/dedx_m11.root", dir) );
601 for (unsigned int i = 0; i < filename.size(); i++)
602 {
603 f_dedx[idx][i] = new TFile(filename[i].c_str(), "READ");
604 const char *name;
605 if (i == 0)
606 name = "d10";
607 else if (i == 1)
608 name = "d11";
609 else if (i == 2)
610 name = "m10";
611 else if (i == 3)
612 name = "m11";
613 else
614 {
615 cout << "Boundary Error! " << endl;
616 exit(1);
617 }
618 h_dedx_p_offset[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_p_offset_%s", name) );
619 h_dedx_p_sigma[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_p_sigma_%s" , name) );
620 h_dedx_m_offset[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_m_offset_%s", name) );
621 h_dedx_m_sigma[idx][i] = (TH2D*)f_dedx[idx][i]->Get( Form("h_dedx_m_sigma_%s" , name) );
622 }
623 //tof_barrel q
624 filename.clear();
625 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_d10.root", dir) );
626 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_d11.root", dir) );
627 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_m10.root", dir) );
628 filename.push_back( path + Form("/share/%s/tof_barrel/tof_q_m11.root", dir) );
629 for (unsigned int i = 0; i < filename.size(); i++)
630 {
631 f_tof_q[idx][i] = new TFile(filename[i].c_str(), "READ");
632 const char *name;
633 if (i == 0)
634 name = "d10";
635 else if (i == 1)
636 name = "d11";
637 else if (i == 2)
638 name = "m10";
639 else if (i == 3)
640 name = "m11";
641 else
642 {
643 cout << "Boundary Error! " << endl;
644 exit(1);
645 }
646 for (unsigned int j = 0; j < 4; j++)
647 {
648 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) );
649 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) );
650 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) );
651 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) );
652 }
653 }
654 //tof_barrel bg&cost
655 filename.clear();
656 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_d10.root", dir) );
657 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_d11.root", dir) );
658 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_m10.root", dir) );
659 filename.push_back( path + Form("/share/%s/tof_barrel/tof_bg_cost_m11.root", dir) );
660 for (unsigned int i = 0; i < filename.size(); i++)
661 {
662 f_tof_bgcost[idx][i] = new TFile(filename[i].c_str(), "READ");
663 const char *name;
664 if (i == 0)
665 name = "d10";
666 else if (i == 1)
667 name = "d11";
668 else if (i == 2)
669 name = "m10";
670 else if (i == 3)
671 name = "m11";
672 else
673 {
674 cout << "Boundary Error! " << endl;
675 exit(1);
676 }
677 for (unsigned int j = 0; j < 4; j++)
678 {
679 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) );
680 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) );
681 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) );
682 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) );
683 }
684 }
685 //tof_barrel wgt
686 filename.clear();
687 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_d10.root", dir) );
688 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_d11.root", dir) );
689 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_m10.root", dir) );
690 filename.push_back( path + Form("/share/%s/tof_barrel/tof_wgt_m11.root", dir) );
691 for (unsigned int i = 0; i < filename.size(); i++)
692 {
693 f_tof_wgt[idx][i] = new TFile(filename[i].c_str(), "READ");
694 const char *name;
695 if (i == 0)
696 name = "d10";
697 else if (i == 1)
698 name = "d11";
699 else if (i == 2)
700 name = "m10";
701 else if (i == 3)
702 name = "m11";
703 else
704 {
705 cout << "Boundary Error! " << endl;
706 exit(1);
707 }
708 for (unsigned int j = 0; j < 15; j++)
709 {
710 for (unsigned int k = 0; k < 5; k++)
711 {
712 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) );
713 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) );
714 }
715 }
716 }
717 //tof_barrel corr
718 filename.clear();
719 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_d10.root", dir) );
720 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_d11.root", dir) );
721 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_m10.root", dir) );
722 filename.push_back( path + Form("/share/%s/tof_barrel/tof_final_m11.root", dir) );
723 for (unsigned int i = 0; i < filename.size(); i++)
724 {
725 f_tof_final[idx][i] = new TFile(filename[i].c_str(), "READ");
726 const char *name;
727 if (i == 0)
728 name = "d10";
729 else if (i == 1)
730 name = "d11";
731 else if (i == 2)
732 name = "m10";
733 else if (i == 3)
734 name = "m11";
735 else
736 {
737 cout << "Boundary Error! " << endl;
738 exit(1);
739 }
740 for (unsigned int j = 0; j < 15; j++)
741 {
742 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) );
743 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) );
744 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) );
745 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) );
746 }
747 }
748 //tof_endcap q
749 filename.clear();
750 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_d10.root", dir) );
751 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_d11.root", dir) );
752 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_m10.root", dir) );
753 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_q_m11.root", dir) );
754 for (unsigned int i = 0; i < filename.size(); i++)
755 {
756 f_tofec_q[idx][i] = new TFile(filename[i].c_str(), "READ");
757 const char *name;
758 if (i == 0)
759 name = "d10";
760 else if (i == 1)
761 name = "d11";
762 else if (i == 2)
763 name = "m10";
764 else if (i == 3)
765 name = "m11";
766 else
767 {
768 cout << "Boundary Error! " << endl;
769 exit(1);
770 }
771 for (unsigned int j = 0; j < 2; j++)
772 {
773 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) );
774 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) );
775 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) );
776 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) );
777 }
778 }
779 //tof_endcap bg
780 filename.clear();
781 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_d10.root", dir) );
782 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_d11.root", dir) );
783 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_m10.root", dir) );
784 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_bg_m11.root", dir) );
785 for (unsigned int i = 0; i < filename.size(); i++)
786 {
787 f_tofec_bg[idx][i] = new TFile(filename[i].c_str(), "READ");
788 const char *name;
789 if (i == 0)
790 name = "d10";
791 else if (i == 1)
792 name = "d11";
793 else if (i == 2)
794 name = "m10";
795 else if (i == 3)
796 name = "m11";
797 else
798 {
799 cout << "Boundary Error! " << endl;
800 exit(1);
801 }
802 for (unsigned int j = 0; j < 2; j++)
803 {
804 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) );
805 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) );
806 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) );
807 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) );
808 }
809 }
810 //tof_endcap cost
811 filename.clear();
812 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_d10.root", dir) );
813 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_d11.root", dir) );
814 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_m10.root", dir) );
815 filename.push_back( path + Form("/share/%s/tof_endcap/tofec_cost_m11.root", dir) );
816 for (unsigned int i = 0; i < filename.size(); i++)
817 {
818 f_tofec_cost[idx][i] = new TFile(filename[i].c_str(), "READ");
819 const char *name;
820 if (i == 0)
821 name = "d10";
822 else if (i == 1)
823 name = "d11";
824 else if (i == 2)
825 name = "m10";
826 else if (i == 3)
827 name = "m11";
828 else
829 {
830 cout << "Boundary Error! " << endl;
831 exit(1);
832 }
833 for (unsigned int j = 0; j < 2; j++)
834 {
835 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) );
836 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) );
837 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) );
838 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) );
839 }
840 }
841 }
842 for (unsigned int idx = 0; idx < 3; idx++)
843 {
844 const char *dir;
845 if (idx == 0)
846 dir = "electron/emc";
847 else if (idx == 1)
848 dir = "kpi/emc_pion";
849 else if (idx == 2)
850 dir = "kpi/emc_kaon";
851 else
852 {
853 cout << "Boundary Error! " << endl;
854 exit(1);
855 }
856 //emc
857 filename.clear();
858 filename.push_back( path + Form("/share/%s/emc_d10.root", dir) );
859 filename.push_back( path + Form("/share/%s/emc_d11.root", dir) );
860 filename.push_back( path + Form("/share/%s/emc_m10.root", dir) );
861 filename.push_back( path + Form("/share/%s/emc_m11.root", dir) );
862 for (unsigned int i = 0; i < filename.size(); i++)
863 {
864 f_emc[idx][i] = new TFile(filename[i].c_str(), "READ");
865 const char *name;
866 if (i == 0)
867 name = "d10";
868 else if (i == 1)
869 name = "d11";
870 else if (i == 2)
871 name = "m10";
872 else if (i == 3)
873 name = "m11";
874 else
875 {
876 cout << "Boundary Error! " << endl;
877 exit(1);
878 }
879 for (unsigned int j = 0; j < 15; j++)
880 {
881 for (unsigned int k = 0; k < 25; k++)
882 {
883 h_emc_ep[idx][i][j][k] = (TH1D*)f_emc[idx][i]->Get( Form("h_ep_%s_%d_%d", name, j, k) );
884 h_emc_e35[idx][i][j][k] = (TH1D*)f_emc[idx][i]->Get( Form("h_e35_%s_%d_%d", name, j, k) );
885 }
886 }
887 }
888 }
889 cout << "Successfully Return from Loading Initializations by package SimplePIDSvc ... " << endl;
890}
891
892int SimplePIDSvc::findBin(double *a, int length, double value)
893{
894 for (int i = 0; i < length; i++)
895 {
896 if (value > a[i] && value <= a[i+1])
897 {
898 return i;
899 }
900 }
901 if (value < a[0])
902 {
903 return 0;
904 }
905 else
906 {
907 return length;
908 }
909}
910
911double SimplePIDSvc::getChi2(int i)
912{
913 return pow(m_dedx_chi[i], 2) + pow(m_tof_chi[i], 2);
914}
915
916void SimplePIDSvc::loadEMCInfo(EvtRecTrack *track)
917{
918 //Initialization
919 for (unsigned int i = 0; i < 5; i++)
920 {
921 m_emc_eop[i] = -99.;
922 m_emc_likelihood[i] = -99.;
923 }
924 m_emc_e = -99.;
925 m_emc_e13 = -99.;
926 m_emc_e35 = -99.;
927 m_emc_sec = -99.;
928 m_emc_lat = -99.;
929 m_lh_electron = -99.;
930
931 if (!track->isEmcShowerValid()) return;
932
933 RecEmcShower* emc_trk = track->emcShower();
934 m_emc_e = emc_trk->energy();
935 for (unsigned int i = 0; i < 5; i++)
936 {
937 m_emc_eop[i] = m_emc_e / fabs(m_p[i]);
938 }
939 double eseed = emc_trk->eSeed();
940 double e3 = emc_trk->e3x3();
941 double e5 = emc_trk->e5x5();
942 m_emc_sec = emc_trk->secondMoment() / 1000.;
943 m_emc_lat = emc_trk->latMoment();
944 if (e3 != 0)
945 {
946 m_emc_e13 = eseed / e3;
947 }
948 if (e5 != 0)
949 {
950 m_emc_e35 = e3 / e5;
951 }
952}
953
954bool SimplePIDSvc::calEMCLikelihood()
955{
956 if (m_emc_eop[0] < 0)
957 return false;
958
959 int idx = getRunIdx(m_run);
960 const Int_t BIN_P = 15;
961 const Int_t BIN_COST = 25;
962 const Int_t BIN_PID = 3;
963 const double EPS = 1e-4;
964 //electron
965 double P[BIN_PID][BIN_P + 1] = {
966 {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},
967 {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},
968 {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}, };
969 double COST[BIN_PID][BIN_COST + 1] = {
970 {-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},
971 {-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},
972 {-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}, };
973
974 double vp, vcost;
975 int pid;
976 int bin_p, bin_cost;
977 for (unsigned int i = 0; i < 4; i++)
978 {
979 //only e, pi ,K
980 if (i == 0)
981 pid = 0;
982 else if (i == 2)
983 pid = 1;
984 else if (i == 3)
985 pid = 2;
986 else
987 continue;
988
989 vp = max(P[pid][0]+EPS, min(P[pid][BIN_P]-EPS, m_p[i]));
990 vcost = max(COST[pid][0]+EPS, min(COST[pid][BIN_COST]-EPS, m_cost[i]));
991 bin_p = findBin(P[pid], BIN_P, vp);
992 bin_cost = findBin(COST[pid], BIN_COST, vcost);
993
994 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);
995 }
996 double a = m_prob[0] > 0 ? m_prob[0] : 0;
997 double b = m_prob[2] > 0 ? m_prob[2] : 0;
998 double c = m_prob[3] > 0 ? m_prob[3] : 0;
999 double sum = a * m_emc_likelihood[0] + b * m_emc_likelihood[2] + c * m_emc_likelihood[3];
1000
1001 if (sum > 0 && m_prob[0] > 0)
1002 {
1003 m_lh_electron = m_prob[0] * m_emc_likelihood[0] / sum;
1004 return true;
1005 }
1006 else
1007 {
1008 return false;
1009 }
1010}
1011
1012bool SimplePIDSvc::ispion()
1013{
1014 if (m_prob[2] > 0.00 && m_prob[2] > m_prob[3])
1015 return true;
1016 else
1017 return false;
1018}
1019
1020bool SimplePIDSvc::iskaon()
1021{
1022 if (m_prob[3] > 0.00 && m_prob[3] > m_prob[2])
1023 return true;
1024 else
1025 return false;
1026}
1027
1028//bool SimplePIDSvc::iselectron_(bool emc)
1029//{
1030// if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
1031// {
1032// if (!emc)
1033// return true;
1034// else if (fabs(m_cost[0]) < 0.7 && m_emc_eop[0] > 0 && m_emc_eop[0] < 0.8)
1035// return false;
1036// 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)
1037// return false;
1038// else if (fabs(m_cost[0]) > 0.85 && m_emc_eop[0] > 0 && m_emc_eop[0] < 0.6)
1039// return false;
1040// else
1041// return true;
1042// }
1043// else
1044// return false;
1045//}
1046
1047bool SimplePIDSvc::iselectron(bool emc)
1048{
1049 if (!emc)
1050 {
1051 if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
1052 return true;
1053 else
1054 return false;
1055 }
1056 else
1057 {
1058 if (calEMCLikelihood())
1059 {
1060 if (m_lh_electron > m_eid_ratio)
1061 return true;
1062 else
1063 return false;
1064 }
1065 else
1066 {
1067 if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
1068 return true;
1069 else
1070 return false;
1071 }
1072 }
1073}
P
double P(RecMdcKalTrack *trk)
Definition: CalibEventSelect.cxx:87
mass
double mass
Definition: CosmicGenerator.cxx:138
sigma
TTree * sigma
Definition: CurConst_dEdXCal.cxx:6
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
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
SimplePIDSvc::SimplePIDSvc
SimplePIDSvc(const std::string &name, ISvcLocator *svcLoc)
Definition: SimplePIDSvc.cxx:20
SimplePIDSvc::~SimplePIDSvc
virtual ~SimplePIDSvc()
Definition: SimplePIDSvc.cxx:29
SimplePIDSvc::getChi2
double getChi2(int i)
Definition: SimplePIDSvc.cxx:911
SimplePIDSvc::preparePID
void preparePID(EvtRecTrack *track)
Definition: SimplePIDSvc.cxx:91
SimplePIDSvc::initialize
virtual StatusCode initialize()
Definition: SimplePIDSvc.cxx:31
SimplePIDSvc::iselectron
bool iselectron(bool emc=true)
Definition: SimplePIDSvc.cxx:1047
SimplePIDSvc::finalize
virtual StatusCode finalize()
Definition: SimplePIDSvc.cxx:45
SimplePIDSvc::queryInterface
virtual StatusCode queryInterface(const InterfaceID &riid, void **ppvIF)
Definition: SimplePIDSvc.cxx:77
TofHitStatus::setStatus
void setStatus(unsigned int status)
Definition: TofHitStatus.cxx:134
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