BOSS 6.6.4.p01
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"
7//#include "EvtRecEvent/EvtRecDTag.h"
8//#include "VertexFit/Helix.h"
10
12#include "TMath.h"
13#include "TFile.h"
14#include "TH1F.h"
15#include <fstream>
16#include "TofCorrection.C"
17
18SimplePIDSvc::SimplePIDSvc(const std::string& name, ISvcLocator* svcLoc)
19 : Service(name, svcLoc)
20{
21 declareProperty("mindedxchi", m_dedxminchi=4);
22 declareProperty("mintofchi", m_tofminchi=2);
23 declareProperty("tofcorrection", m_tofcorrec=true);
24}
25
27{
28}
29
31{
32 MsgStream log(messageService(), name());
33 log << MSG::INFO << "@initialize()" << endreq;
34
35 StatusCode sc = Service::initialize();
36
37 sc = serviceLocator()->service("EventDataSvc", eventSvc_, true);
38 //sc = serviceLocator()->service("VertexDbSvc", vtxSvc_, true);
39
40 //m_dedxminchi=4;
41 //m_tofminchi=2;
42
43 //get eop files
44 string filename;
45 filename = string(getenv("SIMPLEPIDSVCROOT"))+"/share/eop.root";
46 TFile* file=new TFile(filename.c_str(),"read");
47 h_ebarlp=(TH1F*)file->Get("barlp");
48 h_ebarhp=(TH1F*)file->Get("barhp");
49 h_eendlp=(TH1F*)file->Get("endlp");
50 h_eendhp=(TH1F*)file->Get("endhp");
51 h_kbar=(TH1F*)file->Get("barkaon");
52 h_kend=(TH1F*)file->Get("endkaon");
53 h_pibar=(TH1F*)file->Get("barpion");
54 h_piend=(TH1F*)file->Get("endpion");
55
56 return sc;
57}
58
60{
61 MsgStream log(messageService(), name());
62 log << MSG::INFO << "@initialize()" << endreq;
63
64 StatusCode sc = Service::finalize();
65
66 delete h_ebarlp;
67 delete h_ebarhp;
68 delete h_eendlp;
69 delete h_eendhp;
70 delete h_kbar;
71 delete h_kend;
72 delete h_pibar;
73 delete h_piend;
74
75
76 return sc;
77}
78
79StatusCode SimplePIDSvc::queryInterface(const InterfaceID& riid, void** ppvIF)
80{
81 if ( ISimplePIDSvc::interfaceID().versionMatch(riid) ) {
82 *ppvIF = dynamic_cast<ISimplePIDSvc*>(this);
83 }
84 else {
85 return Service::queryInterface(riid, ppvIF);
86 }
87 addRef();
88 return StatusCode::SUCCESS;
89}
90
92 string filename;
93 //filename = string(getenv("SIMPLEPIDSVCROOT"))+"/share/constant/"+string(getenv("BES_RELEASE"))+"_const.txt";
94 filename = string(getenv("SIMPLEPIDSVCROOT"))+"/share/constant/const.txt";
95 //cout<<"const file is:"<<filename<<endl;
96
97 m_datatof.clear();
98 m_mctof.clear();
99 run = fabs(run);
100 std::ifstream fs(filename.c_str());
101 double tag, num;
102 int run_begin, run_end, method;
103 bool isdata=false, ismc=false;
104 while(fs.good() && !fs.eof() ){
105
106 fs >> num;
107 if(num==999){
108 isdata=false;
109 ismc=false;
110 fs >> tag;
111 fs >> run_begin;
112 fs >> run_end;
113 fs >> method;
114 if(tag==99 && run>=run_begin && run<=run_end){
115 isdata=true;
116 ismc=false;
117 }
118 if(tag==-99 && run>=run_begin && run<=run_end){
119 isdata=false;
120 ismc=true;
121 }
122 fs >> num;
123 }
124
125 if(isdata)
126 m_datatof.push_back(num);
127 if(ismc)
128 m_mctof.push_back(num);
129
130 }//end of while
131
132}
134
135
136 SmartDataPtr<Event::EventHeader> eventHeaderpid(eventSvc_,"/Event/EventHeader");
137 m_run=eventHeaderpid->runNumber();
138
139 getConst(m_run);
140 //cout<<"size of constant:"<<m_mctof.size()<<endl;
141 //for(int i=0; i<m_mctof.size();i++)
142 // cout<<m_mctof[i]<<endl;
143 if(m_datatof.size()==0 || m_mctof.size()==0)
144 m_tofcorrec=false;
145
146 //cout<<"track id :"<<track->trackId()<<endl;
147 //track quantity
148 if( track->isMdcKalTrackValid()) {
149 RecMdcKalTrack *mdckalTrk = track->mdcKalTrack();
150
151 for(int pid=0; pid<5; ++pid){
152 HepVector zhelix;
153 if(pid==0){
155 zhelix=mdckalTrk->getZHelixE();
156 }
157 else if(pid==1){
159 zhelix=mdckalTrk->getZHelixMu();
160 }
161 else if(pid==2){
163 zhelix=mdckalTrk->getZHelix();
164 }
165 else if(pid==3){
167 zhelix=mdckalTrk->getZHelixK();
168 }
169 else{
171 zhelix=mdckalTrk->getZHelixP();
172 }
173
174 double kappa=zhelix[2];
175 double theta=CLHEP::pi/2.0-atan(zhelix[4]);
176 m_p[pid]=1.0/fabs(kappa)/sin(theta);
177 m_costh[pid]=cos(theta);
178 }
179 }
180 else{
181 for(int i=0; i<5; ++i){
182 m_p[i]=-99;
183 m_costh[i]=-99;
184 }
185
186 }
187
188
189 //mass
190 m_mass[0]=0.000511;
191 m_mass[1]=0.105658;
192 m_mass[2]=0.14;
193 m_mass[3]=0.494;
194 m_mass[4]=0.94;
195
196 //tof correction factors
197 for(int i=0; i<5; ++i){
198 if(m_tofcorrec){
199 if(m_run>0){
200 m_tofscale1[i]=toflayer1scale(m_p[i]/m_mass[i],fabs(m_costh[i]),m_datatof);
201 m_tofscale2[i]=toflayer2scale(m_p[i]/m_mass[i],fabs(m_costh[i]),m_datatof);
202 }
203 else{
204 m_tofscale1[i]=mctoflayer1scale(m_p[i]/m_mass[i],fabs(m_costh[i]),m_mctof);
205 m_tofscale2[i]=mctoflayer2scale(m_p[i]/m_mass[i],fabs(m_costh[i]),m_mctof);
206 }
207 }
208 else{
209 m_tofscale1[i]=1.0;
210 m_tofscale2[i]=1.0;
211 }
212 }
213
214
215 //eop for electrons
216 m_eop=-99;
217 if( track->isEmcShowerValid() ){
218 RecEmcShower* emcTrk = track->emcShower();
219 double energy= emcTrk->energy();
220 if(fabs(m_p[0])>0)
221 m_eop=energy/fabs(m_p[0]);
222 }
223
224
225 //dE/dx PID
226
227 if( track->isMdcDedxValid()){
228 RecMdcDedx* dedxTrk = track->mdcDedx();
229 m_dedxchi[0]= dedxTrk->chi(0);
230 m_dedxchi[1]= dedxTrk->chi(1);
231 m_dedxchi[2]= dedxTrk->chi(2);
232 m_dedxchi[3]= dedxTrk->chi(3);
233 m_dedxchi[4]= dedxTrk->chi(4);
234 }
235 else{
236 m_dedxchi[0]=-99;
237 m_dedxchi[1]=-99;
238 m_dedxchi[2]=-99;
239 m_dedxchi[3]=-99;
240 m_dedxchi[4]=-99;
241 }
242
243
244 //TOF PID
245
246 m_tofchi[0]=-99;
247 m_tofchi[1]=-99;
248 m_tofchi[2]=-99;
249 m_tofchi[3]=-99;
250 m_tofchi[4]=-99;
251
252 m_tofdt[0]=-99;
253 m_tofdt[1]=-99;
254 m_tofdt[2]=-99;
255 m_tofdt[3]=-99;
256 m_tofdt[4]=-99;
257
258
259 m_tofdt1[0]=-99;
260 m_tofdt1[1]=-99;
261 m_tofdt1[2]=-99;
262 m_tofdt1[3]=-99;
263 m_tofdt1[4]=-99;
264
265 m_tofdt2[0]=-99;
266 m_tofdt2[1]=-99;
267 m_tofdt2[2]=-99;
268 m_tofdt2[3]=-99;
269 m_tofdt2[4]=-99;
270
271 int layer1size=0;
272 int layer2size=0;
273
274
275 if(track->isTofTrackValid() ){
276 SmartRefVector<RecTofTrack> tofTrkCol=track->tofTrack();
277 SmartRefVector<RecTofTrack>::iterator iter_tof=tofTrkCol.begin();
278
279 for(;iter_tof!=tofTrkCol.end();iter_tof++){
280 TofHitStatus* status =new TofHitStatus;
281 status->setStatus( (*iter_tof)->status() );
282
283 /*
284 if( status->is_cluster()){
285 double time=(*iter_tof)->tof();
286 double sigma=(*iter_tof)->sigma(0);
287 double expE_tof=(*iter_tof)->texpElectron();
288 double expMu_tof=(*iter_tof)->texpMuon();
289 double expPi_tof=(*iter_tof)->texpPion();
290 double expK_tof=(*iter_tof)->texpKaon();
291 double expP_tof=(*iter_tof)->texpProton();
292
293 m_tofchi[0]= (time- expE_tof)/sigma ;
294 m_tofchi[1]= (time- expMu_tof)/sigma ;
295 m_tofchi[2]= (time- expPi_tof)/sigma ;
296 m_tofchi[3]= (time- expK_tof)/sigma ;
297 m_tofchi[4]= (time- expP_tof)/sigma ;
298
299 m_tofdt[0]= (time- expE_tof);
300 m_tofdt[1]= (time- expMu_tof);
301 m_tofdt[2]= (time- expPi_tof);
302 m_tofdt[3]= (time- expK_tof);
303 m_tofdt[4]= (time- expP_tof);
304
305 }
306 */
307
308 //two layers
309
310 if(!(status->is_barrel())){//endcap
311
312 if( status->is_counter() && status->layer()==1 ) {
313 layer1size++;
314 for(int i=0; i<5; i++)
315 m_tofdt1[i]=(*iter_tof)->tof() - m_tofscale1[i]*(*iter_tof)->texp(i);
316 m_sigma1=(*iter_tof)->sigma(0);;
317 }
318 }
319 else if( status->is_barrel()) { //barrel
320 if( status->is_counter() ){
321 if(status->layer()==1){ //layer1
322 layer1size++;
323
324 for(int i=0; i<5; i++){
325 m_tofdt1[i]=(*iter_tof)->tof() - m_tofscale1[i]*(*iter_tof)->texp(i);
326 }
327 m_sigma1=(*iter_tof)->sigma(0);;
328 }
329 else if(status->layer()==2){//layer2
330 layer2size++;
331
332 for(int i=0; i<5; i++){
333 m_tofdt2[i]=(*iter_tof)->tof() - m_tofscale2[i]*(*iter_tof)->texp(i);
334 }
335 m_sigma2=(*iter_tof)->sigma(0);;
336 }
337
338 }
339
340 }//end of barrel
341
342
343
344
345 delete status;
346 }
347
348 }//end of tof
349
350 if(layer1size>1){
351 m_tofdt1[0]=-99;
352 m_tofdt1[1]=-99;
353 m_tofdt1[2]=-99;
354 m_tofdt1[3]=-99;
355 m_tofdt1[4]=-99;
356 }
357 if(layer2size>1){
358 m_tofdt2[0]=-99;
359 m_tofdt2[1]=-99;
360 m_tofdt2[2]=-99;
361 m_tofdt2[3]=-99;
362 m_tofdt2[4]=-99;
363 }
364
365 //combine two layers
366
367 combineTOF();
368 calprob();
369
370}
371
373
374
375 for(int i=0; i<5 ; ++i){
376
377 double dt1= m_tofdt1[i];
378 double dt2= m_tofdt2[i];
379 double sigma1= m_sigma1;
380 double sigma2= m_sigma2;
381
382 double sqrtcov=0.041;
383 if(m_tofcorrec){
384 sqrtcov=tofsqrcov(m_p[i]/m_mass[i],m_datatof);
385 if(m_run<0)
386 sqrtcov=mctofsqrcov(m_p[i]/m_mass[i],m_mctof);
387 }
388
389 double weight1=(sigma2*sigma2-sqrtcov*sqrtcov)/(sigma1*sigma1+sigma2*sigma2-2*sqrtcov*sqrtcov);
390 double weight2= (sigma1*sigma1-sqrtcov*sqrtcov)/(sigma1*sigma1+sigma2*sigma2-2*sqrtcov*sqrtcov);
391
392 double dt_comb=-99;
393 double sigma_comb=-99;
394
395 // if(sigma1==0 || sigma2==0) cout<<"zero sigma"<<endl;
396
397 if( dt1!=-99 && dt2!=-99){
398
399 dt_comb= weight1*dt1+weight2*dt2;
400 sigma_comb= sqrt((sigma1*sigma1*sigma2*sigma2 - sqrtcov*sqrtcov*sqrtcov*sqrtcov)/(sigma1*sigma1+sigma2*sigma2-2*sqrtcov*sqrtcov));
401 }
402 else if(dt1!=-99 && dt2==-99){
403 dt_comb = dt1;
404 sigma_comb=sigma1;
405 }
406 else if ( dt1==-99 && dt2!=-99){
407 dt_comb = dt2;
408 sigma_comb=sigma2;
409 }
410
411 if(sigma_comb>0){
412 double sigmascale=1.0;
413 if(m_tofcorrec){
414 sigmascale=tofsigmascale(m_p[i]/m_mass[i],m_datatof);
415 if(m_run<0)
416 sigmascale=mctofsigmascale(m_p[i]/m_mass[i],m_mctof);
417 }
418
419 m_tofchi[i]= dt_comb/sigma_comb/sigmascale;
420
421 }
422 else
423 m_tofchi[i]=-99;
424 }
425
426
427
428}
429
431
432 bool usededx=false;
433 bool usetof=false;
434
435 for(int i=0; i<5 ;i++){
436
437 if(!usededx && fabs(m_dedxchi[i])<m_dedxminchi)
438 usededx=true;
439 if(!usetof && fabs(m_tofchi[i])<m_tofminchi)
440 usetof=true;
441
442 m_dedxonly[i]=false;
443 }
444
445 if(!usededx){
446
447 for(int i=0; i<5; i++)
448 m_dedxchi[i]=-99;
449 }
450
451 if(!usetof){
452
453 for(int i=0; i<5; i++)
454 m_tofchi[i]=-99;
455 }
456
457 /*
458 if( fabs(m_p[2])<0.3)
459 usetof=false;
460 */
461
462 for(int i=0; i<5; ++i){
463 m_prob[i]=-99;
464 double chi2=0;
465 int ndf=0;
466
467
468 if( usededx && usetof ){
469 chi2= pow(m_dedxchi[i],2)+pow(m_tofchi[i],2);
470 ndf=2;
471 }
472 else if( usededx && !usetof ){
473 chi2= pow(m_dedxchi[i],2);
474 ndf=1;
475 m_dedxonly[i]=true;
476 }
477 else if( !usededx && usetof ){
478 chi2= pow(m_tofchi[i],2);
479 ndf=1;
480 }
481 /*
482 if( fabs(m_dedxchi[i])<10 && fabs(m_tofchi[i])<10 ){
483 chi2= pow(m_dedxchi[i],2)+pow(m_tofchi[i],2);
484 ndf=2;
485 }
486 else if( fabs(m_dedxchi[i])<10 && fabs(m_tofchi[i])>10 ){
487 chi2= pow(m_dedxchi[i],2);
488 ndf=1;
489 m_dedxonly[i]=true;
490 }
491 else if( fabs(m_dedxchi[i])>10 && fabs(m_tofchi[i])<10 ){
492 chi2= pow(m_tofchi[i],2);
493 ndf=1;
494 }
495
496 */
497
498 if(ndf>0 && chi2>0)
499 m_prob[i]=TMath::Prob(chi2, ndf);
500
501 }
502
503
504
505}
506
507
509
510
511 if(m_prob[0]>0 && m_prob[0]>m_prob[2] && m_prob[0]>m_prob[3]){
512
513 if(!eop)
514 return true;
515 else{
516
517 if( fabs(m_costh[0])<0.7 && m_eop>0 && m_eop<0.8)
518 return false;
519
520 if( fabs(m_costh[0])>=0.7 && fabs(m_costh[0])<0.8 && m_eop>0 && m_eop<-7.5*fabs(m_costh[0])+6.05)
521 return false;
522
523 if( fabs(m_costh[0])>0.85 && m_eop>0 && m_eop<0.6)
524 return false;
525
526 return true;
527
528 }
529
530 }
531 return false;
532
533
534 /*
535 double pe=-99, pk=-99, ppi=-99;
536
537 if( fabs(m_dedxchi[0])>5 )
538 return false;
539
540 if(m_eop>0.5 && m_eop<1.5){
541 if( fabs(m_costh[0])<0.81 ){
542 pk=h_kbar->GetBinContent(h_kbar->FindBin(m_eop));
543 ppi=h_pibar->GetBinContent(h_pibar->FindBin(m_eop));
544 if(fabs(m_p[0])<1.0)
545 pe=h_ebarlp->GetBinContent(h_ebarlp->FindBin(m_eop));
546 else if(fabs(m_p[0])>=1.0)
547 pe=h_ebarhp->GetBinContent(h_ebarhp->FindBin(m_eop));
548
549 }
550 else if(fabs(m_costh[0])>0.85 ){
551 pk=h_kend->GetBinContent(h_kend->FindBin(m_eop));
552 ppi=h_piend->GetBinContent(h_piend->FindBin(m_eop));
553 if(fabs(m_p[0])<1.0)
554 pe=h_eendlp->GetBinContent(h_eendlp->FindBin(m_eop));
555 else if(fabs(m_p[0])>=1.0)
556 pe=h_eendhp->GetBinContent(h_eendhp->FindBin(m_eop));
557 }
558 }
559
560 double lr=0;
561
562 if( pe>0 ){
563
564 if( fabs(m_dedxchi[0])<5 && fabs(m_tofchi[0])<5){
565
566 lr= pe*TMath::Gaus(m_dedxchi[0])*TMath::Gaus(m_tofchi[0])/(pe*TMath::Gaus(m_dedxchi[0])*TMath::Gaus(m_tofchi[0])+
567 0.9*ppi*TMath::Gaus(m_dedxchi[2])*TMath::Gaus(m_tofchi[2])+
568 0.1*pk*TMath::Gaus(m_dedxchi[3])*TMath::Gaus(m_tofchi[3]));
569 }
570
571 else if( fabs(m_dedxchi[0])<5 && fabs(m_tofchi[0])>=5){
572
573 lr= pe*TMath::Gaus(m_dedxchi[0])/(pe*TMath::Gaus(m_dedxchi[0])+
574 0.9*ppi*TMath::Gaus(m_dedxchi[2])+
575 0.1*pk*TMath::Gaus(m_dedxchi[3]));
576 }
577
578 else if( fabs(m_tofchi[0])<5 && fabs(m_dedxchi[0])>=5){
579
580 lr= pe*TMath::Gaus(m_tofchi[0])/(pe*TMath::Gaus(m_tofchi[0])+
581 0.9*ppi*TMath::Gaus(m_tofchi[2])+
582 0.1*pk*TMath::Gaus(m_tofchi[3]));
583 }
584
585
586 }
587 else{
588 if( fabs(m_dedxchi[0])<5 && fabs(m_tofchi[0])<5){
589
590 lr= TMath::Gaus(m_dedxchi[0])*TMath::Gaus(m_tofchi[0])/(TMath::Gaus(m_dedxchi[0])*TMath::Gaus(m_tofchi[0])+
591 0.9*TMath::Gaus(m_dedxchi[2])*TMath::Gaus(m_tofchi[2])+
592 0.1*TMath::Gaus(m_dedxchi[3])*TMath::Gaus(m_tofchi[3]));
593 }
594
595 else if( fabs(m_dedxchi[0])<5 && fabs(m_tofchi[0])>=5){
596
597 lr= TMath::Gaus(m_dedxchi[0])/(TMath::Gaus(m_dedxchi[0])+
598 0.9*TMath::Gaus(m_dedxchi[2])+
599 0.1*TMath::Gaus(m_dedxchi[3]));
600 }
601
602 else if( fabs(m_tofchi[0])<5 && fabs(m_dedxchi[0])>=5){
603
604 lr= TMath::Gaus(m_tofchi[0])/(TMath::Gaus(m_tofchi[0])+
605 0.9*TMath::Gaus(m_tofchi[2])+
606 0.1*TMath::Gaus(m_tofchi[3]));
607 }
608
609 }
610
611 if( lr>0.8 )
612 return true;
613 return false;
614 */
615
616}
617
619
620 if(m_prob[2]>=0.00 && m_prob[2]>=m_prob[3])
621 return true;
622
623 /*
624 if(m_dedxonly[2] && m_prob[2]>0.10){
625 return true;
626 }
627 */
628 return false;
629}
630
632
633 if(m_prob[3]>=0.00 && m_prob[3]>=m_prob[2]){
634 return true;
635 }
636 /*
637 if(m_dedxonly[3] && m_prob[3]>0.10){
638 return true;
639 }
640 */
641
642 return false;
643}
644
645
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
************Class m_ypar INTEGER m_KeyWgt INTEGER m_nphot INTEGER m_KeyGPS INTEGER m_IsBeamPolarized INTEGER m_EvtGenInterface DOUBLE PRECISION m_Emin DOUBLE PRECISION m_sphot DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_q2 DOUBLE PRECISION m_PolBeam2 DOUBLE PRECISION m_xErrPb *COMMON c_KK2f $ !CMS energy average $ !Spin Polarization vector first beam $ !Spin Polarization vector second beam $ !Beam energy spread[GeV] $ !minimum hadronization energy[GeV] $ !input READ never touch them !$ !debug facility $ !maximum weight $ !inverse alfaQED $ !minimum real photon energy
Definition: KK2f.h:50
double energy() const
Definition: DstEmcShower.h:45
double chi(int i) const
Definition: DstMdcDedx.h:58
static void setPidType(PidType pidType)
bool isMdcDedxValid()
Definition: EvtRecTrack.h:45
RecMdcDedx * mdcDedx()
Definition: EvtRecTrack.h:55
bool isMdcKalTrackValid()
Definition: EvtRecTrack.h:44
SmartRefVector< RecTofTrack > tofTrack()
Definition: EvtRecTrack.h:57
bool isTofTrackValid()
Definition: EvtRecTrack.h:46
RecEmcShower * emcShower()
Definition: EvtRecTrack.h:58
bool isEmcShowerValid()
Definition: EvtRecTrack.h:47
RecMdcKalTrack * mdcKalTrack()
Definition: EvtRecTrack.h:54
static const InterfaceID & interfaceID()
Definition: ISimplePIDSvc.h:17
const HepVector & getZHelix() const
HepVector & getZHelixP()
HepVector & getZHelixK()
HepVector & getZHelixE()
HepVector & getZHelixMu()
SimplePIDSvc(const std::string &name, ISvcLocator *svcLoc)
void getConst(int run)
bool iselectron(bool eop=false)
virtual ~SimplePIDSvc()
void preparePID(EvtRecTrack *track)
virtual StatusCode initialize()
virtual StatusCode finalize()
virtual StatusCode queryInterface(const InterfaceID &riid, void **ppvIF)
bool is_barrel() const
Definition: TofHitStatus.h:26
unsigned int layer() const
Definition: TofHitStatus.h:28
void setStatus(unsigned int status)
bool is_counter() const
Definition: TofHitStatus.h:24
int num[96]
Definition: ranlxd.c:373