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