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"
9#include "DstEvent/TofHitStatus.h"
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);
33 MsgStream log(messageService(), name());
34 log << MSG::INFO <<
"in SimplePIDSvc initialize()" << endreq;
36 StatusCode sc = Service::initialize();
38 sc = serviceLocator()->service(
"EventDataSvc", eventSvc_,
true);
47 MsgStream log(messageService(), name());
48 log << MSG::INFO <<
"in SimplePIDSvc finalize()" << endreq;
50 StatusCode sc = Service::finalize();
52 for (
unsigned int i = 0; i < 2; i++)
54 for (
unsigned int j = 0; j < 4; j++)
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();
66 for (
unsigned int i = 0; i < 3; i++)
68 for (
unsigned int j = 0; j < 4; j++)
85 return Service::queryInterface(riid, ppvIF);
88 return StatusCode::SUCCESS;
94 SmartDataPtr<Event::EventHeader> eventHeaderpid(eventSvc_,
"/Event/EventHeader");
95 m_run = eventHeaderpid->runNumber();
114 for(
unsigned int pid = 0; pid < 5; 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());
125 for(
unsigned int i = 0; i < 5; i++)
128 m_betagamma[i] = -99;
145 if (m_tof_barrel == 1)
147 tofBarrelCorrection();
149 else if (m_tof_barrel == 0)
150 tofEndcapCorrection();
158void SimplePIDSvc::calprob()
160 bool usededx =
false;
163 for (
unsigned int i = 0; i < 5 ;i++)
165 if (!usededx && fabs(m_dedx_chi[i]) < m_dedx_chi_cut)
167 if (!usetof && fabs(m_tof_chi[i]) < m_tof_chi_cut)
170 m_dedx_only[i] =
false;
174 for(
unsigned int i = 0; i < 5; i++)
179 for(
unsigned int i = 0; i < 5; i++)
183 for (
unsigned int i = 0; i < 5; i++)
189 if (usededx && usetof)
191 chi2 = pow(m_dedx_chi[i], 2) + pow(m_tof_chi[i], 2);
194 else if (usededx && !usetof)
196 chi2 = pow(m_dedx_chi[i], 2);
198 m_dedx_only[i] =
true;
200 else if (!usededx && usetof)
202 chi2 = pow(m_tof_chi[i],2);
205 if (ndf > 0 && chi2 > 0)
206 m_prob[i] = TMath::Prob(chi2, ndf);
212int SimplePIDSvc::getRunIdx(
int run_no)
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)
229 else if (run_no >= RUN_BEGIN_DATA_11 && run_no <= RUN_END_DATA_11)
231 else if (run_no >= RUN_BEGIN_MC_10 && run_no <= RUN_END_MC_10)
233 else if (run_no >= RUN_BEGIN_MC_11 && run_no <= RUN_END_MC_11)
242 for (
unsigned int i = 0; i < 8; i++)
244 for (
unsigned int j = 0; j < 5; j++)
245 m_tof_dt[i][j] = -99.;
248 for (
unsigned int i = 0; i < 2; i++)
250 m_tof_zr[i] = -9999.;
251 m_tof_counter[i] = -1;
253 for (
unsigned int i = 0; i < 5; i++)
261 SmartRefVector<RecTofTrack> tofTrk = track->
tofTrack();
262 SmartRefVector<RecTofTrack>::iterator it;
270 for (it = tofTrk.begin(); it != tofTrk.end(); it++)
272 unsigned int st = (*it)->status();
274 if (hitst->
is_raw())
continue;
279 int layer = hitst->
layer();
280 double tof = (*it)->tof();
281 double ph = (*it)->ph();
282 m_tof_counter[layer-1] = (*it)->tofID();
293 m_tof_zr[0] = zrhit[0];
294 m_tof_zr[1] = zrhit[1];
300 idx = ((st & 0xC0) >> 5) + (((st ^ 0x20) & 0x20) >> 5) - 2;
312 if (idx == -1)
continue;
314 for (
unsigned int i = 0; i < 5; i++)
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;
329void SimplePIDSvc::tofBarrelCorrection()
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);
350 for (
unsigned int i = 0; i < 4; i++)
353 int bin_bg, bin_cost, bin_wgt;
357 bg = max(P_LOW+
EPS, min(P_UP-
EPS, m_p[i]));
358 bin_bg = findBin(
P, BIN_P,
bg);
361 else if (i == 2 || i == 3)
363 bg = max(BG_LOW+
EPS, min(BG_UP-
EPS, m_betagamma[i]));
364 bin_bg = findBin(BG, BIN_BG,
bg);
371 double cost = m_cost[i];
372 int charge = m_charge[i];
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;
382 for (
unsigned int j = 0; j < 4; j++)
384 t[j] = m_tof_dt[j][i];
385 if (fabs(
t[j] + 99.) <
EPS)
390 q = max(Q_LOW+
EPS, min(Q_UP-
EPS,
q));
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;
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;
404 bin_wgt = flag[0]*8 + flag[1]*4 + flag[2]*2 + flag[3] - 1;
405 if (bin_wgt == -1)
continue;
407 for (
unsigned int j = 0; j < 4; j++)
410 t[4] +=
t[j] * h_tof_p_wgt[pid][idx][bin_wgt][j]->GetBinContent( bin_bg+1, bin_cost+1 );
412 t[4] +=
t[j] * h_tof_m_wgt[pid][idx][bin_wgt][j]->GetBinContent( bin_bg+1, bin_cost+1 );
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;
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;
432void SimplePIDSvc::tofEndcapCorrection()
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;
446 int idx = getRunIdx(m_run);
450 for (
unsigned int i = 0; i < 4; i++)
456 bg = max(P_LOW+
EPS, min(P_UP-
EPS, m_p[i]));
459 else if (i == 2 || i == 3)
461 bg = max(BG_LOW+
EPS, min(BG_UP-
EPS, m_betagamma[i]));
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;
479 cost = max(COST_EAST_LOW+
EPS, min(COST_EAST_UP-
EPS, cost));
484 cost = max(COST_WEST_LOW+
EPS, min(COST_WEST_UP-
EPS, cost));
486 q = max(Q_LOW+
EPS, min(Q_UP-
EPS,
q));
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;
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;
518 for (
unsigned int i = 0; i < 5; i++)
519 m_dedx_chi[i] = dedx_trk->
chi(i);
523 for (
unsigned int i = 0; i < 5; i++)
528void SimplePIDSvc::dedxCorrection()
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;
540 double offset, sigma;
541 for (
unsigned int i = 0; i < 4; i++)
547 bg = max(P_LOW+
EPS, min(P_UP-
EPS, m_p[i]));
550 else if (i == 2 || i == 3)
552 bg = max(BG_LOW+
EPS, min(BG_UP-
EPS, m_betagamma[i]));
559 double cost = m_cost[i];
560 double charge = m_charge[i];
561 cost = max(COST_LOW+
EPS, min(COST_UP-
EPS, cost));
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;
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;
578void SimplePIDSvc::loadHistogram()
580 string path = getenv(
"SIMPLEPIDSVCROOT");
581 vector<string> filename;
582 for (
unsigned int idx = 0; idx < 2; idx++)
591 cout <<
"Boundary Error! " << endl;
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++)
603 f_dedx[idx][i] =
new TFile(filename[i].
c_str(),
"READ");
615 cout <<
"Boundary Error! " << endl;
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) );
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++)
631 f_tof_q[idx][i] =
new TFile(filename[i].
c_str(),
"READ");
643 cout <<
"Boundary Error! " << endl;
646 for (
unsigned int j = 0; j < 4; j++)
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) );
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++)
662 f_tof_bgcost[idx][i] =
new TFile(filename[i].
c_str(),
"READ");
674 cout <<
"Boundary Error! " << endl;
677 for (
unsigned int j = 0; j < 4; j++)
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) );
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++)
693 f_tof_wgt[idx][i] =
new TFile(filename[i].
c_str(),
"READ");
705 cout <<
"Boundary Error! " << endl;
708 for (
unsigned int j = 0; j < 15; j++)
710 for (
unsigned int k = 0; k < 5; k++)
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) );
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++)
725 f_tof_final[idx][i] =
new TFile(filename[i].
c_str(),
"READ");
737 cout <<
"Boundary Error! " << endl;
740 for (
unsigned int j = 0; j < 15; j++)
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) );
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++)
756 f_tofec_q[idx][i] =
new TFile(filename[i].
c_str(),
"READ");
768 cout <<
"Boundary Error! " << endl;
771 for (
unsigned int j = 0; j < 2; j++)
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) );
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++)
787 f_tofec_bg[idx][i] =
new TFile(filename[i].
c_str(),
"READ");
799 cout <<
"Boundary Error! " << endl;
802 for (
unsigned int j = 0; j < 2; j++)
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) );
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++)
818 f_tofec_cost[idx][i] =
new TFile(filename[i].
c_str(),
"READ");
830 cout <<
"Boundary Error! " << endl;
833 for (
unsigned int j = 0; j < 2; j++)
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) );
842 for (
unsigned int idx = 0; idx < 3; idx++)
846 dir =
"electron/emc";
848 dir =
"kpi/emc_pion";
850 dir =
"kpi/emc_kaon";
853 cout <<
"Boundary Error! " << endl;
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++)
864 f_emc[idx][i] =
new TFile(filename[i].
c_str(),
"READ");
876 cout <<
"Boundary Error! " << endl;
879 for (
unsigned int j = 0; j < 15; j++)
881 for (
unsigned int k = 0; k < 25; k++)
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) );
889 cout <<
"Successfully Return from Loading Initializations by package SimplePIDSvc ... " << endl;
892int SimplePIDSvc::findBin(
double *a,
int length,
double value)
894 for (
int i = 0; i < length; i++)
896 if (value > a[i] && value <= a[i+1])
913 return pow(m_dedx_chi[i], 2) + pow(m_tof_chi[i], 2);
919 for (
unsigned int i = 0; i < 5; i++)
922 m_emc_likelihood[i] = -99.;
929 m_lh_electron = -99.;
934 m_emc_e = emc_trk->
energy();
935 for (
unsigned int i = 0; i < 5; i++)
937 m_emc_eop[i] = m_emc_e / fabs(m_p[i]);
939 double eseed = emc_trk->
eSeed();
940 double e3 = emc_trk->
e3x3();
941 double e5 = emc_trk->
e5x5();
946 m_emc_e13 = eseed / e3;
954bool SimplePIDSvc::calEMCLikelihood()
956 if (m_emc_eop[0] < 0)
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;
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}, };
977 for (
unsigned int i = 0; i < 4; i++)
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);
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);
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];
1001 if (sum > 0 && m_prob[0] > 0)
1003 m_lh_electron = m_prob[0] * m_emc_likelihood[0] / sum;
1014 if (m_prob[2] > 0.00 && m_prob[2] > m_prob[3])
1022 if (m_prob[3] > 0.00 && m_prob[3] > m_prob[2])
1051 if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
1058 if (calEMCLikelihood())
1060 if (m_lh_electron > m_eid_ratio)
1067 if (m_prob[0] > 0 && m_prob[0] > m_prob[2] && m_prob[0] > m_prob[3])
double P(RecMdcKalTrack *trk)
double cos(const BesAngle a)
****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
double secondMoment() const
const Hep3Vector tof1Position() const
const Hep3Vector tof2Position() const
const double theta() const
static void setPidType(PidType pidType)
bool isMdcKalTrackValid()
SmartRefVector< RecTofTrack > tofTrack()
RecEmcShower * emcShower()
RecMdcKalTrack * mdcKalTrack()
static const InterfaceID & interfaceID()
SimplePIDSvc(const std::string &name, ISvcLocator *svcLoc)
void preparePID(EvtRecTrack *track)
virtual StatusCode initialize()
bool iselectron(bool emc=true)
virtual StatusCode finalize()
virtual StatusCode queryInterface(const InterfaceID &riid, void **ppvIF)
unsigned int layer() const
void setStatus(unsigned int status)