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