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