StRoot: StEventUtilities/StuProbabilityPidAlgorithm.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: StuProbabilityPidAlgorithm.cxx,v 1.34 2004/09/19 00:07:28 perev Exp $
00004  *
00005  * Author:Aihong Tang, Richard Witt(FORTRAN version). Kent State University
00006  *        Send questions to aihong@cnr.physics.kent.edu 
00007  ***************************************************************************
00008  *
00009  * Description: A functor that do PID base on Probability (Amplitude) info.
00010  *
00011  ***************************************************************************
00012  *
00013  * $Log: StuProbabilityPidAlgorithm.cxx,v $
00014  * Revision 1.34  2004/09/19 00:07:28  perev
00015  * Small Walgrind leak fixed
00016  *
00017  * Revision 1.33  2004/04/09 15:46:16  aihong
00018  * add isPIDTableRead()
00019  *
00020  * Revision 1.32  2003/09/02 17:58:09  perev
00021  * gcc 3.2 updates + WarnOff
00022  *
00023  * Revision 1.31  2003/06/24 02:53:14  aihong
00024  * update for dAu PIDtable
00025  *
00026  * Revision 1.30  2003/05/02 21:32:08  aihong
00027  * destroy myBandBGFcn in destructor
00028  *
00029  * Revision 1.29  2003/04/30 20:37:12  perev
00030  * Warnings cleanup. Modified lines marked VP
00031  *
00032  * Revision 1.28  2002/12/20 20:26:44  aihong
00033  * let it handle P02gd table
00034  *
00035  * Revision 1.27  2002/12/11 15:35:53  aihong
00036  * put in fabs in processPIDAsFuntion()
00037  *
00038  * Revision 1.26  2002/04/19 16:36:41  perev
00039  * bug fix,init to zero
00040  *
00041  * Revision 1.25  2002/01/17 03:25:27  aihong
00042  * add production Tag to take care of different centrality def. between different productions
00043  *
00044  * Revision 1.24  2001/03/21 18:16:13  aihong
00045  * constructor without StEvent added
00046  *
00047  * Revision 1.23  2001/03/21 17:54:30  aihong
00048  * add processPIDAsFunction()
00049  *
00050  * Revision 1.22  2000/12/28 21:00:57  aihong
00051  * remove temp. fix
00052  *
00053  * Revision 1.21  2000/12/26 23:27:05  aihong
00054  * temperory fix for e amp tail
00055  *
00056  * Revision 1.20  2000/12/20 16:55:16  aihong
00057  * let it survive when no support PIDTable is present
00058  *
00059  * Revision 1.19  2000/12/18 23:51:36  aihong
00060  * mExtrap related bug fixed
00061  *
00062  * Revision 1.10  2000/08/16 12:46:07  aihong
00063  * bug killed
00064  *
00065  * Revision 1.9  2000/08/15 23:04:18  aihong
00066  * speed it up by looking up table
00067  *
00068  * Revision 1.8  2000/07/22 22:45:27  aihong
00069  * change include path
00070  *
00071  * Revision 1.7  2000/07/12 16:29:38  aihong
00072  * change to avoid possible name confliction from StarClassLibary
00073  *
00074  * Revision 1.6  2000/05/24 14:35:41  ullrich
00075  * Added 'const' to compile on Sun CC5.
00076  *
00077  * Revision 1.5  2000/05/05 19:25:39  aihong
00078  * modified ctor
00079  *
00080  * Revision 1.2  2000/03/09 20:45:04  aihong
00081  * add head for Log
00082  *
00083  **************************************************************************/
00084 #include <float.h>
00085 #include <Stsstream.h>
00086 #include "Stiostream.h" //this line should be deleted later
00087 #include "TFile.h"
00088 #include "TF1.h"
00089 
00090 #include "StMessMgr.h"
00091 #include "StPhysicalHelixD.hh"
00092 #include "PhysicalConstants.h"
00093 
00094 #include "StuProbabilityPidAlgorithm.h"
00095 
00096 #include "StEventUtilities/BetheBlochFunction.hh"
00097 #include "StEventUtilities/MaxllBoltz.hh"
00098 #include "StEventUtilities/Linear.hh"
00099 
00100 #include "StEventUtilities/StuFtpcRefMult.hh"
00101 #include "StEventUtilities/StuRefMult.hh"
00102 
00103 //TMap::FindObject goes wild!! TMap::GetValue works.
00104 
00105 //-------------------------------
00106 StuProbabilityPidAlgorithm::StuProbabilityPidAlgorithm(StEvent& ev): 
00107                mPionMinusProb(0.),
00108                mElectronProb(0.),
00109                mKaonMinusProb(0.),
00110                mAntiProtonProb(0.),
00111                mPionPlusProb(0.),
00112                mPositronProb(0.),
00113                mKaonPlusProb(0.),
00114                mProtonProb(0.){
00115 
00116 
00117       PID[0]=-1;//should be sth.standard say unIdentified.
00118       PID[1]=-1;     
00119       PID[2]=-1;
00120       PID[3]=-1;
00121 
00122       mProb[0]=0;
00123       mProb[1]=0;
00124       mProb[2]=0;
00125       mProb[3]=0;
00126 
00127      table = StParticleTable::instance();
00128 
00129      mExtrap=false;
00130      mEvent=&ev;
00131 
00132         //init funtions
00133   myBandBGFcn
00134      =new TF1("myBandBGFcn",BetheBlochFunction, 0.,5., 7); 
00135 
00136   Double_t myPars[7]
00137         ={ 1.072, 0.3199, 1.66032e-07, 1, 1, 2.71172e-07, 0.0005 };
00138 
00139    myBandBGFcn->SetParameters(&myPars[0]);
00140 
00141 }
00142 
00143 //-------------------------------
00144 StuProbabilityPidAlgorithm::StuProbabilityPidAlgorithm(): 
00145                mPionMinusProb(0.),
00146                mElectronProb(0.),
00147                mKaonMinusProb(0.),
00148                mAntiProtonProb(0.),
00149                mPionPlusProb(0.),
00150                mPositronProb(0.),
00151                mKaonPlusProb(0.),
00152                mProtonProb(0.){
00153 
00154 
00155       PID[0]=-1;//should be sth.standard say unIdentified.
00156       PID[1]=-1;     
00157       PID[2]=-1;
00158       PID[3]=-1;
00159 
00160       mProb[0]=0;
00161       mProb[1]=0;
00162       mProb[2]=0;
00163       mProb[3]=0;
00164 
00165      table = StParticleTable::instance();
00166 
00167      mExtrap=false;
00168 
00169         //init funtions
00170   myBandBGFcn
00171      =new TF1("myBandBGFcn",BetheBlochFunction, 0.,5., 7); 
00172 
00173   Double_t myPars[7]
00174         ={ 1.072, 0.3199, 1.66032e-07, 1, 1, 2.71172e-07, 0.0005 };
00175 
00176    myBandBGFcn->SetParameters(&myPars[0]);
00177 
00178 }
00179 
00180 //-------------------------------
00181 StuProbabilityPidAlgorithm::~StuProbabilityPidAlgorithm(){
00182   /* no op */
00183   delete myBandBGFcn;
00184 }
00185 //-------------------------------
00186 void StuProbabilityPidAlgorithm::setDedxMethod(StDedxMethod method){
00187       StuProbabilityPidAlgorithm::mDedxMethod=method;
00188 }
00189 
00190 //-------------------------------
00191 bool StuProbabilityPidAlgorithm::isPIDTableRead(){
00192       return StuProbabilityPidAlgorithm::mPIDTableRead;
00193 }
00194 
00195 
00196 //-------------------------------
00197 StParticleDefinition* StuProbabilityPidAlgorithm::mostLikelihoodParticle(){
00198   
00199        
00200       return table->findParticleByGeantId(PID[0]);
00201 }
00202 
00203 //-------------------------------
00204 StParticleDefinition* StuProbabilityPidAlgorithm::secondLikelihoodParticle(){
00205   
00206        
00207       return table->findParticleByGeantId(PID[1]);
00208 }
00209 
00210 //-------------------------------
00211 StParticleDefinition* StuProbabilityPidAlgorithm::thirdLikelihoodParticle(){
00212   
00213        
00214       return table->findParticleByGeantId(PID[2]);
00215 }
00216 //-------------------------------
00217 StParticleDefinition* StuProbabilityPidAlgorithm::getParticle(int i){
00218 
00219    if (i>=0 && i<4){
00220            return table->findParticleByGeantId(PID[i]);
00221    }   else { 
00222 
00223         gMessMgr->Error()<<"StuProbabilityPidAlgorithm::getParticle(int i), i must be 0,1,2,3 only. "<<endm;
00224 
00225      return 0;
00226    }
00227  }
00228 //-------------------------------
00229 double StuProbabilityPidAlgorithm::getProbability(int i){
00230    if (i>=0 && i<4){
00231            return mProb[i];
00232    }   else { 
00233 
00234         gMessMgr->Error()<<"StuProbabilityPidAlgorithm::getProbability(int i), i must be 0,1,2,3 only. "<<endm;
00235 
00236      return 0.0;
00237  }
00238    
00239 } 
00240 
00241 //-------------------------------
00242 double StuProbabilityPidAlgorithm::mostLikelihoodProbability(){
00243            return mProb[0];
00244  }
00245 //-------------------------------
00246 double StuProbabilityPidAlgorithm::secondLikelihoodProbability(){
00247            return mProb[1];
00248  }
00249 //-------------------------------
00250 double StuProbabilityPidAlgorithm::thirdLikelihoodProbability(){
00251            return mProb[2];
00252  }
00253 
00254 //-------------------------------
00255 bool StuProbabilityPidAlgorithm::isExtrap(){
00256            return mExtrap;
00257 }
00258 
00259 
00260 //-------------------------------
00261 StParticleDefinition* StuProbabilityPidAlgorithm::operator() (const StTrack& theTrack, const StSPtrVecTrackPidTraits& traits){
00262 
00263 
00264 
00265       PID[0]=-1;//should be sth.standard say unIdentified.
00266       PID[1]=-1;     
00267       PID[2]=-1;
00268       PID[3]=-1;
00269 
00270      mProb[0]=0;
00271      mProb[1]=0;
00272      mProb[2]=0;
00273      mProb[3]=0;
00274 
00275      mExtrap=false;
00276 
00277 
00278      mPionMinusProb=0.;
00279      mElectronProb=0.;
00280      mKaonMinusProb=0.;
00281      mAntiProtonProb=0.;
00282      mPionPlusProb=0.;
00283      mPositronProb=0.;
00284      mKaonPlusProb=0.;
00285      mProtonProb=0.;
00286 
00287      if (mPIDTableRead) {
00288 
00289           double rig    =0.0;
00290           double dedx   =0.0;
00291           double dca    =0.0; //in units of cm.
00292           int    nhits  =0;
00293           int    charge =0;
00294           double eta    =0.; 
00295           double cent   =0.; // % central
00296 
00297 
00298            StPrimaryVertex* primaryVtx=mEvent->primaryVertex();
00299     const StPhysicalHelixD& helix=theTrack.geometry()->helix();
00300            dca=helix.distance(primaryVtx->position());
00301 
00302   //cent in cross section
00303 
00304   if (mProductionTag){
00305 
00306   if ( (mProductionTag->GetString()).Contains("P01gl")
00307        || (mProductionTag->GetString()).Contains("P02gd") ){
00308         cent = getCentrality(uncorrectedNumberOfNegativePrimaries(*mEvent));  
00309   }  else if ( (mProductionTag->GetString()).Contains("P03ia_dAu") ){
00310         cent = getCentrality(uncorrectedNumberOfFtpcEastPrimaries(*mEvent)); 
00311   }  else {
00312         gMessMgr->Error()<<"Production tag "<<mProductionTag->GetString().Data()<<" in PIDTable is filled but its name is not recognized ! "<<endm;
00313   }
00314 
00315 
00316   } else { //the first PID table has no production tag
00317         cent = getCentrality(uncorrectedNumberOfNegativePrimaries(*mEvent));
00318   }
00319 
00320 
00321    
00322           const StDedxPidTraits* dedxPidTr=0;
00323 
00324 
00325        charge=(theTrack.geometry())->charge();
00326 
00327        for (int itrait = 0; itrait < int(traits.size()); itrait++){
00328 
00329            dedxPidTr = 0;
00330            if (traits[itrait]->detector() == kTpcId) {
00331              //
00332              // tpc pid trait
00333              //
00334               const StTrackPidTraits* thisTrait = traits[itrait];
00335              //
00336              // perform cast to make the pid trait a dedx trait
00337              //
00338               dedxPidTr = dynamic_cast<const StDedxPidTraits*>(thisTrait);
00339              
00340            }
00341            if (dedxPidTr &&  dedxPidTr->method() == mDedxMethod) break;
00342 
00343        }
00344 
00345            if (dedxPidTr) {
00346              dedx=dedxPidTr->mean();
00347              nhits=dedxPidTr->numberOfPoints();
00348            }
00349 
00350 
00351        if (dedx!=0.0 && nhits>=0 //dedx ==0.0 no sense 
00352            && thisPEnd > 0. && thisEtaEnd > 0. // *End ==0, no PIDTable read.
00353            && thisNHitsEnd > 0. ){
00354 
00355     const StThreeVectorF& p=theTrack.geometry()->momentum();
00356     rig=double(p.mag()/charge);
00357 
00358 
00359 
00360     if (mProductionTag){ //for AuAu, +/- eta were folded together when building PID
00361       //table, for dAu, + and - eta were treated differently.
00362 
00363   if ( (mProductionTag->GetString()).Contains("P01gl")
00364        || (mProductionTag->GetString()).Contains("P02gd") ){
00365     eta=fabs(p.pseudoRapidity());
00366   }  else if ( (mProductionTag->GetString()).Contains("P03ia_dAu") ){
00367     eta=p.pseudoRapidity();
00368   }  else {
00369         gMessMgr->Error()<<"Production tag "<<mProductionTag->GetString().Data()<<" in PIDTable is filled but its name is not recognized ! "<<endm;
00370   }
00371 
00372   } else { //the first PID table has no production tag
00373     eta = fabs(p.pseudoRapidity());
00374   }
00375 
00376     rig   = fabs(rig);
00377     dedx  = (dedx>thisDedxStart) ? dedx : thisDedxStart;
00378     rig   = (rig >thisPStart)  ? rig  : thisPStart;
00379     rig   = (rig <thisPEnd  )  ? rig  : thisPEnd*0.9999;   
00380     eta   = (eta  >thisEtaStart)   ? eta   : thisEtaStart;
00381     eta   = (eta  <thisEtaEnd  )   ? eta   : thisEtaEnd*0.9999;
00382     nhits = (nhits > int(thisNHitsStart)) ? nhits : int(thisNHitsStart);
00383     nhits = (nhits < int(thisNHitsEnd) ) ? nhits : int(thisNHitsEnd-1);
00384 
00385     //----------------get all info. I want for a track. now do PID
00386 
00387     setCalibrations(eta, nhits);
00388 
00389    if (dedx<thisDedxEnd){
00390 
00391    fillPIDByLookUpTable(cent, dca, charge,rig, eta, nhits,dedx);
00392 
00393    } else { lowRigPID(rig,dedx,charge);}
00394  
00395        } else if (dedx==0.0){ fillAsUnknown();}
00396 
00397      // do not do deuteron or higher
00398       myBandBGFcn->SetParameter(3,1);
00399       myBandBGFcn->SetParameter(4,1.45);
00400       if (dedx>myBandBGFcn->Eval(rig,0,0)) fillAsUnknown();
00401 
00402      } else fillAsUnknown();
00403 
00404        fillPIDHypothis();
00405 
00406        return table->findParticleByGeantId(PID[0]);
00407 
00408  }
00409 
00410 //-------------------------------
00411 void StuProbabilityPidAlgorithm::lowRigPID(double rig,double dedx,int theCharge){
00412 
00413        double upper;
00414        double lower;
00415        double rigidity=fabs(rig);      
00416        double mdedx=dedx; 
00417       
00418        double fakeMass=0.;
00419 
00420       //pion
00421       fakeMass=0.32075026;
00422       myBandBGFcn->SetParameter(3,1);
00423       myBandBGFcn->SetParameter(4,0.32075026 );
00424         
00425 
00426       lower =0.;
00427        upper =myBandBGFcn->Eval(rigidity,0,0);
00428 
00429        if (mdedx>lower && mdedx<upper){
00430          PID[0]=(theCharge>0.0)? 8 : 9;  //pi+/-
00431          mProb[0]=1.0;
00432          mProb[1]=0.0;
00433          mProb[2]=0.0;
00434          mProb[3]=0.0;
00435        }
00436 
00437            lower = upper;
00438 
00439        //kaon
00440            fakeMass=0.709707;
00441       myBandBGFcn->SetParameter(3,1);
00442       myBandBGFcn->SetParameter(4,fakeMass);
00443         
00444       upper =myBandBGFcn->Eval(rigidity,0,0);
00445 
00446        if (mdedx>lower && mdedx<upper){
00447          PID[0]=(theCharge>0.0)? 11:12;  //k+/-
00448          mProb[0]=1.0;
00449          mProb[1]=0.0;
00450          mProb[2]=0.0;
00451          mProb[3]=0.0;
00452        }
00453 
00454 
00455        lower = upper;
00456 
00457        //proton/pBar
00458            fakeMass=1.45;
00459       myBandBGFcn->SetParameter(3,1);
00460       myBandBGFcn->SetParameter(4,fakeMass);
00461         
00462       upper =myBandBGFcn->Eval(rigidity,0,0);
00463 
00464        if (mdedx>lower && mdedx<upper){
00465          PID[0]=(theCharge>0.0)? 14:15;  //proton/antiproton
00466          mProb[0]=1.0;
00467          mProb[1]=0.0;
00468          mProb[2]=0.0;
00469          mProb[3]=0.0;
00470        } 
00471 
00472        /*                
00473        lower = upper;
00474 
00475        //deuteron
00476            fakeMass=2.4;
00477       myBandBGFcn->SetParameter(3,1);
00478       myBandBGFcn->SetParameter(4,fakeMass);
00479 
00480       upper =myBandBGFcn->Eval(rigidity,0,0);
00481 
00482        if (mdedx>lower && mdedx<upper){
00483          PID[0]=45;  //deuteron
00484          mProb[0]=1.0;
00485          mProb[1]=0.0;
00486          mProb[2]=0.0;
00487          mProb[3]=0.0;
00488        }
00489 
00490        lower = upper;
00491 
00492        //triton
00493        m = -1.5374;  //New slope needed for deuterons and tritons.
00494        a = 1.8121e-5;
00495        upper = a*::pow(rigidity,m);
00496        if (mdedx>lower && mdedx<upper){
00497          PID[0]=46;  //triton
00498          mProb[0]=1.0;
00499          mProb[1]=0.0;
00500          mProb[2]=0.0;
00501          mProb[3]=0.0;
00502        }
00503 
00504        */
00505 
00506 
00507 
00508 }  
00509 
00510 
00511 
00512 //-------------------------------
00513 void StuProbabilityPidAlgorithm::fillAsUnknown(){
00514 
00515   for (int i=0; i<4; i++) {
00516       PID[i]=-1; mProb[i]=-1;
00517   }
00518 }
00519 
00520 
00521 //-------------------------------
00522 void StuProbabilityPidAlgorithm::readParametersFromFile(TString fileName){
00523 
00524 
00525 
00526       TFile f(fileName,"READ");
00527 
00528       if (f.IsOpen()){
00529 
00530         delete StuProbabilityPidAlgorithm::mEAmp     ;
00531         delete StuProbabilityPidAlgorithm::mECenter  ;
00532         delete StuProbabilityPidAlgorithm::mESigma   ;
00533         delete StuProbabilityPidAlgorithm::mPiAmp    ;
00534         delete StuProbabilityPidAlgorithm::mPiCenter ;
00535         delete StuProbabilityPidAlgorithm::mPiSigma  ;
00536         delete StuProbabilityPidAlgorithm::mKAmp     ;
00537         delete StuProbabilityPidAlgorithm::mKCenter  ;
00538         delete StuProbabilityPidAlgorithm::mKSigma   ;
00539         delete StuProbabilityPidAlgorithm::mPAmp     ;
00540         delete StuProbabilityPidAlgorithm::mPCenter  ;
00541         delete StuProbabilityPidAlgorithm::mPSigma   ;
00542         delete StuProbabilityPidAlgorithm::mEqualyDividableRangeStartSet;
00543         delete StuProbabilityPidAlgorithm::mEqualyDividableRangeEndSet;
00544         delete StuProbabilityPidAlgorithm::mEqualyDividableRangeNBinsSet;
00545         delete StuProbabilityPidAlgorithm::mNoEqualyDividableRangeNBinsSet;
00546         delete StuProbabilityPidAlgorithm::mMultiBinEdgeSet ;
00547         delete StuProbabilityPidAlgorithm::mDcaBinEdgeSet   ;
00548         delete StuProbabilityPidAlgorithm::mBBPrePar;
00549         delete StuProbabilityPidAlgorithm::mBBTurnOver;
00550         delete StuProbabilityPidAlgorithm::mBBOffSet;
00551         delete StuProbabilityPidAlgorithm::mBBScale;
00552         delete StuProbabilityPidAlgorithm::mBBSaturate;
00553         delete StuProbabilityPidAlgorithm::mProductionTag;
00554 
00555         StuProbabilityPidAlgorithm::mEAmp     =(TVectorD* )f.Get("eAmp");
00556         StuProbabilityPidAlgorithm::mECenter  =(TVectorD* )f.Get("eCenter");
00557         StuProbabilityPidAlgorithm::mESigma   =(TVectorD* )f.Get("eSigma");
00558 
00559         StuProbabilityPidAlgorithm::mPiAmp    =(TVectorD* )f.Get("piAmp");
00560         StuProbabilityPidAlgorithm::mPiCenter =(TVectorD* )f.Get("piCenter");
00561         StuProbabilityPidAlgorithm::mPiSigma  =(TVectorD* )f.Get("piSigma");
00562 
00563         StuProbabilityPidAlgorithm::mKAmp     =(TVectorD* )f.Get("kAmp");
00564         StuProbabilityPidAlgorithm::mKCenter  =(TVectorD* )f.Get("kCenter");
00565         StuProbabilityPidAlgorithm::mKSigma   =(TVectorD* )f.Get("kSigma");
00566 
00567         StuProbabilityPidAlgorithm::mPAmp     =(TVectorD* )f.Get("pAmp");
00568         StuProbabilityPidAlgorithm::mPCenter  =(TVectorD* )f.Get("pCenter");
00569         StuProbabilityPidAlgorithm::mPSigma   =(TVectorD* )f.Get("pSigma");
00570 
00571         StuProbabilityPidAlgorithm::mEqualyDividableRangeStartSet
00572                         =(TVectorD* )f.Get("EqualyDividableRangeStartSet");
00573         StuProbabilityPidAlgorithm::mEqualyDividableRangeEndSet
00574                         =(TVectorD* )f.Get("EqualyDividableRangeEndSet");
00575         StuProbabilityPidAlgorithm::mEqualyDividableRangeNBinsSet
00576                         =(TVectorD* )f.Get("EqualyDividableRangeNBinsSet");
00577         StuProbabilityPidAlgorithm::mNoEqualyDividableRangeNBinsSet
00578                         =(TVectorD* )f.Get("NoEqualyDividableRangeNBinsSet");
00579 
00580         StuProbabilityPidAlgorithm::mMultiBinEdgeSet 
00581                         =(TVectorD* )f.Get("MultiBinEdgeSet");
00582         StuProbabilityPidAlgorithm::mDcaBinEdgeSet   
00583                         =(TVectorD* )f.Get("DcaBinEdgeSet");
00584 
00585 
00586         StuProbabilityPidAlgorithm::mBBPrePar
00587                           =(TVectorD* )f.Get("BBPrePar");
00588         StuProbabilityPidAlgorithm::mBBTurnOver
00589                           =(TVectorD* )f.Get("BBTurnOver");
00590         StuProbabilityPidAlgorithm::mBBOffSet
00591                           =(TVectorD* )f.Get("BBOffSet");
00592         StuProbabilityPidAlgorithm::mBBScale
00593                           =(TVectorD* )f.Get("BBScale");
00594         StuProbabilityPidAlgorithm::mBBSaturate
00595                           =(TVectorD* )f.Get("BBSaturate");
00596 
00597         StuProbabilityPidAlgorithm::mProductionTag
00598                           =(TObjString* )f.Get("productionTag");
00599 
00600 
00601 
00602         StuProbabilityPidAlgorithm::thisMultBins
00603                         =int((*mNoEqualyDividableRangeNBinsSet)(0));
00604         StuProbabilityPidAlgorithm::thisDcaBins
00605                         =int((*mNoEqualyDividableRangeNBinsSet)(1));
00606         StuProbabilityPidAlgorithm::thisChargeBins
00607                         =int((*mNoEqualyDividableRangeNBinsSet)(2));
00608 
00609 
00610         StuProbabilityPidAlgorithm::thisPBins
00611                          =int((*mEqualyDividableRangeNBinsSet)(1));
00612         StuProbabilityPidAlgorithm::thisEtaBins
00613                          =int((*mEqualyDividableRangeNBinsSet)(2));
00614         StuProbabilityPidAlgorithm::thisNHitsBins
00615                          =int((*mEqualyDividableRangeNBinsSet)(3));
00616 
00617 
00618         StuProbabilityPidAlgorithm::thisDedxStart
00619                          =(*mEqualyDividableRangeStartSet)(0);
00620         StuProbabilityPidAlgorithm::thisPStart
00621                          =(*mEqualyDividableRangeStartSet)(1);
00622         StuProbabilityPidAlgorithm::thisEtaStart
00623                          =(*mEqualyDividableRangeStartSet)(2);
00624         StuProbabilityPidAlgorithm::thisNHitsStart
00625                          =(*mEqualyDividableRangeStartSet)(3);
00626 
00627         StuProbabilityPidAlgorithm::thisDedxEnd
00628                           =(*mEqualyDividableRangeEndSet)(0);
00629         StuProbabilityPidAlgorithm::thisPEnd
00630                           =(*mEqualyDividableRangeEndSet)(1);
00631         StuProbabilityPidAlgorithm::thisEtaEnd
00632                           =(*mEqualyDividableRangeEndSet)(2);
00633         StuProbabilityPidAlgorithm::thisNHitsEnd
00634                           =(*mEqualyDividableRangeEndSet)(3);
00635 
00636 
00637 
00638 
00639         StuProbabilityPidAlgorithm::mPIDTableRead=true;
00640 
00641       } else if (!f.IsOpen()) {
00642    
00643         gMessMgr->Error()<<"Data file "<<fileName<<" open failed "<<endm;
00644         return;
00645       }
00646 
00647 }
00648 
00649 /* uncomment this method when data base is available and stable.
00650 //-------------------------------
00651 void StuProbabilityPidAlgorithm::readParametersFromTable(St_Table* tb){
00652 
00653 
00654   if (mDataTable!=tb) {//database changed.
00655      StuProbabilityPidAlgorithm::refreshParameters(tb);
00656      mDataTable=tb;
00657   }else if (mDataTable==tb) {return;}//database no change. no need refreshing.
00658 
00659 }
00660 */
00661 
00662 /*uncomment this method when data base is available and stable.
00663 //-------------------------------
00664 void StuProbabilityPidAlgorithm::refreshParameters(St_Table* theTable){
00665    int i;
00666 
00667   if (mDataSet.GetEntries()>0) mDataSet.Delete();
00668 
00669   St_tpcDedxPidAmpDb* temp=(St_tpcDedxPidAmpDb *)theTable;
00670   tpcDedxPidAmpDb_st* pars=(tpcDedxPidAmpDb_st *)temp->GetTable();
00671   
00672   StPidAmpChannelInfoOut* theInfoOut=new StPidAmpChannelInfoOut(0,45,0,FLT_MAX);
00673   
00674   TObjArray* theOnlyChannel=new TObjArray();
00675   theOnlyChannel->Add(theInfoOut);
00676 
00677   StuProbabilityPidAlgorithm::readAType(StElectron::instance(),pars->eMeanPar, pars->eAmpPar, pars->eSigPar,pars->gasCalib, theOnlyChannel);
00678 
00679 
00680 
00681      mDataSet.Add(theOnlyChannel);
00682 
00683 
00684 }
00685 
00686 */
00687 
00688 
00689 
00690 //-------------------------------
00691 int StuProbabilityPidAlgorithm::thisMultBins=0;
00692 int StuProbabilityPidAlgorithm::thisDcaBins=0;
00693 int StuProbabilityPidAlgorithm::thisChargeBins=0;
00694 
00695 int StuProbabilityPidAlgorithm::thisPBins=0;
00696 int StuProbabilityPidAlgorithm::thisEtaBins=0;
00697 int StuProbabilityPidAlgorithm::thisNHitsBins=0;
00698 
00699 double StuProbabilityPidAlgorithm::thisDedxStart=0.;
00700 double StuProbabilityPidAlgorithm::thisDedxEnd=0;
00701 double StuProbabilityPidAlgorithm::thisPStart=0;
00702 double StuProbabilityPidAlgorithm::thisPEnd=0;
00703 double StuProbabilityPidAlgorithm::thisEtaStart=0;
00704 double StuProbabilityPidAlgorithm::thisEtaEnd=0;
00705 double StuProbabilityPidAlgorithm::thisNHitsStart=0;
00706 double StuProbabilityPidAlgorithm::thisNHitsEnd=0;
00707 
00708 
00709 bool StuProbabilityPidAlgorithm::mPIDTableRead=false;
00710 
00711 StDedxMethod  StuProbabilityPidAlgorithm::mDedxMethod=kTruncatedMeanId;
00712 
00713 TVectorD* StuProbabilityPidAlgorithm::mEAmp    = new TVectorD();
00714 TVectorD* StuProbabilityPidAlgorithm::mECenter = new TVectorD();
00715 TVectorD* StuProbabilityPidAlgorithm::mESigma  = new TVectorD();
00716 
00717 TVectorD* StuProbabilityPidAlgorithm::mPiAmp    = new TVectorD();
00718 TVectorD* StuProbabilityPidAlgorithm::mPiCenter = new TVectorD();
00719 TVectorD* StuProbabilityPidAlgorithm::mPiSigma  = new TVectorD();
00720 
00721 TVectorD* StuProbabilityPidAlgorithm::mKAmp    = new TVectorD();
00722 TVectorD* StuProbabilityPidAlgorithm::mKCenter = new TVectorD();
00723 TVectorD* StuProbabilityPidAlgorithm::mKSigma  = new TVectorD();
00724 
00725 TVectorD* StuProbabilityPidAlgorithm::mPAmp    = new TVectorD();
00726 TVectorD* StuProbabilityPidAlgorithm::mPCenter = new TVectorD();
00727 TVectorD* StuProbabilityPidAlgorithm::mPSigma  = new TVectorD();
00728 
00729 
00730 TVectorD* StuProbabilityPidAlgorithm::mEqualyDividableRangeStartSet   = new TVectorD();
00731 TVectorD* StuProbabilityPidAlgorithm::mEqualyDividableRangeEndSet     = new TVectorD();
00732 TVectorD* StuProbabilityPidAlgorithm::mEqualyDividableRangeNBinsSet   = new TVectorD();
00733 TVectorD* StuProbabilityPidAlgorithm::mNoEqualyDividableRangeNBinsSet = new TVectorD();
00734 
00735 TVectorD* StuProbabilityPidAlgorithm::mMultiBinEdgeSet = new TVectorD();
00736 TVectorD* StuProbabilityPidAlgorithm::mDcaBinEdgeSet   = new TVectorD();  
00737 
00738 TVectorD* StuProbabilityPidAlgorithm::mBBPrePar   = new TVectorD(); 
00739 TVectorD* StuProbabilityPidAlgorithm::mBBTurnOver = new TVectorD(); 
00740 TVectorD* StuProbabilityPidAlgorithm::mBBOffSet   = new TVectorD();
00741 TVectorD* StuProbabilityPidAlgorithm::mBBScale    = new TVectorD(); 
00742 TVectorD* StuProbabilityPidAlgorithm::mBBSaturate = new TVectorD(); 
00743 
00744 TObjString* StuProbabilityPidAlgorithm::mProductionTag = new TObjString();
00745 
00746 
00747 
00748 
00749 //-------------------------------
00750 //St_Table* StuProbabilityPidAlgorithm::mDataTable=0;
00751 
00752 
00753 
00754 
00755 //-------------------------------
00756 void StuProbabilityPidAlgorithm::fillPIDByLookUpTable(double myCentrality, double myDca, int myCharge, double myRig, double myEta, int myNhits, double myDedx)
00757 {
00758 
00759 
00760   //assume bound has been checked before they enter this func.
00761   int theMultBin=getCentralityBin(myCentrality);
00762   int theDcaBin = (myDca>(*mDcaBinEdgeSet)(1)) ? 1 : 0;
00763   int theChargeBin=(myCharge>0) ? 1 : 0;
00764   int thePBin=int(thisPBins*myRig/(thisPEnd-thisPStart));
00765   int theEtaBin=int(thisEtaBins*(myEta-thisEtaStart)/(thisEtaEnd-thisEtaStart));
00766   int theNHitsBin=int(thisNHitsBins*float(myNhits)/(thisNHitsEnd-thisNHitsStart));
00767 
00768 
00769     int totalEntry
00770       = thisMultBins*thisDcaBins*thisChargeBins*thisPBins*thisEtaBins*thisNHitsBins;
00771     
00772     int positionPointer=0;
00773     
00774     totalEntry=totalEntry/thisMultBins;
00775     positionPointer=positionPointer+totalEntry*theMultBin;
00776     
00777     totalEntry=totalEntry/thisDcaBins;
00778     positionPointer=positionPointer+totalEntry*theDcaBin;
00779     
00780     totalEntry=totalEntry/thisChargeBins;
00781     positionPointer=positionPointer+totalEntry*theChargeBin;
00782     
00783     totalEntry=totalEntry/thisPBins;
00784     positionPointer=positionPointer+totalEntry*thePBin;
00785     
00786     totalEntry=totalEntry/thisEtaBins;
00787     positionPointer=positionPointer+totalEntry*theEtaBin;
00788     
00789     totalEntry=totalEntry/thisNHitsBins;
00790     positionPointer=positionPointer+totalEntry*theNHitsBin;
00791 
00792 
00793 
00794     //now calculate prob.
00795     TF1 eGaus("eGaus","gaus",thisDedxStart,thisDedxStart);
00796     TF1 piGaus("piGaus","gaus",thisDedxStart,thisDedxStart);
00797     TF1 kGaus("kGaus","gaus",thisDedxStart,thisDedxStart);
00798     TF1 pGaus("pGaus","gaus",thisDedxStart,thisDedxStart);
00799 
00800     
00801     eGaus.SetParameter(0, (*mEAmp)(positionPointer));
00802     eGaus.SetParameter(1, (*mECenter)(positionPointer));
00803     eGaus.SetParameter(2, (*mESigma)(positionPointer));
00804 
00805     piGaus.SetParameter(0, (*mPiAmp)(positionPointer));
00806     piGaus.SetParameter(1, (*mPiCenter)(positionPointer));
00807     piGaus.SetParameter(2, (*mPiSigma)(positionPointer));
00808 
00809     kGaus.SetParameter(0, (*mKAmp)(positionPointer));
00810     kGaus.SetParameter(1, (*mKCenter)(positionPointer));
00811     kGaus.SetParameter(2, (*mKSigma)(positionPointer));
00812 
00813     pGaus.SetParameter(0, (*mPAmp)(positionPointer));
00814     pGaus.SetParameter(1, (*mPCenter)(positionPointer));
00815     pGaus.SetParameter(2, (*mPSigma)(positionPointer));
00816 
00817 
00818     double eContribution=0;
00819    if (mEAmp && mECenter && mESigma) eContribution = eGaus.Eval(myDedx,0.,0.);
00820    //   if (fabs(myRig)>0.95) eContribution=0; //will deal with it later.
00821 
00822 
00823     double piContribution=0;
00824    if (mPiAmp && mPiCenter && mPiSigma) piContribution = piGaus.Eval(myDedx,0.,0.);
00825 
00826     double kContribution=0;
00827    if (mKAmp && mKCenter && mKSigma) kContribution = kGaus.Eval(myDedx,0.,0.);
00828 
00829     double pContribution=0;
00830    if (mPAmp && mPCenter && mPSigma) pContribution = pGaus.Eval(myDedx,0.,0.);
00831 
00832 
00833    double total = eContribution+piContribution+kContribution+pContribution;
00834 
00835    double eProb=0; double piProb=0; double kProb=0; double pProb=0;
00836 
00837    if (total>0) {
00838      eProb=eContribution/total;
00839      piProb=piContribution/total;
00840      kProb=kContribution/total;
00841      pProb=pContribution/total;
00842    }
00843 
00844 
00845             fill(eProb,  int((myCharge>0) ? 2 : 3   ));
00846             fill(piProb, int((myCharge>0) ? 8 : 9   ));
00847             fill(kProb,  int((myCharge>0) ? 11 : 12 ));
00848             fill(pProb,  int((myCharge>0) ? 14 : 15 ));
00849 
00850             int nn=-1;
00851             float PathHeight=1.0e-7;
00852             float halfHeight=(11-2.0)*PathHeight/2.0;
00853 
00854      if (fabs(myRig) > 0.3) {
00855       if ((myDedx<(((*mPiCenter)(positionPointer))+halfHeight-nn*PathHeight)) 
00856  && (myDedx > ( ((*mKCenter)(positionPointer))-halfHeight+nn*PathHeight) )) 
00857         mExtrap=true;
00858      }
00859 }
00860 //-------------------------------
00861 void StuProbabilityPidAlgorithm::fillPIDHypothis(){
00862   
00863   for (int i=0; i<4; i++){
00864 
00865     switch (PID[i]) {
00866 
00867     case 2 : mPositronProb=mProb[i];
00868                break;
00869     case 3 : mElectronProb=mProb[i];
00870                break;
00871     case 8 : mPionPlusProb=mProb[i];
00872                break;
00873     case 9 : mPionMinusProb=mProb[i];
00874                break;
00875     case 11 : mKaonPlusProb=mProb[i];
00876                break;
00877     case 12 : mKaonMinusProb=mProb[i];
00878                break;
00879     case 14 : mProtonProb=mProb[i];
00880                break;
00881     case 15 : mAntiProtonProb=mProb[i];
00882                break;
00883     }
00884   }
00885 
00886 }
00887 //-------------------------------
00888 int StuProbabilityPidAlgorithm::getCentralityBin(double theCent){
00889   int theBin=0; //input % centrality. output the corres. bin.
00890   
00891   for (int mm = (thisMultBins-1); mm>0; mm--) {
00892     if (theCent< (*mMultiBinEdgeSet)(mm) ) theBin=mm;
00893   }
00894 
00895   return theBin;
00896 
00897 }
00898 
00899 //-------------------------------
00900 
00901 
00902 double StuProbabilityPidAlgorithm::getCentrality(int theMult){
00903 
00904 
00905   if (mProductionTag){
00906 
00907 
00908   if ( (mProductionTag->GetString()).Contains("P01gl")
00909        || (mProductionTag->GetString()).Contains("P02gd") ){
00910   return  getCentrality_P01gl(theMult);
00911   } else if ((mProductionTag->GetString()).Contains("P03ia_dAu")){
00912   return  getCentrality_P03ia_dAu(theMult);
00913   }  else {
00914   gMessMgr->Error()<<"Production tag "<<mProductionTag->GetString().Data()<<" in PIDTable is filled but its name is not recognized ! "<<endm;
00915   }
00916 
00917 
00918   }  else {//the first PID table does not have a productionTag
00919 
00920      // limits 
00921     // For Cut Set 1       
00922     // 225 ~ 3%            
00923     // 215 ~ 5%            
00924     // 200 ~ 7%            
00925     // 180 ~ 10%           
00926     // 140 ~ 18%           
00927     // 130 ~ 20%           
00928     // 120 ~ 23%
00929     // 115 ~ 24%
00930     // 100 ~ 28%
00931 
00932     //  highMedLimit = 180;  // 10%
00933     //  medLowLimit = 108;   //~26%
00934      
00935   if (theMult > 225 ) return 0.03;
00936   else if (theMult > 215 ) return 0.05;
00937   else if (theMult > 200 ) return 0.07;
00938   else if (theMult > 180 ) return 0.10;
00939   else if (theMult > 140 ) return 0.18;
00940   else if (theMult > 130 ) return 0.20; 
00941   else if (theMult > 120 ) return 0.23;
00942   else if (theMult > 115 ) return 0.24;
00943   else if (theMult > 100 ) return 0.28;
00944   else return 0.99;
00945 
00946   }
00947   return 0.99;
00948 
00949 
00950 }
00951 //-------------------------------
00952 
00953 
00954 double StuProbabilityPidAlgorithm::getCentrality_P01gl(int theMult){
00955      
00956   //from Zhangbu's study     
00957   //5% 474 
00958   //10% 409 
00959   //20% 302 
00960   //30% 216 
00961   //40% 149 
00962   //50% 98 
00963   //60% 59 
00964   //70% 32 
00965   //80% 14 
00966 
00967   if (theMult*2. > 474 )      return 0.03; //use Nch instead of hminus. so *2.
00968   else if (theMult*2. > 409 ) return 0.10;
00969   else if (theMult*2. > 302 ) return 0.20;
00970   else if (theMult*2. > 216 ) return 0.30;
00971   else if (theMult*2. > 149 ) return 0.40;
00972   else if (theMult*2. > 98  ) return 0.50; 
00973   else if (theMult*2. > 59  ) return 0.60;
00974   else if (theMult*2. > 32  ) return 0.70;
00975   else if (theMult*2. > 14  ) return 0.80;
00976   else return 0.99;
00977 }
00978 
00979 //-------------------------------
00980 
00981 
00982 double StuProbabilityPidAlgorithm::getCentrality_P03ia_dAu(int theMult){
00983      
00984   //from Joern's study    
00985   // * centrality bin    multiplicity FTPC east    percentOfEvents
00986   // *        1                   <=11                 100-40
00987   // *        2                   <=17                 40-20
00988   // *        3                   >=18                 20-0
00989 
00990   if (theMult >= 18 )      return 0.19; //do not need to be exact. 
00991   else if (theMult >= 12 ) return 0.39; //just for gettting bin # correctly.
00992   else if (theMult > 0 )   return 0.8;
00993   else return 0.99;
00994 }
00995 
00996 //-------------------------------
00997 void StuProbabilityPidAlgorithm::fill(double prob, int geantId){
00998 
00999               if (prob>mProb[0]) {
01000          mProb[3]=mProb[2];
01001          mProb[2]=mProb[1];
01002          mProb[1]=mProb[0];
01003          mProb[0]=prob;
01004 
01005          PID[3]=PID[2];    
01006          PID[2]=PID[1];
01007          PID[1]=PID[0];
01008          PID[0]=geantId;
01009             }
01010             else if (prob>mProb[1]){
01011           mProb[3]=mProb[2];
01012           mProb[2]=mProb[1];
01013           mProb[1]=prob;
01014 
01015           PID[3]=PID[2];
01016           PID[2]=PID[1];
01017           PID[1] =geantId;
01018             }
01019             else if (prob>mProb[2]){
01020             mProb[3]=mProb[2];
01021             mProb[2]=prob;
01022             PID[3]=PID[2];
01023             PID[2]=geantId;
01024             }
01025               else if (prob>=mProb[3]){
01026                 mProb[3]=prob;
01027                 PID[3]=geantId;
01028               }
01029 
01030 }  
01031 //------------------------------------------------
01032 int StuProbabilityPidAlgorithm::getCalibPosition(double theEta, int theNHits){
01033 
01034   int theEtaBin=int(thisEtaBins*(theEta-thisEtaStart)/(thisEtaEnd-thisEtaStart));
01035   int theNHitsBin=int(thisNHitsBins*float(theNHits)/(thisNHitsEnd-thisNHitsStart));
01036 
01037     int totalEntry
01038       = thisEtaBins*thisNHitsBins;
01039     
01040     int positionPointer=0;
01041     
01042     totalEntry=totalEntry/thisEtaBins;
01043     positionPointer=positionPointer+totalEntry*theEtaBin;
01044     
01045     totalEntry=totalEntry/thisNHitsBins;
01046     positionPointer=positionPointer+totalEntry*theNHitsBin;
01047 
01048 
01049 
01050     return positionPointer;
01051 }
01052 //--------------------------------------------
01053 void StuProbabilityPidAlgorithm::setCalibrations(double theEta, int theNhits){
01054    
01055   int thePosition=getCalibPosition(theEta,theNhits);
01056 
01057    myBandBGFcn->SetParameter(5,(*mBBScale)(thePosition));
01058    myBandBGFcn->SetParameter(2,(*mBBOffSet)(thePosition));
01059 
01060 
01061  if (mProductionTag){
01062 
01063    if ( (mProductionTag->GetString()).Contains("P01gl") 
01064      || (mProductionTag->GetString()).Contains("P02gd")){
01065      
01066    myBandBGFcn->SetParameter(0,(*mBBPrePar)(thePosition));
01067    myBandBGFcn->SetParameter(1,(*mBBTurnOver)(thePosition));
01068    myBandBGFcn->SetParameter(6,(*mBBSaturate)(thePosition));
01069 
01070    }
01071  }
01072 
01073 
01074 }
01075 //-------------------------------------cent,dca,charge,rig,eta,nhits,dedx
01076 void StuProbabilityPidAlgorithm::processPIDAsFunction (double theCent, double theDca, int theCharge, double theRig, double theEta, int theNhits, double theDedx){
01077 
01078 
01079 
01080       PID[0]=-1;//should be sth.standard say unIdentified.
01081       PID[1]=-1;     
01082       PID[2]=-1;
01083       PID[3]=-1;
01084 
01085      mProb[0]=0;
01086      mProb[1]=0;
01087      mProb[2]=0;
01088      mProb[3]=0;
01089 
01090      mExtrap=false;
01091 
01092 
01093      mPionMinusProb=0.;
01094      mElectronProb=0.;
01095      mKaonMinusProb=0.;
01096      mAntiProtonProb=0.;
01097      mPionPlusProb=0.;
01098      mPositronProb=0.;
01099      mKaonPlusProb=0.;
01100      mProtonProb=0.;
01101 
01102      if (mPIDTableRead) {
01103 
01104           double rig    =fabs(theRig);
01105           double dedx   =theDedx;
01106           double dca    =theDca; //in units of cm.
01107           int    nhits  =theNhits;
01108           int    charge =theCharge;
01109           double eta    =0.; 
01110           double cent   =theCent; // % central
01111 
01112 
01113     if (mProductionTag){ //for AuAu, +/- eta were folded together when building PID
01114       //table, for dAu, + and - eta were treated differently.
01115 
01116   if ( (mProductionTag->GetString()).Contains("P01gl")
01117        || (mProductionTag->GetString()).Contains("P02gd") ){
01118     eta=fabs(theEta);
01119   }  else if ( (mProductionTag->GetString()).Contains("P03ia_dAu") ){
01120     eta=theEta;
01121   }  else {
01122         gMessMgr->Error()<<"Production tag "<<mProductionTag->GetString().Data()<<" in PIDTable is filled but its name is not recognized ! "<<endm;
01123   }
01124 
01125   } else { //the first PID table has no production tag
01126     eta = fabs(theEta);
01127   }
01128 
01129 
01130 
01131 
01132 
01133        if (dedx!=0.0 && nhits>=0 //dedx ==0.0 no sense 
01134            && thisPEnd > 0. && thisEtaEnd > 0. // *End ==0, no PIDTable read.
01135            && thisNHitsEnd > 0. ){
01136 
01137     rig   = fabs(rig);
01138     dedx  = (dedx>thisDedxStart) ? dedx : thisDedxStart;
01139     rig   = (rig >thisPStart)  ? rig  : thisPStart;
01140     rig   = (rig <thisPEnd  )  ? rig  : thisPEnd*0.9999;   
01141     eta   = (eta  >thisEtaStart)   ? eta   : thisEtaStart;
01142     eta   = (eta  <thisEtaEnd  )   ? eta   : thisEtaEnd*0.9999;
01143     nhits = (nhits > int(thisNHitsStart)) ? nhits : int(thisNHitsStart);
01144     nhits = (nhits < int(thisNHitsEnd) ) ? nhits : int(thisNHitsEnd-1);
01145 
01146     //----------------get all info. I want for a track. now do PID
01147 
01148     setCalibrations(eta, nhits);
01149 
01150    if (dedx<thisDedxEnd){
01151 
01152    fillPIDByLookUpTable(cent, dca, charge,rig, eta, nhits,dedx);
01153 
01154    } else { lowRigPID(rig,dedx,charge);}
01155  
01156        } else if (dedx==0.0){ fillAsUnknown();}
01157 
01158      // do not do deuteron or higher
01159       myBandBGFcn->SetParameter(3,1);
01160       myBandBGFcn->SetParameter(4,1.45);
01161       if (dedx>myBandBGFcn->Eval(rig,0,0)) fillAsUnknown();
01162 
01163      } else fillAsUnknown();
01164 
01165        fillPIDHypothis();
01166 
01167 
01168 
01169  }