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   fNumCorrel(0), fPoolDepth(0), fCorrelType(NULL), fGenQA(kFALSE),
  31   minTriggPt(0), maxTriggPt(0), bwTriggPt(0), minAssocPt(0), maxAssocPt(0), bwAssocPt(0),
  32   fITSRefit(kFALSE), fTPCRefit(kFALSE), fTRDRefit(kFALSE), fRejectKinkChild(kFALSE),
  33   fMaxNsigmaVtx(0), fMaxITSChi2(0), fMaxTPCChi2(0), fMinNClusITS(0), fMinNClusTPC(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(fBinning[centr]) delete [] fBinning[centr];
  44   fNumBins[centr] = 0;
  45   if(fBinning[zvert]) delete [] fBinning[zvert];
  46   fNumBins[zvert] = 0;
  47 }
  48
  49 void AliJetCorrelSelector::SetCorrelTypes(UInt_t s, UInt_t * const v){
  50   // fills the array of correlation types
  51   if(s<1){std::cerr<<"AliJetCorrelSelector::SetCorrelTypes - empty array"<<std::endl; exit(-1);}
  52   if(s>kMAXNUMCORREL){std::cerr<<"AliJetCorrelSelector: correlation array too big!"<<std::endl; exit(-1);}
  53   fNumCorrel = s;
  54   fCorrelType = new UInt_t[fNumCorrel];
  55   for(UInt_t k=0; k<fNumCorrel; k++){
  56     if(v[k]>99){
  57       std::cerr<<"AliJetCorrelSelector::SetCorrelTypes - read error? val["<<k<<"]="<<v[k]<<std::endl;
  58       exit(-1);
  59     }
  60     else fCorrelType[k] = v[k];
  61   }
  62 }
  63
  64 void AliJetCorrelSelector::SetBinningCentr(UInt_t s, Float_t * const v){
  65   // fills array of centrality bins
  66   if(s<1){std::cerr<<"AliJetCorrelSelector::SetBinningCentr - empty array"<<std::endl; exit(-1);}
  67   if(s>kMAXCENTBIN){std::cerr<<"AliJetCorrelSelector: centrality array too big!"<<std::endl; exit(-1);}
  68   fNumBins[centr] = s;
  69   fBinning[centr] = new Float_t[fNumBins[centr]];
  70   for(UInt_t k=0; k<fNumBins[centr]; k++){
  71     if(TMath::Abs(v[k])>999.){
  72       std::cerr<<"AliJetCorrelSelector::SetBinningCentr - read error? val["<<k<<"]="<<v[k]<<std::endl;
  73       exit(-1);
  74     }
  75     else fBinning[centr][k] = v[k];
  76   }
  77 }
  78
  79 void AliJetCorrelSelector::SetBinningZvert(UInt_t s, Float_t * const v){
  80   // fills array of vertex bins
  81   if(s<1){std::cerr<<"AliJetCorrelSelector::SetBinningZvert - empty array"<<std::endl; exit(-1);}
  82   if(s>kMAXVERTBIN){std::cerr<<"AliJetCorrelSelector: vertex array too big!"<<std::endl; exit(-1);}
  83   fNumBins[zvert] = s;
  84   fBinning[zvert] = new Float_t[fNumBins[zvert]];
  85   for(UInt_t k=0; k<fNumBins[zvert]; k++){
  86     if(TMath::Abs(v[k])>999.){
  87       std::cerr<<"AliJetCorrelSelector::SetBinningZvert - read error? val["<<k<<"]="<<v[k]<<std::endl;
  88       exit(-1);
  89     }
  90     else fBinning[zvert][k] = v[k];
  91   }
  92 }
  93
  94 void AliJetCorrelSelector::SetBinningTrigg(Float_t min, Float_t max, Float_t bw){
  95   // sets trigger Pt binning
  96   minTriggPt = min;
  97   maxTriggPt = max;
  98   bwTriggPt  = bw;
  99 }
 100
 101 void AliJetCorrelSelector::SetBinningAssoc(Float_t min, Float_t max, Float_t bw){
 102   // sets associated Pt binning
 103   minAssocPt = min;
 104   maxAssocPt = max;
 105   bwAssocPt  = bw;
 106 }
 107
 108 void AliJetCorrelSelector::Show(){
 109   // print out all user selections
 110   std::cout<<"Generic selections: "<<std::endl<<" PoolDepth="<<fPoolDepth;
 111   std::cout<<std::endl<<" Correlation Types: ";
 112   for(UInt_t k=0; k<fNumCorrel; k++) std::cout<<fCorrelType[k]<<" ";
 113   std::cout<<std::endl<<" Centrality/Multiplicity binning: ";
 114   for(UInt_t k=0; k<fNumBins[centr]; k++) std::cout<<fBinning[centr][k]<<" ";
 115   std::cout<<std::endl<<" Vertex binning: ";
 116   for(UInt_t k=0; k<fNumBins[zvert]; k++) std::cout<<fBinning[zvert][k]<<" ";
 117   std::cout<<std::endl<<" Trigg binning:"<<minTriggPt<<"->"<<maxTriggPt<<"/"<<bwTriggPt;
 118   std::cout<<std::endl<<" Assoc binning:"<<minAssocPt<<"->"<<maxAssocPt<<"/"<<bwAssocPt;
 119   std::cout<<std::endl<<"Track selections: "<<std::endl<<" MaxNsigmaVtx="<<fMaxNsigmaVtx<<std::endl
 120            <<" MaxITSChi2="<<fMaxITSChi2<<" MaxTPCChi2="<<fMaxTPCChi2<<std::endl
 121            <<" MinNClusITS="<<fMinNClusITS<<" MinNClusTPC="<<fMinNClusTPC<<std::endl
 122            <<" ITSRefit="<<fITSRefit<<" TPCRefit="<<fTPCRefit<<" TRDRefit="<<fTRDRefit<<std::endl
 123            <<" RejectKinkChild="<<fRejectKinkChild<<std::endl;
 124 }
 125
 126 Float_t AliJetCorrelSelector::BinBorder(BinType_t cType, UInt_t k){
 127   // returns bin margins
 128   if(k<=NoOfBins(cType)) return fBinning[cType][k];
 129   else {std::cerr<<"BinBorder Error: bin of type "<<cType<<" outside range "<<k<<std::endl;  exit(0);}
 130 }
 131
 132 Int_t AliJetCorrelSelector::GetBin(BinType_t cType, Float_t val){
 133   // returns bin number
 134   Int_t iBin=-1; UInt_t nBins=NoOfBins(cType);
 135   for(UInt_t i=0; i<nBins; i++)
 136     if(BinBorder(cType,i)<=val && val<BinBorder(cType,i+1)) iBin=i;
 137   return iBin;
 138 }
 139
 140 //////////////////////////////////////////////////////////
 141 // Cutting Methods
 142 /////////////////////////////////////////////////////////
 143
 144 Bool_t AliJetCorrelSelector::GoodTrackPair(CorrelTrack_t *t1, CorrelTrack_t *t2){
 145   // meant for two-track cuts (like TPC entrance); but hopes are that single-track cuts,
 146   // like a high no of TPC clusters, will avoid double counting of split tracks...
 147   t1->Show(); t2->Show();
 148   return kTRUE;
 149 }
 150
 151 Bool_t AliJetCorrelSelector::LowQualityTrack(AliESDtrack* track){
 152   // selects low quality tracks
 153   UInt_t status = track->GetStatus();
 154   if(fITSRefit && !(status & AliESDtrack::kITSrefit)) return kTRUE;
 155   if(fTPCRefit && !(status & AliESDtrack::kTPCrefit)) return kTRUE;
 156
 157   UInt_t nClusITS = track->GetITSclusters(0);
 158   UInt_t nClusTPC = track->GetTPCclusters(0); // or track->GetTPCNcls() ?
 159   if(nClusITS<fMinNClusITS) return kTRUE;
 160   if(nClusTPC<fMinNClusTPC) return kTRUE;
 161
 162   Float_t chi2ITS=-1., chi2TPC=-1.;
 163   if(nClusITS!=0) chi2ITS = track->GetITSchi2()/Float_t(nClusITS);
 164   if(nClusTPC!=0) chi2TPC = track->GetTPCchi2()/Float_t(nClusTPC);
 165   if(chi2ITS<0 || chi2ITS>fMaxITSChi2) return kTRUE;
 166   if(chi2TPC<0 || chi2TPC>fMaxTPCChi2) return kTRUE;
 167
 168   if(fRejectKinkChild && track->GetKinkIndex(0)>0) return kTRUE;
 169   if(GetSigmaToVertex(track)>fMaxNsigmaVtx) return kTRUE;
 170
 171   return kFALSE;
 172 }
 173
 174 Bool_t AliJetCorrelSelector::PassPID(AliESDtrack* track, PartType_t PartType){
 175   // checks if a track has the required ID
 176   Bool_t hasReqPID = kFALSE;
 177   Stat_t fPid;
 178   Stat_t fWeight;
 179   GetPID(track, fPid, fWeight);
 180
 181   switch(PartType){
 182     case hadron:
 183       // if(fPID!=0) hasReqPID = kTRUE;
 184       hasReqPID = kTRUE;
 185       break;
 186     case electron:
 187       // if(fTRDRefit && !(status & AliESDtrack::kTRDrefit)) hasReqPID = kFALSE;
 188       // if(fPID!=0) hasReqPID = kFALSE;
 189       hasReqPID = kTRUE;
 190       break;
 191     case pion:
 192       // if(fPID!=2) hasReqPID = kFALSE;
 193       hasReqPID = kTRUE;
 194       break;
 195     case kaon:
 196       // if(fPID!=3) hasReqPID = kFALSE;
 197       hasReqPID = kTRUE;
 198       break;
 199     case proton:
 200       // if(fPID!=4) hasReqPID = kFALSE;
 201       hasReqPID = kTRUE;
 202       break;
 203     default:
 204       std::cerr<<"AliJetCorrelSelector::PassPID() - ERROR: wrong type!"<<std::endl;
 205       exit(-1);
 206   }
 207   return hasReqPID;
 208 }
 209
 210 Float_t AliJetCorrelSelector::GetSigmaToVertex(AliESDtrack* track){
 211   // Calculates the number of sigma to the vertex; from ANALYSIS/AliESDtrackCuts
 212   Float_t b[2], bRes[2], bCov[3];
 213   track->GetImpactParameters(b,bCov);
 214   if(bCov[0]<=0 || bCov[2]<=0) return -1.;
 215   bRes[0] = TMath::Sqrt(bCov[0]);
 216   bRes[1] = TMath::Sqrt(bCov[2]);
 217
 218   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 219
 220   if(TMath::Exp(-d*d/2)<1e-10) return 1000;
 221
 222   Float_t nSigma = TMath::ErfInverse(1-TMath::Exp(-d*d/2))*TMath::Sqrt(2);
 223   return nSigma;
 224 }
 225
 226 void AliJetCorrelSelector::GetPID(AliESDtrack* track, Stat_t& fpid, Stat_t& fweight){
 227   // Finds most probable particle: 0=Electron, 1=Muon, 2=Pion, 3=Kaon, 4=Proton
 228   fpid = -1;
 229   fweight = -1;
 230
 231   Double_t wpart[5], wpartbayes[5];
 232   track->GetESDpid(wpart); // probability of the different particle types
 233   Double_t c[5]={1., 1., 1., 1., 1.}; // Tentative particle type "concentrations"
 234   //  Double_t c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
 235
 236   //Bayes formula
 237   Double_t rcc = 0.;
 238   for(Int_t i=0; i<5; i++) {rcc += c[i]*wpart[i];}
 239   if(TMath::Abs(rcc)<1e-10) return;
 240   for(Int_t i=0; i<5; i++) {wpartbayes[i] = c[i]*wpart[i]/rcc;}
 241
 242   Int_t ipid=-1;
 243   Float_t max=0.;
 244   for(Int_t i=0; i<5; i++) {
 245     if(wpartbayes[i]>max) {
 246       ipid = i;
 247       max = wpartbayes[i];
 248     }
 249   }
 250
 251   fpid = ipid;
 252   fweight = max;
 253 }