StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StuProbabilityPidAlgorithm.cxx
1 /***************************************************************************
2  *
3  * $Id: StuProbabilityPidAlgorithm.cxx,v 1.34 2004/09/19 00:07:28 perev Exp $
4  *
5  * Author:Aihong Tang, Richard Witt(FORTRAN version). Kent State University
6  * Send questions to aihong@cnr.physics.kent.edu
7  ***************************************************************************
8  *
9  * Description: A functor that do PID base on Probability (Amplitude) info.
10  *
11  ***************************************************************************
12  *
13  * $Log: StuProbabilityPidAlgorithm.cxx,v $
14  * Revision 1.34 2004/09/19 00:07:28 perev
15  * Small Walgrind leak fixed
16  *
17  * Revision 1.33 2004/04/09 15:46:16 aihong
18  * add isPIDTableRead()
19  *
20  * Revision 1.32 2003/09/02 17:58:09 perev
21  * gcc 3.2 updates + WarnOff
22  *
23  * Revision 1.31 2003/06/24 02:53:14 aihong
24  * update for dAu PIDtable
25  *
26  * Revision 1.30 2003/05/02 21:32:08 aihong
27  * destroy myBandBGFcn in destructor
28  *
29  * Revision 1.29 2003/04/30 20:37:12 perev
30  * Warnings cleanup. Modified lines marked VP
31  *
32  * Revision 1.28 2002/12/20 20:26:44 aihong
33  * let it handle P02gd table
34  *
35  * Revision 1.27 2002/12/11 15:35:53 aihong
36  * put in fabs in processPIDAsFuntion()
37  *
38  * Revision 1.26 2002/04/19 16:36:41 perev
39  * bug fix,init to zero
40  *
41  * Revision 1.25 2002/01/17 03:25:27 aihong
42  * add production Tag to take care of different centrality def. between different productions
43  *
44  * Revision 1.24 2001/03/21 18:16:13 aihong
45  * constructor without StEvent added
46  *
47  * Revision 1.23 2001/03/21 17:54:30 aihong
48  * add processPIDAsFunction()
49  *
50  * Revision 1.22 2000/12/28 21:00:57 aihong
51  * remove temp. fix
52  *
53  * Revision 1.21 2000/12/26 23:27:05 aihong
54  * temperory fix for e amp tail
55  *
56  * Revision 1.20 2000/12/20 16:55:16 aihong
57  * let it survive when no support PIDTable is present
58  *
59  * Revision 1.19 2000/12/18 23:51:36 aihong
60  * mExtrap related bug fixed
61  *
62  * Revision 1.10 2000/08/16 12:46:07 aihong
63  * bug killed
64  *
65  * Revision 1.9 2000/08/15 23:04:18 aihong
66  * speed it up by looking up table
67  *
68  * Revision 1.8 2000/07/22 22:45:27 aihong
69  * change include path
70  *
71  * Revision 1.7 2000/07/12 16:29:38 aihong
72  * change to avoid possible name confliction from StarClassLibary
73  *
74  * Revision 1.6 2000/05/24 14:35:41 ullrich
75  * Added 'const' to compile on Sun CC5.
76  *
77  * Revision 1.5 2000/05/05 19:25:39 aihong
78  * modified ctor
79  *
80  * Revision 1.2 2000/03/09 20:45:04 aihong
81  * add head for Log
82  *
83  **************************************************************************/
84 #include <float.h>
85 #include <Stsstream.h>
86 #include "Stiostream.h" //this line should be deleted later
87 #include "TFile.h"
88 #include "TF1.h"
89 
90 #include "StMessMgr.h"
91 #include "StPhysicalHelixD.hh"
92 #include "PhysicalConstants.h"
93 
94 #include "StuProbabilityPidAlgorithm.h"
95 
96 #include "StEventUtilities/BetheBlochFunction.hh"
97 #include "StEventUtilities/MaxllBoltz.hh"
98 #include "StEventUtilities/Linear.hh"
99 
100 #include "StEventUtilities/StuFtpcRefMult.hh"
101 #include "StEventUtilities/StuRefMult.hh"
102 
103 //TMap::FindObject goes wild!! TMap::GetValue works.
104 
105 //-------------------------------
106 StuProbabilityPidAlgorithm::StuProbabilityPidAlgorithm(StEvent& ev):
107  mPionMinusProb(0.),
108  mElectronProb(0.),
109  mKaonMinusProb(0.),
110  mAntiProtonProb(0.),
111  mPionPlusProb(0.),
112  mPositronProb(0.),
113  mKaonPlusProb(0.),
114  mProtonProb(0.){
115 
116 
117  PID[0]=-1;//should be sth.standard say unIdentified.
118  PID[1]=-1;
119  PID[2]=-1;
120  PID[3]=-1;
121 
122  mProb[0]=0;
123  mProb[1]=0;
124  mProb[2]=0;
125  mProb[3]=0;
126 
127  table = StParticleTable::instance();
128 
129  mExtrap=false;
130  mEvent=&ev;
131 
132  //init funtions
133  myBandBGFcn
134  =new TF1("myBandBGFcn",BetheBlochFunction, 0.,5., 7);
135 
136  Double_t myPars[7]
137  ={ 1.072, 0.3199, 1.66032e-07, 1, 1, 2.71172e-07, 0.0005 };
138 
139  myBandBGFcn->SetParameters(&myPars[0]);
140 
141 }
142 
143 //-------------------------------
144 StuProbabilityPidAlgorithm::StuProbabilityPidAlgorithm():
145  mPionMinusProb(0.),
146  mElectronProb(0.),
147  mKaonMinusProb(0.),
148  mAntiProtonProb(0.),
149  mPionPlusProb(0.),
150  mPositronProb(0.),
151  mKaonPlusProb(0.),
152  mProtonProb(0.){
153 
154 
155  PID[0]=-1;//should be sth.standard say unIdentified.
156  PID[1]=-1;
157  PID[2]=-1;
158  PID[3]=-1;
159 
160  mProb[0]=0;
161  mProb[1]=0;
162  mProb[2]=0;
163  mProb[3]=0;
164 
165  table = StParticleTable::instance();
166 
167  mExtrap=false;
168 
169  //init funtions
170  myBandBGFcn
171  =new TF1("myBandBGFcn",BetheBlochFunction, 0.,5., 7);
172 
173  Double_t myPars[7]
174  ={ 1.072, 0.3199, 1.66032e-07, 1, 1, 2.71172e-07, 0.0005 };
175 
176  myBandBGFcn->SetParameters(&myPars[0]);
177 
178 }
179 
180 //-------------------------------
181 StuProbabilityPidAlgorithm::~StuProbabilityPidAlgorithm(){
182  /* no op */
183  delete myBandBGFcn;
184 }
185 //-------------------------------
186 void StuProbabilityPidAlgorithm::setDedxMethod(StDedxMethod method){
187  StuProbabilityPidAlgorithm::mDedxMethod=method;
188 }
189 
190 //-------------------------------
191 bool StuProbabilityPidAlgorithm::isPIDTableRead(){
192  return StuProbabilityPidAlgorithm::mPIDTableRead;
193 }
194 
195 
196 //-------------------------------
197 StParticleDefinition* StuProbabilityPidAlgorithm::mostLikelihoodParticle(){
198 
199 
200  return table->findParticleByGeantId(PID[0]);
201 }
202 
203 //-------------------------------
204 StParticleDefinition* StuProbabilityPidAlgorithm::secondLikelihoodParticle(){
205 
206 
207  return table->findParticleByGeantId(PID[1]);
208 }
209 
210 //-------------------------------
211 StParticleDefinition* StuProbabilityPidAlgorithm::thirdLikelihoodParticle(){
212 
213 
214  return table->findParticleByGeantId(PID[2]);
215 }
216 //-------------------------------
217 StParticleDefinition* StuProbabilityPidAlgorithm::getParticle(int i){
218 
219  if (i>=0 && i<4){
220  return table->findParticleByGeantId(PID[i]);
221  } else {
222 
223  gMessMgr->Error()<<"StuProbabilityPidAlgorithm::getParticle(int i), i must be 0,1,2,3 only. "<<endm;
224 
225  return 0;
226  }
227  }
228 //-------------------------------
229 double StuProbabilityPidAlgorithm::getProbability(int i){
230  if (i>=0 && i<4){
231  return mProb[i];
232  } else {
233 
234  gMessMgr->Error()<<"StuProbabilityPidAlgorithm::getProbability(int i), i must be 0,1,2,3 only. "<<endm;
235 
236  return 0.0;
237  }
238 
239 }
240 
241 //-------------------------------
242 double StuProbabilityPidAlgorithm::mostLikelihoodProbability(){
243  return mProb[0];
244  }
245 //-------------------------------
246 double StuProbabilityPidAlgorithm::secondLikelihoodProbability(){
247  return mProb[1];
248  }
249 //-------------------------------
250 double StuProbabilityPidAlgorithm::thirdLikelihoodProbability(){
251  return mProb[2];
252  }
253 
254 //-------------------------------
255 bool StuProbabilityPidAlgorithm::isExtrap(){
256  return mExtrap;
257 }
258 
259 
260 //-------------------------------
261 StParticleDefinition* StuProbabilityPidAlgorithm::operator() (const StTrack& theTrack, const StSPtrVecTrackPidTraits& traits){
262 
263 
264 
265  PID[0]=-1;//should be sth.standard say unIdentified.
266  PID[1]=-1;
267  PID[2]=-1;
268  PID[3]=-1;
269 
270  mProb[0]=0;
271  mProb[1]=0;
272  mProb[2]=0;
273  mProb[3]=0;
274 
275  mExtrap=false;
276 
277 
278  mPionMinusProb=0.;
279  mElectronProb=0.;
280  mKaonMinusProb=0.;
281  mAntiProtonProb=0.;
282  mPionPlusProb=0.;
283  mPositronProb=0.;
284  mKaonPlusProb=0.;
285  mProtonProb=0.;
286 
287  if (mPIDTableRead) {
288 
289  double rig =0.0;
290  double dedx =0.0;
291  double dca =0.0; //in units of cm.
292  int nhits =0;
293  int charge =0;
294  double eta =0.;
295  double cent =0.; // % central
296 
297 
298  StPrimaryVertex* primaryVtx=mEvent->primaryVertex();
299  const StPhysicalHelixD& helix=theTrack.geometry()->helix();
300  dca=helix.distance(primaryVtx->position());
301 
302  //cent in cross section
303 
304  if (mProductionTag){
305 
306  if ( (mProductionTag->GetString()).Contains("P01gl")
307  || (mProductionTag->GetString()).Contains("P02gd") ){
308  cent = getCentrality(uncorrectedNumberOfNegativePrimaries(*mEvent));
309  } else if ( (mProductionTag->GetString()).Contains("P03ia_dAu") ){
310  cent = getCentrality(uncorrectedNumberOfFtpcEastPrimaries(*mEvent));
311  } else {
312  gMessMgr->Error()<<"Production tag "<<mProductionTag->GetString().Data()<<" in PIDTable is filled but its name is not recognized ! "<<endm;
313  }
314 
315 
316  } else { //the first PID table has no production tag
317  cent = getCentrality(uncorrectedNumberOfNegativePrimaries(*mEvent));
318  }
319 
320 
321 
322  const StDedxPidTraits* dedxPidTr=0;
323 
324 
325  charge=(theTrack.geometry())->charge();
326 
327  for (int itrait = 0; itrait < int(traits.size()); itrait++){
328 
329  dedxPidTr = 0;
330  if (traits[itrait]->detector() == kTpcId) {
331  //
332  // tpc pid trait
333  //
334  const StTrackPidTraits* thisTrait = traits[itrait];
335  //
336  // perform cast to make the pid trait a dedx trait
337  //
338  dedxPidTr = dynamic_cast<const StDedxPidTraits*>(thisTrait);
339 
340  }
341  if (dedxPidTr && dedxPidTr->method() == mDedxMethod) break;
342 
343  }
344 
345  if (dedxPidTr) {
346  dedx=dedxPidTr->mean();
347  nhits=dedxPidTr->numberOfPoints();
348  }
349 
350 
351  if (dedx!=0.0 && nhits>=0 //dedx ==0.0 no sense
352  && thisPEnd > 0. && thisEtaEnd > 0. // *End ==0, no PIDTable read.
353  && thisNHitsEnd > 0. ){
354 
355  const StThreeVectorF& p=theTrack.geometry()->momentum();
356  rig=double(p.mag()/charge);
357 
358 
359 
360  if (mProductionTag){ //for AuAu, +/- eta were folded together when building PID
361  //table, for dAu, + and - eta were treated differently.
362 
363  if ( (mProductionTag->GetString()).Contains("P01gl")
364  || (mProductionTag->GetString()).Contains("P02gd") ){
365  eta=fabs(p.pseudoRapidity());
366  } else if ( (mProductionTag->GetString()).Contains("P03ia_dAu") ){
367  eta=p.pseudoRapidity();
368  } else {
369  gMessMgr->Error()<<"Production tag "<<mProductionTag->GetString().Data()<<" in PIDTable is filled but its name is not recognized ! "<<endm;
370  }
371 
372  } else { //the first PID table has no production tag
373  eta = fabs(p.pseudoRapidity());
374  }
375 
376  rig = fabs(rig);
377  dedx = (dedx>thisDedxStart) ? dedx : thisDedxStart;
378  rig = (rig >thisPStart) ? rig : thisPStart;
379  rig = (rig <thisPEnd ) ? rig : thisPEnd*0.9999;
380  eta = (eta >thisEtaStart) ? eta : thisEtaStart;
381  eta = (eta <thisEtaEnd ) ? eta : thisEtaEnd*0.9999;
382  nhits = (nhits > int(thisNHitsStart)) ? nhits : int(thisNHitsStart);
383  nhits = (nhits < int(thisNHitsEnd) ) ? nhits : int(thisNHitsEnd-1);
384 
385  //----------------get all info. I want for a track. now do PID
386 
387  setCalibrations(eta, nhits);
388 
389  if (dedx<thisDedxEnd){
390 
391  fillPIDByLookUpTable(cent, dca, charge,rig, eta, nhits,dedx);
392 
393  } else { lowRigPID(rig,dedx,charge);}
394 
395  } else if (dedx==0.0){ fillAsUnknown();}
396 
397  // do not do deuteron or higher
398  myBandBGFcn->SetParameter(3,1);
399  myBandBGFcn->SetParameter(4,1.45);
400  if (dedx>myBandBGFcn->Eval(rig,0,0)) fillAsUnknown();
401 
402  } else fillAsUnknown();
403 
404  fillPIDHypothis();
405 
406  return table->findParticleByGeantId(PID[0]);
407 
408  }
409 
410 //-------------------------------
411 void StuProbabilityPidAlgorithm::lowRigPID(double rig,double dedx,int theCharge){
412 
413  double upper;
414  double lower;
415  double rigidity=fabs(rig);
416  double mdedx=dedx;
417 
418  double fakeMass=0.;
419 
420  //pion
421  fakeMass=0.32075026;
422  myBandBGFcn->SetParameter(3,1);
423  myBandBGFcn->SetParameter(4,0.32075026 );
424 
425 
426  lower =0.;
427  upper =myBandBGFcn->Eval(rigidity,0,0);
428 
429  if (mdedx>lower && mdedx<upper){
430  PID[0]=(theCharge>0.0)? 8 : 9; //pi+/-
431  mProb[0]=1.0;
432  mProb[1]=0.0;
433  mProb[2]=0.0;
434  mProb[3]=0.0;
435  }
436 
437  lower = upper;
438 
439  //kaon
440  fakeMass=0.709707;
441  myBandBGFcn->SetParameter(3,1);
442  myBandBGFcn->SetParameter(4,fakeMass);
443 
444  upper =myBandBGFcn->Eval(rigidity,0,0);
445 
446  if (mdedx>lower && mdedx<upper){
447  PID[0]=(theCharge>0.0)? 11:12; //k+/-
448  mProb[0]=1.0;
449  mProb[1]=0.0;
450  mProb[2]=0.0;
451  mProb[3]=0.0;
452  }
453 
454 
455  lower = upper;
456 
457  //proton/pBar
458  fakeMass=1.45;
459  myBandBGFcn->SetParameter(3,1);
460  myBandBGFcn->SetParameter(4,fakeMass);
461 
462  upper =myBandBGFcn->Eval(rigidity,0,0);
463 
464  if (mdedx>lower && mdedx<upper){
465  PID[0]=(theCharge>0.0)? 14:15; //proton/antiproton
466  mProb[0]=1.0;
467  mProb[1]=0.0;
468  mProb[2]=0.0;
469  mProb[3]=0.0;
470  }
471 
472  /*
473  lower = upper;
474 
475  //deuteron
476  fakeMass=2.4;
477  myBandBGFcn->SetParameter(3,1);
478  myBandBGFcn->SetParameter(4,fakeMass);
479 
480  upper =myBandBGFcn->Eval(rigidity,0,0);
481 
482  if (mdedx>lower && mdedx<upper){
483  PID[0]=45; //deuteron
484  mProb[0]=1.0;
485  mProb[1]=0.0;
486  mProb[2]=0.0;
487  mProb[3]=0.0;
488  }
489 
490  lower = upper;
491 
492  //triton
493  m = -1.5374; //New slope needed for deuterons and tritons.
494  a = 1.8121e-5;
495  upper = a*::pow(rigidity,m);
496  if (mdedx>lower && mdedx<upper){
497  PID[0]=46; //triton
498  mProb[0]=1.0;
499  mProb[1]=0.0;
500  mProb[2]=0.0;
501  mProb[3]=0.0;
502  }
503 
504  */
505 
506 
507 
508 }
509 
510 
511 
512 //-------------------------------
513 void StuProbabilityPidAlgorithm::fillAsUnknown(){
514 
515  for (int i=0; i<4; i++) {
516  PID[i]=-1; mProb[i]=-1;
517  }
518 }
519 
520 
521 //-------------------------------
522 void StuProbabilityPidAlgorithm::readParametersFromFile(TString fileName){
523 
524 
525 
526  TFile f(fileName,"READ");
527 
528  if (f.IsOpen()){
529 
530  delete StuProbabilityPidAlgorithm::mEAmp ;
531  delete StuProbabilityPidAlgorithm::mECenter ;
532  delete StuProbabilityPidAlgorithm::mESigma ;
533  delete StuProbabilityPidAlgorithm::mPiAmp ;
534  delete StuProbabilityPidAlgorithm::mPiCenter ;
535  delete StuProbabilityPidAlgorithm::mPiSigma ;
536  delete StuProbabilityPidAlgorithm::mKAmp ;
537  delete StuProbabilityPidAlgorithm::mKCenter ;
538  delete StuProbabilityPidAlgorithm::mKSigma ;
539  delete StuProbabilityPidAlgorithm::mPAmp ;
540  delete StuProbabilityPidAlgorithm::mPCenter ;
541  delete StuProbabilityPidAlgorithm::mPSigma ;
542  delete StuProbabilityPidAlgorithm::mEqualyDividableRangeStartSet;
543  delete StuProbabilityPidAlgorithm::mEqualyDividableRangeEndSet;
544  delete StuProbabilityPidAlgorithm::mEqualyDividableRangeNBinsSet;
545  delete StuProbabilityPidAlgorithm::mNoEqualyDividableRangeNBinsSet;
546  delete StuProbabilityPidAlgorithm::mMultiBinEdgeSet ;
547  delete StuProbabilityPidAlgorithm::mDcaBinEdgeSet ;
548  delete StuProbabilityPidAlgorithm::mBBPrePar;
549  delete StuProbabilityPidAlgorithm::mBBTurnOver;
550  delete StuProbabilityPidAlgorithm::mBBOffSet;
551  delete StuProbabilityPidAlgorithm::mBBScale;
552  delete StuProbabilityPidAlgorithm::mBBSaturate;
553  delete StuProbabilityPidAlgorithm::mProductionTag;
554 
555  StuProbabilityPidAlgorithm::mEAmp =(TVectorD* )f.Get("eAmp");
556  StuProbabilityPidAlgorithm::mECenter =(TVectorD* )f.Get("eCenter");
557  StuProbabilityPidAlgorithm::mESigma =(TVectorD* )f.Get("eSigma");
558 
559  StuProbabilityPidAlgorithm::mPiAmp =(TVectorD* )f.Get("piAmp");
560  StuProbabilityPidAlgorithm::mPiCenter =(TVectorD* )f.Get("piCenter");
561  StuProbabilityPidAlgorithm::mPiSigma =(TVectorD* )f.Get("piSigma");
562 
563  StuProbabilityPidAlgorithm::mKAmp =(TVectorD* )f.Get("kAmp");
564  StuProbabilityPidAlgorithm::mKCenter =(TVectorD* )f.Get("kCenter");
565  StuProbabilityPidAlgorithm::mKSigma =(TVectorD* )f.Get("kSigma");
566 
567  StuProbabilityPidAlgorithm::mPAmp =(TVectorD* )f.Get("pAmp");
568  StuProbabilityPidAlgorithm::mPCenter =(TVectorD* )f.Get("pCenter");
569  StuProbabilityPidAlgorithm::mPSigma =(TVectorD* )f.Get("pSigma");
570 
571  StuProbabilityPidAlgorithm::mEqualyDividableRangeStartSet
572  =(TVectorD* )f.Get("EqualyDividableRangeStartSet");
573  StuProbabilityPidAlgorithm::mEqualyDividableRangeEndSet
574  =(TVectorD* )f.Get("EqualyDividableRangeEndSet");
575  StuProbabilityPidAlgorithm::mEqualyDividableRangeNBinsSet
576  =(TVectorD* )f.Get("EqualyDividableRangeNBinsSet");
577  StuProbabilityPidAlgorithm::mNoEqualyDividableRangeNBinsSet
578  =(TVectorD* )f.Get("NoEqualyDividableRangeNBinsSet");
579 
580  StuProbabilityPidAlgorithm::mMultiBinEdgeSet
581  =(TVectorD* )f.Get("MultiBinEdgeSet");
582  StuProbabilityPidAlgorithm::mDcaBinEdgeSet
583  =(TVectorD* )f.Get("DcaBinEdgeSet");
584 
585 
586  StuProbabilityPidAlgorithm::mBBPrePar
587  =(TVectorD* )f.Get("BBPrePar");
588  StuProbabilityPidAlgorithm::mBBTurnOver
589  =(TVectorD* )f.Get("BBTurnOver");
590  StuProbabilityPidAlgorithm::mBBOffSet
591  =(TVectorD* )f.Get("BBOffSet");
592  StuProbabilityPidAlgorithm::mBBScale
593  =(TVectorD* )f.Get("BBScale");
594  StuProbabilityPidAlgorithm::mBBSaturate
595  =(TVectorD* )f.Get("BBSaturate");
596 
597  StuProbabilityPidAlgorithm::mProductionTag
598  =(TObjString* )f.Get("productionTag");
599 
600 
601 
602  StuProbabilityPidAlgorithm::thisMultBins
603  =int((*mNoEqualyDividableRangeNBinsSet)(0));
604  StuProbabilityPidAlgorithm::thisDcaBins
605  =int((*mNoEqualyDividableRangeNBinsSet)(1));
606  StuProbabilityPidAlgorithm::thisChargeBins
607  =int((*mNoEqualyDividableRangeNBinsSet)(2));
608 
609 
610  StuProbabilityPidAlgorithm::thisPBins
611  =int((*mEqualyDividableRangeNBinsSet)(1));
612  StuProbabilityPidAlgorithm::thisEtaBins
613  =int((*mEqualyDividableRangeNBinsSet)(2));
614  StuProbabilityPidAlgorithm::thisNHitsBins
615  =int((*mEqualyDividableRangeNBinsSet)(3));
616 
617 
618  StuProbabilityPidAlgorithm::thisDedxStart
619  =(*mEqualyDividableRangeStartSet)(0);
620  StuProbabilityPidAlgorithm::thisPStart
621  =(*mEqualyDividableRangeStartSet)(1);
622  StuProbabilityPidAlgorithm::thisEtaStart
623  =(*mEqualyDividableRangeStartSet)(2);
624  StuProbabilityPidAlgorithm::thisNHitsStart
625  =(*mEqualyDividableRangeStartSet)(3);
626 
627  StuProbabilityPidAlgorithm::thisDedxEnd
628  =(*mEqualyDividableRangeEndSet)(0);
629  StuProbabilityPidAlgorithm::thisPEnd
630  =(*mEqualyDividableRangeEndSet)(1);
631  StuProbabilityPidAlgorithm::thisEtaEnd
632  =(*mEqualyDividableRangeEndSet)(2);
633  StuProbabilityPidAlgorithm::thisNHitsEnd
634  =(*mEqualyDividableRangeEndSet)(3);
635 
636 
637 
638 
639  StuProbabilityPidAlgorithm::mPIDTableRead=true;
640 
641  } else if (!f.IsOpen()) {
642 
643  gMessMgr->Error()<<"Data file "<<fileName<<" open failed "<<endm;
644  return;
645  }
646 
647 }
648 
649 /* uncomment this method when data base is available and stable.
650 //-------------------------------
651 void StuProbabilityPidAlgorithm::readParametersFromTable(St_Table* tb){
652 
653 
654  if (mDataTable!=tb) {//database changed.
655  StuProbabilityPidAlgorithm::refreshParameters(tb);
656  mDataTable=tb;
657  }else if (mDataTable==tb) {return;}//database no change. no need refreshing.
658 
659 }
660 */
661 
662 /*uncomment this method when data base is available and stable.
663 //-------------------------------
664 void StuProbabilityPidAlgorithm::refreshParameters(St_Table* theTable){
665  int i;
666 
667  if (mDataSet.GetEntries()>0) mDataSet.Delete();
668 
669  St_tpcDedxPidAmpDb* temp=(St_tpcDedxPidAmpDb *)theTable;
670  tpcDedxPidAmpDb_st* pars=(tpcDedxPidAmpDb_st *)temp->GetTable();
671 
672  StPidAmpChannelInfoOut* theInfoOut=new StPidAmpChannelInfoOut(0,45,0,FLT_MAX);
673 
674  TObjArray* theOnlyChannel=new TObjArray();
675  theOnlyChannel->Add(theInfoOut);
676 
677  StuProbabilityPidAlgorithm::readAType(StElectron::instance(),pars->eMeanPar, pars->eAmpPar, pars->eSigPar,pars->gasCalib, theOnlyChannel);
678 
679 
680 
681  mDataSet.Add(theOnlyChannel);
682 
683 
684 }
685 
686 */
687 
688 
689 
690 //-------------------------------
691 int StuProbabilityPidAlgorithm::thisMultBins=0;
692 int StuProbabilityPidAlgorithm::thisDcaBins=0;
693 int StuProbabilityPidAlgorithm::thisChargeBins=0;
694 
695 int StuProbabilityPidAlgorithm::thisPBins=0;
696 int StuProbabilityPidAlgorithm::thisEtaBins=0;
697 int StuProbabilityPidAlgorithm::thisNHitsBins=0;
698 
699 double StuProbabilityPidAlgorithm::thisDedxStart=0.;
700 double StuProbabilityPidAlgorithm::thisDedxEnd=0;
701 double StuProbabilityPidAlgorithm::thisPStart=0;
702 double StuProbabilityPidAlgorithm::thisPEnd=0;
703 double StuProbabilityPidAlgorithm::thisEtaStart=0;
704 double StuProbabilityPidAlgorithm::thisEtaEnd=0;
705 double StuProbabilityPidAlgorithm::thisNHitsStart=0;
706 double StuProbabilityPidAlgorithm::thisNHitsEnd=0;
707 
708 
709 bool StuProbabilityPidAlgorithm::mPIDTableRead=false;
710 
711 StDedxMethod StuProbabilityPidAlgorithm::mDedxMethod=kTruncatedMeanId;
712 
713 TVectorD* StuProbabilityPidAlgorithm::mEAmp = new TVectorD();
714 TVectorD* StuProbabilityPidAlgorithm::mECenter = new TVectorD();
715 TVectorD* StuProbabilityPidAlgorithm::mESigma = new TVectorD();
716 
717 TVectorD* StuProbabilityPidAlgorithm::mPiAmp = new TVectorD();
718 TVectorD* StuProbabilityPidAlgorithm::mPiCenter = new TVectorD();
719 TVectorD* StuProbabilityPidAlgorithm::mPiSigma = new TVectorD();
720 
721 TVectorD* StuProbabilityPidAlgorithm::mKAmp = new TVectorD();
722 TVectorD* StuProbabilityPidAlgorithm::mKCenter = new TVectorD();
723 TVectorD* StuProbabilityPidAlgorithm::mKSigma = new TVectorD();
724 
725 TVectorD* StuProbabilityPidAlgorithm::mPAmp = new TVectorD();
726 TVectorD* StuProbabilityPidAlgorithm::mPCenter = new TVectorD();
727 TVectorD* StuProbabilityPidAlgorithm::mPSigma = new TVectorD();
728 
729 
730 TVectorD* StuProbabilityPidAlgorithm::mEqualyDividableRangeStartSet = new TVectorD();
731 TVectorD* StuProbabilityPidAlgorithm::mEqualyDividableRangeEndSet = new TVectorD();
732 TVectorD* StuProbabilityPidAlgorithm::mEqualyDividableRangeNBinsSet = new TVectorD();
733 TVectorD* StuProbabilityPidAlgorithm::mNoEqualyDividableRangeNBinsSet = new TVectorD();
734 
735 TVectorD* StuProbabilityPidAlgorithm::mMultiBinEdgeSet = new TVectorD();
736 TVectorD* StuProbabilityPidAlgorithm::mDcaBinEdgeSet = new TVectorD();
737 
738 TVectorD* StuProbabilityPidAlgorithm::mBBPrePar = new TVectorD();
739 TVectorD* StuProbabilityPidAlgorithm::mBBTurnOver = new TVectorD();
740 TVectorD* StuProbabilityPidAlgorithm::mBBOffSet = new TVectorD();
741 TVectorD* StuProbabilityPidAlgorithm::mBBScale = new TVectorD();
742 TVectorD* StuProbabilityPidAlgorithm::mBBSaturate = new TVectorD();
743 
744 TObjString* StuProbabilityPidAlgorithm::mProductionTag = new TObjString();
745 
746 
747 
748 
749 //-------------------------------
750 //St_Table* StuProbabilityPidAlgorithm::mDataTable=0;
751 
752 
753 
754 
755 //-------------------------------
756 void StuProbabilityPidAlgorithm::fillPIDByLookUpTable(double myCentrality, double myDca, int myCharge, double myRig, double myEta, int myNhits, double myDedx)
757 {
758 
759 
760  //assume bound has been checked before they enter this func.
761  int theMultBin=getCentralityBin(myCentrality);
762  int theDcaBin = (myDca>(*mDcaBinEdgeSet)(1)) ? 1 : 0;
763  int theChargeBin=(myCharge>0) ? 1 : 0;
764  int thePBin=int(thisPBins*myRig/(thisPEnd-thisPStart));
765  int theEtaBin=int(thisEtaBins*(myEta-thisEtaStart)/(thisEtaEnd-thisEtaStart));
766  int theNHitsBin=int(thisNHitsBins*float(myNhits)/(thisNHitsEnd-thisNHitsStart));
767 
768 
769  int totalEntry
770  = thisMultBins*thisDcaBins*thisChargeBins*thisPBins*thisEtaBins*thisNHitsBins;
771 
772  int positionPointer=0;
773 
774  totalEntry=totalEntry/thisMultBins;
775  positionPointer=positionPointer+totalEntry*theMultBin;
776 
777  totalEntry=totalEntry/thisDcaBins;
778  positionPointer=positionPointer+totalEntry*theDcaBin;
779 
780  totalEntry=totalEntry/thisChargeBins;
781  positionPointer=positionPointer+totalEntry*theChargeBin;
782 
783  totalEntry=totalEntry/thisPBins;
784  positionPointer=positionPointer+totalEntry*thePBin;
785 
786  totalEntry=totalEntry/thisEtaBins;
787  positionPointer=positionPointer+totalEntry*theEtaBin;
788 
789  totalEntry=totalEntry/thisNHitsBins;
790  positionPointer=positionPointer+totalEntry*theNHitsBin;
791 
792 
793 
794  //now calculate prob.
795  TF1 eGaus("eGaus","gaus",thisDedxStart,thisDedxStart);
796  TF1 piGaus("piGaus","gaus",thisDedxStart,thisDedxStart);
797  TF1 kGaus("kGaus","gaus",thisDedxStart,thisDedxStart);
798  TF1 pGaus("pGaus","gaus",thisDedxStart,thisDedxStart);
799 
800 
801  eGaus.SetParameter(0, (*mEAmp)(positionPointer));
802  eGaus.SetParameter(1, (*mECenter)(positionPointer));
803  eGaus.SetParameter(2, (*mESigma)(positionPointer));
804 
805  piGaus.SetParameter(0, (*mPiAmp)(positionPointer));
806  piGaus.SetParameter(1, (*mPiCenter)(positionPointer));
807  piGaus.SetParameter(2, (*mPiSigma)(positionPointer));
808 
809  kGaus.SetParameter(0, (*mKAmp)(positionPointer));
810  kGaus.SetParameter(1, (*mKCenter)(positionPointer));
811  kGaus.SetParameter(2, (*mKSigma)(positionPointer));
812 
813  pGaus.SetParameter(0, (*mPAmp)(positionPointer));
814  pGaus.SetParameter(1, (*mPCenter)(positionPointer));
815  pGaus.SetParameter(2, (*mPSigma)(positionPointer));
816 
817 
818  double eContribution=0;
819  if (mEAmp && mECenter && mESigma) eContribution = eGaus.Eval(myDedx,0.,0.);
820  // if (fabs(myRig)>0.95) eContribution=0; //will deal with it later.
821 
822 
823  double piContribution=0;
824  if (mPiAmp && mPiCenter && mPiSigma) piContribution = piGaus.Eval(myDedx,0.,0.);
825 
826  double kContribution=0;
827  if (mKAmp && mKCenter && mKSigma) kContribution = kGaus.Eval(myDedx,0.,0.);
828 
829  double pContribution=0;
830  if (mPAmp && mPCenter && mPSigma) pContribution = pGaus.Eval(myDedx,0.,0.);
831 
832 
833  double total = eContribution+piContribution+kContribution+pContribution;
834 
835  double eProb=0; double piProb=0; double kProb=0; double pProb=0;
836 
837  if (total>0) {
838  eProb=eContribution/total;
839  piProb=piContribution/total;
840  kProb=kContribution/total;
841  pProb=pContribution/total;
842  }
843 
844 
845  fill(eProb, int((myCharge>0) ? 2 : 3 ));
846  fill(piProb, int((myCharge>0) ? 8 : 9 ));
847  fill(kProb, int((myCharge>0) ? 11 : 12 ));
848  fill(pProb, int((myCharge>0) ? 14 : 15 ));
849 
850  int nn=-1;
851  float PathHeight=1.0e-7;
852  float halfHeight=(11-2.0)*PathHeight/2.0;
853 
854  if (fabs(myRig) > 0.3) {
855  if ((myDedx<(((*mPiCenter)(positionPointer))+halfHeight-nn*PathHeight))
856  && (myDedx > ( ((*mKCenter)(positionPointer))-halfHeight+nn*PathHeight) ))
857  mExtrap=true;
858  }
859 }
860 //-------------------------------
861 void StuProbabilityPidAlgorithm::fillPIDHypothis(){
862 
863  for (int i=0; i<4; i++){
864 
865  switch (PID[i]) {
866 
867  case 2 : mPositronProb=mProb[i];
868  break;
869  case 3 : mElectronProb=mProb[i];
870  break;
871  case 8 : mPionPlusProb=mProb[i];
872  break;
873  case 9 : mPionMinusProb=mProb[i];
874  break;
875  case 11 : mKaonPlusProb=mProb[i];
876  break;
877  case 12 : mKaonMinusProb=mProb[i];
878  break;
879  case 14 : mProtonProb=mProb[i];
880  break;
881  case 15 : mAntiProtonProb=mProb[i];
882  break;
883  }
884  }
885 
886 }
887 //-------------------------------
888 int StuProbabilityPidAlgorithm::getCentralityBin(double theCent){
889  int theBin=0; //input % centrality. output the corres. bin.
890 
891  for (int mm = (thisMultBins-1); mm>0; mm--) {
892  if (theCent< (*mMultiBinEdgeSet)(mm) ) theBin=mm;
893  }
894 
895  return theBin;
896 
897 }
898 
899 //-------------------------------
900 
901 
902 double StuProbabilityPidAlgorithm::getCentrality(int theMult){
903 
904 
905  if (mProductionTag){
906 
907 
908  if ( (mProductionTag->GetString()).Contains("P01gl")
909  || (mProductionTag->GetString()).Contains("P02gd") ){
910  return getCentrality_P01gl(theMult);
911  } else if ((mProductionTag->GetString()).Contains("P03ia_dAu")){
912  return getCentrality_P03ia_dAu(theMult);
913  } else {
914  gMessMgr->Error()<<"Production tag "<<mProductionTag->GetString().Data()<<" in PIDTable is filled but its name is not recognized ! "<<endm;
915  }
916 
917 
918  } else {//the first PID table does not have a productionTag
919 
920  // limits
921  // For Cut Set 1
922  // 225 ~ 3%
923  // 215 ~ 5%
924  // 200 ~ 7%
925  // 180 ~ 10%
926  // 140 ~ 18%
927  // 130 ~ 20%
928  // 120 ~ 23%
929  // 115 ~ 24%
930  // 100 ~ 28%
931 
932  // highMedLimit = 180; // 10%
933  // medLowLimit = 108; //~26%
934 
935  if (theMult > 225 ) return 0.03;
936  else if (theMult > 215 ) return 0.05;
937  else if (theMult > 200 ) return 0.07;
938  else if (theMult > 180 ) return 0.10;
939  else if (theMult > 140 ) return 0.18;
940  else if (theMult > 130 ) return 0.20;
941  else if (theMult > 120 ) return 0.23;
942  else if (theMult > 115 ) return 0.24;
943  else if (theMult > 100 ) return 0.28;
944  else return 0.99;
945 
946  }
947  return 0.99;
948 
949 
950 }
951 //-------------------------------
952 
953 
954 double StuProbabilityPidAlgorithm::getCentrality_P01gl(int theMult){
955 
956  //from Zhangbu's study
957  //5% 474
958  //10% 409
959  //20% 302
960  //30% 216
961  //40% 149
962  //50% 98
963  //60% 59
964  //70% 32
965  //80% 14
966 
967  if (theMult*2. > 474 ) return 0.03; //use Nch instead of hminus. so *2.
968  else if (theMult*2. > 409 ) return 0.10;
969  else if (theMult*2. > 302 ) return 0.20;
970  else if (theMult*2. > 216 ) return 0.30;
971  else if (theMult*2. > 149 ) return 0.40;
972  else if (theMult*2. > 98 ) return 0.50;
973  else if (theMult*2. > 59 ) return 0.60;
974  else if (theMult*2. > 32 ) return 0.70;
975  else if (theMult*2. > 14 ) return 0.80;
976  else return 0.99;
977 }
978 
979 //-------------------------------
980 
981 
982 double StuProbabilityPidAlgorithm::getCentrality_P03ia_dAu(int theMult){
983 
984  //from Joern's study
985  // * centrality bin multiplicity FTPC east percentOfEvents
986  // * 1 <=11 100-40
987  // * 2 <=17 40-20
988  // * 3 >=18 20-0
989 
990  if (theMult >= 18 ) return 0.19; //do not need to be exact.
991  else if (theMult >= 12 ) return 0.39; //just for gettting bin # correctly.
992  else if (theMult > 0 ) return 0.8;
993  else return 0.99;
994 }
995 
996 //-------------------------------
997 void StuProbabilityPidAlgorithm::fill(double prob, int geantId){
998 
999  if (prob>mProb[0]) {
1000  mProb[3]=mProb[2];
1001  mProb[2]=mProb[1];
1002  mProb[1]=mProb[0];
1003  mProb[0]=prob;
1004 
1005  PID[3]=PID[2];
1006  PID[2]=PID[1];
1007  PID[1]=PID[0];
1008  PID[0]=geantId;
1009  }
1010  else if (prob>mProb[1]){
1011  mProb[3]=mProb[2];
1012  mProb[2]=mProb[1];
1013  mProb[1]=prob;
1014 
1015  PID[3]=PID[2];
1016  PID[2]=PID[1];
1017  PID[1] =geantId;
1018  }
1019  else if (prob>mProb[2]){
1020  mProb[3]=mProb[2];
1021  mProb[2]=prob;
1022  PID[3]=PID[2];
1023  PID[2]=geantId;
1024  }
1025  else if (prob>=mProb[3]){
1026  mProb[3]=prob;
1027  PID[3]=geantId;
1028  }
1029 
1030 }
1031 //------------------------------------------------
1032 int StuProbabilityPidAlgorithm::getCalibPosition(double theEta, int theNHits){
1033 
1034  int theEtaBin=int(thisEtaBins*(theEta-thisEtaStart)/(thisEtaEnd-thisEtaStart));
1035  int theNHitsBin=int(thisNHitsBins*float(theNHits)/(thisNHitsEnd-thisNHitsStart));
1036 
1037  int totalEntry
1038  = thisEtaBins*thisNHitsBins;
1039 
1040  int positionPointer=0;
1041 
1042  totalEntry=totalEntry/thisEtaBins;
1043  positionPointer=positionPointer+totalEntry*theEtaBin;
1044 
1045  totalEntry=totalEntry/thisNHitsBins;
1046  positionPointer=positionPointer+totalEntry*theNHitsBin;
1047 
1048 
1049 
1050  return positionPointer;
1051 }
1052 //--------------------------------------------
1053 void StuProbabilityPidAlgorithm::setCalibrations(double theEta, int theNhits){
1054 
1055  int thePosition=getCalibPosition(theEta,theNhits);
1056 
1057  myBandBGFcn->SetParameter(5,(*mBBScale)(thePosition));
1058  myBandBGFcn->SetParameter(2,(*mBBOffSet)(thePosition));
1059 
1060 
1061  if (mProductionTag){
1062 
1063  if ( (mProductionTag->GetString()).Contains("P01gl")
1064  || (mProductionTag->GetString()).Contains("P02gd")){
1065 
1066  myBandBGFcn->SetParameter(0,(*mBBPrePar)(thePosition));
1067  myBandBGFcn->SetParameter(1,(*mBBTurnOver)(thePosition));
1068  myBandBGFcn->SetParameter(6,(*mBBSaturate)(thePosition));
1069 
1070  }
1071  }
1072 
1073 
1074 }
1075 //-------------------------------------cent,dca,charge,rig,eta,nhits,dedx
1076 void StuProbabilityPidAlgorithm::processPIDAsFunction (double theCent, double theDca, int theCharge, double theRig, double theEta, int theNhits, double theDedx){
1077 
1078 
1079 
1080  PID[0]=-1;//should be sth.standard say unIdentified.
1081  PID[1]=-1;
1082  PID[2]=-1;
1083  PID[3]=-1;
1084 
1085  mProb[0]=0;
1086  mProb[1]=0;
1087  mProb[2]=0;
1088  mProb[3]=0;
1089 
1090  mExtrap=false;
1091 
1092 
1093  mPionMinusProb=0.;
1094  mElectronProb=0.;
1095  mKaonMinusProb=0.;
1096  mAntiProtonProb=0.;
1097  mPionPlusProb=0.;
1098  mPositronProb=0.;
1099  mKaonPlusProb=0.;
1100  mProtonProb=0.;
1101 
1102  if (mPIDTableRead) {
1103 
1104  double rig =fabs(theRig);
1105  double dedx =theDedx;
1106  double dca =theDca; //in units of cm.
1107  int nhits =theNhits;
1108  int charge =theCharge;
1109  double eta =0.;
1110  double cent =theCent; // % central
1111 
1112 
1113  if (mProductionTag){ //for AuAu, +/- eta were folded together when building PID
1114  //table, for dAu, + and - eta were treated differently.
1115 
1116  if ( (mProductionTag->GetString()).Contains("P01gl")
1117  || (mProductionTag->GetString()).Contains("P02gd") ){
1118  eta=fabs(theEta);
1119  } else if ( (mProductionTag->GetString()).Contains("P03ia_dAu") ){
1120  eta=theEta;
1121  } else {
1122  gMessMgr->Error()<<"Production tag "<<mProductionTag->GetString().Data()<<" in PIDTable is filled but its name is not recognized ! "<<endm;
1123  }
1124 
1125  } else { //the first PID table has no production tag
1126  eta = fabs(theEta);
1127  }
1128 
1129 
1130 
1131 
1132 
1133  if (dedx!=0.0 && nhits>=0 //dedx ==0.0 no sense
1134  && thisPEnd > 0. && thisEtaEnd > 0. // *End ==0, no PIDTable read.
1135  && thisNHitsEnd > 0. ){
1136 
1137  rig = fabs(rig);
1138  dedx = (dedx>thisDedxStart) ? dedx : thisDedxStart;
1139  rig = (rig >thisPStart) ? rig : thisPStart;
1140  rig = (rig <thisPEnd ) ? rig : thisPEnd*0.9999;
1141  eta = (eta >thisEtaStart) ? eta : thisEtaStart;
1142  eta = (eta <thisEtaEnd ) ? eta : thisEtaEnd*0.9999;
1143  nhits = (nhits > int(thisNHitsStart)) ? nhits : int(thisNHitsStart);
1144  nhits = (nhits < int(thisNHitsEnd) ) ? nhits : int(thisNHitsEnd-1);
1145 
1146  //----------------get all info. I want for a track. now do PID
1147 
1148  setCalibrations(eta, nhits);
1149 
1150  if (dedx<thisDedxEnd){
1151 
1152  fillPIDByLookUpTable(cent, dca, charge,rig, eta, nhits,dedx);
1153 
1154  } else { lowRigPID(rig,dedx,charge);}
1155 
1156  } else if (dedx==0.0){ fillAsUnknown();}
1157 
1158  // do not do deuteron or higher
1159  myBandBGFcn->SetParameter(3,1);
1160  myBandBGFcn->SetParameter(4,1.45);
1161  if (dedx>myBandBGFcn->Eval(rig,0,0)) fillAsUnknown();
1162 
1163  } else fillAsUnknown();
1164 
1165  fillPIDHypothis();
1166 
1167 
1168 
1169  }
StHelix::distance
double distance(const StThreeVector< double > &p, bool scanPeriods=true) const
minimal distance between point and helix
Definition: StHelix.cc:240
Generated by   doxygen 1.8.5