PWG2/SPECTRA/AliProtonAnalysis.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$ */
  15
  16 //-----------------------------------------------------------------
  17 //                 AliProtonAnalysis class
  18 //   This is the class to deal with the proton analysis
  19 //   Origin: Panos Christakoglou, UOA-CERN, Panos.Christakoglou@cern.ch
  20 //-----------------------------------------------------------------
  21 #include <Riostream.h>
  22 #include <TFile.h>
  23 #include <TSystem.h>
  24 #include <TF1.h>
  25 #include <TH2F.h>
  26 #include <TH1D.h>
  27 #include <TH1I.h>
  28 #include <TParticle.h>
  29
  30 #include "AliProtonAnalysis.h"
  31
  32 #include <AliExternalTrackParam.h>
  33 #include <AliAODEvent.h>
  34 #include <AliESDEvent.h>
  35 #include <AliLog.h>
  36 #include <AliPID.h>
  37
  38 #include <AliStack.h>
  39
  40 ClassImp(AliProtonAnalysis)
  41
  42 //____________________________________________________________________//
  43 AliProtonAnalysis::AliProtonAnalysis() :
  44   TObject(),
  45   fNBinsY(0), fMinY(0), fMaxY(0),
  46   fNBinsPt(0), fMinPt(0), fMaxPt(0),
  47   fMinTPCClusters(0), fMinITSClusters(0),
  48   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
  49   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
  50   fMaxSigmaToVertex(0),
  51   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  52   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  53   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
  54   fMaxSigmaToVertexFlag(kFALSE),
  55   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
  56   fFunctionProbabilityFlag(kFALSE),
  57   fElectronFunction(0), fMuonFunction(0),
  58   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
  59   fUseTPCOnly(kFALSE), fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0) {
  60   //Default constructor
  61   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
  62 }
  63
  64 //____________________________________________________________________//
  65 AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) :
  66   TObject(),
  67   fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
  68   fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
  69   fMinTPCClusters(0), fMinITSClusters(0),
  70   fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
  71   fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
  72   fMaxSigmaToVertex(0),
  73   fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
  74   fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
  75   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
  76   fMaxSigmaToVertexFlag(kFALSE),
  77   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
  78   fFunctionProbabilityFlag(kFALSE),
  79   fElectronFunction(0), fMuonFunction(0),
  80   fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
  81   fUseTPCOnly(kFALSE), fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0) {
  82   //Default constructor
  83
  84   fHistEvents = new TH1I("fHistEvents","Analyzed events",1,0,1);
  85
  86   fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  87   fHistYPtProtons->SetStats(kTRUE);
  88   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  89   fHistYPtProtons->GetXaxis()->SetTitle("y");
  90   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
  91
  92   fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
  93   fHistYPtAntiProtons->SetStats(kTRUE);
  94   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
  95   fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
  96   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
  97 }
  98
  99 //____________________________________________________________________//
 100 AliProtonAnalysis::~AliProtonAnalysis() {
 101   //Default destructor
 102
 103 }
 104
 105 //____________________________________________________________________//
 106 void AliProtonAnalysis::InitHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) {
 107   fNBinsY = nbinsY;
 108   fMinY = fLowY;
 109   fMaxY = fHighY;
 110   fNBinsPt = nbinsPt;
 111   fMinPt = fLowPt;
 112   fMaxPt = fHighPt;
 113
 114   fHistEvents = new TH1I("fHistEvents","Anallyzed events",1,0,1);
 115
 116   fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 117   fHistYPtProtons->SetStats(kTRUE);
 118   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
 119   fHistYPtProtons->GetXaxis()->SetTitle("y");
 120   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
 121
 122   fHistYPtAntiProtons = new TH2F("fHistYPtAntiProtons","y-Pt Antiprotons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
 123   fHistYPtAntiProtons->SetStats(kTRUE);
 124   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
 125   fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
 126   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
 127 }
 128
 129 //____________________________________________________________________//
 130 Bool_t AliProtonAnalysis::ReadFromFile(const char* filename) {
 131   Bool_t status = kTRUE;
 132
 133   TFile *file = TFile::Open(filename);
 134   if(!file) {
 135     cout<<"Could not find the input file "<<filename<<endl;
 136     status = kFALSE;
 137   }
 138
 139   TList *list = (TList *)file->Get("outputList1");
 140   if(list) {
 141     cout<<"Retrieving objects from the list "<<list->GetName()<<"..."<<endl;
 142     fHistYPtProtons = (TH2F *)list->At(0);
 143     fHistYPtAntiProtons = (TH2F *)list->At(1);
 144     fHistEvents = (TH1I *)list->At(2);
 145   }
 146   else if(!list) {
 147     cout<<"Retrieving objects from the file... "<<endl;
 148     fHistYPtProtons = (TH2F *)file->Get("fHistYPtProtons");
 149     fHistYPtAntiProtons = (TH2F *)file->Get("fHistYPtAntiProtons");
 150     fHistEvents = (TH1I *)file->Get("fHistEvents");
 151   }
 152   if((!fHistYPtProtons)||(!fHistYPtAntiProtons)||(!fHistEvents)) {
 153     cout<<"Input containers were not found!!!"<<endl;
 154     status = kFALSE;
 155   }
 156   else {
 157     fHistYPtProtons->Sumw2();
 158     fHistYPtAntiProtons->Sumw2();
 159   }
 160
 161   return status;
 162 }
 163
 164 //____________________________________________________________________//
 165 TH1D *AliProtonAnalysis::GetProtonYHistogram() {
 166   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 167   TH1D *fYProtons = (TH1D *)fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"e");
 168   fYProtons->SetStats(kFALSE);
 169   fYProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
 170   fYProtons->SetTitle("dN/dy protons");
 171   fYProtons->SetMarkerStyle(kFullCircle);
 172   fYProtons->SetMarkerColor(4);
 173   if(nAnalyzedEvents > 0)
 174     fYProtons->Scale(1./nAnalyzedEvents);
 175
 176   return fYProtons;
 177 }
 178
 179 //____________________________________________________________________//
 180 TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
 181   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 182   TH1D *fYAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"e");
 183   fYAntiProtons->SetStats(kFALSE);
 184   fYAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
 185   fYAntiProtons->SetTitle("dN/dy antiprotons");
 186   fYAntiProtons->SetMarkerStyle(kFullCircle);
 187   fYAntiProtons->SetMarkerColor(4);
 188   if(nAnalyzedEvents > 0)
 189     fYAntiProtons->Scale(1./nAnalyzedEvents);
 190
 191   return fYAntiProtons;
 192 }
 193
 194 //____________________________________________________________________//
 195 TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
 196   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 197   TH1D *fPtProtons = (TH1D *)fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
 198   fPtProtons->SetStats(kFALSE);
 199   fPtProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
 200   fPtProtons->SetTitle("dN/dPt protons");
 201   fPtProtons->SetMarkerStyle(kFullCircle);
 202   fPtProtons->SetMarkerColor(4);
 203   if(nAnalyzedEvents > 0)
 204     fPtProtons->Scale(1./nAnalyzedEvents);
 205
 206   return fPtProtons;
 207 }
 208
 209 //____________________________________________________________________//
 210 TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
 211   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 212   TH1D *fPtAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e");
 213   fPtAntiProtons->SetStats(kFALSE);
 214   fPtAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
 215   fPtAntiProtons->SetTitle("dN/dPt antiprotons");
 216   fPtAntiProtons->SetMarkerStyle(kFullCircle);
 217   fPtAntiProtons->SetMarkerColor(4);
 218   if(nAnalyzedEvents > 0)
 219     fPtAntiProtons->Scale(1./nAnalyzedEvents);
 220
 221   return fPtAntiProtons;
 222 }
 223
 224 //____________________________________________________________________//
 225 TH1D *AliProtonAnalysis::GetYRatioHistogram() {
 226   TH1D *fYProtons = GetProtonYHistogram();
 227   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 228
 229   TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 230   hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
 231   hRatioY->SetMarkerStyle(kFullCircle);
 232   hRatioY->SetMarkerColor(4);
 233   hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
 234   hRatioY->GetYaxis()->SetTitleOffset(1.4);
 235   hRatioY->GetXaxis()->SetTitle("y");
 236   hRatioY->GetXaxis()->SetTitleColor(1);
 237   hRatioY->SetStats(kFALSE);
 238
 239   return hRatioY;
 240 }
 241
 242 //____________________________________________________________________//
 243 TH1D *AliProtonAnalysis::GetPtRatioHistogram() {
 244   TH1D *fPtProtons = GetProtonPtHistogram();
 245   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 246
 247   TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 248   hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
 249   hRatioPt->SetMarkerStyle(kFullCircle);
 250   hRatioPt->SetMarkerColor(4);
 251   hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
 252   hRatioPt->GetYaxis()->SetTitleOffset(1.4);
 253   hRatioPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 254   hRatioPt->GetXaxis()->SetTitleColor(1);
 255   hRatioPt->SetStats(kFALSE);
 256
 257   return hRatioPt;
 258 }
 259
 260 //____________________________________________________________________//
 261 TH1D *AliProtonAnalysis::GetYAsymmetryHistogram() {
 262   TH1D *fYProtons = GetProtonYHistogram();
 263   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 264
 265   TH1D *hsum = new TH1D("hsumY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 266   hsum->Add(fYProtons,fYAntiProtons,1.0,1.0);
 267
 268   TH1D *hdiff = new TH1D("hdiffY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 269   hdiff->Add(fYProtons,fYAntiProtons,1.0,-1.0);
 270
 271   TH1D *hAsymmetryY = new TH1D("hAsymmetryY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 272   hAsymmetryY->Divide(hdiff,hsum,2.0,1.);
 273   hAsymmetryY->SetMarkerStyle(kFullCircle);
 274   hAsymmetryY->SetMarkerColor(4);
 275   hAsymmetryY->GetYaxis()->SetTitle("A_{p}");
 276   hAsymmetryY->GetYaxis()->SetTitleOffset(1.4);
 277   hAsymmetryY->GetXaxis()->SetTitle("y");
 278   hAsymmetryY->GetXaxis()->SetTitleColor(1);
 279   hAsymmetryY->SetStats(kFALSE);
 280
 281   return hAsymmetryY;
 282 }
 283
 284 //____________________________________________________________________//
 285 TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
 286   TH1D *fPtProtons = GetProtonPtHistogram();
 287   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 288
 289   TH1D *hsum = new TH1D("hsumPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 290   hsum->Add(fPtProtons,fPtAntiProtons,1.0,1.0);
 291
 292   TH1D *hdiff = new TH1D("hdiffPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 293   hdiff->Add(fPtProtons,fPtAntiProtons,1.0,-1.0);
 294
 295   TH1D *hAsymmetryPt = new TH1D("hAsymmetryPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 296   hAsymmetryPt->Divide(hdiff,hsum,2.0,1.);
 297   hAsymmetryPt->SetMarkerStyle(kFullCircle);
 298   hAsymmetryPt->SetMarkerColor(4);
 299   hAsymmetryPt->GetYaxis()->SetTitle("A_{p}");
 300   hAsymmetryPt->GetYaxis()->SetTitleOffset(1.4);
 301   hAsymmetryPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 302   hAsymmetryPt->GetXaxis()->SetTitleColor(1);
 303   hAsymmetryPt->SetStats(kFALSE);
 304
 305   return hAsymmetryPt;
 306 }
 307
 308 //____________________________________________________________________//
 309 Double_t AliProtonAnalysis::GetParticleFraction(Int_t i, Double_t p) {
 310   Double_t partFrac=0;
 311   if(fFunctionProbabilityFlag) {
 312     if(i == 0) partFrac = fElectronFunction->Eval(p);
 313     if(i == 1) partFrac = fMuonFunction->Eval(p);
 314     if(i == 2) partFrac = fPionFunction->Eval(p);
 315     if(i == 3) partFrac = fKaonFunction->Eval(p);
 316     if(i == 4) partFrac = fProtonFunction->Eval(p);
 317   }
 318   else partFrac = fPartFrac[i];
 319
 320   return partFrac;
 321 }
 322
 323 //____________________________________________________________________//
 324 void AliProtonAnalysis::Analyze(AliESDEvent* fESD) {
 325   //Main analysis part - ESD
 326   fHistEvents->Fill(0); //number of analyzed events
 327   Double_t Pt = 0.0, P = 0.0;
 328   Int_t nGoodTracks = fESD->GetNumberOfTracks();
 329   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
 330     AliESDtrack* track = fESD->GetTrack(iTracks);
 331     Double_t probability[5];
 332
 333     if(IsAccepted(track)) {
 334       if(fUseTPCOnly) {
 335         AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 336         if(!tpcTrack) continue;
 337         Pt = tpcTrack->Pt();
 338         P = tpcTrack->P();
 339
 340         //pid
 341         track->GetTPCpid(probability);
 342         Double_t rcc = 0.0;
 343         for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 344           rcc += probability[i]*GetParticleFraction(i,P);
 345         if(rcc == 0.0) continue;
 346         Double_t w[5];
 347         for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 348           w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
 349         Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
 350         if(fParticleType == 4) {
 351           if(tpcTrack->Charge() > 0)
 352             fHistYPtProtons->Fill(Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz()),Pt);
 353           else if(tpcTrack->Charge() < 0)
 354             fHistYPtAntiProtons->Fill(Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz()),Pt);
 355         }//proton check
 356       }//TPC only tracks
 357       else if(!fUseTPCOnly) {
 358         Pt = track->Pt();
 359         P = track->P();
 360
 361         //pid
 362         track->GetESDpid(probability);
 363         Double_t rcc = 0.0;
 364         for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 365           rcc += probability[i]*GetParticleFraction(i,P);
 366         if(rcc == 0.0) continue;
 367         Double_t w[5];
 368         for(Int_t i = 0; i < AliPID::kSPECIES; i++)
 369           w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
 370         Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
 371         if(fParticleType == 4) {
 372           //cout<<"(Anti)protons found..."<<endl;
 373           if(track->Charge() > 0)
 374             fHistYPtProtons->Fill(Rapidity(track->Px(),track->Py(),track->Pz()),Pt);
 375           else if(track->Charge() < 0)
 376             fHistYPtAntiProtons->Fill(Rapidity(track->Px(),track->Py(),track->Pz()),Pt);
 377         }//proton check
 378       }//combined tracking
 379     }//cuts
 380   }//track loop
 381 }
 382
 383 //____________________________________________________________________//
 384 void AliProtonAnalysis::Analyze(AliAODEvent* fAOD) {
 385   //Main analysis part - AOD
 386   fHistEvents->Fill(0); //number of analyzed events
 387   Int_t nGoodTracks = fAOD->GetNumberOfTracks();
 388   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
 389     AliAODTrack* track = fAOD->GetTrack(iTracks);
 390     Double_t Pt = track->Pt();
 391     Double_t P = track->P();
 392
 393     //pid
 394     Double_t probability[10];
 395     track->GetPID(probability);
 396     Double_t rcc = 0.0;
 397     for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*GetParticleFraction(i,P);
 398     if(rcc == 0.0) continue;
 399     Double_t w[10];
 400     for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
 401     Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIESN,w);
 402     if(fParticleType == 4) {
 403       if(track->Charge() > 0) fHistYPtProtons->Fill(track->Y(fParticleType),Pt);
 404       else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(track->Y(fParticleType),Pt);
 405     }//proton check
 406   }//track loop
 407 }
 408
 409 //____________________________________________________________________//
 410 void AliProtonAnalysis::Analyze(AliStack* stack) {
 411   //Main analysis part - MC
 412   fHistEvents->Fill(0); //number of analyzed events
 413   for(Int_t i = 0; i < stack->GetNprimary(); i++) {
 414     TParticle *particle = stack->Particle(i);
 415     if(particle->Pt() < 0.1) continue;
 416     if(TMath::Abs(particle->Eta()) > 1.0) continue;
 417     Int_t pdgcode = particle->GetPdgCode();
 418     if(pdgcode == 2212) fHistYPtProtons->Fill(Rapidity(particle->Px(),
 419                                                        particle->Py(),
 420                                                        particle->Pz()),
 421                                               particle->Pt());
 422     if(pdgcode == -2212) fHistYPtAntiProtons->Fill(Rapidity(particle->Px(),
 423                                                             particle->Py(),
 424                                                             particle->Pz()),
 425                                                    particle->Pt());
 426   }//particle loop
 427 }
 428
 429 //____________________________________________________________________//
 430 Bool_t AliProtonAnalysis::IsAccepted(AliESDtrack* track) {
 431   // Checks if the track is excluded from the cuts
 432   Double_t Pt = 0.0, Px = 0.0, Py = 0.0, Pz = 0.0;
 433   if(fUseTPCOnly) {
 434     AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 435     if(!tpcTrack) {
 436       Pt = 0.0; Px = 0.0; Py = 0.0; Pz = 0.0;
 437     }
 438     else {
 439       Pt = tpcTrack->Pt();
 440       Px = tpcTrack->Px();
 441       Py = tpcTrack->Py();
 442       Pz = tpcTrack->Pz();
 443     }
 444   }
 445   else{
 446     Pt = track->Pt();
 447     Px = track->Px();
 448     Py = track->Py();
 449     Pz = track->Pz();
 450   }
 451
 452   Int_t  fIdxInt[200];
 453   Int_t nClustersITS = track->GetITSclusters(fIdxInt);
 454   Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
 455
 456   Float_t chi2PerClusterITS = -1;
 457   Float_t chi2PerClusterTPC = -1;
 458   if (nClustersTPC!=0)
 459     chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);
 460
 461   Double_t extCov[15];
 462   track->GetExternalCovariance(extCov);
 463
 464   if(fMinITSClustersFlag)
 465     if(nClustersITS < fMinITSClusters) return kFALSE;
 466   if(fMinTPCClustersFlag)
 467     if(nClustersTPC < fMinTPCClusters) return kFALSE;
 468   if(fMaxChi2PerTPCClusterFlag)
 469     if(chi2PerClusterTPC > fMaxChi2PerTPCCluster) return kFALSE;
 470   if(fMaxChi2PerITSClusterFlag)
 471     if(chi2PerClusterITS > fMaxChi2PerITSCluster) return kFALSE;
 472   if(fMaxCov11Flag)
 473     if(extCov[0] > fMaxCov11) return kFALSE;
 474   if(fMaxCov22Flag)
 475     if(extCov[2] > fMaxCov22) return kFALSE;
 476   if(fMaxCov33Flag)
 477     if(extCov[5] > fMaxCov33) return kFALSE;
 478   if(fMaxCov44Flag)
 479     if(extCov[9] > fMaxCov44) return kFALSE;
 480   if(fMaxCov55Flag)
 481     if(extCov[14] > fMaxCov55) return kFALSE;
 482   if(fMaxSigmaToVertexFlag)
 483     if(GetSigmaToVertex(track) > fMaxSigmaToVertex) return kFALSE;
 484   if(fITSRefitFlag)
 485     if ((track->GetStatus() & AliESDtrack::kITSrefit) == 0) return kFALSE;
 486   if(fTPCRefitFlag)
 487     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
 488
 489   if((Pt < fMinPt) || (Pt > fMaxPt)) return kFALSE;
 490   if((Rapidity(Px,Py,Pz) < fMinY) || (Rapidity(Px,Py,Pz) > fMaxY)) return kFALSE;
 491
 492   return kTRUE;
 493 }
 494
 495 //____________________________________________________________________//
 496 Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) {
 497   // Calculates the number of sigma to the vertex.
 498
 499   Float_t b[2];
 500   Float_t bRes[2];
 501   Float_t bCov[3];
 502   esdTrack->GetImpactParameters(b,bCov);
 503   if (bCov[0]<=0 || bCov[2]<=0) {
 504     //AliDebug(1, "Estimated b resolution lower or equal zero!");
 505     bCov[0]=0; bCov[2]=0;
 506   }
 507   bRes[0] = TMath::Sqrt(bCov[0]);
 508   bRes[1] = TMath::Sqrt(bCov[2]);
 509
 510   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 511
 512   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 513
 514   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 515
 516   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 517
 518   return d;
 519 }
 520
 521 //____________________________________________________________________//
 522 Double_t AliProtonAnalysis::Rapidity(Double_t Px, Double_t Py, Double_t Pz) {
 523   //returns the rapidity of the proton - to be removed
 524   Double_t fMass = 9.38270000000000048e-01;
 525
 526   Double_t P = TMath::Sqrt(TMath::Power(Px,2) +
 527                            TMath::Power(Py,2) +
 528                            TMath::Power(Pz,2));
 529   Double_t energy = TMath::Sqrt(P*P + fMass*fMass);
 530   Double_t y = -999;
 531   if(energy != Pz)
 532     y = 0.5*TMath::Log((energy + Pz)/(energy - Pz));
 533
 534   return y;
 535 }
 536
 537 //____________________________________________________________________//
 538 Bool_t AliProtonAnalysis::PrintMean(TH1 *hist, Double_t edge) {
 539   //calculates the mean value of the ratio/asymmetry within \pm edge
 540   Double_t sum = 0.0;
 541   Int_t nentries = 0;
 542   //calculate the mean
 543   for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
 544     Double_t x = hist->GetBinCenter(i+1);
 545     Double_t y = hist->GetBinContent(i+1);
 546     if(TMath::Abs(x) < edge) {
 547       sum += y;
 548       nentries += 1;
 549     }
 550   }
 551   Double_t mean = 0.0;
 552   if(nentries != 0)
 553     mean = sum/nentries;
 554
 555   //calculate the error
 556   for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
 557     Double_t x = hist->GetBinCenter(i+1);
 558     Double_t y = hist->GetBinContent(i+1);
 559     if(TMath::Abs(x) < edge) {
 560       sum += TMath::Power((mean - y),2);
 561       nentries += 1;
 562     }
 563   }
 564
 565   Double_t error = 0.0;
 566   if(nentries != 0)
 567     error =  TMath::Sqrt(sum)/nentries;
 568
 569   cout<<"========================================="<<endl;
 570   cout<<"Input distribution: "<<hist->GetName()<<endl;
 571   cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
 572   cout<<"Mean value :"<<mean<<endl;
 573   cout<<"Error: "<<error<<endl;
 574   cout<<"========================================="<<endl;
 575
 576   return 0;
 577 }
 578
 579 //____________________________________________________________________//
 580 Bool_t AliProtonAnalysis::PrintYields(TH1 *hist, Double_t edge) {
 581   //calculates the (anti)proton yields within the \pm edge
 582   Double_t sum = 0.0, sumerror = 0.0;
 583   Double_t error = 0.0;
 584   for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
 585     Double_t x = hist->GetBinCenter(i+1);
 586     Double_t y = hist->GetBinContent(i+1);
 587     if(TMath::Abs(x) < edge) {
 588       sum += y;
 589       sumerror += TMath::Power(hist->GetBinError(i+1),2);
 590     }
 591   }
 592
 593   error = TMath::Sqrt(sumerror);
 594
 595   cout<<"========================================="<<endl;
 596   cout<<"Input distribution: "<<hist->GetName()<<endl;
 597   cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
 598   cout<<"Yields :"<<sum<<endl;
 599   cout<<"Error: "<<error<<endl;
 600   cout<<"========================================="<<endl;
 601
 602   return 0;
 603 }
 604
 605
 606