PWG4/JetCorrel/AliJetCorrelSelector.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id: $ */
  16
  17 //__________________________________________
  18 // Class for user selections. Object is created by ConfigJetCorrel.C macro
  19 // and passed to AddTaskJetCorrel.C running macro.
  20 //-- Author: Paul Constantin
  21
  22 #include "AliJetCorrelSelector.h"
  23
  24 using namespace std;
  25 using namespace JetCorrelHD;
  26
  27 ClassImp(AliJetCorrelSelector)
  28
  29 AliJetCorrelSelector::AliJetCorrelSelector() :
  30   fGenQA(kFALSE), fNumCorrel(0), nEvtTriggs(0), fPoolDepth(0), fCorrelType(NULL), fEvtTriggs(NULL),
  31   minTriggPt(0), maxTriggPt(0), bwTriggPt(0), minAssocPt(0), maxAssocPt(0), bwAssocPt(0),
  32   fITSRefit(kFALSE), fTPCRefit(kFALSE), fTRDRefit(kFALSE), fRejectKinkChild(kFALSE), fMaxEta(0),
  33   fMaxNsigmaVtx(0), fMaxTrkVtx(0), fMaxITSChi2(0), fMaxTPCChi2(0), fMinNClusITS(0), fMinNClusTPC(0), trkMinProx(0) {
  34   // (default) constructor
  35   fNumBins[centr] = 0; fBinning[centr] = NULL;
  36   fNumBins[zvert] = 0; fBinning[zvert] = NULL;
  37 }
  38
  39 AliJetCorrelSelector::~AliJetCorrelSelector(){
  40   // destructor
  41   if(fCorrelType) delete [] fCorrelType;
  42   fNumCorrel = 0;
  43   if(fEvtTriggs) delete [] fEvtTriggs;
  44   nEvtTriggs = 0;
  45   if(fBinning[centr]) delete [] fBinning[centr];
  46   fNumBins[centr] = 0;
  47   if(fBinning[zvert]) delete [] fBinning[zvert];
  48   fNumBins[zvert] = 0;
  49 }
  50
  51 void AliJetCorrelSelector::SetCorrelTypes(UInt_t s, UInt_t * const v){
  52   // fills the array of correlation types
  53   if(s<1){std::cerr<<"AliJetCorrelSelector::SetCorrelTypes - empty array"<<std::endl; exit(-1);}
  54   if(s>kMAXNUMCORREL){std::cerr<<"AliJetCorrelSelector: correlation array too big!"<<std::endl; exit(-1);}
  55   fNumCorrel = s;
  56   fCorrelType = new UInt_t[fNumCorrel];
  57   for(UInt_t k=0; k<fNumCorrel; k++){
  58     if(v[k]>99){
  59       std::cerr<<"AliJetCorrelSelector::SetCorrelTypes - read error? val["<<k<<"]="<<v[k]<<std::endl;
  60       exit(-1);
  61     }
  62     else fCorrelType[k] = v[k];
  63   }
  64 }
  65
  66 void AliJetCorrelSelector::SetTriggers(UInt_t s, TString * const v){
  67   // fills the array of event triggers
  68   if(s<1){std::cerr<<"AliJetCorrelSelector::SetTriggers - empty array"<<std::endl; exit(-1);}
  69   if(s>9){std::cerr<<"AliJetCorrelSelector: event trigger array too big!"<<std::endl; exit(-1);}
  70   nEvtTriggs = s;
  71   fEvtTriggs = new TString[nEvtTriggs];
  72   for(UInt_t k=0; k<nEvtTriggs; k++){
  73     if(!v[k].IsAscii()){
  74       std::cerr<<"AliJetCorrelSelector::SetTriggers - read error? val["<<k<<"]="<<v[k]<<std::endl;
  75       exit(-1);
  76     }
  77     else fEvtTriggs[k] = ToUpper(v[k]);
  78   }
  79 }
  80
  81
  82 void AliJetCorrelSelector::SetBinningCentr(UInt_t s, Float_t * const v){
  83   // fills array of centrality bins
  84   if(s<1){std::cerr<<"AliJetCorrelSelector::SetBinningCentr - empty array"<<std::endl; exit(-1);}
  85   if(s>kMAXCENTBIN){std::cerr<<"AliJetCorrelSelector: centrality array too big!"<<std::endl; exit(-1);}
  86   fNumBins[centr] = s;
  87   fBinning[centr] = new Float_t[fNumBins[centr]];
  88   for(UInt_t k=0; k<fNumBins[centr]; k++){
  89     if(TMath::Abs(v[k])>999.){
  90       std::cerr<<"AliJetCorrelSelector::SetBinningCentr - read error? val["<<k<<"]="<<v[k]<<std::endl;
  91       exit(-1);
  92     }
  93     else fBinning[centr][k] = v[k];
  94   }
  95 }
  96
  97 void AliJetCorrelSelector::SetBinningZvert(UInt_t s, Float_t * const v){
  98   // fills array of vertex bins
  99   if(s<1){std::cerr<<"AliJetCorrelSelector::SetBinningZvert - empty array"<<std::endl; exit(-1);}
 100   if(s>kMAXVERTBIN){std::cerr<<"AliJetCorrelSelector: vertex array too big!"<<std::endl; exit(-1);}
 101   fNumBins[zvert] = s;
 102   fBinning[zvert] = new Float_t[fNumBins[zvert]];
 103   for(UInt_t k=0; k<fNumBins[zvert]; k++){
 104     if(TMath::Abs(v[k])>999.){
 105       std::cerr<<"AliJetCorrelSelector::SetBinningZvert - read error? val["<<k<<"]="<<v[k]<<std::endl;
 106       exit(-1);
 107     }
 108     else fBinning[zvert][k] = v[k];
 109   }
 110 }
 111
 112 void AliJetCorrelSelector::SetBinningTrigg(Float_t min, Float_t max, Float_t bw){
 113   // sets trigger Pt binning
 114   minTriggPt = min;
 115   maxTriggPt = max;
 116   bwTriggPt  = bw;
 117 }
 118
 119 void AliJetCorrelSelector::SetBinningAssoc(Float_t min, Float_t max, Float_t bw){
 120   // sets associated Pt binning
 121   minAssocPt = min;
 122   maxAssocPt = max;
 123   bwAssocPt  = bw;
 124 }
 125
 126 void AliJetCorrelSelector::Show(){
 127   // print out all user selections
 128   std::cout<<"Generic selections: "<<std::endl<<" PoolDepth="<<fPoolDepth;
 129   std::cout<<std::endl<<" Correlation Types: ";
 130   for(UInt_t k=0; k<fNumCorrel; k++) std::cout<<fCorrelType[k]<<" ";
 131   std::cout<<std::endl<<" Event Triggers: ";
 132   for(UInt_t k=0; k<nEvtTriggs; k++) std::cout<<fEvtTriggs[k]<<" ";
 133   std::cout<<std::endl<<" Centrality/Multiplicity binning: ";
 134   for(UInt_t k=0; k<fNumBins[centr]; k++) std::cout<<fBinning[centr][k]<<" ";
 135   std::cout<<std::endl<<" Vertex binning: ";
 136   for(UInt_t k=0; k<fNumBins[zvert]; k++) std::cout<<fBinning[zvert][k]<<" ";
 137   std::cout<<std::endl<<" Trigg binning:"<<minTriggPt<<"->"<<maxTriggPt<<"/"<<bwTriggPt;
 138   std::cout<<std::endl<<" Assoc binning:"<<minAssocPt<<"->"<<maxAssocPt<<"/"<<bwAssocPt;
 139   std::cout<<std::endl<<"Track selections: "<<std::endl
 140            <<" MaxEta="<<fMaxEta<<" MaxTrkVtx="<<fMaxTrkVtx<<" MaxNsigmaVtx="<<fMaxNsigmaVtx<<std::endl
 141            <<" MaxITSChi2="<<fMaxITSChi2<<" MaxTPCChi2="<<fMaxTPCChi2<<std::endl
 142            <<" MinNClusITS="<<fMinNClusITS<<" MinNClusTPC="<<fMinNClusTPC<<std::endl
 143            <<" ITSRefit="<<fITSRefit<<" TPCRefit="<<fTPCRefit<<" TRDRefit="<<fTRDRefit<<std::endl
 144            <<" RejectKinkChild="<<fRejectKinkChild<<" minTrackPairTPCDist="<<trkMinProx<<std::endl;
 145 }
 146
 147 Float_t AliJetCorrelSelector::BinBorder(BinType_t cType, UInt_t k){
 148   // returns bin margins
 149   if(k<=NoOfBins(cType)) return fBinning[cType][k];
 150   else {std::cerr<<"BinBorder Error: bin of type "<<cType<<" outside range "<<k<<std::endl; exit(0);}
 151 }
 152
 153 Int_t AliJetCorrelSelector::GetBin(BinType_t cType, Float_t val){
 154   // returns bin number
 155   Int_t iBin=-1; UInt_t nBins=NoOfBins(cType);
 156   for(UInt_t i=0; i<nBins; i++)
 157     if(BinBorder(cType,i)<=val && val<BinBorder(cType,i+1)) iBin=i;
 158   return iBin;
 159 }
 160
 161 //////////////////////////////////////////////////////////
 162 // Cutting Methods
 163 /////////////////////////////////////////////////////////
 164
 165 Bool_t AliJetCorrelSelector::SelectedEvtTrigger(AliVEvent *fEVT){
 166   // matches the event trigger classes with the user trigger classes
 167   if(fEVT->InheritsFrom("AliESDEvent")){
 168     const AliESDEvent *esd = (AliESDEvent*)fEVT;
 169     TString trigClass = esd->GetFiredTriggerClasses();
 170     if(nEvtTriggs==1 && fEvtTriggs[0].Contains("ALL")) return kTRUE;
 171     for(UInt_t k=0; k<nEvtTriggs; k++)
 172       if(trigClass.Contains(fEvtTriggs[k])) return kTRUE;
 173     return kFALSE;
 174   } else {std::cerr<<"AliJetCorrelSelector::SelectedEvtTrigger ERROR: not an ESD event!"<<std::endl; exit(0);}
 175 }
 176
 177 Bool_t AliJetCorrelSelector::CloseTrackPair(Float_t dist){
 178   // applies two-track cut (dist at TPC entrance); it is possible that single-track cuts,
 179   // like fraction of shared TPC clusters, will avoid inclusion of split tracks...
 180   if(dist>trkMinProx) return kFALSE;
 181   return kTRUE;
 182 }
 183
 184 Bool_t AliJetCorrelSelector::LowQualityTrack(AliESDtrack* track){
 185   // selects low quality tracks
 186   if(track->Eta()>fMaxEta) return kTRUE;
 187   UInt_t status = track->GetStatus();
 188   if(fITSRefit && !(status & AliESDtrack::kITSrefit)) return kTRUE;
 189   if(fTPCRefit && !(status & AliESDtrack::kTPCrefit)) return kTRUE;
 190
 191 //   UInt_t nClusITS = track->GetITSclusters(0);
 192 //   if(nClusITS<fMinNClusITS) return kTRUE;
 193 //   Float_t chi2ITS=-1.;
 194 //   if(nClusITS!=0) chi2ITS = track->GetITSchi2()/Float_t(nClusITS);
 195 //   if(chi2ITS<0 || chi2ITS>fMaxITSChi2) return kTRUE;
 196
 197   UInt_t nClusTPC = track->GetTPCclusters(0); // or track->GetTPCNcls() ?
 198   if(nClusTPC<fMinNClusTPC) return kTRUE;
 199   Float_t chi2TPC=-1.;
 200   if(nClusTPC!=0) chi2TPC = track->GetTPCchi2()/Float_t(nClusTPC);
 201   if(chi2TPC<0 || chi2TPC>fMaxTPCChi2) return kTRUE;
 202
 203   if(fRejectKinkChild && track->GetKinkIndex(0)>0) return kTRUE;
 204 //  Float_t sigTrkVtx = GetSigmaToVertex(track);
 205 //  if(sigTrkVtx<0 || sigTrkVtx>fMaxNsigmaVtx) return kTRUE;
 206   // instead of track-vertex DCA sigma cut, apply value-cut:
 207   Float_t b[2], bCov[3];
 208   track->GetImpactParameters(b,bCov);
 209   if((b[0]*b[0]+b[1]*b[1])>(fMaxTrkVtx*fMaxTrkVtx)) return kTRUE;
 210
 211   return kFALSE;
 212 }
 213
 214 Bool_t AliJetCorrelSelector::PassPID(AliESDtrack* track, PartType_t PartType){
 215   // checks if a track has the required ID
 216   Bool_t hasReqPID = kFALSE;
 217   Stat_t fPid;
 218   Stat_t fWeight;
 219   GetPID(track, fPid, fWeight);
 220
 221   switch(PartType){
 222     case hadron:
 223       // if(fPID!=0) hasReqPID = kTRUE;
 224       hasReqPID = kTRUE;
 225       break;
 226     case electron:
 227       // if(fTRDRefit && !(status & AliESDtrack::kTRDrefit)) hasReqPID = kFALSE;
 228       // if(fPID!=0) hasReqPID = kFALSE;
 229       hasReqPID = kTRUE;
 230       break;
 231     case pion:
 232       // if(fPID!=2) hasReqPID = kFALSE;
 233       hasReqPID = kTRUE;
 234       break;
 235     case kaon:
 236       // if(fPID!=3) hasReqPID = kFALSE;
 237       hasReqPID = kTRUE;
 238       break;
 239     case proton:
 240       // if(fPID!=4) hasReqPID = kFALSE;
 241       hasReqPID = kTRUE;
 242       break;
 243     default:
 244       std::cerr<<"AliJetCorrelSelector::PassPID() - ERROR: wrong type!"<<std::endl;
 245       exit(-1);
 246   }
 247   return hasReqPID;
 248 }
 249
 250 Float_t AliJetCorrelSelector::GetSigmaToVertex(AliESDtrack* track){
 251   // Calculates the number of sigma to the vertex; from ANALYSIS/AliESDtrackCuts
 252   Float_t b[2], bRes[2], bCov[3];
 253   track->GetImpactParameters(b,bCov);
 254   if(bCov[0]<=0 || bCov[2]<=0) return -1.;
 255   bRes[0] = TMath::Sqrt(bCov[0]);
 256   bRes[1] = TMath::Sqrt(bCov[2]);
 257
 258   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 259
 260   if(TMath::Exp(-d*d/2)<1e-10) return 1000;
 261
 262   Float_t nSigma = TMath::ErfInverse(1-TMath::Exp(-d*d/2))*TMath::Sqrt(2);
 263   return nSigma;
 264 }
 265
 266 void AliJetCorrelSelector::GetPID(AliESDtrack* track, Stat_t& fpid, Stat_t& fweight){
 267   // Finds most probable particle: 0=Electron, 1=Muon, 2=Pion, 3=Kaon, 4=Proton
 268   fpid = -1;
 269   fweight = -1;
 270
 271   Double_t wpart[5], wpartbayes[5];
 272   track->GetESDpid(wpart); // probability of the different particle types
 273   Double_t c[5]={1., 1., 1., 1., 1.}; // Tentative particle type "concentrations"
 274   //  Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
 275
 276   //Bayes formula
 277   Double_t rcc = 0.;
 278   for(Int_t i=0; i<5; i++) {rcc += c[i]*wpart[i];}
 279   if(TMath::Abs(rcc)<1e-10) return;
 280   for(Int_t i=0; i<5; i++) {wpartbayes[i] = c[i]*wpart[i]/rcc;}
 281
 282   Int_t ipid=-1;
 283   Float_t max=0.;
 284   for(Int_t i=0; i<5; i++) {
 285     if(wpartbayes[i]>max) {
 286       ipid = i;
 287       max = wpartbayes[i];
 288     }
 289   }
 290
 291   fpid = ipid;
 292   fweight = max;
 293 }