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