PWG2/SPECTRA/AliProtonAnalysisBase.cxx

   1 /**************************************************************************
   2  * Author: Panos Christakoglou.                                           *
   3  * Contributors are mentioned in the code where appropriate.              *
   4  *                                                                        *
   5  * Permission to use, copy, modify and distribute this software and its   *
   6  * documentation strictly for non-commercial purposes is hereby granted   *
   7  * without fee, provided that the above copyright notice appears in all   *
   8  * copies and that both the copyright notice and this permission notice   *
   9  * appear in the supporting documentation. The authors make no claims     *
  10  * about the suitability of this software for any purpose. It is          *
  11  * provided "as is" without express or implied warranty.                  *
  12  **************************************************************************/
  13
  14 /* $Id: AliProtonAnalysisBase.cxx 31056 2009-02-16 14:31:41Z pchrist $ */
  15
  16 //-----------------------------------------------------------------
  17 //                 AliProtonAnalysisBase class
  18 //   This is the class to deal with the proton analysis
  19 //   Origin: Panos Christakoglou | Panos.Christakoglou@cern.ch
  20 //-----------------------------------------------------------------
  21 #include <Riostream.h>
  22 #include <TCanvas.h>
  23 #include <TLatex.h>
  24 #include <TF1.h>
  25 #include <TList.h>
  26 #include <TH1F.h>
  27
  28 #include <AliExternalTrackParam.h>
  29 #include <AliESDEvent.h>
  30 #include <AliPID.h>
  31 #include <AliVertexerTracks.h>
  32 #include <AliTPCpidESD.h>
  33 class AliLog;
  34 class AliESDVertex;
  35
  36 #include "AliProtonAnalysisBase.h"
  37
  38 ClassImp(AliProtonAnalysisBase)
  39
  40 //____________________________________________________________________//
  41 AliProtonAnalysisBase::AliProtonAnalysisBase() :
  42   TObject(),  fProtonAnalysisLevel("ESD"),
  43   fTriggerMode(kMB2), fProtonAnalysisMode(kTPC), fProtonPIDMode(kBayesian),
  44   fAnalysisEtaMode(kFALSE),
  45   fVxMax(100.), fVyMax(100.), fVzMax(100.),
  46   fNBinsX(0), fMinX(0), fMaxX(0),
  47   fNBinsY(0), fMinY(0), fMaxY(0),
  48   fMinTPCClusters(0), fMinITSClusters(0),
  49   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
  50   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
  51   fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
  52   fMaxDCAXY(0), fMaxDCAXYTPC(0),
  53   fMaxDCAZ(0), fMaxDCAZTPC(0),
  54   fMaxDCA3D(0), fMaxDCA3DTPC(0),
  55   fMaxConstrainChi2(0), fMinTPCdEdxPoints(0),
  56   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  57   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  58   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE),
  59   fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
  60   fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
  61   fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
  62   fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
  63   fMaxDCA3DFlag(kFALSE), fMaxDCA3DTPCFlag(kFALSE),
  64   fMaxConstrainChi2Flag(kFALSE),
  65   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
  66   fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE), fTOFpidFlag(kFALSE),
  67   fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
  68   fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
  69   fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
  70   fMinTPCdEdxPointsFlag(kFALSE),
  71   fFunctionProbabilityFlag(kFALSE),
  72   fNSigma(0),
  73   fElectronFunction(0), fMuonFunction(0),
  74   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
  75   fDebugMode(kFALSE) {
  76   //Default constructor
  77   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
  78   for(Int_t i = 0; i < 24; i++) {
  79     fdEdxMean[i] = 0.0;
  80     fdEdxSigma[i] = 0.0;
  81   }
  82   fListVertexQA = new TList();
  83   fListVertexQA->SetName("fListVertexQA");
  84   TH1F *gHistVx = new TH1F("gHistVx",
  85                            "Vx distribution;V_{x} [cm];Entries",
  86                            100,-5.,5.);
  87   gHistVx->SetFillColor(kRed-2);
  88   fListVertexQA->Add(gHistVx);
  89   TH1F *gHistVxAccepted = new TH1F("gHistVxaccepted",
  90                                    "Vx distribution;V_{x} [cm];Entries",
  91                                    100,-5.,5.);
  92   fListVertexQA->Add(gHistVxAccepted);
  93   TH1F *gHistVy = new TH1F("gHistVy",
  94                            "Vy distribution;V_{y} [cm];Entries",
  95                            100,-5.,5.);
  96   gHistVy->SetFillColor(kRed-2);
  97   fListVertexQA->Add(gHistVy);
  98   TH1F *gHistVyAccepted = new TH1F("gHistVyaccepted",
  99                                    "Vy distribution;V_{y} [cm];Entries",
 100                                    100,-5.,5.);
 101   fListVertexQA->Add(gHistVyAccepted);
 102   TH1F *gHistVz = new TH1F("gHistVz",
 103                            "Vz distribution;V_{z} [cm];Entries",
 104                            100,-25.,25.);
 105   gHistVz->SetFillColor(kRed-2);
 106   fListVertexQA->Add(gHistVz);
 107   TH1F *gHistVzAccepted = new TH1F("gHistVzaccepted",
 108                                    "Vz distribution;V_{z} [cm];Entries",
 109                                    100,-25.,25.);
 110   fListVertexQA->Add(gHistVzAccepted);
 111
 112
 113 }
 114
 115 //____________________________________________________________________//
 116 AliProtonAnalysisBase::~AliProtonAnalysisBase() {
 117   //Default destructor
 118   if(fElectronFunction) delete fElectronFunction;
 119   if(fMuonFunction) delete fMuonFunction;
 120   if(fPionFunction) delete fPionFunction;
 121   if(fKaonFunction) delete fKaonFunction;
 122   if(fProtonFunction) delete fProtonFunction;
 123   if(fListVertexQA) delete fListVertexQA;
 124 }
 125
 126 //____________________________________________________________________//
 127 Double_t AliProtonAnalysisBase::GetParticleFraction(Int_t i, Double_t p) {
 128   //Return the a priori probs
 129   Double_t partFrac=0;
 130   if(fFunctionProbabilityFlag) {
 131     if(i == 0) partFrac = fElectronFunction->Eval(p);
 132     if(i == 1) partFrac = fMuonFunction->Eval(p);
 133     if(i == 2) partFrac = fPionFunction->Eval(p);
 134     if(i == 3) partFrac = fKaonFunction->Eval(p);
 135     if(i == 4) partFrac = fProtonFunction->Eval(p);
 136   }
 137   else partFrac = fPartFrac[i];
 138
 139   return partFrac;
 140 }
 141
 142 //____________________________________________________________________//
 143 Bool_t AliProtonAnalysisBase::IsInPhaseSpace(AliESDtrack* const track) {
 144   // Checks if the track is outside the analyzed y-Pt phase space
 145   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 146   Double_t eta = 0.0;
 147
 148   if((fProtonAnalysisMode == kTPC) || (fProtonAnalysisMode == kHybrid)) {
 149     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 150     if(!tpcTrack) {
 151       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
 152     }
 153     else {
 154       gPt = tpcTrack->Pt();
 155       gPx = tpcTrack->Px();
 156       gPy = tpcTrack->Py();
 157       gPz = tpcTrack->Pz();
 158       eta = tpcTrack->Eta();
 159     }
 160   }//standalone TPC or Hybrid TPC approaches
 161   else {
 162     gPt = track->Pt();
 163     gPx = track->Px();
 164     gPy = track->Py();
 165     gPz = track->Pz();
 166     eta = track->Eta();
 167   }
 168
 169   if((gPt < fMinY) || (gPt > fMaxY)) {
 170       if(fDebugMode)
 171         Printf("IsInPhaseSpace: Track rejected because it has a Pt value of %lf (accepted interval: %lf - %lf)",gPt,fMinY,fMaxY);
 172       return kFALSE;
 173   }
 174   if(fAnalysisEtaMode) {
 175     if((eta < fMinX) || (eta > fMaxX)) {
 176       if(fDebugMode)
 177         Printf("IsInPhaseSpace: Track rejected because it has an eta value of %lf (accepted interval: %lf - %lf)",eta,fMinX,fMaxX);
 178       return kFALSE;
 179     }
 180   }
 181   else {
 182     if((Rapidity(gPx,gPy,gPz) < fMinX) || (Rapidity(gPx,gPy,gPz) > fMaxX)) {
 183       if(fDebugMode)
 184         Printf("IsInPhaseSpace: Track rejected because it has a y value of %lf (accepted interval: %lf - %lf)",Rapidity(gPx,gPy,gPz),fMinX,fMaxX);
 185       return kFALSE;
 186     }
 187   }
 188
 189   return kTRUE;
 190 }
 191
 192 //____________________________________________________________________//
 193 Bool_t AliProtonAnalysisBase::IsAccepted(AliESDEvent *esd,
 194                                          const AliESDVertex *vertex,
 195                                          AliESDtrack* track) {
 196   // Checks if the track is excluded from the cuts
 197   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 198   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
 199   Double_t dca3D = 0.0;
 200
 201   if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
 202     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 203     if(!tpcTrack) {
 204       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
 205       dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
 206       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
 207     }
 208     else {
 209       gPt = tpcTrack->Pt();
 210       gPx = tpcTrack->Px();
 211       gPy = tpcTrack->Py();
 212       gPz = tpcTrack->Pz();
 213       tpcTrack->PropagateToDCA(vertex,
 214                                esd->GetMagneticField(),
 215                                100.,dca,cov);
 216     }
 217   }//standalone TPC or hybrid TPC approaches
 218   else {
 219     gPt = track->Pt();
 220     gPx = track->Px();
 221     gPy = track->Py();
 222     gPz = track->Pz();
 223     track->PropagateToDCA(vertex,
 224                           esd->GetMagneticField(),
 225                           100.,dca,cov);
 226   }
 227   dca3D = TMath::Sqrt(TMath::Power(dca[0],2) +
 228                       TMath::Power(dca[1],2));
 229
 230   Int_t  fIdxInt[200];
 231   Int_t nClustersITS = track->GetITSclusters(fIdxInt);
 232   Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
 233
 234   Float_t chi2PerClusterITS = -1;
 235   if (nClustersITS!=0)
 236     chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
 237   Float_t chi2PerClusterTPC = -1;
 238   if (nClustersTPC!=0)
 239     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
 240
 241   Double_t extCov[15];
 242   track->GetExternalCovariance(extCov);
 243
 244   if(fPointOnITSLayer1Flag) {
 245     if(!track->HasPointOnITSLayer(0)) {
 246       if(fDebugMode)
 247         Printf("IsAccepted: Track rejected because it doesn't have a point on the 1st ITS layer");
 248       return kFALSE;
 249     }
 250   }
 251   if(fPointOnITSLayer2Flag) {
 252     if(!track->HasPointOnITSLayer(1)) {
 253       if(fDebugMode)
 254         Printf("IsAccepted: Track rejected because it doesn't have a point on the 2nd ITS layer");
 255       return kFALSE;
 256     }
 257   }
 258   if(fPointOnITSLayer3Flag) {
 259     if(!track->HasPointOnITSLayer(2)) {
 260       if(fDebugMode)
 261         Printf("IsAccepted: Track rejected because it doesn't have a point on the 3rd ITS layer");
 262       return kFALSE;
 263     }
 264   }
 265   if(fPointOnITSLayer4Flag) {
 266     if(!track->HasPointOnITSLayer(3)) {
 267       if(fDebugMode)
 268         Printf("IsAccepted: Track rejected because it doesn't have a point on the 4th ITS layer");
 269       return kFALSE;
 270     }
 271   }
 272   if(fPointOnITSLayer5Flag) {
 273     if(!track->HasPointOnITSLayer(4)) {
 274       if(fDebugMode)
 275         Printf("IsAccepted: Track rejected because it doesn't have a point on the 5th ITS layer");
 276       return kFALSE;
 277     }
 278   }
 279   if(fPointOnITSLayer6Flag) {
 280     if(!track->HasPointOnITSLayer(5)) {
 281       if(fDebugMode)
 282         Printf("IsAccepted: Track rejected because it doesn't have a point on the 6th ITS layer");
 283       return kFALSE;
 284     }
 285   }
 286   if(fMinITSClustersFlag) {
 287     if(nClustersITS < fMinITSClusters) {
 288       if(fDebugMode)
 289         Printf("IsAccepted: Track rejected because it has %d ITS points (min. requested: %d)",nClustersITS,fMinITSClusters);
 290       return kFALSE;
 291     }
 292   }
 293   if(fMaxChi2PerITSClusterFlag) {
 294     if(chi2PerClusterITS > fMaxChi2PerITSCluster) {
 295       if(fDebugMode)
 296         Printf("IsAccepted: Track rejected because it has a chi2 per ITS cluster %lf (max. requested: %lf)",chi2PerClusterITS,fMaxChi2PerITSCluster);
 297       return kFALSE;
 298     }
 299   }
 300   if(fMinTPCClustersFlag) {
 301     if(nClustersTPC < fMinTPCClusters) {
 302       if(fDebugMode)
 303         Printf("IsAccepted: Track rejected because it has %d TPC clusters (min. requested: %d)",nClustersTPC,fMinTPCClusters);
 304       return kFALSE;
 305     }
 306   }
 307   if(fMaxChi2PerTPCClusterFlag) {
 308     if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) {
 309       if(fDebugMode)
 310         Printf("IsAccepted: Track rejected because it has a chi2 per TPC cluster %lf (max. requested: %lf)",chi2PerClusterTPC,fMaxChi2PerTPCCluster);
 311       return kFALSE;
 312     }
 313   }
 314   if(fMaxCov11Flag) {
 315     if(extCov[0] > fMaxCov11) {
 316       if(fDebugMode)
 317         Printf("IsAccepted: Track rejected because it has a cov11 value of %lf (max. requested: %lf)",extCov[0],fMaxCov11);
 318       return kFALSE;
 319     }
 320   }
 321   if(fMaxCov22Flag) {
 322     if(extCov[2] > fMaxCov22) {
 323       if(fDebugMode)
 324         Printf("IsAccepted: Track rejected because it has a cov22 value of %lf (max. requested: %lf)",extCov[2],fMaxCov22);
 325       return kFALSE;
 326     }
 327   }
 328   if(fMaxCov33Flag) {
 329     if(extCov[5] > fMaxCov33) {
 330       if(fDebugMode)
 331         Printf("IsAccepted: Track rejected because it has a cov33 value of %lf (max. requested: %lf)",extCov[5],fMaxCov33);
 332       return kFALSE;
 333     }
 334   }
 335   if(fMaxCov44Flag) {
 336     if(extCov[9] > fMaxCov44) {
 337       if(fDebugMode)
 338         Printf("IsAccepted: Track rejected because it has a cov44 value of %lf (max. requested: %lf)",extCov[9],fMaxCov44);
 339       return kFALSE;
 340     }
 341   }
 342   if(fMaxCov55Flag) {
 343     if(extCov[14] > fMaxCov55) {
 344       if(fDebugMode)
 345         Printf("IsAccepted: Track rejected because it has a cov55 value of %lf (max. requested: %lf)",extCov[14],fMaxCov55);
 346       return kFALSE;
 347     }
 348   }
 349   if(fMaxSigmaToVertexFlag) {
 350     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) {
 351       if(fDebugMode)
 352         Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertex);
 353       return kFALSE;
 354     }
 355   }
 356   if(fMaxSigmaToVertexTPCFlag) {
 357     if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) {
 358       if(fDebugMode)
 359         Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertexTPC);
 360       return kFALSE;
 361     }
 362   }
 363   if(fMaxDCAXYFlag) {
 364     if(TMath::Abs(dca[0]) > fMaxDCAXY) {
 365       if(fDebugMode)
 366         Printf("IsAccepted: Track rejected because it has a value of dca(xy) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXY);
 367       return kFALSE;
 368     }
 369   }
 370   if(fMaxDCAXYTPCFlag) {
 371     if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) {
 372       if(fDebugMode)
 373         Printf("IsAccepted: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXYTPC);
 374       return kFALSE;
 375     }
 376   }
 377   if(fMaxDCAZFlag) {
 378     if(TMath::Abs(dca[1]) > fMaxDCAZ) {
 379       if(fDebugMode)
 380         Printf("IsAccepted: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZ);
 381       return kFALSE;
 382     }
 383   }
 384   if(fMaxDCAZTPCFlag) {
 385     if(TMath::Abs(dca[1]) > fMaxDCAZTPC) {
 386       if(fDebugMode)
 387         Printf("IsAccepted: Track rejected because it has a value of dca(z) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZTPC);
 388       return kFALSE;
 389     }
 390   }
 391   if(fMaxDCA3DFlag) {
 392     if(TMath::Abs(dca3D) > fMaxDCA3D) {
 393       if(fDebugMode)
 394         Printf("IsAccepted: Track rejected because it has a value of dca(3D) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3D);
 395       return kFALSE;
 396     }
 397   }
 398   if(fMaxDCA3DTPCFlag) {
 399     if(TMath::Abs(dca3D) > fMaxDCA3DTPC)  {
 400       if(fDebugMode)
 401         Printf("IsAccepted: Track rejected because it has a value of dca(3D) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3DTPC);
 402       return kFALSE;
 403     }
 404   }
 405   if(fMaxConstrainChi2Flag) {
 406     if(track->GetConstrainedChi2() > 0)
 407       if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2)  {
 408       if(fDebugMode)
 409         Printf("IsAccepted: Track rejected because it has a value of the constrained chi2 to the vertex of %lf (max. requested: %lf)",TMath::Log(track->GetConstrainedChi2()),fMaxConstrainChi2);
 410       return kFALSE;
 411       }
 412   }
 413   if(fMinTPCdEdxPointsFlag) {
 414     if(track->GetTPCsignalN() < fMinTPCdEdxPoints) {
 415       if(fDebugMode)
 416         Printf("IsAccepted: Track rejected because it has %d TPC points for the calculation of the energy loss (min. requested: %d)",track->GetTPCsignalN(),fMinTPCdEdxPoints);
 417       return kFALSE;
 418     }
 419   }
 420   if(fITSRefitFlag) {
 421     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
 422       if(fDebugMode)
 423         Printf("IsAccepted: Track rejected because it has no ITS refit flag");
 424       return kFALSE;
 425     }
 426   }
 427   if(fTPCRefitFlag) {
 428     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
 429       if(fDebugMode)
 430         Printf("IsAccepted: Track rejected because it has no TPC refit flag");
 431       return kFALSE;
 432     }
 433   }
 434   if(fESDpidFlag) {
 435     if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) {
 436       if(fDebugMode)
 437         Printf("IsAccepted: Track rejected because it has no ESD pid flag");
 438       return kFALSE;
 439     }
 440   }
 441   if(fTPCpidFlag) {
 442     if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) {
 443       if(fDebugMode)
 444         Printf("IsAccepted: Track rejected because it has no TPC pid flag");
 445       return kFALSE;
 446     }
 447   }
 448   if(fTOFpidFlag) {
 449     if ((track->GetStatus() & AliESDtrack::kTOFpid) == 0) {
 450       if(fDebugMode)
 451         Printf("IsAccepted: Track rejected because it has no TOF pid flag");
 452       return kFALSE;
 453     }
 454   }
 455
 456   return kTRUE;
 457 }
 458
 459 //____________________________________________________________________//
 460 Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
 461   // Calculates the number of sigma to the vertex.
 462
 463   Float_t b[2];
 464   Float_t bRes[2];
 465   Float_t bCov[3];
 466   if((fProtonAnalysisMode == kTPC)&&(fProtonAnalysisMode != kHybrid))
 467     esdTrack->GetImpactParametersTPC(b,bCov);
 468   else
 469     esdTrack->GetImpactParameters(b,bCov);
 470
 471   if (bCov[0]<=0 || bCov[2]<=0) {
 472     //AliDebug(1, "Estimated b resolution lower or equal zero!");
 473     bCov[0]=0; bCov[2]=0;
 474   }
 475   bRes[0] = TMath::Sqrt(bCov[0]);
 476   bRes[1] = TMath::Sqrt(bCov[2]);
 477
 478   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 479
 480   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 481
 482   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 483
 484   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 485
 486   return d;
 487 }
 488
 489 //____________________________________________________________________//
 490 Double_t AliProtonAnalysisBase::Rapidity(Double_t gPx,
 491                                          Double_t gPy,
 492                                          Double_t gPz) const {
 493   //returns the rapidity of the proton - to be removed
 494   Double_t fMass = 9.38270000000000048e-01;
 495
 496   Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) +
 497                            TMath::Power(gPy,2) +
 498                            TMath::Power(gPz,2));
 499   Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
 500   Double_t y = -999;
 501   if(energy != gPz)
 502     y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
 503
 504   return y;
 505 }
 506
 507 //________________________________________________________________________
 508 const AliESDVertex* AliProtonAnalysisBase::GetVertex(AliESDEvent* esd,
 509                                                      AnalysisMode mode,
 510                                                      Double_t gVxMax,
 511                                                      Double_t gVyMax,
 512                                                      Double_t gVzMax) {
 513   // Get the vertex from the ESD and returns it if the vertex is valid
 514   // Second argument decides which vertex is used (this selects
 515   // also the quality criteria that are applied)
 516   const AliESDVertex* vertex = 0;
 517   if(mode == kHybrid)
 518     vertex = esd->GetPrimaryVertexSPD();
 519   else if(mode == kTPC){
 520     Double_t kBz = esd->GetMagneticField();
 521     AliVertexerTracks vertexer(kBz);
 522     vertexer.SetTPCMode();
 523     AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
 524     esd->SetPrimaryVertexTPC(vTPC);
 525     for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
 526       AliESDtrack *t = esd->GetTrack(i);
 527       t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
 528     }
 529     delete vTPC;
 530     vertex = esd->GetPrimaryVertexTPC();
 531   }
 532   else if(mode == kGlobal)
 533     vertex = esd->GetPrimaryVertex();
 534   else
 535     Printf("GetVertex: ERROR: Invalid second argument %d", mode);
 536
 537   if(!vertex) {
 538     if(fDebugMode)
 539       Printf("GetVertex: Event rejected because there is no valid vertex object");
 540     return 0;
 541   }
 542
 543   // check Ncontributors
 544   if(vertex->GetNContributors() <= 0) {
 545     if(fDebugMode)
 546       Printf("GetVertex: Event rejected because the number of contributors for the vertex determination is <= 0");
 547     return 0;
 548   }
 549
 550   // check resolution
 551   Double_t zRes = vertex->GetZRes();
 552   if(zRes == 0) {
 553     if(fDebugMode)
 554       Printf("GetVertex: Event rejected because the value of the vertex resolution in z is 0");
 555     return 0;
 556   }
 557   ((TH1F *)(fListVertexQA->At(0)))->Fill(vertex->GetXv());
 558   ((TH1F *)(fListVertexQA->At(2)))->Fill(vertex->GetYv());
 559   ((TH1F *)(fListVertexQA->At(4)))->Fill(vertex->GetZv());
 560   //check position
 561   if(TMath::Abs(vertex->GetXv()) > gVxMax) {
 562     if(fDebugMode)
 563       Printf("GetVertex: Event rejected because it has a Vx value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetXv()),gVxMax,gVxMax);
 564     return 0;
 565   }
 566   if(TMath::Abs(vertex->GetYv()) > gVyMax)  {
 567     if(fDebugMode)
 568       Printf("GetVertex: Event rejected because it has a Vy value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetYv()),gVyMax,gVyMax);
 569     return 0;
 570   }
 571   if(TMath::Abs(vertex->GetZv()) > gVzMax)  {
 572     if(fDebugMode)
 573       Printf("GetVertex: Event rejected because it has a Vz value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetZv()),gVzMax,gVzMax);
 574     return 0;
 575   }
 576   ((TH1F *)(fListVertexQA->At(1)))->Fill(vertex->GetXv());
 577   ((TH1F *)(fListVertexQA->At(3)))->Fill(vertex->GetYv());
 578   ((TH1F *)(fListVertexQA->At(5)))->Fill(vertex->GetZv());
 579
 580   return vertex;
 581 }
 582
 583 //________________________________________________________________________
 584 Bool_t AliProtonAnalysisBase::IsEventTriggered(const AliESDEvent *esd,
 585                                                TriggerMode trigger) {
 586   // check if the event was triggered
 587   ULong64_t triggerMask = esd->GetTriggerMask();
 588
 589   // definitions from p-p.cfg
 590   ULong64_t spdFO = (1 << 14);
 591   ULong64_t v0left = (1 << 11);
 592   ULong64_t v0right = (1 << 12);
 593
 594   switch (trigger) {
 595   case kMB1: {
 596     if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
 597       return kTRUE;
 598     break;
 599   }
 600   case kMB2: {
 601     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
 602       return kTRUE;
 603     break;
 604   }
 605   case kSPDFASTOR: {
 606     if (triggerMask & spdFO)
 607       return kTRUE;
 608     break;
 609   }
 610   }//switch
 611
 612   return kFALSE;
 613 }
 614
 615 //________________________________________________________________________
 616 TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
 617   // return the list of cuts and their cut values for reference
 618   TLatex l;
 619   l.SetTextAlign(12);
 620   l.SetTextSize(0.04);
 621
 622   TString listOfCuts;
 623
 624   TCanvas *c = new TCanvas("cListOfCuts","List of cuts",0,0,900,600);
 625   c->SetFillColor(10); c->SetHighLightColor(41);
 626   c->Divide(3,2);
 627
 628   c->cd(1)->SetFillColor(10);
 629   l.DrawLatex(0.3,0.9,"Analysis details: List of cuts\n\n");
 630
 631   listOfCuts = "Analysis level: "; listOfCuts += fProtonAnalysisLevel;
 632   l.DrawLatex(0.1,0.82,listOfCuts.Data());
 633   listOfCuts = "Analysis mode: ";
 634   if(fProtonAnalysisMode == kTPC) listOfCuts += "TPC standalone";
 635   if(fProtonAnalysisMode == kHybrid) listOfCuts += "Hybrid TPC";
 636   if(fProtonAnalysisMode == kGlobal) listOfCuts += "Global tracking";
 637   l.DrawLatex(0.1,0.74,listOfCuts.Data());
 638   listOfCuts = "Trigger mode: ";
 639   if(fTriggerMode == kMB1) listOfCuts += "Minimum bias 1";
 640   if(fTriggerMode == kMB2) listOfCuts += "Minimum bias 2";
 641   if(fTriggerMode == kSPDFASTOR) listOfCuts += "FastOR (SPD)";
 642   l.DrawLatex(0.1,0.66,listOfCuts.Data());
 643   listOfCuts = "PID mode: ";
 644   if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
 645   if(fProtonPIDMode == kRatio) listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.})";
 646   if(fProtonPIDMode == kSigma1) {
 647     listOfCuts += "N_{#sigma}(1) area: "; listOfCuts += fNSigma;
 648     listOfCuts += " #sigma";
 649   }
 650   if(fProtonPIDMode == kSigma2) {
 651     listOfCuts += "N_{#sigma}(2) area: "; listOfCuts += fNSigma;
 652     listOfCuts += " #sigma";
 653   }
 654   l.DrawLatex(0.1,0.58,listOfCuts.Data());
 655   listOfCuts = "Accepted vertex diamond: ";
 656   l.DrawLatex(0.1,0.5,listOfCuts.Data());
 657   listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
 658   l.DrawLatex(0.6,0.5,listOfCuts.Data());
 659   listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
 660   l.DrawLatex(0.6,0.4,listOfCuts.Data());
 661   listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
 662   l.DrawLatex(0.6,0.3,listOfCuts.Data());
 663   listOfCuts = "Phase space: ";
 664   l.DrawLatex(0.1,0.2,listOfCuts.Data());
 665   if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";
 666   else listOfCuts = "|y| < ";
 667   listOfCuts += TMath::Abs(fMinX);
 668   l.DrawLatex(0.35,0.2,listOfCuts.Data());
 669   listOfCuts = "N_{bins} = ";
 670   listOfCuts += fNBinsX; listOfCuts += " (binning: ";
 671   listOfCuts += (fMaxX-fMinX)/fNBinsX; listOfCuts += ")";
 672   l.DrawLatex(0.35,0.15,listOfCuts.Data());
 673   listOfCuts = ""; listOfCuts += fMinY; listOfCuts += " < P_{T} < ";
 674   listOfCuts += fMaxY; listOfCuts += "GeV/c";
 675   l.DrawLatex(0.35,0.1,listOfCuts.Data());
 676   listOfCuts = "N_{bins} = ";
 677   listOfCuts += fNBinsY; listOfCuts += " (binning: ";
 678   listOfCuts += (fMaxY-fMinY)/fNBinsY; listOfCuts += ")";
 679   l.DrawLatex(0.35,0.05,listOfCuts.Data());
 680
 681   c->cd(2)->SetFillColor(2);
 682   l.DrawLatex(0.3,0.9,"ITS related cuts");
 683   listOfCuts = "Request a cluster on SPD1: ";
 684   listOfCuts += fPointOnITSLayer1Flag;
 685   l.DrawLatex(0.1,0.8,listOfCuts.Data());
 686   listOfCuts = "Request a cluster on SPD2: ";
 687   listOfCuts += fPointOnITSLayer2Flag;
 688   l.DrawLatex(0.1,0.7,listOfCuts.Data());
 689   listOfCuts = "Request a cluster on SDD1: ";
 690   listOfCuts += fPointOnITSLayer3Flag;
 691   l.DrawLatex(0.1,0.6,listOfCuts.Data());
 692   listOfCuts = "Request a cluster on SDD2: ";
 693   listOfCuts += fPointOnITSLayer4Flag;
 694   l.DrawLatex(0.1,0.5,listOfCuts.Data());
 695   listOfCuts = "Request a cluster on SSD1: ";
 696   listOfCuts += fPointOnITSLayer5Flag;
 697   l.DrawLatex(0.1,0.4,listOfCuts.Data());
 698   listOfCuts = "Request a cluster on SSD2: ";
 699   listOfCuts += fPointOnITSLayer6Flag;
 700   l.DrawLatex(0.1,0.3,listOfCuts.Data());
 701   listOfCuts = "Minimum number of ITS clusters: ";
 702   if(fMinITSClustersFlag) listOfCuts += fMinITSClusters;
 703   else listOfCuts += "Not used";
 704   l.DrawLatex(0.1,0.2,listOfCuts.Data());
 705   listOfCuts = "Maximum #chi^{2} per ITS cluster: ";
 706   if(fMaxChi2PerITSClusterFlag) listOfCuts += fMaxChi2PerITSCluster;
 707   else listOfCuts += "Not used";
 708   l.DrawLatex(0.1,0.1,listOfCuts.Data());
 709
 710   c->cd(3)->SetFillColor(3);
 711   l.DrawLatex(0.3,0.9,"TPC related cuts");
 712   listOfCuts = "Minimum number of TPC clusters: ";
 713   if(fMinTPCClustersFlag) listOfCuts += fMinTPCClusters;
 714   else listOfCuts += "Not used";
 715   l.DrawLatex(0.1,0.7,listOfCuts.Data());
 716   listOfCuts = "Maximum #chi^{2} per TPC cluster: ";
 717   if(fMaxChi2PerTPCClusterFlag) listOfCuts += fMaxChi2PerTPCCluster;
 718   else listOfCuts += "Not used";
 719   l.DrawLatex(0.1,0.5,listOfCuts.Data());
 720   listOfCuts = "Minimum number of TPC points for the dE/dx: ";
 721   if(fMinTPCdEdxPointsFlag) listOfCuts += fMinTPCdEdxPoints;
 722   else listOfCuts += "Not used";
 723   l.DrawLatex(0.1,0.3,listOfCuts.Data());
 724
 725   c->cd(4)->SetFillColor(4);
 726   l.DrawLatex(0.3,0.9,"Tracking related cuts");
 727   listOfCuts = "Maximum cov11: ";
 728   if(fMaxCov11Flag) listOfCuts += fMaxCov11;
 729   else listOfCuts += "Not used";
 730   l.DrawLatex(0.1,0.75,listOfCuts.Data());
 731   listOfCuts = "Maximum cov22: ";
 732   if(fMaxCov22Flag) listOfCuts += fMaxCov22;
 733   else listOfCuts += "Not used";
 734   l.DrawLatex(0.1,0.6,listOfCuts.Data());
 735   listOfCuts = "Maximum cov33: ";
 736   if(fMaxCov33Flag) listOfCuts += fMaxCov33;
 737   else listOfCuts += "Not used";
 738   l.DrawLatex(0.1,0.45,listOfCuts.Data());
 739   listOfCuts = "Maximum cov44: ";
 740   if(fMaxCov44Flag) listOfCuts += fMaxCov44;
 741   else listOfCuts += "Not used";
 742   l.DrawLatex(0.1,0.3,listOfCuts.Data());
 743   listOfCuts = "Maximum cov55: ";
 744   if(fMaxCov55Flag) listOfCuts += fMaxCov55;
 745   else listOfCuts += "Not used";
 746   l.DrawLatex(0.1,0.15,listOfCuts.Data());
 747
 748   c->cd(5)->SetFillColor(5);
 749   l.DrawLatex(0.3,0.9,"DCA related cuts");
 750   listOfCuts = "Maximum sigma to vertex: ";
 751   if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
 752   else listOfCuts += "Not used";
 753   l.DrawLatex(0.1,0.8,listOfCuts.Data());
 754   listOfCuts = "Maximum sigma to vertex (TPC): ";
 755   if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
 756   else listOfCuts += "Not used";
 757   l.DrawLatex(0.1,0.72,listOfCuts.Data());
 758   listOfCuts = "Maximum DCA in xy: ";
 759   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
 760   else listOfCuts += "Not used";
 761   l.DrawLatex(0.1,0.64,listOfCuts.Data());
 762   listOfCuts = "Maximum DCA in z: ";
 763   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
 764   else listOfCuts += "Not used";
 765   l.DrawLatex(0.1,0.56,listOfCuts.Data());
 766   listOfCuts = "Maximum DCA in 3D: ";
 767   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
 768   else listOfCuts += "Not used";
 769   l.DrawLatex(0.1,0.48,listOfCuts.Data());
 770   listOfCuts = "Maximum DCA in xy (TPC): ";
 771   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
 772   else listOfCuts += "Not used";
 773   l.DrawLatex(0.1,0.4,listOfCuts.Data());
 774   listOfCuts = "Maximum DCA in z (TPC): ";
 775   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
 776   else listOfCuts += "Not used";
 777   l.DrawLatex(0.1,0.32,listOfCuts.Data());
 778   listOfCuts = "Maximum DCA in 3D (TPC): ";
 779   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
 780   else listOfCuts += "Not used";
 781   l.DrawLatex(0.1,0.24,listOfCuts.Data());
 782   listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
 783   if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
 784   else listOfCuts += "Not used";
 785   l.DrawLatex(0.1,0.16,listOfCuts.Data());
 786
 787   c->cd(6)->SetFillColor(6);
 788   l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
 789   listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
 790   l.DrawLatex(0.1,0.7,listOfCuts.Data());
 791   listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag;
 792   l.DrawLatex(0.1,0.5,listOfCuts.Data());
 793   listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
 794   l.DrawLatex(0.1,0.3,listOfCuts.Data());
 795   listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
 796   l.DrawLatex(0.1,0.1,listOfCuts.Data());
 797
 798   return c;
 799 }
 800
 801 //________________________________________________________________________
 802 Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
 803   //Function that checks if a track is a proton
 804   Double_t probability[5];
 805   Double_t gPt = 0.0, gP = 0.0, gEta = 0.0;
 806   Long64_t fParticleType = 0;
 807
 808   //Bayesian approach for the PID
 809   if(fProtonPIDMode == kBayesian) {
 810     if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
 811       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 812       if(tpcTrack) {
 813         gPt = tpcTrack->Pt();
 814         gP = tpcTrack->P();
 815         track->GetTPCpid(probability);
 816       }
 817     }//TPC standalone or Hybrid TPC
 818     else if(fProtonAnalysisMode == kGlobal) {
 819       gPt = track->Pt();
 820       gP = track->P();
 821       track->GetESDpid(probability);
 822     }//Global tracking
 823
 824     Double_t rcc = 0.0;
 825     for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 826       rcc += probability[i]*GetParticleFraction(i,gP);
 827     if(rcc != 0.0) {
 828       Double_t w[5];
 829       for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 830         w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
 831       fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
 832     }
 833     if(fParticleType == 4)
 834       return kTRUE;
 835   }
 836   //Ratio of the measured over the theoretical dE/dx a la STAR
 837   else if(fProtonPIDMode == kRatio) {
 838     Printf("The kRatio mode is not implemented yet!!!");
 839     return kFALSE;
 840   }
 841   //Definition of an N-sigma area around the dE/dx vs P band
 842   else if(fProtonPIDMode == kSigma1) {
 843     Double_t fAlephParameters[5];
 844     fAlephParameters[0] = 4.23232575531564326e+00; //50*0.76176e-1;
 845     fAlephParameters[1] = 8.68482806165147636e+00; //10.632;
 846     fAlephParameters[2] = 1.34000000000000005e-05; //0.13279e-4;
 847     fAlephParameters[3] = 2.30445734159456084e+00; //1.8631;
 848     fAlephParameters[4] = 2.25624744086878559e+00; //1.9479;
 849
 850     AliTPCpidESD *fPidObject = new AliTPCpidESD();
 851     fPidObject->SetBetheBlochParameters(fAlephParameters[0]/50.,fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
 852
 853     Double_t nsigma = TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kProton));
 854     if(nsigma <= fNSigma)
 855       return kTRUE;
 856   }//kSigma1 PID method
 857   //Another definition of an N-sigma area around the dE/dx vs P band
 858   else if(fProtonPIDMode == kSigma2) {
 859     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 860     if(tpcTrack) {
 861       gPt = tpcTrack->Pt();
 862       gP = tpcTrack->P();
 863       gEta = tpcTrack->Eta();
 864     }
 865     //We start the P slices at >0.3GeV/c with a bining of 50MeV/c ==> Int_t(0.3001/0.05) = 6
 866     Int_t nbinP = Int_t(gP/0.05) - 6;
 867     Double_t tpcSignal = track->GetTPCsignal();
 868     Double_t dEdxTheory = Bethe(gP/9.38270000000000048e-01);
 869     Double_t dEdxSigma = fdEdxSigma[nbinP];
 870     Double_t nsigma = TMath::Abs(tpcSignal - dEdxTheory)/(tpcSignal*(dEdxSigma/TMath::Sqrt(track->GetTPCsignalN())));
 871     if(nsigma <= fNSigma)
 872       return kTRUE;
 873   }
 874
 875   return kFALSE;
 876 }
 877
 878 //________________________________________________________________________
 879 void AliProtonAnalysisBase::SetdEdxBandInfo(const char *filename) {
 880   // This function is used in case the kSigma1 or kSigma2 PID mode is selected
 881   // It takes as an argument the name of the ascii file (for the time being)
 882   // that is generated as a prior process.
 883   // This ascii file has three columns: The min. P value (bins of 50MeV/c)
 884   // the mean and the sigma of the dE/dx distributions for protons coming
 885   // from a gaussian fit.
 886   ifstream in;
 887   in.open(filename);
 888
 889   Double_t gPtMin = 0.0;
 890   Int_t iCounter = 0;
 891   while(in.good()) {
 892     in >> gPtMin >> fdEdxMean[iCounter] >> fdEdxSigma[iCounter];
 893     if(fDebugMode)
 894       Printf("Momentum bin: %d - Min momentum: %lf - mean(dE/dx): %lf - sigma(dE/dx): %lf",iCounter+1,gPtMin,fdEdxMean[iCounter],fdEdxSigma[iCounter]);
 895     iCounter += 1;
 896   }
 897 }
 898
 899 //________________________________________________________________________
 900 Double_t AliProtonAnalysisBase::Bethe(Double_t bg) const {
 901   // This is the Bethe-Bloch function normalised to 1 at the minimum
 902   // We renormalize it based on the MC information
 903   // WARNING: To be revised soon!!!
 904   // This is just a temporary fix
 905   Double_t normalization = 49.2;
 906   Double_t bg2=bg*bg;
 907   Double_t beta2 = bg2/(1.+ bg2);
 908   Double_t bb = 8.62702e-2*(9.14550 - beta2 - TMath::Log(3.51000e-5 + 1./bg2))/beta2;
 909   //
 910   const Float_t kmeanCorrection =0.1;
 911   Double_t meanCorrection =(1+(bb-1)*kmeanCorrection);
 912   bb *= meanCorrection;
 913
 914   return normalization*bb;
 915 }
 916
 917