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