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 | Panos.Christakoglou@cern.ch
  20 //-----------------------------------------------------------------
  21 #include <Riostream.h>
  22 #include <TFile.h>
  23 #include <TSystem.h>
  24 #include <TF1.h>
  25 #include <TH2D.h>
  26 #include <TH1D.h>
  27 #include <TH1I.h>
  28 #include <TParticle.h>
  29
  30 #include <AliExternalTrackParam.h>
  31 #include <AliAODEvent.h>
  32 #include <AliESDEvent.h>
  33 //#include <AliLog.h>
  34 #include <AliPID.h>
  35 #include <AliStack.h>
  36 #include <AliCFContainer.h>
  37 #include <AliCFEffGrid.h>
  38 #include <AliCFDataGrid.h>
  39 //#include <AliESDVertex.h>
  40 class AliLog;
  41 class AliESDVertex;
  42
  43 #include "AliProtonAnalysis.h"
  44 #include "AliProtonAnalysisBase.h"
  45
  46 ClassImp(AliProtonAnalysis)
  47
  48 //____________________________________________________________________//
  49 AliProtonAnalysis::AliProtonAnalysis() :
  50   TObject(), fProtonAnalysisBase(0),
  51   fNBinsY(0), fMinY(0), fMaxY(0),
  52   fNBinsPt(0), fMinPt(0), fMaxPt(0),
  53   fProtonContainer(0), fAntiProtonContainer(0),
  54   fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
  55   fEffGridListProtons(0), fCorrectionListProtons2D(0),
  56   fEfficiencyListProtons1D(0), fCorrectionListProtons1D(0),
  57   fEffGridListAntiProtons(0), fCorrectionListAntiProtons2D(0),
  58   fEfficiencyListAntiProtons1D(0), fCorrectionListAntiProtons1D(0),
  59   fCorrectProtons(0), fCorrectAntiProtons(0) {
  60   //Default constructor
  61 }
  62
  63 //____________________________________________________________________//
  64 AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY,
  65                                      Float_t fLowY, Float_t fHighY,
  66                                      Int_t nbinsPt,
  67                                      Float_t fLowPt, Float_t fHighPt) :
  68   TObject(), fProtonAnalysisBase(0),
  69   fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
  70   fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
  71   fProtonContainer(0), fAntiProtonContainer(0),
  72   fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
  73   fEffGridListProtons(0), fCorrectionListProtons2D(0),
  74   fEfficiencyListProtons1D(0), fCorrectionListProtons1D(0),
  75   fEffGridListAntiProtons(0), fCorrectionListAntiProtons2D(0),
  76   fEfficiencyListAntiProtons1D(0), fCorrectionListAntiProtons1D(0),
  77   fCorrectProtons(0), fCorrectAntiProtons(0){
  78   //Default constructor
  79   fHistEvents = new TH1I("fHistEvents","Analyzed events",1,0,1);
  80
  81   fHistYPtProtons = new TH2D("fHistYPtProtons","Protons",
  82                              fNBinsY,fMinY,fMaxY,
  83                              fNBinsPt,fMinPt,fMaxPt);
  84   fHistYPtProtons->SetStats(kTRUE);
  85   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV/c]");
  86   if(fProtonAnalysisBase->GetEtaMode())
  87     fHistYPtProtons->GetXaxis()->SetTitle("#eta");
  88   else
  89     fHistYPtProtons->GetXaxis()->SetTitle("y");
  90   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
  91
  92   fHistYPtAntiProtons = new TH2D("fHistYPtAntiProtons","Antiprotons",
  93                                  fNBinsY,fMinY,fMaxY,
  94                                  fNBinsPt,fMinPt,fMaxPt);
  95   fHistYPtAntiProtons->SetStats(kTRUE);
  96   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV/c]");
  97   if(fProtonAnalysisBase->GetEtaMode())
  98     fHistYPtAntiProtons->GetXaxis()->SetTitle("#eta");
  99   else
 100     fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
 101   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
 102
 103   //setting up the containers
 104   Int_t iBin[2];
 105   iBin[0] = nbinsY;
 106   iBin[1] = nbinsPt;
 107   Double_t *binLimY = new Double_t[nbinsY+1];
 108   Double_t *binLimPt = new Double_t[nbinsPt+1];
 109   //values for bin lower bounds
 110   for(Int_t i = 0; i <= nbinsY; i++)
 111     binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;
 112   for(Int_t i = 0; i <= nbinsPt; i++)
 113     binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;
 114
 115   fProtonContainer = new AliCFContainer("containerProtons",
 116                                         "container for protons",
 117                                         1,2,iBin);
 118   fProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
 119   fProtonContainer->SetBinLimits(1,binLimPt); //pT
 120   fAntiProtonContainer = new AliCFContainer("containerAntiProtons",
 121                                             "container for antiprotons",
 122                                             1,2,iBin);
 123   fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta
 124   fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
 125 }
 126
 127 //____________________________________________________________________//
 128 AliProtonAnalysis::~AliProtonAnalysis() {
 129   //Default destructor
 130   if(fProtonAnalysisBase) delete fProtonAnalysisBase;
 131
 132   if(fHistEvents) delete fHistEvents;
 133   if(fHistYPtProtons) delete fHistYPtProtons;
 134   if(fHistYPtAntiProtons) delete fHistYPtAntiProtons;
 135   if(fProtonContainer) delete fProtonContainer;
 136   if(fAntiProtonContainer) delete fAntiProtonContainer;
 137
 138   if(fEffGridListProtons) delete fEffGridListProtons;
 139   if(fCorrectionListProtons2D) delete fCorrectionListProtons2D;
 140   if(fEfficiencyListProtons1D) delete fEfficiencyListProtons1D;
 141   if(fCorrectionListProtons1D) delete fCorrectionListProtons1D;
 142   if(fEffGridListAntiProtons) delete fEffGridListAntiProtons;
 143   if(fCorrectionListAntiProtons2D) delete fCorrectionListAntiProtons2D;
 144   if(fEfficiencyListAntiProtons1D) delete fEfficiencyListAntiProtons1D;
 145   if(fCorrectionListAntiProtons1D) delete fCorrectionListAntiProtons1D;
 146   if(fCorrectProtons) delete fCorrectProtons;
 147   if(fCorrectAntiProtons) delete fCorrectAntiProtons;
 148 }
 149
 150 //____________________________________________________________________//
 151 void AliProtonAnalysis::InitAnalysisHistograms(Int_t nbinsY,
 152                                                Float_t fLowY, Float_t fHighY,
 153                                                Int_t nbinsPt,
 154                                                Float_t fLowPt, Float_t fHighPt) {
 155   //Initializes the histograms
 156   fNBinsY = nbinsY;
 157   fMinY = fLowY;
 158   fMaxY = fHighY;
 159   fNBinsPt = nbinsPt;
 160   fMinPt = fLowPt;
 161   fMaxPt = fHighPt;
 162
 163   fHistEvents = new TH1I("fHistEvents","Analyzed events",1,0,1);
 164
 165   fHistYPtProtons = new TH2D("fHistYPtProtons","Protons",
 166                              fNBinsY,fMinY,fMaxY,
 167                              fNBinsPt,fMinPt,fMaxPt);
 168   fHistYPtProtons->SetStats(kTRUE);
 169   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV/c]");
 170   if(fProtonAnalysisBase->GetEtaMode())
 171     fHistYPtProtons->GetXaxis()->SetTitle("#eta");
 172   else
 173     fHistYPtProtons->GetXaxis()->SetTitle("y");
 174   fHistYPtProtons->GetXaxis()->SetTitleColor(1);
 175
 176   fHistYPtAntiProtons = new TH2D("fHistYPtAntiProtons","Antiprotons",
 177                                  fNBinsY,fMinY,fMaxY,
 178                                  fNBinsPt,fMinPt,fMaxPt);
 179   fHistYPtAntiProtons->SetStats(kTRUE);
 180   fHistYPtAntiProtons->GetYaxis()->SetTitle("P_{T} [GeV/c]");
 181   if(fProtonAnalysisBase->GetEtaMode())
 182     fHistYPtAntiProtons->GetXaxis()->SetTitle("#eta");
 183   else
 184     fHistYPtAntiProtons->GetXaxis()->SetTitle("y");
 185   fHistYPtAntiProtons->GetXaxis()->SetTitleColor(1);
 186
 187   //setting up the containers
 188   Int_t iBin[2];
 189   iBin[0] = nbinsY;
 190   iBin[1] = nbinsPt;
 191   Double_t *binLimY = new Double_t[nbinsY+1];
 192   Double_t *binLimPt = new Double_t[nbinsPt+1];
 193   //values for bin lower bounds
 194   for(Int_t i = 0; i <= nbinsY; i++)
 195     binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;
 196   for(Int_t i = 0; i <= nbinsPt; i++)
 197     binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;
 198
 199   fProtonContainer = new AliCFContainer("containerProtons",
 200                                         "container for protons",
 201                                         1,2,iBin);
 202   fProtonContainer->SetBinLimits(0,binLimY); //rapidity
 203   fProtonContainer->SetBinLimits(1,binLimPt); //pT
 204   fAntiProtonContainer = new AliCFContainer("containerAntiProtons",
 205                                             "container for antiprotons",
 206                                             1,2,iBin);
 207   fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity
 208   fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT
 209 }
 210
 211 //____________________________________________________________________//
 212 Bool_t AliProtonAnalysis::ReadFromFile(const char* filename) {
 213   //Read the containers from the existing file
 214   Bool_t status = kTRUE;
 215
 216   TFile *file = TFile::Open(filename);
 217   if(!file) {
 218     cout<<"Could not find the input file "<<filename<<endl;
 219     status = kFALSE;
 220   }
 221
 222   TList *list = (TList *)file->Get("outputList");
 223   if(list) {
 224     cout<<"Retrieving objects from the list "<<list->GetName()<<"..."<<endl;
 225     fHistYPtProtons = (TH2D *)list->At(0);
 226     fHistYPtAntiProtons = (TH2D *)list->At(1);
 227     fHistEvents = (TH1I *)list->At(2);
 228     fProtonContainer = (AliCFContainer *)list->At(3);
 229     fAntiProtonContainer = (AliCFContainer *)list->At(4);
 230   }
 231   else if(!list) {
 232     cout<<"Retrieving objects from the file... "<<endl;
 233     fHistYPtProtons = (TH2D *)file->Get("fHistYPtProtons");
 234     fHistYPtAntiProtons = (TH2D *)file->Get("fHistYPtAntiProtons");
 235     fHistEvents = (TH1I *)file->Get("fHistEvents");
 236     fProtonContainer = (AliCFContainer *)file->Get("containerProtons");
 237     fAntiProtonContainer = (AliCFContainer *)file->Get("containerAntiProtons");
 238   }
 239   if((!fHistYPtProtons)||(!fHistYPtAntiProtons)||(!fHistEvents)
 240      ||(!fProtonContainer)||(!fAntiProtonContainer)) {
 241     cout<<"Input containers were not found!!!"<<endl;
 242     status = kFALSE;
 243   }
 244   else {
 245     //fHistYPtProtons = fProtonContainer->ShowProjection(0,1,0);
 246     //fHistYPtAntiProtons = fAntiProtonContainer->ShowProjection(0,1,0);
 247     fHistYPtProtons->Sumw2();
 248     fHistYPtAntiProtons->Sumw2();
 249   }
 250
 251   return status;
 252 }
 253
 254 //____________________________________________________________________//
 255 TH1D *AliProtonAnalysis::GetProtonYHistogram() {
 256   //Get the y histogram for protons
 257   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 258
 259   TH1D *fYProtons = (TH1D *)fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"");
 260   //TH1D *fYProtons = fProtonContainer->ShowProjection(0,0); //variable-step
 261
 262   fYProtons->SetStats(kFALSE);
 263   fYProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
 264   fYProtons->SetTitle("dN/dy protons");
 265   fYProtons->SetMarkerStyle(kFullCircle);
 266   fYProtons->SetMarkerColor(4);
 267   if(nAnalyzedEvents > 0)
 268   fYProtons->Scale(1./nAnalyzedEvents);
 269
 270   return fYProtons;
 271 }
 272
 273 //____________________________________________________________________//
 274 TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
 275   //Get the y histogram for antiprotons
 276   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 277
 278   TH1D *fYAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"");
 279   //TH1D *fYAntiProtons = fAntiProtonContainer->ShowProjection(0,0);//variable-step
 280
 281   fYAntiProtons->SetStats(kFALSE);
 282   fYAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
 283   fYAntiProtons->SetTitle("dN/dy antiprotons");
 284   fYAntiProtons->SetMarkerStyle(kFullCircle);
 285   fYAntiProtons->SetMarkerColor(4);
 286   if(nAnalyzedEvents > 0)
 287     fYAntiProtons->Scale(1./nAnalyzedEvents);
 288
 289   return fYAntiProtons;
 290 }
 291
 292 //____________________________________________________________________//
 293 TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
 294   //Get the Pt histogram for protons
 295   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 296
 297   TH1D *fPtProtons = (TH1D *)fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"");
 298   //TH1D *fPtProtons = fProtonContainer->ShowProjection(1,0); //variable-step
 299
 300   fPtProtons->SetStats(kFALSE);
 301   fPtProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
 302   fPtProtons->SetTitle("dN/dPt protons");
 303   fPtProtons->SetMarkerStyle(kFullCircle);
 304   fPtProtons->SetMarkerColor(4);
 305   if(nAnalyzedEvents > 0)
 306     fPtProtons->Scale(1./nAnalyzedEvents);
 307
 308   return fPtProtons;
 309 }
 310
 311 //____________________________________________________________________//
 312 TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
 313   //Get the Pt histogram for antiprotons
 314   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 315
 316   TH1D *fPtAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"");
 317   //TH1D *fPtAntiProtons = fAntiProtonContainer->ShowProjection(1,0); //variable-step
 318
 319   fPtAntiProtons->SetStats(kFALSE);
 320   fPtAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
 321   fPtAntiProtons->SetTitle("dN/dPt antiprotons");
 322   fPtAntiProtons->SetMarkerStyle(kFullCircle);
 323   fPtAntiProtons->SetMarkerColor(4);
 324   if(nAnalyzedEvents > 0)
 325     fPtAntiProtons->Scale(1./nAnalyzedEvents);
 326
 327   return fPtAntiProtons;
 328 }
 329
 330 //____________________________________________________________________//
 331 TH1D *AliProtonAnalysis::GetProtonCorrectedYHistogram() {
 332   //Get the corrected y histogram for protons
 333   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 334
 335   TH1D *fYProtons = fCorrectProtons->Project(0); //0: rapidity
 336
 337   fYProtons->SetStats(kFALSE);
 338   fYProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
 339   fYProtons->GetXaxis()->SetTitle("y");
 340   fYProtons->SetTitle("dN/dy protons");
 341   fYProtons->SetMarkerStyle(kFullCircle);
 342   fYProtons->SetMarkerColor(4);
 343   if(nAnalyzedEvents > 0)
 344     fYProtons->Scale(1./nAnalyzedEvents);
 345
 346   return fYProtons;
 347 }
 348
 349 //____________________________________________________________________//
 350 TH1D *AliProtonAnalysis::GetAntiProtonCorrectedYHistogram() {
 351   //Get the corrected y histogram for antiprotons
 352   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 353
 354   TH1D *fYAntiProtons = fCorrectAntiProtons->Project(0); //0: rapidity
 355
 356   fYAntiProtons->SetStats(kFALSE);
 357   fYAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
 358   fYAntiProtons->GetXaxis()->SetTitle("y");
 359   fYAntiProtons->SetTitle("dN/dy protons");
 360   fYAntiProtons->SetMarkerStyle(kFullCircle);
 361   fYAntiProtons->SetMarkerColor(4);
 362   if(nAnalyzedEvents > 0)
 363     fYAntiProtons->Scale(1./nAnalyzedEvents);
 364
 365   return fYAntiProtons;
 366 }
 367
 368 //____________________________________________________________________//
 369 TH1D *AliProtonAnalysis::GetProtonCorrectedPtHistogram() {
 370   //Get the corrected Pt histogram for protons
 371   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 372
 373   TH1D *fPtProtons = fCorrectProtons->Project(0); //0: rapidity
 374
 375   fPtProtons->SetStats(kFALSE);
 376   fPtProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
 377   fPtProtons->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 378   fPtProtons->SetTitle("dN/dPt protons");
 379   fPtProtons->SetMarkerStyle(kFullCircle);
 380   fPtProtons->SetMarkerColor(4);
 381   if(nAnalyzedEvents > 0)
 382     fPtProtons->Scale(1./nAnalyzedEvents);
 383
 384   return fPtProtons;
 385 }
 386
 387 //____________________________________________________________________//
 388 TH1D *AliProtonAnalysis::GetAntiProtonCorrectedPtHistogram() {
 389   //Get the corrected Pt histogram for antiprotons
 390   Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
 391
 392   TH1D *fPtAntiProtons = fCorrectAntiProtons->Project(0); //0: rapidity
 393
 394   fPtAntiProtons->SetStats(kFALSE);
 395   fPtAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
 396   fPtAntiProtons->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 397   fPtAntiProtons->SetTitle("dN/dPt antiprotons");
 398   fPtAntiProtons->SetMarkerStyle(kFullCircle);
 399   fPtAntiProtons->SetMarkerColor(4);
 400   if(nAnalyzedEvents > 0)
 401     fPtAntiProtons->Scale(1./nAnalyzedEvents);
 402
 403   return fPtAntiProtons;
 404 }
 405
 406 //____________________________________________________________________//
 407 TH1D *AliProtonAnalysis::GetYRatioHistogram() {
 408   //Returns the rapidity dependence of the ratio (uncorrected)
 409   TH1D *fYProtons = GetProtonYHistogram();
 410   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 411
 412   TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 413   hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
 414   hRatioY->SetMarkerStyle(kFullCircle);
 415   hRatioY->SetMarkerColor(4);
 416   hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
 417   hRatioY->GetYaxis()->SetTitleOffset(1.4);
 418   hRatioY->GetXaxis()->SetTitle("y");
 419   hRatioY->GetXaxis()->SetTitleColor(1);
 420   hRatioY->SetStats(kFALSE);
 421
 422   return hRatioY;
 423 }
 424
 425 //____________________________________________________________________//
 426 TH1D *AliProtonAnalysis::GetYRatioCorrectedHistogram(TH2D *gCorrectionProtons,
 427                                                      TH2D *gCorrectionAntiProtons) {
 428   //Returns the rapidity dependence of the ratio (corrected)
 429   fHistYPtProtons->Multiply(gCorrectionProtons);
 430   TH1D *fYProtons = GetProtonYHistogram();
 431   fHistYPtAntiProtons->Multiply(gCorrectionAntiProtons);
 432   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 433
 434   TH1D *hRatioY = new TH1D("hRatioY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 435   hRatioY->Divide(fYAntiProtons,fYProtons,1.0,1.0);
 436   hRatioY->SetMarkerStyle(kFullCircle);
 437   hRatioY->SetMarkerColor(4);
 438   hRatioY->GetYaxis()->SetTitle("#bar{p}/p");
 439   hRatioY->GetYaxis()->SetTitleOffset(1.4);
 440   hRatioY->GetXaxis()->SetTitle("y");
 441   hRatioY->GetXaxis()->SetTitleColor(1);
 442   hRatioY->SetStats(kFALSE);
 443
 444   return hRatioY;
 445 }
 446
 447 //____________________________________________________________________//
 448 TH1D *AliProtonAnalysis::GetPtRatioHistogram() {
 449   //Returns the pT dependence of the ratio (uncorrected)
 450   TH1D *fPtProtons = GetProtonPtHistogram();
 451   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 452
 453   TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 454   hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
 455   hRatioPt->SetMarkerStyle(kFullCircle);
 456   hRatioPt->SetMarkerColor(4);
 457   hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
 458   hRatioPt->GetYaxis()->SetTitleOffset(1.4);
 459   hRatioPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 460   hRatioPt->GetXaxis()->SetTitleColor(1);
 461   hRatioPt->SetStats(kFALSE);
 462
 463   return hRatioPt;
 464 }
 465
 466 //____________________________________________________________________//
 467 TH1D *AliProtonAnalysis::GetPtRatioCorrectedHistogram(TH2D *gCorrectionProtons,
 468                                                       TH2D *gCorrectionAntiProtons) {
 469   //Returns the Pt dependence of the ratio (corrected)
 470   fHistYPtProtons->Multiply(gCorrectionProtons);
 471   TH1D *fPtProtons = GetProtonPtHistogram();
 472   fHistYPtAntiProtons->Multiply(gCorrectionAntiProtons);
 473   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 474
 475   TH1D *hRatioPt = new TH1D("hRatioPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 476   hRatioPt->Divide(fPtAntiProtons,fPtProtons,1.0,1.0);
 477   hRatioPt->SetMarkerStyle(kFullCircle);
 478   hRatioPt->SetMarkerColor(4);
 479   hRatioPt->GetYaxis()->SetTitle("#bar{p}/p");
 480   hRatioPt->GetYaxis()->SetTitleOffset(1.4);
 481   hRatioPt->GetXaxis()->SetTitle("y");
 482   hRatioPt->GetXaxis()->SetTitleColor(1);
 483   hRatioPt->SetStats(kFALSE);
 484
 485   return hRatioPt;
 486 }
 487
 488 //____________________________________________________________________//
 489 TH1D *AliProtonAnalysis::GetYAsymmetryHistogram() {
 490   //Returns the rapidity dependence of the asymmetry (uncorrected)
 491   TH1D *fYProtons = GetProtonYHistogram();
 492   TH1D *fYAntiProtons = GetAntiProtonYHistogram();
 493
 494   TH1D *hsum = new TH1D("hsumY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 495   hsum->Add(fYProtons,fYAntiProtons,1.0,1.0);
 496
 497   TH1D *hdiff = new TH1D("hdiffY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 498   hdiff->Add(fYProtons,fYAntiProtons,1.0,-1.0);
 499
 500   TH1D *hAsymmetryY = new TH1D("hAsymmetryY","",fYProtons->GetNbinsX(),fYProtons->GetXaxis()->GetXmin(),fYProtons->GetXaxis()->GetXmax());
 501   hAsymmetryY->Divide(hdiff,hsum,2.0,1.);
 502   hAsymmetryY->SetMarkerStyle(kFullCircle);
 503   hAsymmetryY->SetMarkerColor(4);
 504   hAsymmetryY->GetYaxis()->SetTitle("A_{p}");
 505   hAsymmetryY->GetYaxis()->SetTitleOffset(1.4);
 506   hAsymmetryY->GetXaxis()->SetTitle("y");
 507   hAsymmetryY->GetXaxis()->SetTitleColor(1);
 508   hAsymmetryY->SetStats(kFALSE);
 509
 510   return hAsymmetryY;
 511 }
 512
 513 //____________________________________________________________________//
 514 TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
 515   //Returns the pT dependence of the asymmetry (uncorrected)
 516   TH1D *fPtProtons = GetProtonPtHistogram();
 517   TH1D *fPtAntiProtons = GetAntiProtonPtHistogram();
 518
 519   TH1D *hsum = new TH1D("hsumPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 520   hsum->Add(fPtProtons,fPtAntiProtons,1.0,1.0);
 521
 522   TH1D *hdiff = new TH1D("hdiffPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 523   hdiff->Add(fPtProtons,fPtAntiProtons,1.0,-1.0);
 524
 525   TH1D *hAsymmetryPt = new TH1D("hAsymmetryPt","",fPtProtons->GetNbinsX(),fPtProtons->GetXaxis()->GetXmin(),fPtProtons->GetXaxis()->GetXmax());
 526   hAsymmetryPt->Divide(hdiff,hsum,2.0,1.);
 527   hAsymmetryPt->SetMarkerStyle(kFullCircle);
 528   hAsymmetryPt->SetMarkerColor(4);
 529   hAsymmetryPt->GetYaxis()->SetTitle("A_{p}");
 530   hAsymmetryPt->GetYaxis()->SetTitleOffset(1.4);
 531   hAsymmetryPt->GetXaxis()->SetTitle("P_{T} [GeV/c]");
 532   hAsymmetryPt->GetXaxis()->SetTitleColor(1);
 533   hAsymmetryPt->SetStats(kFALSE);
 534
 535   return hAsymmetryPt;
 536 }
 537
 538 //____________________________________________________________________//
 539 void AliProtonAnalysis::Analyze(AliESDEvent* esd,
 540                                 const AliESDVertex *vertex) {
 541   //Main analysis part - ESD
 542   Int_t nTracks = 0;
 543   Int_t nIdentifiedProtons = 0, nIdentifiedAntiProtons = 0;
 544   Int_t nSurvivedProtons = 0, nSurvivedAntiProtons = 0;
 545
 546   fHistEvents->Fill(0); //number of analyzed events
 547   Double_t containerInput[2] ;
 548   Double_t gPt = 0.0, gP = 0.0;
 549   nTracks = esd->GetNumberOfTracks();
 550   for(Int_t iTracks = 0; iTracks < nTracks; iTracks++) {
 551     AliESDtrack* track = esd->GetTrack(iTracks);
 552     AliESDtrack trackTPC;
 553
 554     //in case it's a TPC only track relate it to the proper vertex
 555     /*if(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC) {
 556       Float_t p[2],cov[3];
 557       track->GetImpactParametersTPC(p,cov);
 558       if (p[0]==0 && p[1]==0)
 559         track->RelateToVertexTPC(((AliESDEvent*)esd)->GetPrimaryVertexTPC(),esd->GetMagneticField(),kVeryBig);
 560       if (!track->FillTPCOnlyTrack(trackTPC)) {
 561         continue;
 562       }
 563       track = &trackTPC ;
 564       }*/
 565
 566     if((fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC)||(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kHybrid)) {
 567       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
 568       if(!tpcTrack) continue;
 569       gPt = tpcTrack->Pt();
 570       gP = tpcTrack->P();
 571
 572       if(fProtonAnalysisBase->IsProton(track)) {
 573           if(tpcTrack->Charge() > 0) {
 574             nIdentifiedProtons += 1;
 575             if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) continue;//track cuts
 576             if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue; //track outside the analyzed y-Pt
 577             nSurvivedProtons += 1;
 578             if(fProtonAnalysisBase->GetEtaMode()) {
 579               fHistYPtProtons->Fill(tpcTrack->Eta(),
 580                                     gPt);
 581               containerInput[0] = tpcTrack->Eta();
 582             }
 583             else {
 584               fHistYPtProtons->Fill(fProtonAnalysisBase->Rapidity(tpcTrack->Px(),
 585                                                                   tpcTrack->Py(),
 586                                                                   tpcTrack->Pz()),
 587                                     gPt);
 588               //fill the container
 589               containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),
 590                                                                 tpcTrack->Py(),
 591                                                                 tpcTrack->Pz());
 592             }
 593             containerInput[1] = gPt;
 594             fProtonContainer->Fill(containerInput,0);
 595           }//protons
 596           else if(tpcTrack->Charge() < 0) {
 597             nIdentifiedAntiProtons += 1;
 598             if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) continue;//track cuts
 599             if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue; //track outside the analyzed y-Pt
 600             nSurvivedAntiProtons += 1;
 601             if(fProtonAnalysisBase->GetEtaMode()) {
 602               fHistYPtAntiProtons->Fill(tpcTrack->Eta(),
 603                                         gPt);
 604               containerInput[0] = tpcTrack->Eta();
 605             }
 606             else {
 607               fHistYPtAntiProtons->Fill(fProtonAnalysisBase->Rapidity(tpcTrack->Px(),
 608                                                                       tpcTrack->Py(),
 609                                                                       tpcTrack->Pz()),
 610                                         gPt);
 611               //fill the container
 612               containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),
 613                                                                 tpcTrack->Py(),
 614                                                                 tpcTrack->Pz());
 615             }
 616             containerInput[1] = gPt;
 617             fAntiProtonContainer->Fill(containerInput,0);
 618           }//antiprotons
 619       }//proton check
 620     }//TPC only tracks
 621     else if(fProtonAnalysisBase->GetAnalysisMode() == AliProtonAnalysisBase::kGlobal) {
 622       gPt = track->Pt();
 623       gP = track->P();
 624
 625       if(fProtonAnalysisBase->IsProton(track)) {
 626         if(track->Charge() > 0) {
 627           nIdentifiedProtons += 1;
 628           if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) continue;//track cuts
 629           if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue; //track outside the analyzed y-Pt
 630           nSurvivedProtons += 1;
 631           if(fProtonAnalysisBase->GetEtaMode()) {
 632             fHistYPtProtons->Fill(track->Eta(),
 633                                   gPt);
 634             containerInput[0] = track->Eta();
 635           }
 636           else {
 637             fHistYPtProtons->Fill(fProtonAnalysisBase->Rapidity(track->Px(),
 638                                                                 track->Py(),
 639                                                                 track->Pz()),
 640                                   gPt);
 641             //fill the container
 642             containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),
 643                                                               track->Py(),
 644                                                               track->Pz());
 645           }
 646           containerInput[1] = gPt;
 647           fProtonContainer->Fill(containerInput,0);
 648         }//protons
 649         else if(track->Charge() < 0) {
 650           nIdentifiedAntiProtons += 1;
 651           if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) continue;//track cuts
 652           if(!fProtonAnalysisBase->IsInPhaseSpace(track)) continue; //track outside the analyzed y-Pt
 653           nSurvivedAntiProtons += 1;
 654           if(fProtonAnalysisBase->GetEtaMode()) {
 655             fHistYPtAntiProtons->Fill(track->Eta(),
 656                                       gPt);
 657             containerInput[0] = track->Eta();
 658           }
 659           else {
 660             fHistYPtAntiProtons->Fill(fProtonAnalysisBase->Rapidity(track->Px(),
 661                                                                     track->Py(),
 662                                                                     track->Pz()),
 663                                       gPt);
 664             //fill the container
 665             containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),
 666                                                               track->Py(),
 667                                                               track->Pz());
 668           }
 669           containerInput[1] = gPt;
 670           fAntiProtonContainer->Fill(containerInput,0);
 671         }//antiprotons
 672       }//proton check
 673     }//combined tracking
 674   }//track loop
 675
 676   if(fProtonAnalysisBase->GetDebugMode())
 677     Printf("Initial number of tracks: %d | Identified (anti)protons: %d - %d | Survived (anti)protons: %d - %d",nTracks,nIdentifiedProtons,nIdentifiedAntiProtons,nSurvivedProtons,nSurvivedAntiProtons);
 678 }
 679
 680 //____________________________________________________________________//
 681 void AliProtonAnalysis::Analyze(AliAODEvent* const fAOD) {
 682   //Main analysis part - AOD
 683   fHistEvents->Fill(0); //number of analyzed events
 684   Int_t nTracks = fAOD->GetNumberOfTracks();
 685   for(Int_t iTracks = 0; iTracks < nTracks; iTracks++) {
 686     AliAODTrack* track = fAOD->GetTrack(iTracks);
 687     Double_t gPt = track->Pt();
 688     Double_t gP = track->P();
 689
 690     //pid
 691     Double_t probability[10];
 692     track->GetPID(probability);
 693     Double_t rcc = 0.0;
 694     for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*fProtonAnalysisBase->GetParticleFraction(i,gP);
 695     if(rcc == 0.0) continue;
 696     Double_t w[10];
 697     for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*fProtonAnalysisBase->GetParticleFraction(i,gP)/rcc;
 698     Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIESN,w);
 699     if(fParticleType == 4) {
 700       if(track->Charge() > 0)
 701         fHistYPtProtons->Fill(track->Y(fParticleType),gPt);
 702       else if(track->Charge() < 0)
 703         fHistYPtAntiProtons->Fill(track->Y(fParticleType),gPt);
 704     }//proton check
 705   }//track loop
 706 }
 707
 708 //____________________________________________________________________//
 709 void AliProtonAnalysis::Analyze(AliStack* const stack,
 710                                 Bool_t iInclusive) {
 711   //Main analysis part - MC
 712   fHistEvents->Fill(0); //number of analyzed events
 713
 714   Int_t nParticles = 0;
 715   //inclusive protons -
 716   if(iInclusive) nParticles = stack->GetNtrack();
 717   else nParticles = stack->GetNprimary();
 718
 719   for(Int_t i = 0; i < nParticles; i++) {
 720     TParticle *particle = stack->Particle(i);
 721     if(!particle) continue;
 722
 723     //in case of inclusive protons reject the secondaries from hadronic inter.
 724     if(particle->GetUniqueID() == 13) continue;
 725
 726     if(TMath::Abs(particle->Eta()) > 1.0) continue;
 727     if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt)) continue;
 728     if((fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz()) > fMaxY)||(fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz()) < fMinY)) continue;
 729
 730     Int_t pdgcode = particle->GetPdgCode();
 731     if(pdgcode == 2212) fHistYPtProtons->Fill(fProtonAnalysisBase->Rapidity(particle->Px(),
 732                                                                             particle->Py(),
 733                                                                             particle->Pz()),
 734                                               particle->Pt());
 735     if(pdgcode == -2212) fHistYPtAntiProtons->Fill(fProtonAnalysisBase->Rapidity(particle->Px(),
 736                                                                                  particle->Py(),
 737                                                                                  particle->Pz()),
 738                                                    particle->Pt());
 739   }//particle loop
 740 }
 741
 742 //____________________________________________________________________//
 743 Bool_t AliProtonAnalysis::PrintMean(TH1 *hist, Double_t edge) {
 744   //calculates the mean value of the ratio/asymmetry within \pm edge
 745   Double_t sum = 0.0;
 746   Int_t nentries = 0;
 747   //calculate the mean
 748   for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
 749     Double_t x = hist->GetBinCenter(i+1);
 750     Double_t y = hist->GetBinContent(i+1);
 751     if(TMath::Abs(x) < edge) {
 752       sum += y;
 753       nentries += 1;
 754     }
 755   }
 756   Double_t mean = 0.0;
 757   if(nentries != 0)
 758     mean = sum/nentries;
 759
 760   //calculate the error
 761   for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
 762     Double_t x = hist->GetBinCenter(i+1);
 763     Double_t y = hist->GetBinContent(i+1);
 764     if(TMath::Abs(x) < edge) {
 765       sum += TMath::Power((mean - y),2);
 766       nentries += 1;
 767     }
 768   }
 769
 770   Double_t error = 0.0;
 771   if(nentries != 0)
 772     error =  TMath::Sqrt(sum)/nentries;
 773
 774   cout<<"========================================="<<endl;
 775   cout<<"Input distribution: "<<hist->GetName()<<endl;
 776   cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
 777   cout<<"Mean value :"<<mean<<endl;
 778   cout<<"Error: "<<error<<endl;
 779   cout<<"========================================="<<endl;
 780
 781   return 0;
 782 }
 783
 784 //____________________________________________________________________//
 785 Bool_t AliProtonAnalysis::PrintYields(TH1 *hist, Double_t edge) {
 786   //calculates the (anti)proton yields within the \pm edge
 787   Double_t sum = 0.0, sumerror = 0.0;
 788   Double_t error = 0.0;
 789   for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
 790     Double_t x = hist->GetBinCenter(i+1);
 791     Double_t y = hist->GetBinContent(i+1);
 792     if(TMath::Abs(x) < edge) {
 793       sum += y;
 794       sumerror += TMath::Power(hist->GetBinError(i+1),2);
 795     }
 796   }
 797
 798   error = TMath::Sqrt(sumerror);
 799
 800   cout<<"========================================="<<endl;
 801   cout<<"Input distribution: "<<hist->GetName()<<endl;
 802   cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
 803   cout<<"Yields :"<<sum<<endl;
 804   cout<<"Error: "<<error<<endl;
 805   cout<<"========================================="<<endl;
 806
 807   return 0;
 808 }
 809
 810 //____________________________________________________________________//
 811 void AliProtonAnalysis::Correct(Int_t step) {
 812   //Applies the correction maps to the initial containers
 813   fCorrectProtons = new AliCFDataGrid("correctProtons",
 814                                       "corrected data",
 815                                       *fProtonContainer);
 816   fCorrectProtons->SetMeasured(0);
 817   fCorrectProtons->ApplyEffCorrection(*(AliCFEffGrid *)fEffGridListProtons->At(step));
 818
 819   fCorrectAntiProtons = new AliCFDataGrid("correctAntiProtons",
 820                                           "corrected data",
 821                                           *fAntiProtonContainer);
 822   fCorrectAntiProtons->SetMeasured(0);
 823   fCorrectAntiProtons->ApplyEffCorrection(*(AliCFEffGrid *)fEffGridListAntiProtons->At(step));
 824 }
 825
 826 //____________________________________________________________________//
 827 Bool_t AliProtonAnalysis::ReadCorrectionContainer(const char* filename) {
 828   // Reads the outout of the correction framework task
 829   // Creates the correction maps
 830   // Puts the results in the different TList objects
 831   Bool_t status = kTRUE;
 832
 833   TFile *file = TFile::Open(filename);
 834   if(!file) {
 835     cout<<"Could not find the input CORRFW file "<<filename<<endl;
 836     status = kFALSE;
 837   }
 838
 839   //________________________________________//
 840   //Protons
 841   fEffGridListProtons = new TList();
 842   fCorrectionListProtons2D = new TList();
 843   fEfficiencyListProtons1D = new TList();
 844   fCorrectionListProtons1D = new TList();
 845
 846   AliCFContainer *corrfwContainerProtons = (AliCFContainer*) (file->Get("containerProtons"));
 847   if(!corrfwContainerProtons) {
 848     cout<<"CORRFW container for protons not found!"<<endl;
 849     status = kFALSE;
 850   }
 851
 852   Int_t nSteps = corrfwContainerProtons->GetNStep();
 853   TH2D *gYPt[4];
 854   //currently the GRID is formed by the y-pT parameters
 855   //Add Vz as a next step
 856   Int_t iRap = 0, iPt = 1;
 857   AliCFEffGrid *effProtonsStep0Step1 = new AliCFEffGrid("eff10",
 858                                          "effProtonsStep0Step1",
 859                                          *corrfwContainerProtons);
 860   effProtonsStep0Step1->CalculateEfficiency(1,0); //eff= step1/step0
 861   fEffGridListProtons->Add(effProtonsStep0Step1);
 862   gYPt[0] = effProtonsStep0Step1->Project(iRap,iPt);
 863   fCorrectionListProtons2D->Add(gYPt[0]);
 864
 865   AliCFEffGrid *effProtonsStep0Step2 = new AliCFEffGrid("eff20",
 866                                          "effProtonsStep0Step2",
 867                                          *corrfwContainerProtons);
 868   effProtonsStep0Step2->CalculateEfficiency(2,0); //eff= step2/step0
 869   fEffGridListProtons->Add(effProtonsStep0Step2);
 870   gYPt[1] = effProtonsStep0Step2->Project(iRap,iPt);
 871   fCorrectionListProtons2D->Add(gYPt[1]);
 872
 873   AliCFEffGrid *effProtonsStep0Step3 = new AliCFEffGrid("eff30",
 874                                          "effProtonsStep0Step3",
 875                                          *corrfwContainerProtons);
 876   effProtonsStep0Step3->CalculateEfficiency(3,0); //eff= step1/step0
 877   fEffGridListProtons->Add(effProtonsStep0Step3);
 878   gYPt[2] = effProtonsStep0Step3->Project(iRap,iPt);
 879   fCorrectionListProtons2D->Add(gYPt[2]);
 880
 881   TH1D *gEfficiency[2][3]; //efficiency as a function of pT and of y (raws-[2])
 882   TH1D *gCorrection[2][3]; //efficiency as a function of pT and of y (raws-[2])
 883   TString gTitle;
 884   //Get the projection of the efficiency maps
 885   for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
 886     gEfficiency[iParameter][0] = effProtonsStep0Step1->Project(iParameter);
 887     gTitle = "ProtonsEfficiency_Parameter"; gTitle += iParameter+1;
 888     gTitle += "_Step0_Step1";
 889     gEfficiency[iParameter][0]->SetName(gTitle.Data());
 890     fEfficiencyListProtons1D->Add(gEfficiency[iParameter][0]);
 891     gTitle = "ProtonsCorrection_Parameter"; gTitle += iParameter+1;
 892     gTitle += "_Step0_Step1";
 893     gCorrection[iParameter][0] = new TH1D(gTitle.Data(),
 894                                           gTitle.Data(),
 895                                           gEfficiency[iParameter][0]->GetNbinsX(),
 896                                           gEfficiency[iParameter][0]->GetXaxis()->GetXmin(),
 897                                           gEfficiency[iParameter][0]->GetXaxis()->GetXmax());
 898     //initialisation of the correction
 899     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][0]->GetNbinsX(); iBin++)
 900       gCorrection[iParameter][0]->SetBinContent(iBin,1.0);
 901
 902     gEfficiency[iParameter][1] = effProtonsStep0Step2->Project(iParameter);
 903     gTitle = "ProtonsEfficiency_Parameter"; gTitle += iParameter+1;
 904     gTitle += "_Step0_Step2";
 905     gEfficiency[iParameter][1]->SetName(gTitle.Data());
 906     fEfficiencyListProtons1D->Add(gEfficiency[iParameter][1]);
 907     gTitle = "ProtonsCorrection_Parameter"; gTitle += iParameter+1;
 908     gTitle += "_Step0_Step2";
 909     gCorrection[iParameter][1] = new TH1D(gTitle.Data(),
 910                                           gTitle.Data(),
 911                                           gEfficiency[iParameter][1]->GetNbinsX(),
 912                                           gEfficiency[iParameter][1]->GetXaxis()->GetXmin(),
 913                                           gEfficiency[iParameter][1]->GetXaxis()->GetXmax());
 914     //initialisation of the correction
 915     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][1]->GetNbinsX(); iBin++)
 916       gCorrection[iParameter][1]->SetBinContent(iBin,1.0);
 917
 918     gEfficiency[iParameter][2] = effProtonsStep0Step3->Project(iParameter);
 919     gTitle = "ProtonsEfficiency_Parameter"; gTitle += iParameter+1;
 920     gTitle += "_Step0_Step3";
 921     gEfficiency[iParameter][2]->SetName(gTitle.Data());
 922     fEfficiencyListProtons1D->Add(gEfficiency[iParameter][2]);
 923     gTitle = "ProtonsCorrection_Parameter"; gTitle += iParameter+1;
 924     gTitle += "_Step0_Step3";
 925     gCorrection[iParameter][2] = new TH1D(gTitle.Data(),
 926                                           gTitle.Data(),
 927                                           gEfficiency[iParameter][2]->GetNbinsX(),
 928                                           gEfficiency[iParameter][2]->GetXaxis()->GetXmin(),
 929                                           gEfficiency[iParameter][2]->GetXaxis()->GetXmax());
 930     //initialisation of the correction
 931     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][2]->GetNbinsX(); iBin++)
 932       gCorrection[iParameter][2]->SetBinContent(iBin,1.0);
 933   }//parameter loop
 934   //Calculate the 1D correction parameters as a function of y and pT
 935   for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
 936     for(Int_t iStep = 1; iStep < nSteps; iStep++) {
 937       gCorrection[iParameter][iStep-1]->Divide(gEfficiency[iParameter][iStep-1]);
 938       fCorrectionListProtons1D->Add(gCorrection[iParameter][iStep-1]);
 939     }
 940   }
 941
 942   //________________________________________//
 943   //AntiProtons
 944   fEffGridListAntiProtons = new TList();
 945   fCorrectionListAntiProtons2D = new TList();
 946   fEfficiencyListAntiProtons1D = new TList();
 947   fCorrectionListAntiProtons1D = new TList();
 948
 949   AliCFContainer *corrfwContainerAntiProtons = (AliCFContainer*) (file->Get("containerAntiProtons"));
 950   if(!corrfwContainerAntiProtons) {
 951     cout<<"CORRFW container for antiprotons not found!"<<endl;
 952     status = kFALSE;
 953   }
 954
 955   nSteps = corrfwContainerAntiProtons->GetNStep();
 956   //currently the GRID is formed by the y-pT parameters
 957   //Add Vz as a next step
 958   AliCFEffGrid *effAntiProtonsStep0Step1 = new AliCFEffGrid("eff10",
 959                                          "effAntiProtonsStep0Step1",
 960                                          *corrfwContainerAntiProtons);
 961   effAntiProtonsStep0Step1->CalculateEfficiency(1,0); //eff= step1/step0
 962   fEffGridListAntiProtons->Add(effAntiProtonsStep0Step1);
 963   gYPt[0] = effAntiProtonsStep0Step1->Project(iRap,iPt);
 964   fCorrectionListAntiProtons2D->Add(gYPt[0]);
 965
 966   AliCFEffGrid *effAntiProtonsStep0Step2 = new AliCFEffGrid("eff20",
 967                                          "effAntiProtonsStep0Step2",
 968                                          *corrfwContainerAntiProtons);
 969   effAntiProtonsStep0Step2->CalculateEfficiency(2,0); //eff= step2/step0
 970   fEffGridListAntiProtons->Add(effAntiProtonsStep0Step2);
 971   gYPt[1] = effAntiProtonsStep0Step2->Project(iRap,iPt);
 972   fCorrectionListAntiProtons2D->Add(gYPt[1]);
 973
 974   AliCFEffGrid *effAntiProtonsStep0Step3 = new AliCFEffGrid("eff30",
 975                                          "effAntiProtonsStep0Step3",
 976                                          *corrfwContainerAntiProtons);
 977   effAntiProtonsStep0Step3->CalculateEfficiency(3,0); //eff= step1/step0
 978   fEffGridListAntiProtons->Add(effAntiProtonsStep0Step3);
 979   gYPt[2] = effAntiProtonsStep0Step3->Project(iRap,iPt);
 980   fCorrectionListAntiProtons2D->Add(gYPt[2]);
 981
 982   //Get the projection of the efficiency maps
 983   for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
 984     gEfficiency[iParameter][0] = effAntiProtonsStep0Step1->Project(iParameter);
 985     gTitle = "AntiProtonsEfficiency_Parameter"; gTitle += iParameter+1;
 986     gTitle += "_Step0_Step1";
 987     gEfficiency[iParameter][0]->SetName(gTitle.Data());
 988     fEfficiencyListAntiProtons1D->Add(gEfficiency[iParameter][0]);
 989     gTitle = "AntiProtonsCorrection_Parameter"; gTitle += iParameter+1;
 990     gTitle += "_Step0_Step1";
 991     gCorrection[iParameter][0] = new TH1D(gTitle.Data(),
 992                                           gTitle.Data(),
 993                                           gEfficiency[iParameter][0]->GetNbinsX(),
 994                                           gEfficiency[iParameter][0]->GetXaxis()->GetXmin(),
 995                                           gEfficiency[iParameter][0]->GetXaxis()->GetXmax());
 996     //initialisation of the correction
 997     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][0]->GetNbinsX(); iBin++)
 998       gCorrection[iParameter][0]->SetBinContent(iBin,1.0);
 999
