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