PWGLF/SPECTRA/AntiprotonToProton/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 <AliESDpid.h>
  33 #include <AliTPCPIDResponse.h>
  34 class AliLog;
  35 class AliESDVertex;
  36
  37 #include "AliProtonAnalysisBase.h"
  38
  39 ClassImp(AliProtonAnalysisBase)
  40
  41 //____________________________________________________________________//
  42 AliProtonAnalysisBase::AliProtonAnalysisBase() :
  43   TObject(),  fProtonAnalysisLevel("ESD"), fAnalysisMC(kFALSE),
  44   fTriggerMode(kMB2), kUseOnlineTrigger(kFALSE), kUseOfflineTrigger(kFALSE),
  45   fPhysicsSelection(0),
  46   fProtonAnalysisMode(kTPC), fProtonPIDMode(kBayesian),
  47   fAnalysisEtaMode(kFALSE),
  48   fRunQAAnalysis(kFALSE),
  49   fVxMax(100.), fVyMax(100.), fVzMax(100.), fMinNumOfContributors(0),
  50   fNBinsX(0), fMinX(0.), fMaxX(0.),
  51   fNBinsY(0), fMinY(0.), fMaxY(0.),
  52   fMinTPCClusters(0), fMinITSClusters(0),
  53   fMaxChi2PerTPCCluster(0.), fMaxChi2PerITSCluster(0.),
  54   fMaxCov11(0.), fMaxCov22(0.), fMaxCov33(0.), fMaxCov44(0.), fMaxCov55(0.),
  55   fMaxSigmaToVertex(0.), fMaxSigmaToVertexTPC(0.),
  56   fMaxDCAXY(0.), fMaxDCAXYTPC(0.),
  57   fMaxDCAZ(0.), fMaxDCAZTPC(0.),
  58   fMaxDCA3D(0.), fMaxDCA3DTPC(0.),
  59   fMaxConstrainChi2(0.), fMinTPCdEdxPoints(0),
  60   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  61   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  62   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE),
  63   fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
  64   fMaxSigmaToVertexFlag(kFALSE), fMaxSigmaToVertexTPCFlag(kFALSE),
  65   fMaxDCAXYFlag(kFALSE), fMaxDCAXYTPCFlag(kFALSE),
  66   fMaxDCAZFlag(kFALSE), fMaxDCAZTPCFlag(kFALSE),
  67   fMaxDCA3DFlag(kFALSE), fMaxDCA3DTPCFlag(kFALSE),
  68   fMaxConstrainChi2Flag(kFALSE),
  69   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
  70   fESDpidFlag(kFALSE), fTPCpidFlag(kFALSE), fTOFpidFlag(kFALSE),
  71   fPointOnSPDLayersFlag(0), fPointOnSDDLayersFlag(0), fPointOnSSDLayersFlag(0),
  72   fPointOnITSLayer1Flag(0), fPointOnITSLayer2Flag(0),
  73   fPointOnITSLayer3Flag(0), fPointOnITSLayer4Flag(0),
  74   fPointOnITSLayer5Flag(0), fPointOnITSLayer6Flag(0),
  75   fMinTPCdEdxPointsFlag(kFALSE),
  76   fPtDependentDcaXY(0), fPtDependentDcaXYFlag(kFALSE), fNSigmaDCAXY(0),
  77   fFunctionProbabilityFlag(kFALSE),
  78   fNSigma(0), fNRatio(0.),
  79   fElectronFunction(0), fMuonFunction(0),
  80   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
  81   fDebugMode(kFALSE), fListVertexQA(0) {
  82   //Default constructor
  83   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
  84   /*for(Int_t i = 0; i < 24; i++) {
  85     fdEdxMean[i] = 0.0;
  86     fdEdxSigma[i] = 0.0;
  87     }*/
  88
  89   fListVertexQA = new TList();
  90   fListVertexQA->SetName("fListVertexQA");
  91   TH1F *gHistVx = new TH1F("gHistVx",
  92                            "Vx distribution;V_{x} [cm];Entries",
  93                            100,-5.,5.);
  94   gHistVx->SetFillColor(kRed-2);
  95   fListVertexQA->Add(gHistVx);
  96   TH1F *gHistVxAccepted = new TH1F("gHistVxaccepted",
  97                                    "Vx distribution;V_{x} [cm];Entries",
  98                                    100,-5.,5.);
  99   fListVertexQA->Add(gHistVxAccepted);
 100   TH1F *gHistVy = new TH1F("gHistVy",
 101                            "Vy distribution;V_{y} [cm];Entries",
 102                            100,-5.,5.);
 103   gHistVy->SetFillColor(kRed-2);
 104   fListVertexQA->Add(gHistVy);
 105   TH1F *gHistVyAccepted = new TH1F("gHistVyaccepted",
 106                                    "Vy distribution;V_{y} [cm];Entries",
 107                                    100,-5.,5.);
 108   fListVertexQA->Add(gHistVyAccepted);
 109   TH1F *gHistVz = new TH1F("gHistVz",
 110                            "Vz distribution;V_{z} [cm];Entries",
 111                            100,-25.,25.);
 112   gHistVz->SetFillColor(kRed-2);
 113   fListVertexQA->Add(gHistVz);
 114   TH1F *gHistVzAccepted = new TH1F("gHistVzaccepted",
 115                                    "Vz distribution;V_{z} [cm];Entries",
 116                                    100,-25.,25.);
 117   fListVertexQA->Add(gHistVzAccepted);
 118
 119   TH1F *gHistNumberOfContributors = new TH1F("gHistNumberOfContributors",
 120                                              "Number of contributors;N_{contr.};Entries",
 121                                              100,0.,100.);
 122   fListVertexQA->Add(gHistNumberOfContributors);
 123
 124
 125 }
 126
 127 //____________________________________________________________________//
 128 AliProtonAnalysisBase::~AliProtonAnalysisBase() {
 129   //Default destructor
 130   if(fElectronFunction) delete fElectronFunction;
 131   if(fMuonFunction) delete fMuonFunction;
 132   if(fPionFunction) delete fPionFunction;
 133   if(fKaonFunction) delete fKaonFunction;
 134   if(fProtonFunction) delete fProtonFunction;
 135   if(fListVertexQA) delete fListVertexQA;
 136   if(fPtDependentDcaXY) delete fPtDependentDcaXY;
 137   if(fPhysicsSelection) delete fPhysicsSelection;
 138 }
 139
 140 //____________________________________________________________________//
 141 Double_t AliProtonAnalysisBase::GetParticleFraction(Int_t i, Double_t p) {
 142   //Return the a priori probs
 143   Double_t partFrac=0;
 144   if(fFunctionProbabilityFlag) {
 145     if(i == 0) partFrac = fElectronFunction->Eval(p);
 146     if(i == 1) partFrac = fMuonFunction->Eval(p);
 147     if(i == 2) partFrac = fPionFunction->Eval(p);
 148     if(i == 3) partFrac = fKaonFunction->Eval(p);
 149     if(i == 4) partFrac = fProtonFunction->Eval(p);
 150   }
 151   else partFrac = fPartFrac[i];
 152
 153   return partFrac;
 154 }
 155
 156 //____________________________________________________________________//
 157 Bool_t AliProtonAnalysisBase::IsInPhaseSpace(AliESDtrack* const track) {
 158   // Checks if the track is outside the analyzed y-Pt phase space
 159   Double_t gP = 0.0, gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 160   Double_t eta = 0.0;
 161
 162   if((fProtonAnalysisMode == kTPC) || (fProtonAnalysisMode == kHybrid) || (fProtonAnalysisMode == kFullHybrid)) {
 163     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 164     if(!tpcTrack) {
 165       gP = 0.0; gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0; eta = -10.0;
 166     }
 167     else {
 168       gP = tpcTrack->P();
 169       gPt = tpcTrack->Pt();
 170       gPx = tpcTrack->Px();
 171       gPy = tpcTrack->Py();
 172       gPz = tpcTrack->Pz();
 173       eta = tpcTrack->Eta();
 174     }
 175   }//standalone TPC or Hybrid TPC approaches
 176   else {
 177     gP = track->P();
 178     gPt = track->Pt();
 179     gPx = track->Px();
 180     gPy = track->Py();
 181     gPz = track->Pz();
 182     eta = track->Eta();
 183   }
 184
 185   if((gPt < fMinY) || (gPt > fMaxY)) {
 186       if(fDebugMode)
 187         Printf("IsInPhaseSpace: Track rejected because it has a Pt value of %lf (accepted interval: %lf - %lf)",gPt,fMinY,fMaxY);
 188       return kFALSE;
 189   }
 190   if((gP < fMinY) || (gP > fMaxY)) {
 191       if(fDebugMode)
 192         Printf("IsInPhaseSpace: Track rejected because it has a P value of %lf (accepted interval: %lf - %lf)",gP,fMinY,fMaxY);
 193       return kFALSE;
 194   }
 195   if(fAnalysisEtaMode) {
 196     if((eta < fMinX) || (eta > fMaxX)) {
 197       if(fDebugMode)
 198         Printf("IsInPhaseSpace: Track rejected because it has an eta value of %lf (accepted interval: %lf - %lf)",eta,fMinX,fMaxX);
 199       return kFALSE;
 200     }
 201   }
 202   else {
 203     if((Rapidity(gPx,gPy,gPz) < fMinX) || (Rapidity(gPx,gPy,gPz) > fMaxX)) {
 204       if(fDebugMode)
 205         Printf("IsInPhaseSpace: Track rejected because it has a y value of %lf (accepted interval: %lf - %lf)",Rapidity(gPx,gPy,gPz),fMinX,fMaxX);
 206       return kFALSE;
 207     }
 208   }
 209
 210   return kTRUE;
 211 }
 212
 213 //____________________________________________________________________//
 214 Bool_t AliProtonAnalysisBase::IsAccepted(AliESDtrack* track) {
 215   // Checks if the track is excluded from the cuts
 216   //Int_t  fIdxInt[200];
 217   //Int_t nClustersITS = track->GetITSclusters(fIdxInt);
 218   //Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
 219   Int_t nClustersITS = track->GetITSclusters(0x0);
 220   Int_t nClustersTPC = track->GetTPCclusters(0x0);
 221
 222   Float_t chi2PerClusterITS = -1;
 223   if (nClustersITS!=0)
 224     chi2PerClusterITS = track->GetITSchi2()/Float_t(nClustersITS);
 225   Float_t chi2PerClusterTPC = -1;
 226   if (nClustersTPC!=0)
 227     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
 228
 229   Double_t extCov[15];
 230   track->GetExternalCovariance(extCov);
 231
 232   if(fPointOnSPDLayersFlag) {
 233     if((!track->HasPointOnITSLayer(0))&&(!track->HasPointOnITSLayer(1))) {
 234       if(fDebugMode)
 235         Printf("IsAccepted: Track rejected because it doesn't have a point on either SPD layers");
 236       return kFALSE;
 237     }
 238   }
 239   if(fPointOnSDDLayersFlag) {
 240     if((!track->HasPointOnITSLayer(2))&&(!track->HasPointOnITSLayer(3))) {
 241       if(fDebugMode)
 242         Printf("IsAccepted: Track rejected because it doesn't have a point on either SDD layers");
 243       return kFALSE;
 244     }
 245   }
 246   if(fPointOnSSDLayersFlag) {
 247     if((!track->HasPointOnITSLayer(4))&&(!track->HasPointOnITSLayer(5))) {
 248       if(fDebugMode)
 249         Printf("IsAccepted: Track rejected because it doesn't have a point on either SSD layers");
 250       return kFALSE;
 251     }
 252   }
 253   if(fPointOnITSLayer1Flag) {
 254     if(!track->HasPointOnITSLayer(0)) {
 255       if(fDebugMode)
 256         Printf("IsAccepted: Track rejected because it doesn't have a point on the 1st ITS layer");
 257       return kFALSE;
 258     }
 259   }
 260   if(fPointOnITSLayer2Flag) {
 261     if(!track->HasPointOnITSLayer(1)) {
 262       if(fDebugMode)
 263         Printf("IsAccepted: Track rejected because it doesn't have a point on the 2nd ITS layer");
 264       return kFALSE;
 265     }
 266   }
 267   if(fPointOnITSLayer3Flag) {
 268     if(!track->HasPointOnITSLayer(2)) {
 269       if(fDebugMode)
 270         Printf("IsAccepted: Track rejected because it doesn't have a point on the 3rd ITS layer");
 271       return kFALSE;
 272     }
 273   }
 274   if(fPointOnITSLayer4Flag) {
 275     if(!track->HasPointOnITSLayer(3)) {
 276       if(fDebugMode)
 277         Printf("IsAccepted: Track rejected because it doesn't have a point on the 4th ITS layer");
 278       return kFALSE;
 279     }
 280   }
 281   if(fPointOnITSLayer5Flag) {
 282     if(!track->HasPointOnITSLayer(4)) {
 283       if(fDebugMode)
 284         Printf("IsAccepted: Track rejected because it doesn't have a point on the 5th ITS layer");
 285       return kFALSE;
 286     }
 287   }
 288   if(fPointOnITSLayer6Flag) {
 289     if(!track->HasPointOnITSLayer(5)) {
 290       if(fDebugMode)
 291         Printf("IsAccepted: Track rejected because it doesn't have a point on the 6th ITS layer");
 292       return kFALSE;
 293     }
 294   }
 295   if(fMinITSClustersFlag) {
 296     if(nClustersITS < fMinITSClusters) {
 297       if(fDebugMode)
 298         Printf("IsAccepted: Track rejected because it has %d ITS points (min. requested: %d)",nClustersITS,fMinITSClusters);
 299       return kFALSE;
 300     }
 301   }
 302   if(fMaxChi2PerITSClusterFlag) {
 303     if(chi2PerClusterITS > fMaxChi2PerITSCluster) {
 304       if(fDebugMode)
 305         Printf("IsAccepted: Track rejected because it has a chi2 per ITS cluster %lf (max. requested: %lf)",chi2PerClusterITS,fMaxChi2PerITSCluster);
 306       return kFALSE;
 307     }
 308   }
 309   if(fMinTPCClustersFlag) {
 310     if(nClustersTPC < fMinTPCClusters) {
 311       if(fDebugMode)
 312         Printf("IsAccepted: Track rejected because it has %d TPC clusters (min. requested: %d)",nClustersTPC,fMinTPCClusters);
 313       return kFALSE;
 314     }
 315   }
 316   if(fMaxChi2PerTPCClusterFlag) {
 317     if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) {
 318       if(fDebugMode)
 319         Printf("IsAccepted: Track rejected because it has a chi2 per TPC cluster %lf (max. requested: %lf)",chi2PerClusterTPC,fMaxChi2PerTPCCluster);
 320       return kFALSE;
 321     }
 322   }
 323   if(fMaxCov11Flag) {
 324     if(extCov[0] > fMaxCov11) {
 325       if(fDebugMode)
 326         Printf("IsAccepted: Track rejected because it has a cov11 value of %lf (max. requested: %lf)",extCov[0],fMaxCov11);
 327       return kFALSE;
 328     }
 329   }
 330   if(fMaxCov22Flag) {
 331     if(extCov[2] > fMaxCov22) {
 332       if(fDebugMode)
 333         Printf("IsAccepted: Track rejected because it has a cov22 value of %lf (max. requested: %lf)",extCov[2],fMaxCov22);
 334       return kFALSE;
 335     }
 336   }
 337   if(fMaxCov33Flag) {
 338     if(extCov[5] > fMaxCov33) {
 339       if(fDebugMode)
 340         Printf("IsAccepted: Track rejected because it has a cov33 value of %lf (max. requested: %lf)",extCov[5],fMaxCov33);
 341       return kFALSE;
 342     }
 343   }
 344   if(fMaxCov44Flag) {
 345     if(extCov[9] > fMaxCov44) {
 346       if(fDebugMode)
 347         Printf("IsAccepted: Track rejected because it has a cov44 value of %lf (max. requested: %lf)",extCov[9],fMaxCov44);
 348       return kFALSE;
 349     }
 350   }
 351   if(fMaxCov55Flag) {
 352     if(extCov[14] > fMaxCov55) {
 353       if(fDebugMode)
 354         Printf("IsAccepted: Track rejected because it has a cov55 value of %lf (max. requested: %lf)",extCov[14],fMaxCov55);
 355       return kFALSE;
 356     }
 357   }
 358   if(fMinTPCdEdxPointsFlag) {
 359     if(track->GetTPCsignalN() < fMinTPCdEdxPoints) {
 360       if(fDebugMode)
 361         Printf("IsAccepted: Track rejected because it has %d TPC points for the calculation of the energy loss (min. requested: %d)",track->GetTPCsignalN(),fMinTPCdEdxPoints);
 362       return kFALSE;
 363     }
 364   }
 365   if(fITSRefitFlag) {
 366     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) {
 367       if(fDebugMode)
 368         Printf("IsAccepted: Track rejected because it has no ITS refit flag");
 369       return kFALSE;
 370     }
 371   }
 372   if(fTPCRefitFlag) {
 373     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
 374       if(fDebugMode)
 375         Printf("IsAccepted: Track rejected because it has no TPC refit flag");
 376       return kFALSE;
 377     }
 378   }
 379   if(fESDpidFlag) {
 380     if ((track->GetStatus() & AliESDtrack::kESDpid) == 0) {
 381       if(fDebugMode)
 382         Printf("IsAccepted: Track rejected because it has no ESD pid flag");
 383       return kFALSE;
 384     }
 385   }
 386   if(fTPCpidFlag) {
 387     if ((track->GetStatus() & AliESDtrack::kTPCpid) == 0) {
 388       if(fDebugMode)
 389         Printf("IsAccepted: Track rejected because it has no TPC pid flag");
 390       return kFALSE;
 391     }
 392   }
 393   if(fTOFpidFlag) {
 394     if ((track->GetStatus() & AliESDtrack::kTOFpid) == 0) {
 395       if(fDebugMode)
 396         Printf("IsAccepted: Track rejected because it has no TOF pid flag");
 397       return kFALSE;
 398     }
 399   }
 400
 401   return kTRUE;
 402 }
 403
 404 //____________________________________________________________________//
 405 void AliProtonAnalysisBase::SetPtDependentDCAxy(Int_t nSigma,
 406                                                 Double_t p0,
 407                                                 Double_t p1,
 408                                                 Double_t p2) {
 409   //Pt dependent dca cut in xy
 410   fNSigmaDCAXY = nSigma;
 411   fPtDependentDcaXY = new TF1("fPtDependentDcaXY","[0]+[1]/x^[2]",0.1,10.1);
 412   fPtDependentDcaXY->SetParameter(0,p0);
 413   fPtDependentDcaXY->SetParameter(1,p1);
 414   fPtDependentDcaXY->SetParameter(2,p2);
 415   fPtDependentDcaXYFlag = kTRUE;
 416 }
 417
 418 //____________________________________________________________________//
 419 Bool_t AliProtonAnalysisBase::IsPrimary(AliESDEvent *esd,
 420                                         const AliESDVertex *vertex,
 421                                         AliESDtrack* track) {
 422   // Checks if the track is a primary-like candidate
 423   const Double_t kMicrometer2Centimeter = 0.0001;
 424   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 425   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
 426   Double_t dca3D = 0.0;
 427   Float_t dcaXY = 0.0, dcaZ = 0.0;
 428
 429   if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)) {
 430     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 431     if(!tpcTrack) {
 432       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
 433       dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
 434       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
 435     }
 436     else {
 437       gPt = tpcTrack->Pt();
 438       gPx = tpcTrack->Px();
 439       gPy = tpcTrack->Py();
 440       gPz = tpcTrack->Pz();
 441       tpcTrack->PropagateToDCA(vertex,
 442                                esd->GetMagneticField(),
 443                                100.,dca,cov);
 444     }
 445   }//standalone TPC or hybrid TPC approaches
 446   if(fProtonAnalysisMode == kFullHybrid) {
 447     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 448     AliExternalTrackParam cParam;
 449     if(!tpcTrack) {
 450       gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
 451       dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
 452       cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
 453     }
 454     else {
 455       gPt = tpcTrack->Pt();
 456       gPx = tpcTrack->Px();
 457       gPy = tpcTrack->Py();
 458       gPz = tpcTrack->Pz();
 459       track->RelateToVertex(vertex,
 460                             esd->GetMagneticField(),
 461                             100.,&cParam);
 462       track->GetImpactParameters(dcaXY,dcaZ);
 463       dca[0] = dcaXY; dca[1] = dcaZ;
 464     }
 465   }//standalone TPC or hybrid TPC approaches
 466   else {
 467     gPt = track->Pt();
 468     gPx = track->Px();
 469     gPy = track->Py();
 470     gPz = track->Pz();
 471     track->PropagateToDCA(vertex,
 472                           esd->GetMagneticField(),
 473                           100.,dca,cov);
 474   }
 475   dca3D = TMath::Sqrt(TMath::Power(dca[0],2) +
 476                       TMath::Power(dca[1],2));
 477
 478   if(fMaxSigmaToVertexFlag) {
 479     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) {
 480       if(fDebugMode)
 481         Printf("IsPrimary: Track rejected because it has a %lf sigmas to vertex (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertex);
 482       return kFALSE;
 483     }
 484   }
 485   if(fMaxSigmaToVertexTPCFlag) {
 486     if(GetSigmaToVertex(track) > fMaxSigmaToVertexTPC) {
 487       if(fDebugMode)
 488         Printf("IsPrimary: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",GetSigmaToVertex(track),fMaxSigmaToVertexTPC);
 489       return kFALSE;
 490     }
 491   }
 492   if(fMaxDCAXYFlag) {
 493     if(TMath::Abs(dca[0]) > fMaxDCAXY) {
 494       if(fDebugMode)
 495         Printf("IsPrimary: Track rejected because it has a value of dca(xy) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXY);
 496       return kFALSE;
 497     }
 498   }
 499   if(fMaxDCAXYTPCFlag) {
 500     if(TMath::Abs(dca[0]) > fMaxDCAXYTPC) {
 501       if(fDebugMode)
 502         Printf("IsPrimary: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[0]),fMaxDCAXYTPC);
 503       return kFALSE;
 504     }
 505   }
 506   if(fMaxDCAZFlag) {
 507     if(TMath::Abs(dca[1]) > fMaxDCAZ) {
 508       if(fDebugMode)
 509         Printf("IsPrimary: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZ);
 510       return kFALSE;
 511     }
 512   }
 513   if(fMaxDCAZTPCFlag) {
 514     if(TMath::Abs(dca[1]) > fMaxDCAZTPC) {
 515       if(fDebugMode)
 516         Printf("IsPrimary: Track rejected because it has a value of dca(z) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca[1]),fMaxDCAZTPC);
 517       return kFALSE;
 518     }
 519   }
 520   if(fMaxDCA3DFlag) {
 521     if(TMath::Abs(dca3D) > fMaxDCA3D) {
 522       if(fDebugMode)
 523         Printf("IsPrimary: Track rejected because it has a value of dca(3D) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3D);
 524       return kFALSE;
 525     }
 526   }
 527   if(fMaxDCA3DTPCFlag) {
 528     if(TMath::Abs(dca3D) > fMaxDCA3DTPC)  {
 529       if(fDebugMode)
 530         Printf("IsPrimary: Track rejected because it has a value of dca(3D) (TPC) of %lf (max. requested: %lf)",TMath::Abs(dca3D),fMaxDCA3DTPC);
 531       return kFALSE;
 532     }
 533   }
 534   if(fMaxConstrainChi2Flag) {
 535     if(track->GetConstrainedChi2() > 0)
 536       if(TMath::Log(track->GetConstrainedChi2()) > fMaxConstrainChi2)  {
 537       if(fDebugMode)
 538         Printf("IsPrimary: Track rejected because it has a value of the constrained chi2 to the vertex of %lf (max. requested: %lf)",TMath::Log(track->GetConstrainedChi2()),fMaxConstrainChi2);
 539       return kFALSE;
 540       }
 541   }
 542   if(fPtDependentDcaXYFlag) {
 543     if(TMath::Abs(dca[0]) > kMicrometer2Centimeter*fNSigmaDCAXY*fPtDependentDcaXY->Eval(gPt)) {
 544       if(fDebugMode)
 545         Printf("IsPrimary: Track rejected because it has a value of the dca(xy) higher than the %d sigma pt dependent cut: %lf (max. requested: %lf)",fNSigmaDCAXY,TMath::Abs(dca[0]),fNSigmaDCAXY*fPtDependentDcaXY->Eval(gPt));
 546       return kFALSE;
 547     }
 548   }
 549
 550   return kTRUE;
 551 }
 552
 553 //____________________________________________________________________//
 554 Float_t AliProtonAnalysisBase::GetSigmaToVertex(AliESDtrack* esdTrack) const {
 555   // Calculates the number of sigma to the vertex.
 556   Float_t b[2];
 557   Float_t bRes[2];
 558   Float_t bCov[3];
 559   if((fProtonAnalysisMode == kTPC)&&(fProtonAnalysisMode != kHybrid)&&(fProtonAnalysisMode != kFullHybrid))
 560     esdTrack->GetImpactParametersTPC(b,bCov);
 561   else
 562     esdTrack->GetImpactParameters(b,bCov);
 563
 564   if (bCov[0]<=0 || bCov[2]<=0) {
 565     //AliDebug(1, "Estimated b resolution lower or equal zero!");
 566     bCov[0]=0; bCov[2]=0;
 567   }
 568   bRes[0] = TMath::Sqrt(bCov[0]);
 569   bRes[1] = TMath::Sqrt(bCov[2]);
 570
 571   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 572
 573   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 574
 575   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 576
 577   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 578
 579   return d;
 580 }
 581
 582 //____________________________________________________________________//
 583 Double_t AliProtonAnalysisBase::Rapidity(Double_t gPx,
 584                                          Double_t gPy,
 585                                          Double_t gPz) const {
 586   //returns the rapidity of the proton - to be removed
 587   Double_t fMass = 9.38270000000000048e-01;
 588
 589   Double_t gP = TMath::Sqrt(TMath::Power(gPx,2) +
 590                            TMath::Power(gPy,2) +
 591                            TMath::Power(gPz,2));
 592   Double_t energy = TMath::Sqrt(gP*gP + fMass*fMass);
 593   Double_t y = -999;
 594   if(energy != gPz)
 595     y = 0.5*TMath::Log((energy + gPz)/(energy - gPz));
 596
 597   return y;
 598 }
 599
 600 //________________________________________________________________________
 601 const AliESDVertex* AliProtonAnalysisBase::GetVertex(AliESDEvent* esd,
 602                                                      AnalysisMode mode,
 603                                                      Double_t gVxMax,
 604                                                      Double_t gVyMax,
 605                                                      Double_t gVzMax) {
 606   // Get the vertex from the ESD and returns it if the vertex is valid
 607   // Second argument decides which vertex is used (this selects
 608   // also the quality criteria that are applied)
 609   const AliESDVertex* vertex = 0;
 610   if((mode == kHybrid)||(mode == kFullHybrid))
 611     vertex = esd->GetPrimaryVertexSPD();
 612   else if(mode == kTPC){
 613     Double_t kBz = esd->GetMagneticField();
 614     AliVertexerTracks vertexer(kBz);
 615     vertexer.SetTPCMode();
 616     AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
 617     esd->SetPrimaryVertexTPC(vTPC);
 618     for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
 619       AliESDtrack *t = esd->GetTrack(i);
 620       t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
 621     }
 622     delete vTPC;
 623     vertex = esd->GetPrimaryVertexTPC();
 624   }
 625   else if(mode == kGlobal)
 626     vertex = esd->GetPrimaryVertex();
 627   else
 628     Printf("GetVertex: ERROR: Invalid second argument %d", mode);
 629
 630   if(!vertex) {
 631     if(fDebugMode)
 632       Printf("GetVertex: Event rejected because there is no valid vertex object");
 633     return 0;
 634   }
 635
 636   // check Ncontributors
 637   if(vertex->GetNContributors() <= 0) {
 638     if(fDebugMode)
 639       Printf("GetVertex: Event rejected because the number of contributors for the vertex determination is <= 0");
 640     return 0;
 641   }
 642
 643   // check resolution
 644   Double_t zRes = vertex->GetZRes();
 645   if(zRes == 0) {
 646     if(fDebugMode)
 647       Printf("GetVertex: Event rejected because the value of the vertex resolution in z is 0");
 648     return 0;
 649   }
 650   ((TH1F *)(fListVertexQA->At(0)))->Fill(vertex->GetXv());
 651   ((TH1F *)(fListVertexQA->At(2)))->Fill(vertex->GetYv());
 652   ((TH1F *)(fListVertexQA->At(4)))->Fill(vertex->GetZv());
 653
 654   //check position
 655   if(TMath::Abs(vertex->GetXv()) > gVxMax) {
 656     if(fDebugMode)
 657       Printf("GetVertex: Event rejected because it has a Vx value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetXv()),gVxMax,gVxMax);
 658     return 0;
 659   }
 660   if(TMath::Abs(vertex->GetYv()) > gVyMax)  {
 661     if(fDebugMode)
 662       Printf("GetVertex: Event rejected because it has a Vy value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetYv()),gVyMax,gVyMax);
 663     return 0;
 664   }
 665   if(TMath::Abs(vertex->GetZv()) > gVzMax)  {
 666     if(fDebugMode)
 667       Printf("GetVertex: Event rejected because it has a Vz value of %lf cm (accepted interval: -%lf - %lf)",TMath::Abs(vertex->GetZv()),gVzMax,gVzMax);
 668     return 0;
 669   }
 670   ((TH1F *)(fListVertexQA->At(1)))->Fill(vertex->GetXv());
 671   ((TH1F *)(fListVertexQA->At(3)))->Fill(vertex->GetYv());
 672   ((TH1F *)(fListVertexQA->At(5)))->Fill(vertex->GetZv());
 673   ((TH1F *)(fListVertexQA->At(6)))->Fill(vertex->GetNContributors());
 674
 675   //check number of contributors
 676   if(fMinNumOfContributors > 0) {
 677     if(fMinNumOfContributors > vertex->GetNContributors()) {
 678       if(fDebugMode)
 679         Printf("GetVertex: Event rejected because it has %d number of contributors (requested minimum: %d)",vertex->GetNContributors(),fMinNumOfContributors);
 680
 681       return 0;
 682     }
 683   }
 684
 685   return vertex;
 686 }
 687
 688 //________________________________________________________________________
 689 Bool_t AliProtonAnalysisBase::IsEventTriggered(const AliESDEvent *esd,
 690                                                TriggerMode trigger) {
 691   // check if the event was triggered
 692   ULong64_t triggerMask = esd->GetTriggerMask();
 693   TString firedTriggerClass = esd->GetFiredTriggerClasses();
 694
 695   if(fAnalysisMC) {
 696     // definitions from p-p.cfg
 697     ULong64_t spdFO = (1 << 14);
 698     ULong64_t v0left = (1 << 11);
 699     ULong64_t v0right = (1 << 12);
 700
 701     switch (trigger) {
 702     case kMB1: {
 703       if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
 704         return kTRUE;
 705       break;
 706     }
 707     case kMB2: {
 708       if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
 709         return kTRUE;
 710       break;
 711     }
 712     case kSPDFASTOR: {
 713       if (triggerMask & spdFO)
 714         return kTRUE;
 715       break;
 716     }
 717     }//switch
 718   }
 719   else {
 720     if(kUseOnlineTrigger) {
 721       if(firedTriggerClass.Contains("CINT1B-ABCE-NOPF-ALL"))
 722         return kTRUE;
 723     }
 724     else if(!kUseOnlineTrigger)
 725       return kTRUE;
 726   }
 727
 728   return kFALSE;
 729 }
 730
 731 //________________________________________________________________________
 732 TCanvas *AliProtonAnalysisBase::GetListOfCuts() {
 733   // return the list of cuts and their cut values for reference
 734   TLatex l;
 735   l.SetTextAlign(12);
 736   l.SetTextSize(0.04);
 737
 738   TString listOfCuts;
 739
 740   TCanvas *c = new TCanvas("cListOfCuts","List of cuts",0,0,900,600);
 741   c->SetFillColor(10); c->SetHighLightColor(41);
 742   c->Divide(3,2);
 743
 744   c->cd(1)->SetFillColor(10);
 745   l.DrawLatex(0.3,0.9,"Analysis details: List of cuts\n\n");
 746
 747   listOfCuts = "Analysis level: "; listOfCuts += fProtonAnalysisLevel;
 748   l.DrawLatex(0.1,0.82,listOfCuts.Data());
 749   listOfCuts = "Analysis mode: ";
 750   if(fProtonAnalysisMode == kTPC) listOfCuts += "TPC standalone";
 751   if(fProtonAnalysisMode == kHybrid) listOfCuts += "Hybrid TPC";
 752   if(fProtonAnalysisMode == kFullHybrid) listOfCuts += "Full Hybrid TPC";
 753   if(fProtonAnalysisMode == kGlobal) listOfCuts += "Global tracking";
 754   l.DrawLatex(0.1,0.74,listOfCuts.Data());
 755   listOfCuts = "Trigger mode: ";
 756   if(fTriggerMode == kMB1) listOfCuts += "Minimum bias 1";
 757   if(fTriggerMode == kMB2) listOfCuts += "Minimum bias 2";
 758   if(fTriggerMode == kSPDFASTOR) listOfCuts += "FastOR (SPD)";
 759   l.DrawLatex(0.1,0.66,listOfCuts.Data());
 760   listOfCuts = "PID mode: ";
 761   if(fProtonPIDMode == kBayesian) listOfCuts += "Bayesian PID";
 762   if(fProtonPIDMode == kRatio) {
 763     listOfCuts += "Z = ln((dE/dx)_{exp.}/(dE/dx)_{theor.}) > ";
 764     listOfCuts += fNRatio;
 765   }
 766   if(fProtonPIDMode == kSigma) {
 767     listOfCuts += "N_{#sigma} area: "; listOfCuts += fNSigma;
 768     listOfCuts += " #sigma";
 769   }
 770   //if(fProtonPIDMode == kSigma2) {
 771   //listOfCuts += "N_{#sigma}(2) area: "; listOfCuts += fNSigma;
 772   //listOfCuts += " #sigma";
 773   //}
 774   l.DrawLatex(0.1,0.58,listOfCuts.Data());
 775   listOfCuts = "Accepted vertex diamond: ";
 776   l.DrawLatex(0.1,0.52,listOfCuts.Data());
 777   listOfCuts = "|V_{x}| < "; listOfCuts += fVxMax; listOfCuts += " [cm]";
 778   l.DrawLatex(0.6,0.52,listOfCuts.Data());
 779   listOfCuts = "|V_{y}| < "; listOfCuts += fVyMax; listOfCuts += " [cm]";
 780   l.DrawLatex(0.6,0.45,listOfCuts.Data());
 781   listOfCuts = "|V_{z}| < "; listOfCuts += fVzMax; listOfCuts += " [cm]";
 782   l.DrawLatex(0.6,0.38,listOfCuts.Data());
 783   listOfCuts = "N_{contributors} > "; listOfCuts += fMinNumOfContributors;
 784   l.DrawLatex(0.6,0.31,listOfCuts.Data());
 785   listOfCuts = "Phase space: ";
 786   l.DrawLatex(0.1,0.2,listOfCuts.Data());
 787   if(fAnalysisEtaMode) listOfCuts = "|#eta| < ";
 788   else listOfCuts = "|y| < ";
 789   listOfCuts += TMath::Abs(fMinX);
 790   l.DrawLatex(0.35,0.2,listOfCuts.Data());
 791   listOfCuts = "N_{bins} = ";
 792   listOfCuts += fNBinsX; listOfCuts += " (binning: ";
 793   listOfCuts += (fMaxX-fMinX)/fNBinsX; listOfCuts += ")";
 794   l.DrawLatex(0.35,0.15,listOfCuts.Data());
 795   listOfCuts = ""; listOfCuts += fMinY; listOfCuts += " < P_{T} < ";
 796   listOfCuts += fMaxY; listOfCuts += "GeV/c";
 797   l.DrawLatex(0.35,0.1,listOfCuts.Data());
 798   listOfCuts = "N_{bins} = ";
 799   listOfCuts += fNBinsY; listOfCuts += " (binning: ";
 800   listOfCuts += (fMaxY-fMinY)/fNBinsY; listOfCuts += ")";
 801   l.DrawLatex(0.35,0.05,listOfCuts.Data());
 802
 803   c->cd(2)->SetFillColor(2);
 804   l.DrawLatex(0.3,0.95,"ITS related cuts");
 805   listOfCuts = "Request a cluster on either of the SPD layers: ";
 806   listOfCuts += fPointOnSPDLayersFlag;
 807   l.DrawLatex(0.1,0.9,listOfCuts.Data());
 808   listOfCuts = "Request a cluster on either of the SDD layers: ";
 809   listOfCuts += fPointOnSDDLayersFlag;
 810   l.DrawLatex(0.1,0.83,listOfCuts.Data());
 811   listOfCuts = "Request a cluster on either of the SSD layers: ";
 812   listOfCuts += fPointOnSSDLayersFlag;
 813   l.DrawLatex(0.1,0.76,listOfCuts.Data());
 814
 815   listOfCuts = "Request a cluster on SPD1: ";
 816   listOfCuts += fPointOnITSLayer1Flag;
 817   l.DrawLatex(0.1,0.69,listOfCuts.Data());
 818   listOfCuts = "Request a cluster on SPD2: ";
 819   listOfCuts += fPointOnITSLayer2Flag;
 820   l.DrawLatex(0.1,0.62,listOfCuts.Data());
 821   listOfCuts = "Request a cluster on SDD1: ";
 822   listOfCuts += fPointOnITSLayer3Flag;
 823   l.DrawLatex(0.1,0.55,listOfCuts.Data());
 824   listOfCuts = "Request a cluster on SDD2: ";
 825   listOfCuts += fPointOnITSLayer4Flag;
 826   l.DrawLatex(0.1,0.48,listOfCuts.Data());
 827   listOfCuts = "Request a cluster on SSD1: ";
 828   listOfCuts += fPointOnITSLayer5Flag;
 829   l.DrawLatex(0.1,0.41,listOfCuts.Data());
 830   listOfCuts = "Request a cluster on SSD2: ";
 831   listOfCuts += fPointOnITSLayer6Flag;
 832   l.DrawLatex(0.1,0.34,listOfCuts.Data());
 833   listOfCuts = "Minimum number of ITS clusters: ";
 834   if(fMinITSClustersFlag) listOfCuts += fMinITSClusters;
 835   else listOfCuts += "Not used";
 836   l.DrawLatex(0.1,0.27,listOfCuts.Data());
 837   listOfCuts = "Maximum #chi^{2} per ITS cluster: ";
 838   if(fMaxChi2PerITSClusterFlag) listOfCuts += fMaxChi2PerITSCluster;
 839   else listOfCuts += "Not used";
 840   l.DrawLatex(0.1,0.2,listOfCuts.Data());
 841
 842   c->cd(3)->SetFillColor(3);
 843   l.DrawLatex(0.3,0.9,"TPC related cuts");
 844   listOfCuts = "Minimum number of TPC clusters: ";
 845   if(fMinTPCClustersFlag) listOfCuts += fMinTPCClusters;
 846   else listOfCuts += "Not used";
 847   l.DrawLatex(0.1,0.7,listOfCuts.Data());
 848   listOfCuts = "Maximum #chi^{2} per TPC cluster: ";
 849   if(fMaxChi2PerTPCClusterFlag) listOfCuts += fMaxChi2PerTPCCluster;
 850   else listOfCuts += "Not used";
 851   l.DrawLatex(0.1,0.5,listOfCuts.Data());
 852   listOfCuts = "Minimum number of TPC points for the dE/dx: ";
 853   if(fMinTPCdEdxPointsFlag) listOfCuts += fMinTPCdEdxPoints;
 854   else listOfCuts += "Not used";
 855   l.DrawLatex(0.1,0.3,listOfCuts.Data());
 856
 857   c->cd(4)->SetFillColor(4);
 858   l.DrawLatex(0.3,0.9,"Tracking related cuts");
 859   listOfCuts = "Maximum cov11: ";
 860   if(fMaxCov11Flag) listOfCuts += fMaxCov11;
 861   else listOfCuts += "Not used";
 862   l.DrawLatex(0.1,0.75,listOfCuts.Data());
 863   listOfCuts = "Maximum cov22: ";
 864   if(fMaxCov22Flag) listOfCuts += fMaxCov22;
 865   else listOfCuts += "Not used";
 866   l.DrawLatex(0.1,0.6,listOfCuts.Data());
 867   listOfCuts = "Maximum cov33: ";
 868   if(fMaxCov33Flag) listOfCuts += fMaxCov33;
 869   else listOfCuts += "Not used";
 870   l.DrawLatex(0.1,0.45,listOfCuts.Data());
 871   listOfCuts = "Maximum cov44: ";
 872   if(fMaxCov44Flag) listOfCuts += fMaxCov44;
 873   else listOfCuts += "Not used";
 874   l.DrawLatex(0.1,0.3,listOfCuts.Data());
 875   listOfCuts = "Maximum cov55: ";
 876   if(fMaxCov55Flag) listOfCuts += fMaxCov55;
 877   else listOfCuts += "Not used";
 878   l.DrawLatex(0.1,0.15,listOfCuts.Data());
 879
 880   c->cd(5)->SetFillColor(5);
 881   l.DrawLatex(0.3,0.9,"DCA related cuts");
 882   listOfCuts = "Maximum sigma to vertex: ";
 883   if(fMaxSigmaToVertexFlag) listOfCuts += fMaxSigmaToVertex;
 884   else listOfCuts += "Not used";
 885   l.DrawLatex(0.1,0.8,listOfCuts.Data());
 886   listOfCuts = "Maximum sigma to vertex (TPC): ";
 887   if(fMaxSigmaToVertexTPCFlag) listOfCuts += fMaxSigmaToVertexTPC;
 888   else listOfCuts += "Not used";
 889   l.DrawLatex(0.1,0.72,listOfCuts.Data());
 890   listOfCuts = "Maximum DCA in xy: ";
 891   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXY;
 892   else listOfCuts += "Not used";
 893   l.DrawLatex(0.1,0.64,listOfCuts.Data());
 894   listOfCuts = "Maximum DCA in z: ";
 895   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZ;
 896   else listOfCuts += "Not used";
 897   l.DrawLatex(0.1,0.56,listOfCuts.Data());
 898   listOfCuts = "Maximum DCA in 3D: ";
 899   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3D;
 900   else listOfCuts += "Not used";
 901   l.DrawLatex(0.1,0.48,listOfCuts.Data());
 902   listOfCuts = "Maximum DCA in xy (TPC): ";
 903   if(fMaxDCAXYFlag) listOfCuts += fMaxDCAXYTPC;
 904   else listOfCuts += "Not used";
 905   l.DrawLatex(0.1,0.4,listOfCuts.Data());
 906   listOfCuts = "Maximum DCA in z (TPC): ";
 907   if(fMaxDCAZFlag) listOfCuts += fMaxDCAZTPC;
 908   else listOfCuts += "Not used";
 909   l.DrawLatex(0.1,0.32,listOfCuts.Data());
 910   listOfCuts = "Maximum DCA in 3D (TPC): ";
 911   if(fMaxDCA3DFlag) listOfCuts += fMaxDCA3DTPC;
 912   else listOfCuts += "Not used";
 913   l.DrawLatex(0.1,0.24,listOfCuts.Data());
 914   listOfCuts = "Maximum constrained #chi^{2} (vertex): ";
 915   if(fMaxConstrainChi2Flag) listOfCuts += fMaxConstrainChi2;
 916   else listOfCuts += "Not used";
 917   l.DrawLatex(0.1,0.16,listOfCuts.Data());
 918
 919   c->cd(6)->SetFillColor(6);
 920   l.DrawLatex(0.3,0.9,"Tracking bits related cuts");
 921   listOfCuts = "Request the ITS refit bit: "; listOfCuts += fITSRefitFlag;
 922   l.DrawLatex(0.1,0.7,listOfCuts.Data());
 923   listOfCuts = "Request the TPC refit bit: "; listOfCuts += fTPCRefitFlag;
 924   l.DrawLatex(0.1,0.5,listOfCuts.Data());
 925   listOfCuts = "Request the TPC pid bit: "; listOfCuts += fTPCpidFlag;
 926   l.DrawLatex(0.1,0.3,listOfCuts.Data());
 927   listOfCuts = "Request the ESD pid bit: "; listOfCuts += fESDpidFlag;
 928   l.DrawLatex(0.1,0.1,listOfCuts.Data());
 929
 930   return c;
 931 }
 932
 933 //________________________________________________________________________
 934 Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
 935   //Function that checks if a track is a proton
 936   Double_t probability[5] = {0.,0.,0.,0.,0.};
 937   Double_t gPt = 0.0, gP = 0.0, gEta = 0.0;
 938   Long64_t fParticleType = 0;
 939
 940   //Bayesian approach for the PID
 941   if(fProtonPIDMode == kBayesian) {
 942     if((fProtonAnalysisMode == kTPC)||(fProtonAnalysisMode == kHybrid)||(fProtonAnalysisMode == kFullHybrid)) {
 943       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 944       if(tpcTrack) {
 945         gPt = tpcTrack->Pt();
 946         gP = tpcTrack->P();
 947         track->GetTPCpid(probability);
 948       }
 949       else return 0;
 950     }//TPC standalone or Hybrid TPC
 951     else if(fProtonAnalysisMode == kGlobal) {
 952       gPt = track->Pt();
 953       gP = track->P();
 954       track->GetESDpid(probability);
 955     }//Global tracking
 956
 957     Double_t rcc = 0.0;
 958     for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 959       rcc += probability[i]*GetParticleFraction(i,gP);
 960     if(rcc != 0.0) {
 961       Double_t w[5];
 962       for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 963         w[i] = probability[i]*GetParticleFraction(i,gP)/rcc;
 964       fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
 965     }
 966     if(fParticleType == 4)
 967       return kTRUE;
 968   }//Bayesian pid
 969   //Ratio of the measured over the theoretical dE/dx a la STAR
 970   else if(fProtonPIDMode == kRatio) {
 971     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 972     if(tpcTrack) {
 973       gPt = tpcTrack->Pt();
 974       gP = track->GetInnerParam()->P();
 975       gEta = tpcTrack->Eta();
 976     }
 977     Double_t fAlephParameters[5] = {0.,0.,0.,0.,0.};
 978     if(fAnalysisMC) {
 979       fAlephParameters[0] = 2.15898e+00/50.;
 980       fAlephParameters[1] = 1.75295e+01;
 981       fAlephParameters[2] = 3.40030e-09;
 982       fAlephParameters[3] = 1.96178e+00;
 983       fAlephParameters[4] = 3.91720e+00;
 984     }
 985     else {
 986       fAlephParameters[0] = 0.0283086;
 987       fAlephParameters[1] = 2.63394e+01;
 988       fAlephParameters[2] = 5.04114e-11;
 989       fAlephParameters[3] = 2.12543e+00;
 990       fAlephParameters[4] = 4.88663e+00;
 991     }
 992
 993     AliTPCPIDResponse *tpcResponse = new AliTPCPIDResponse();
 994     tpcResponse->SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
 995
 996     Double_t normalizeddEdx = -10.;
 997     if((track->GetTPCsignal() > 0.0) && (tpcResponse->GetExpectedSignal(gP,AliPID::kProton) > 0.0))
 998       normalizeddEdx = TMath::Log(track->GetTPCsignal()/tpcResponse->GetExpectedSignal(gP,AliPID::kProton));
 999
1000     delete tpcResponse;
1001
1002     if(normalizeddEdx >= fNRatio)
1003       return kTRUE;
1004   }//kRatio PID mode
1005   //Definition of an N-sigma area around the dE/dx vs P band
1006   else if(fProtonPIDMode == kSigma) {
1007     Double_t fAlephParameters[5];
1008     if(fAnalysisMC) {
1009       fAlephParameters[0] = 2.15898e+00/50.;
1010       fAlephParameters[1] = 1.75295e+01;
1011       fAlephParameters[2] = 3.40030e-09;
1012       fAlephParameters[3] = 1.96178e+00;
1013       fAlephParameters[4] = 3.91720e+00;
1014     }
1015     else {
1016       fAlephParameters[0] = 0.0283086;
1017       fAlephParameters[1] = 2.63394e+01;
1018       fAlephParameters[2] = 5.04114e-11;
1019       fAlephParameters[3] = 2.12543e+00;
1020       fAlephParameters[4] = 4.88663e+00;
1021     }
1022
1023     Double_t nsigma = 100.0;
1024     AliTPCPIDResponse *tpcResponse = new AliTPCPIDResponse();
1025     tpcResponse->SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
1026
1027     Double_t mom = track->GetP();
1028     const AliExternalTrackParam *in = track->GetInnerParam();
1029     if (in)
1030       mom = in->GetP();
1031
1032     nsigma = TMath::Abs(tpcResponse->GetNumberOfSigmas(mom,track->GetTPCsignal(),track->GetTPCsignalN(),AliPID::kProton));
1033
1034     delete tpcResponse;
1035     if(nsigma <= fNSigma)
1036       return kTRUE;
1037   }//kSigma PID method
1038
1039   return kFALSE;
1040 }
1041
1042
1043