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