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