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 <TFile.h>
  23 #include <TSystem.h>
  24 #include <TF1.h>
  25 #include <TH2D.h>
  26 #include <TH1D.h>
  27 #include <TH1I.h>
  28 #include <TParticle.h>
  29 #include <TCanvas.h>
  30 #include <TLatex.h>
  31
  32 #include "AliProtonAnalysisBase.h"
  33
  34 #include <AliExternalTrackParam.h>
  35 #include <AliESDEvent.h>
  36 #include <AliLog.h>
  37 #include <AliPID.h>
  38 #include <AliESDVertex.h>
  39 #include <AliVertexerTracks.h>
  40
  41 ClassImp(AliProtonAnalysisBase)
  42
  43 //____________________________________________________________________//
  44 AliProtonAnalysisBase::AliProtonAnalysisBase() :
  45   TObject(),  fProtonAnalysisLevel("ESD"),
  46   fTriggerMode(kMB2), fProtonAnalysisMode(kTPC), fProtonPIDMode(kBayesian),
  47   fAnalysisEtaMode(kFALSE),
  48   fVxMax(100.), fVyMax(100.), fVzMax(100.),
  49   fNBinsX(0), fMinX(0), fMaxX(0),
  50   fNBinsY(0), fMinY(0), fMaxY(0),
  51   fMinTPCClusters(0), fMinITSClusters(0),
  52   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
  53   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
  54   fMaxSigmaToVertex(0), fMaxSigmaToVertexTPC(0),
  55   fMaxDCAXY(0), fMaxDCAXYTPC(0),
  56   fMaxDCAZ(0), fMaxDCAZTPC(0),
  57   fMaxDCA3D(0), fMaxDCA3DTPC(0),
  58   fMaxConstrainChi2(0),
  59   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  60   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  61   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE),
  62   fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
  63   fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
  64   fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
  65   fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
  66   fMaxDCA3DFlag(kFALSE), fMaxDCA3DTPCFlag(kFALSE),
  67   fMaxConstrainChi2Flag(kFALSE),
  68   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
  69   fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE),
  70   fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
  71   fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
  72   fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
  73   fFunctionProbabilityFlag(kFALSE),
  74   fElectronFunction(0), fMuonFunction(0),
  75   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
  76   fDebugMode(kFALSE) {
  77   //Default constructor
  78   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
  79 }
  80
  81 //____________________________________________________________________//
  82 AliProtonAnalysisBase::~AliProtonAnalysisBase() {
  83   //Default destructor
  84   if(fElectronFunction) delete fElectronFunction;
  85   if(fMuonFunction) delete fMuonFunction;
  86   if(fPionFunction) delete fPionFunction;
  87   if(fKaonFunction) delete fKaonFunction;
  88   if(fProtonFunction) delete fProtonFunction;
  89 }
  90
  91 //____________________________________________________________________//
  92 Double_t AliProtonAnalysisBase::GetParticleFraction(Int_t i, Double_t p) {
  93   //Return the a priori probs
  94   Double_t partFrac=0;
  95   if(fFunctionProbabilityFlag) {
  96     if(i == 0) partFrac = fElectronFunction->Eval(p);
  97     if(i == 1) partFrac = fMuonFunction->Eval(p);
  98     if(i == 2) partFrac = fPionFunction->Eval(p);
  99     if(i == 3) partFrac = fKaonFunction->Eval(p);
 100     if(i == 4) partFrac = fProtonFunction->Eval(p);
 101   }
 102   else partFrac = fPartFrac[i];
 103
 104   return partFrac;
 105 }
 106
 107 //____________________________________________________________________//
 108 Bool_t AliProtonAnalysisBase::IsInPhaseSpace(AliESDtrack* const track) {
 109   // Checks if the track is outside the analyzed y-Pt phase space
 110   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 111   Double_t eta = 0.0;
 112
 113   if(kTPC) {
 114     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 115     if(!tpcTrack) {
 116       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
 117     }
 118     else {
 119       gPt = tpcTrack->Pt();
 120       gPx = tpcTrack->Px();
 121       gPy = tpcTrack->Py();
 122       gPz = tpcTrack->Pz();
 123       eta = tpcTrack->Eta();
 124     }
 125   }
 126   else {
 127     gPt = track->Pt();
 128     gPx = track->Px();
 129     gPy = track->Py();
 130     gPz = track->Pz();
 131     eta = track->Eta();
 132   }
 133
 134   if((gPt < fMinY) || (gPt > fMaxY)) {
 135       if(fDebugMode)
 136         Printf("IsInPhaseSpace: Track rejected because it has a Pt value of %lf (accepted interval: %lf - %lf)",gPt,fMinY,fMaxY);
 137       return kFALSE;
 138   }
 139   if(fAnalysisEtaMode) {
 140     if((eta < fMinX) || (eta > fMaxX)) {
 141       if(fDebugMode)
 142         Printf("IsInPhaseSpace: Track rejected because it has an eta value of %lf (accepted interval: %lf - %lf)",eta,fMinX,fMaxX);
 143       return kFALSE;
 144     }
 145   }
 146   else {
 147     if((Rapidity(gPx,gPy,gPz) < fMinX) || (Rapidity(gPx,gPy,gPz) > fMaxX)) {
 148       if(fDebugMode)
 149         Printf("IsInPhaseSpace: Track rejected because it has a y value of %lf (accepted interval: %lf - %lf)",Rapidity(gPx,gPy,gPz),fMinX,fMaxX);
 150       return kFALSE;
 151     }
 152   }
 153
 154   return kTRUE;
 155 }
 156
 157 //____________________________________________________________________//
 158 Bool_t AliProtonAnalysisBase::IsAccepted(AliESDEvent *esd,
 159                                          const AliESDVertex *vertex,
 160                                          AliESDtrack* track) {
 161   // Checks if the track is excluded from the cuts
 162   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 163   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
 164   Double_t dca3D = 0.0;
 165
 166   if((kTPC)&&(!kHybrid)) {
 167     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 168     if(!tpcTrack) {
 169       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
 170       dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
 171       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
 172     }
 173     else {
 174       gPt = tpcTrack->Pt();
 175       gPx = tpcTrack->Px();
 176       gPy = tpcTrack->Py();
 177       gPz = tpcTrack->Pz();
 178       tpcTrack->PropagateToDCA(vertex,
 179                                esd->GetMagneticField(),
 180                                100.,dca,cov);
 181     }
 182   }
 183   else if(kHybrid) {
 184      AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 185     if(!tpcTrack) {
 186       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
 187       dca[0] = -100.; dca[1] = -100.;
 188       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
 189     }
 190     else {
 191       gPt = tpcTrack->Pt();
 192       gPx = tpcTrack->Px();
 193       gPy = tpcTrack->Py();
 194       gPz = tpcTrack->Pz();
 195       tpcTrack->PropagateToDCA(vertex,
 196                                esd->GetMagneticField(),
 197                                100.,dca,cov);
 198     }
 199   }
 200   else{
 201     gPt = track->Pt();
 202     gPx = track->Px();
 203     gPy = track->Py();
 204     gPz = track->Pz();
 205     track->PropagateToDCA(vertex,
 206                           esd->GetMagneticField(),
 207                           100.,dca,cov);
 208   }
 209   dca3D = TMath::Sqrt(TMath::Power(dca[0],2) +
 210                       TMath::Power(dca[1],2));
 211
 212   Int_t  fIdxInt[200];
 213   Int_t nClustersITS = track->GetITSclusters(fIdxInt);
 214   Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
 215
 216   Float_t chi2PerClusterITS = -1;
 217   if (nClustersITS!=0)
 218     chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
 219   Float_t chi2PerClusterTPC = -1;
 220   if (nClustersTPC!=0)
 221     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
 222
 223   Double_t extCov[15];
 224   track->GetExternalCovariance(extCov);
 225
 226   if(fPointOnITSLayer1Flag) {
 227     if(!track->HasPointOnITSLayer(0)) {
 228       if(fDebugMode)
 229         Printf("IsAccepted: Track rejected because it doesn't have a point on the 1st ITS layer");
 230       return kFALSE;
 231     }
 232   }
 233   if(fPointOnITSLayer2Flag) {
 234     if(!track->HasPointOnITSLayer(1)) {
 235       if(fDebugMode)
 236         Printf("IsAccepted: Track rejected because it doesn't have a point on the 2nd ITS layer");
 237       return kFALSE;
 238     }
 239   }
 240   if(fPointOnITSLayer3Flag) {
 241     if(!track->HasPointOnITSLayer(2)) {
 242       if(fDebugMode)
 243         Printf("IsAccepted: Track rejected because it doesn't have a point on the 3rd ITS layer");
 244       return kFALSE;
 245     }
 246   }
 247   if(fPointOnITSLayer4Flag) {
 248     if(!track->HasPointOnITSLayer(3)) {
 249       if(fDebugMode)
 250         Printf("IsAccepted: Track rejected because it doesn't have a point on the 4th ITS layer");
 251       return kFALSE;
 252     }
 253   }
 254   if(fPointOnITSLayer5Flag) {
 255     if(!track->HasPointOnITSLayer(4)) {
 256       if(fDebugMode)
 257         Printf("IsAccepted: Track rejected because it doesn't have a point on the 5th ITS layer");
 258       return kFALSE;
 259     }
 260   }
 261   if(fPointOnITSLayer6Flag) {
 262     if(!track->HasPointOnITSLayer(5)) {
 263       if(fDebugMode)
 264         Printf("IsAccepted: Track rejected because it doesn't have a point on the 6th ITS layer");
 265       return kFALSE;
 266     }
 267   }
 268   if(fMinITSClustersFlag) {
 269     if(nClustersITS < fMinITSClusters) {
 270       if(fDebugMode)
 271         Printf("IsAccepted: Track rejected because it has %d ITS points (min. requested: %d)",nClustersITS,fMinITSClusters);
 272       return kFALSE;
 273     }
 274   }
 275   if(fMaxChi2PerITSClusterFlag) {
 276     if(chi2PerClusterITS > fMaxChi2PerITSCluster) {
 277       if(fDebugMode)
 278         Printf("IsAccepted: Track rejected because it has a chi2 per ITS cluster %lf (max. requested: %lf)",chi2PerClusterITS,fMaxChi2PerITSCluster);
 279       return kFALSE;
 280     }
 281   }
 282   if(fMinTPCClustersFlag) {
 283     if(nClustersTPC < fMinTPCClusters) {
 284       if(fDebugMode)
 285         Printf("IsAccepted: Track rejected because it has %d TPC clusters (min. requested: %d)",nClustersTPC,fMinTPCClusters);
 286       return kFALSE;
 287     }
 288   }
 289   if(fMaxChi2PerTPCClusterFlag) {
 290     if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) {
 291       if(fDebugMode)
 292         Printf("IsAccepted: Track rejected because it has a chi2 per TPC cluster %lf (max. requested: %lf)",chi2PerClusterTPC,fMaxChi2PerTPCCluster);
 293       return kFALSE;
 294     }
 295   }
 296   if(fMaxCov11Flag) {
 297     if(extCov[0] > fMaxCov11) {
 298       if(fDebugMode)
 299         Printf("IsAccepted: Track rejected because it has a cov11 value of %lf (max. requested: %lf)",extCov[0],fMaxCov11);
 300       return kFALSE;
 301     }
 302   }
 303   if(fMaxCov22Flag) {
 304     if(extCov[2] > fMaxCov22) {
 305       if(fDebugMode)
 306         Printf("IsAccepted: Track rejected because it has a cov22 value of %lf (max. requested: %lf)",extCov[2],fMaxCov22);
 307       return kFALSE;
 308     }
 309   }
 310   if(fMaxCov33Flag) {
 311     if(extCov[5] > fMaxCov33) {
 312       if(fDebugMode)
 313         Printf("IsAccepted: Track rejected because it has a cov33 value of %lf (max. requested: %lf)",extCov[5],fMaxCov33);
 314       return kFALSE;
 315     }
 316   }
 317   if(fMaxCov44Flag) {
 318     if(extCov[9] > fMaxCov44) {
 319       if(fDebugMode)
 320         Printf("IsAccepted: Track rejected because it has a cov44 value of %lf (max. requested: %lf)",extCov[9],fMaxCov44);
 321       return kFALSE;
 322     }
 323   }
 324   if(fMaxCov55Flag) {
 325     if(extCov[14] > fMaxCov55) {
 326       if(fDebugMode)
 327         Printf("IsAccepted: Track rejected because it has a cov55 value of %lf (max. requested: %lf)",extCov[14],fMaxCov55);
 328       return kFALSE;
 329     }
 330   }
 331   if(fMaxSigmaToVertexFlag) {
 332     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) {
 333       if(fDebugMode)
 334         Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertex);
 335       return kFALSE;
 336     }
 337   }
 338   if(fMaxSigmaToVertexTPCFlag) {
 339     if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) {
 340       if(fDebugMode)
 341         Printf("IsAccepted: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertexTPC);
 342       return kFALSE;
 343     }
 344   }
 345   if(fMaxDCAXYFlag) {
 346     if(TMath::Abs(dca[0]) > fMaxDCAXY) {
 347       if(fDebugMode)
 348         Printf("IsAccepted: Track rejected because it has a value of dca(xy) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXY);
 349       return kFALSE;
 350     }
 351   }
 352   if(fMaxDCAXYTPCFlag) {
 353     if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) {
 354       if(fDebugMode)
 355         Printf("IsAccepted: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXYTPC);
 356       return kFALSE;
 357     }
 358   }
 359   if(fMaxDCAZFlag) {
 360     if(TMath::Abs(dca[1]) > fMaxDCAZ) {
 361       if(fDebugMode)
 362         Printf("IsAccepted: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZ);
 363       return kFALSE;
 364     }
 365   }
 366   if(fMaxDCAZTPCFlag) {
 367     if(TMath::Abs(dca[1]) > fMaxDCAZTPC) {
 368       if(fDebugMode)
 369         Printf("IsAccepted: Track rejected because it has a value of dca(z) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZTPC);
 370       return kFALSE;
 371     }
 372   }
 373   if(fMaxDCA3DFlag) {
 374     if(TMath::Abs(dca3D) > fMaxDCA3D) {
 375       if(fDebugMode)
 376         Printf("IsAccepted: Track rejected because it has a value of dca(3D) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3D);
 377       return kFALSE;
 378     }
 379   }
 380   if(fMaxDCA3DTPCFlag) {
 381     if(TMath::Abs(dca3D) > fMaxDCA3DTPC)  {
 382       if(fDebugMode)
 383         Printf("IsAccepted: Track rejected because it has a value of dca(3D) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3DTPC);
 384       return kFALSE;
 385     }
 386   }
 387   if(fMaxConstrainChi2Flag) {
 388     if(track->GetConstrainedChi2() > 0)
 389       if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2)  {
 390       if(fDebugMode)
 391         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);
 392       return kFALSE;
 393       }
 394   }
 395   if(fITSRefitFlag) {
 396     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
 397       if(fDebugMode)
 398         Printf("IsAccepted: Track rejected because it has no ITS refit flag");
 399       return kFALSE;
 400     }
 401   }
 402   if(fTPCRefitFlag) {
 403     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
 404       if(fDebugMode)
 405         Printf("IsAccepted: Track rejected because it has no TPC refit flag");
 406       return kFALSE;
 407     }
 408   }
 409   if(fESDpidFlag) {
 410     if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) {
 411       if(fDebugMode)
 412         Printf("IsAccepted: Track rejected because it has no ESD pid flag");
 413       return kFALSE;
 414     }
 415   }
 416   if(fTPCpidFlag) {
 417     if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) {
 418       if(fDebugMode)
 419         Printf("IsAccepted: Track rejected because it has no TPC pid flag");
 420       return kFALSE;
 421     }
 422   }
 423
 424   return kTRUE;
 425 }
 426
 427 //____________________________________________________________________//
 428 Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
 429   // Calculates the number of sigma to the vertex.
 430
 431   Float_t b[2];
 432   Float_t bRes[2];
 433   Float_t bCov[3];
 434   if((kTPC)&&(!kHybrid))
 435     esdTrack->GetImpactParametersTPC(b,bCov);
 436   else
 437     esdTrack->GetImpactParameters(b,bCov);
 438
 439   if (bCov[0]<=0 || bCov[2]<=0) {
 440     //AliDebug(1, "Estimated b resolution lower or equal zero!");
 441     bCov[0]=0; bCov[2]=0;
 442   }
 443   bRes[0] = TMath::Sqrt(bCov[0]);
 444   bRes[1] = TMath::Sqrt(bCov[2]);
 445
 446   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 447
 448   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 449
 450   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 451
 452   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 453
 454   return d;
 455 }
 456
 457 //____________________________________________________________________//
 458 Double_t AliProtonAnalysisBase::Rapidity(Double_t gPx,
 459                                          Double_t gPy,
 460                                          Double_t gPz) const {
 461   //returns the rapidity of the proton - to be removed
 462   Double_t fMass = 9.38270000000000048e-01;
 463
 464   Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) +
 465                            TMath::Power(gPy,2) +
 466                            TMath::Power(gPz,2));
 467   Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
 468   Double_t y = -999;
 469   if(energy != gPz)
 470     y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
 471
 472   return y;
 473 }
 474
 475 //________________________________________________________________________
 476 const AliESDVertex* AliProtonAnalysisBase::GetVertex(AliESDEvent* esd,
 477                                                      AnalysisMode mode,
 478                                                      Double_t gVxMax,
 479                                                      Double_t gVyMax,
 480                                                      Double_t gVzMax) {
 481   // Get the vertex from the ESD and returns it if the vertex is valid
 482   // Second argument decides which vertex is used (this selects
 483   // also the quality criteria that are applied)
 484   const AliESDVertex* vertex = 0;
 485   if(mode == kHybrid)
 486     vertex = esd->GetPrimaryVertexSPD();
 487   else if(mode == kTPC){
 488     Double_t kBz = esd->GetMagneticField();
 489     AliVertexerTracks vertexer(kBz);
 490     vertexer.SetTPCMode();
 491     AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
 492     esd->SetPrimaryVertexTPC(vTPC);
 493     for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
 494       AliESDtrack *t = esd->GetTrack(i);
 495       t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
 496     }
 497     delete vTPC;
 498     vertex = esd->GetPrimaryVertexTPC();
 499   }
 500   else if(mode == kGlobal)
 501     vertex = esd->GetPrimaryVertex();
 502   else
 503     Printf("GetVertex: ERROR: Invalid second argument %d", mode);
 504
 505   if(!vertex) {
 506     if(fDebugMode)
 507       Printf("GetVertex: Event rejected because there is no valid vertex object");
 508     return 0;
 509   }
 510
 511   // check Ncontributors
 512   if(vertex->GetNContributors() <= 0) {
 513     if(fDebugMode)
 514       Printf("GetVertex: Event rejected because the number of contributors for the vertex determination is <= 0");
 515     return 0;
 516   }
 517
 518   // check resolution
 519   Double_t zRes = vertex->GetZRes();
 520   if(zRes == 0) {
 521     if(fDebugMode)
 522       Printf("GetVertex: Event rejected because the value of the vertex resolution in z is 0");
 523     return 0;
 524   }
 525
 526   //check position
 527   if(TMath::Abs(vertex->GetXv()) > gVxMax) {
 528     if(fDebugMode)
 529       Printf("GetVertex: Event rejected because it has a Vx value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetXv()),gVxMax,gVxMax);
 530     return 0;
 531   }
 532   if(TMath::Abs(vertex->GetYv()) > gVyMax)  {
 533     if(fDebugMode)
 534       Printf("GetVertex: Event rejected because it has a Vy value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetYv()),gVyMax,gVyMax);
 535     return 0;
 536   }
 537   if(TMath::Abs(vertex->GetZv()) > gVzMax)  {
 538     if(fDebugMode)
 539       Printf("GetVertex: Event rejected because it has a Vz value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetZv()),gVzMax,gVzMax);
 540     return 0;
 541   }
 542
 543   return vertex;
 544 }
 545
 546 //________________________________________________________________________
 547 Bool_t AliProtonAnalysisBase::IsEventTriggered(const AliESDEvent *esd,
 548                                                TriggerMode trigger) {
 549   // check if the event was triggered
 550   ULong64_t triggerMask = esd->GetTriggerMask();
 551
 552   // definitions from p-p.cfg
 553   ULong64_t spdFO = (1 << 14);
 554   ULong64_t v0left = (1 << 11);
 555   ULong64_t v0right = (1 << 12);
 556
 557   switch (trigger) {
 558   case kMB1: {
 559     if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
 560       return kTRUE;
 561     break;
 562   }
 563   case kMB2: {
 564     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
 565       return kTRUE;
 566     break;
 567   }
 568   case kSPDFASTOR: {
 569     if (triggerMask & spdFO)
 570       return kTRUE;
 571     break;
 572   }
 573   }//switch
 574
 575   return kFALSE;
 576 }
 577
 578 //________________________________________________________________________
 579 TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
 580   // return the list of cuts and their cut values for reference
 581   TLatex l;
 582   l.SetTextAlign(12);
 583   l.SetTextSize(0.04);
 584
 585   TString listOfCuts;
 586
 587   TCanvas *c = new TCanvas("cListOfCuts","List of cuts",0,0,900,600);
 588   c->SetFillColor(10); c->SetHighLightColor(41);
 589   c->Divide(3,2);
 590
 591   c->cd(1)->SetFillColor(10);
 592   l.DrawLatex(0.3,0.9,"Analysis details: List of cuts\n\n");
 593
 594   listOfCuts = "Analysis level: "; listOfCuts += fProtonAnalysisLevel;
 595   l.DrawLatex(0.1,0.82,listOfCuts.Data());
 596   listOfCuts = "Analysis mode: ";
 597   if(fProtonAnalysisMode == kTPC) listOfCuts += "TPC standalone";
 598   if(fProtonAnalysisMode == kHybrid) listOfCuts += "Hybrid TPC";
 599   if(fProtonAnalysisMode == kGlobal) listOfCuts += "Global tracking";
 600   l.DrawLatex(0.1,0.74,listOfCuts.Data());
 601   listOfCuts = "Trigger mode: ";
 602   if(fTriggerMode == kMB1) listOfCuts += "Minimum bias 1";
 603   if(fTriggerMode == kMB2) listOfCuts += "Minimum bias 2";
 604   if(fTriggerMode == kSPDFASTOR) listOfCuts += "FastOR (SPD)";
 605   l.DrawLatex(0.1,0.66,listOfCuts.Data());
 606   listOfCuts = "PID mode: ";
 607   if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
 608   if(fProtonPIDMode == kRatio) listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.})";
 609   if(fProtonPIDMode == kSigma) listOfCuts += "N__{#sigma} area";
 610   l.DrawLatex(0.1,0.58,listOfCuts.Data());
 611   listOfCuts = "Accepted vertex diamond: ";
 612   l.DrawLatex(0.1,0.5,listOfCuts.Data());
 613   listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
 614   l.DrawLatex(0.6,0.5,listOfCuts.Data());
 615   listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
 616   l.DrawLatex(0.6,0.4,listOfCuts.Data());
 617   listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
 618   l.DrawLatex(0.6,0.3,listOfCuts.Data());
 619   listOfCuts = "Phase space: ";
 620   l.DrawLatex(0.1,0.2,listOfCuts.Data());
 621   if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";
 622   else listOfCuts = "|y| < ";
 623   listOfCuts += TMath::Abs(fMinX);
 624   l.DrawLatex(0.35,0.2,listOfCuts.Data());
 625   listOfCuts = "N_{bins} = ";
 626   listOfCuts += fNBinsX; listOfCuts += " (binning: ";
 627   listOfCuts += (fMaxX-fMinX)/fNBinsX; listOfCuts += ")";
 628   l.DrawLatex(0.35,0.15,listOfCuts.Data());
 629   listOfCuts = ""; listOfCuts += fMinY; listOfCuts += " < P_{T} < ";
 630   listOfCuts += fMaxY; listOfCuts += "GeV/c";
 631   l.DrawLatex(0.35,0.1,listOfCuts.Data());
 632   listOfCuts = "N_{bins} = ";
 633   listOfCuts += fNBinsY; listOfCuts += " (binning: ";
 634   listOfCuts += (fMaxY-fMinY)/fNBinsY; listOfCuts += ")";
 635   l.DrawLatex(0.35,0.05,listOfCuts.Data());
 636
 637   c->cd(2)->SetFillColor(2);
 638   l.DrawLatex(0.3,0.9,"ITS related cuts");
 639   listOfCuts = "Request a cluster on SPD1: ";
 640   listOfCuts += fPointOnITSLayer1Flag;
 641   l.DrawLatex(0.1,0.8,listOfCuts.Data());
 642   listOfCuts = "Request a cluster on SPD2: ";
 643   listOfCuts += fPointOnITSLayer2Flag;
 644   l.DrawLatex(0.1,0.7,listOfCuts.Data());
 645   listOfCuts = "Request a cluster on SDD1: ";
 646   listOfCuts += fPointOnITSLayer3Flag;
 647   l.DrawLatex(0.1,0.6,listOfCuts.Data());
 648   listOfCuts = "Request a cluster on SDD2: ";
 649   listOfCuts += fPointOnITSLayer4Flag;
 650   l.DrawLatex(0.1,0.5,listOfCuts.Data());
 651   listOfCuts = "Request a cluster on SSD1: ";
 652   listOfCuts += fPointOnITSLayer5Flag;
 653   l.DrawLatex(0.1,0.4,listOfCuts.Data());
 654   listOfCuts = "Request a cluster on SSD2: ";
 655   listOfCuts += fPointOnITSLayer6Flag;
 656   l.DrawLatex(0.1,0.3,listOfCuts.Data());
 657   listOfCuts = "Minimum number of ITS clusters: ";
 658   if(fMinITSClustersFlag) listOfCuts += fMinITSClusters;
 659   else listOfCuts += "Not used";
 660   l.DrawLatex(0.1,0.2,listOfCuts.Data());
 661   listOfCuts = "Maximum #chi^{2} per ITS cluster: ";
 662   if(fMaxChi2PerITSClusterFlag) listOfCuts += fMaxChi2PerITSCluster;
 663   else listOfCuts += "Not used";
 664   l.DrawLatex(0.1,0.1,listOfCuts.Data());
 665
 666   c->cd(3)->SetFillColor(3);
 667   l.DrawLatex(0.3,0.9,"TPC related cuts");
 668   listOfCuts = "Minimum number of TPC clusters: ";
 669   if(fMinTPCClustersFlag) listOfCuts += fMinTPCClusters;
 670   else listOfCuts += "Not used";
 671   l.DrawLatex(0.1,0.7,listOfCuts.Data());
 672   listOfCuts = "Maximum #chi^{2} per TPC cluster: ";
 673   if(fMaxChi2PerTPCClusterFlag) listOfCuts += fMaxChi2PerTPCCluster;
 674   else listOfCuts += "Not used";
 675   l.DrawLatex(0.1,0.5,listOfCuts.Data());
 676
 677   c->cd(4)->SetFillColor(4);
 678   l.DrawLatex(0.3,0.9,"Tracking related cuts");
 679   listOfCuts = "Maximum cov11: ";
 680   if(fMaxCov11Flag) listOfCuts += fMaxCov11;
 681   else listOfCuts += "Not used";
 682   l.DrawLatex(0.1,0.75,listOfCuts.Data());
 683   listOfCuts = "Maximum cov22: ";
 684   if(fMaxCov22Flag) listOfCuts += fMaxCov22;
 685   else listOfCuts += "Not used";
 686   l.DrawLatex(0.1,0.6,listOfCuts.Data());
 687   listOfCuts = "Maximum cov33: ";
 688   if(fMaxCov33Flag) listOfCuts += fMaxCov33;
 689   else listOfCuts += "Not used";
 690   l.DrawLatex(0.1,0.45,listOfCuts.Data());
 691   listOfCuts = "Maximum cov44: ";
 692   if(fMaxCov44Flag) listOfCuts += fMaxCov44;
 693   else listOfCuts += "Not used";
 694   l.DrawLatex(0.1,0.3,listOfCuts.Data());
 695   listOfCuts = "Maximum cov55: ";
 696   if(fMaxCov55Flag) listOfCuts += fMaxCov55;
 697   else listOfCuts += "Not used";
 698   l.DrawLatex(0.1,0.15,listOfCuts.Data());
 699
 700   c->cd(5)->SetFillColor(5);
 701   l.DrawLatex(0.3,0.9,"DCA related cuts");
 702   listOfCuts = "Maximum sigma to vertex: ";
 703   if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
 704   else listOfCuts += "Not used";
 705   l.DrawLatex(0.1,0.8,listOfCuts.Data());
 706   listOfCuts = "Maximum sigma to vertex (TPC): ";
 707   if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
 708   else listOfCuts += "Not used";
 709   l.DrawLatex(0.1,0.72,listOfCuts.Data());
 710   listOfCuts = "Maximum DCA in xy: ";
 711   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
 712   else listOfCuts += "Not used";
 713   l.DrawLatex(0.1,0.64,listOfCuts.Data());
 714   listOfCuts = "Maximum DCA in z: ";
 715   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
 716   else listOfCuts += "Not used";
 717   l.DrawLatex(0.1,0.56,listOfCuts.Data());
 718   listOfCuts = "Maximum DCA in 3D: ";
 719   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
 720   else listOfCuts += "Not used";
 721   l.DrawLatex(0.1,0.48,listOfCuts.Data());
 722   listOfCuts = "Maximum DCA in xy (TPC): ";
 723   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
 724   else listOfCuts += "Not used";
 725   l.DrawLatex(0.1,0.4,listOfCuts.Data());
 726   listOfCuts = "Maximum DCA in z (TPC): ";
 727   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
 728   else listOfCuts += "Not used";
 729   l.DrawLatex(0.1,0.32,listOfCuts.Data());
 730   listOfCuts = "Maximum DCA in 3D (TPC): ";
 731   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
 732   else listOfCuts += "Not used";
 733   l.DrawLatex(0.1,0.24,listOfCuts.Data());
 734   listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
 735   if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
 736   else listOfCuts += "Not used";
 737   l.DrawLatex(0.1,0.16,listOfCuts.Data());
 738
 739   c->cd(6)->SetFillColor(6);
 740   l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
 741   listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
 742   l.DrawLatex(0.1,0.7,listOfCuts.Data());
 743   listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag;
 744   l.DrawLatex(0.1,0.5,listOfCuts.Data());
 745   listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
 746   l.DrawLatex(0.1,0.3,listOfCuts.Data());
 747   listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
 748   l.DrawLatex(0.1,0.1,listOfCuts.Data());
 749
 750   return c;
 751 }
 752
 753 //________________________________________________________________________
 754 Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
 755   //Function that checks if a track is a proton
 756   Double_t probability[5];
 757   Double_t gPt = 0.0, gP = 0.0;
 758   Long64_t fParticleType = 0;
 759
 760   //Bayesian approach for the PID
 761   if(fProtonPIDMode == kBayesian) {
 762     if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
 763       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 764       if(tpcTrack) {
 765         gPt = tpcTrack->Pt();
 766         gP = tpcTrack->P();
 767         track->GetTPCpid(probability);
 768       }
 769     }//TPC standalone or Hybrid TPC
 770     else if(fProtonAnalysisMode == kGlobal) {
 771       gPt = track->Pt();
 772       gP = track->P();
 773       track->GetESDpid(probability);
 774     }//Global tracking
 775
 776     Double_t rcc = 0.0;
 777     for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 778       rcc += probability[i]*GetParticleFraction(i,gP);
 779     if(rcc != 0.0) {
 780       Double_t w[5];
 781       for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 782         w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
 783       fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
 784     }
 785     if(fParticleType == 4)
 786       return kTRUE;
 787   }
 788   //Ratio of the measured over the theoretical dE/dx a la STAR
 789   else if(fProtonPIDMode == kRatio) {
 790     Printf("The kRatio mode is not implemented yet!!!");
 791     return kFALSE;
 792   }
 793   //Definition of an N-sigma area around the dE/dx vs P band
 794   else if(fProtonPIDMode == kSigma) {
 795     Printf("The kSigma mode is not implemented yet!!!");
 796     return kFALSE;
 797   }
 798
 799   return kFALSE;
 800 }
 801
 802
 803