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