1000     gEfficiency[iParameter][1] = effAntiProtonsStep0Step2->Project(iParameter);
1001     gTitle = "AntiProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1002     gTitle += "_Step0_Step2";
1003     gEfficiency[iParameter][1]->SetName(gTitle.Data());
1004     fEfficiencyListAntiProtons1D->Add(gEfficiency[iParameter][1]);
1005     gTitle = "AntiProtonsCorrection_Parameter"; gTitle += iParameter+1;
1006     gTitle += "_Step0_Step2";
1007     gCorrection[iParameter][1] = new TH1D(gTitle.Data(),
1008                                           gTitle.Data(),
1009                                           gEfficiency[iParameter][1]->GetNbinsX(),
1010                                           gEfficiency[iParameter][1]->GetXaxis()->GetXmin(),
1011                                           gEfficiency[iParameter][1]->GetXaxis()->GetXmax());
1012     //initialisation of the correction
1013     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][1]->GetNbinsX(); iBin++)
1014       gCorrection[iParameter][1]->SetBinContent(iBin,1.0);
1015
1016     gEfficiency[iParameter][2] = effAntiProtonsStep0Step3->Project(iParameter);
1017     gTitle = "AntiProtonsEfficiency_Parameter"; gTitle += iParameter+1;
1018     gTitle += "_Step0_Step3";
1019     gEfficiency[iParameter][2]->SetName(gTitle.Data());
1020     fEfficiencyListAntiProtons1D->Add(gEfficiency[iParameter][2]);
1021     gTitle = "AntiProtonsCorrection_Parameter"; gTitle += iParameter+1;
1022     gTitle += "_Step0_Step3";
1023     gCorrection[iParameter][2] = new TH1D(gTitle.Data(),
1024                                           gTitle.Data(),
1025                                           gEfficiency[iParameter][2]->GetNbinsX(),
1026                                           gEfficiency[iParameter][2]->GetXaxis()->GetXmin(),
1027                                           gEfficiency[iParameter][2]->GetXaxis()->GetXmax());
1028     //initialisation of the correction
1029     for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][2]->GetNbinsX(); iBin++)
1030       gCorrection[iParameter][2]->SetBinContent(iBin,1.0);
1031   }//parameter loop
1032   //Calculate the 1D correction parameters as a function of y and pT
1033   for(Int_t iParameter = 0; iParameter < 2; iParameter++) {
1034     for(Int_t iStep = 1; iStep < nSteps; iStep++) {
1035       gCorrection[iParameter][iStep-1]->Divide(gEfficiency[iParameter][iStep-1]);
1036       fCorrectionListAntiProtons1D->Add(gCorrection[iParameter][iStep-1]);
1037     }
1038   }
1039
1040   return status;
1041 }
1042
1043
1044
1045
1046
1047
1048
1049
1050