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