PWG2/SPECTRA/AliProtonAnalysisBase.cxx

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