PWG2/SPECTRA/AliProtonAnalysis.cxx

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