]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/SPECTRA/AliProtonAnalysis.cxx
First attempt to provide an interface to the correction framework
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliProtonAnalysis.cxx
index 50504ff4d2bc8d05e1f401bbea04a75b8acf13e9..752fdbc56f99added58a1fe1ec1d80a2c0eae185 100644 (file)
 #include <Riostream.h>
 #include <TFile.h>
 #include <TSystem.h>
+#include <TF1.h>
 #include <TH2F.h>
 #include <TH1D.h>
+#include <TH1I.h>
+#include <TParticle.h>
 
 #include "AliProtonAnalysis.h"
 
+#include <AliExternalTrackParam.h>
+#include <AliAODEvent.h>
 #include <AliESDEvent.h>
 #include <AliLog.h>
+#include <AliPID.h>
+#include <AliStack.h>
+#include <AliCFContainer.h>
+#include <AliCFEffGrid.h>
+#include <AliCFEffGrid.h>
 
 ClassImp(AliProtonAnalysis)
 
@@ -45,9 +55,16 @@ AliProtonAnalysis::AliProtonAnalysis() :
   fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
   fMaxSigmaToVertexFlag(kFALSE),
   fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
-  fHistYPtProtons(0), fHistYPtAntiProtons(0) {
+  fFunctionProbabilityFlag(kFALSE), 
+  fElectronFunction(0), fMuonFunction(0),
+  fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
+  fUseTPCOnly(kFALSE), fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
+  fCorrectionList2D(0), fEfficiencyList1D(0), fCorrectionList1D(0) {
   //Default constructor
   for(Int_t i = 0; i < 5; i++) fPartFrac[i] = 0.0;
+  fCorrectionList2D = new TList(); 
+  fEfficiencyList1D = new TList(); 
+  fCorrectionList1D = new TList();
 }
 
 //____________________________________________________________________//
@@ -55,9 +72,24 @@ AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY
   TObject(),
   fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),
   fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),
-  fHistYPtProtons(0), fHistYPtAntiProtons(0) {
+  fMinTPCClusters(0), fMinITSClusters(0),
+  fMaxChi2PerTPCCluster(0), fMaxChi2PerITSCluster(0),
+  fMaxCov11(0), fMaxCov22(0), fMaxCov33(0), fMaxCov44(0), fMaxCov55(0),
+  fMaxSigmaToVertex(0),
+  fMinTPCClustersFlag(kFALSE), fMinITSClustersFlag(kFALSE),
+  fMaxChi2PerTPCClusterFlag(kFALSE), fMaxChi2PerITSClusterFlag(kFALSE),
+  fMaxCov11Flag(kFALSE), fMaxCov22Flag(kFALSE), fMaxCov33Flag(kFALSE), fMaxCov44Flag(kFALSE), fMaxCov55Flag(kFALSE),
+  fMaxSigmaToVertexFlag(kFALSE),
+  fITSRefitFlag(kFALSE), fTPCRefitFlag(kFALSE),
+  fFunctionProbabilityFlag(kFALSE), 
+  fElectronFunction(0), fMuonFunction(0),
+  fPionFunction(0), fKaonFunction(0), fProtonFunction(0),
+  fUseTPCOnly(kFALSE), fHistEvents(0), fHistYPtProtons(0), fHistYPtAntiProtons(0),
+  fCorrectionList2D(0), fEfficiencyList1D(0), fCorrectionList1D(0) {
   //Default constructor
 
+  fHistEvents = new TH1I("fHistEvents","Analyzed events",1,0,1);
+
   fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
   fHistYPtProtons->SetStats(kTRUE);
   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
@@ -74,7 +106,12 @@ AliProtonAnalysis::AliProtonAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY
 //____________________________________________________________________//
 AliProtonAnalysis::~AliProtonAnalysis() {
   //Default destructor
-  
+  if(fHistEvents) delete fHistEvents;
+  if(fHistYPtProtons) delete fHistYPtProtons;
+  if(fHistYPtAntiProtons) delete fHistYPtAntiProtons;
+  if(fCorrectionList2D) delete fCorrectionList2D;
+  if(fEfficiencyList1D) delete fEfficiencyList1D;
+  if(fCorrectionList1D) delete fCorrectionList1D;
 }
 
 //____________________________________________________________________//
@@ -86,6 +123,8 @@ void AliProtonAnalysis::InitHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHig
   fMinPt = fLowPt;
   fMaxPt = fHighPt;
 
+  fHistEvents = new TH1I("fHistEvents","Anallyzed events",1,0,1);
+
   fHistYPtProtons = new TH2F("fHistYPtProtons","y-Pt Protons",fNBinsY,fMinY,fMaxY,fNBinsPt,fMinPt,fMaxPt);
   fHistYPtProtons->SetStats(kTRUE);
   fHistYPtProtons->GetYaxis()->SetTitle("P_{T} [GeV]");
@@ -100,58 +139,96 @@ void AliProtonAnalysis::InitHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHig
 }
 
 //____________________________________________________________________//
-void AliProtonAnalysis::ReadFromFile(const char* filename) {
+Bool_t AliProtonAnalysis::ReadFromFile(const char* filename) {
+  Bool_t status = kTRUE;
+
   TFile *file = TFile::Open(filename);
-  fHistYPtProtons = (TH2F *)file->Get("fHistYPtProtons");
-  fHistYPtAntiProtons = (TH2F *)file->Get("fHistYPtAntiProtons");
-  fHistYPtProtons->Sumw2();
-  fHistYPtAntiProtons->Sumw2();
+  if(!file) {
+    cout<<"Could not find the input file "<<filename<<endl;
+    status = kFALSE;
+  }
+
+  TList *list = (TList *)file->Get("outputList1");
+  if(list) {
+    cout<<"Retrieving objects from the list "<<list->GetName()<<"..."<<endl; 
+    fHistYPtProtons = (TH2F *)list->At(0);
+    fHistYPtAntiProtons = (TH2F *)list->At(1);
+    fHistEvents = (TH1I *)list->At(2);
+  }
+  else if(!list) {
+    cout<<"Retrieving objects from the file... "<<endl;
+    fHistYPtProtons = (TH2F *)file->Get("fHistYPtProtons");
+    fHistYPtAntiProtons = (TH2F *)file->Get("fHistYPtAntiProtons");
+    fHistEvents = (TH1I *)file->Get("fHistEvents");
+  }
+  if((!fHistYPtProtons)||(!fHistYPtAntiProtons)||(!fHistEvents)) {
+    cout<<"Input containers were not found!!!"<<endl;
+    status = kFALSE;
+  }
+  else {
+    fHistYPtProtons->Sumw2();
+    fHistYPtAntiProtons->Sumw2();
+  }
+
+  return status;
 }
 
 //____________________________________________________________________//
 TH1D *AliProtonAnalysis::GetProtonYHistogram() {
+  Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
   TH1D *fYProtons = (TH1D *)fHistYPtProtons->ProjectionX("fYProtons",0,fHistYPtProtons->GetYaxis()->GetNbins(),"e"); 
   fYProtons->SetStats(kFALSE);
-  fYProtons->GetYaxis()->SetTitle("dN/dy");
+  fYProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
   fYProtons->SetTitle("dN/dy protons");
   fYProtons->SetMarkerStyle(kFullCircle);
   fYProtons->SetMarkerColor(4);
+  if(nAnalyzedEvents > 0)
+    fYProtons->Scale(1./nAnalyzedEvents);
 
   return fYProtons;
 }
 
 //____________________________________________________________________//
 TH1D *AliProtonAnalysis::GetAntiProtonYHistogram() {
+  Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
   TH1D *fYAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionX("fYAntiProtons",0,fHistYPtAntiProtons->GetYaxis()->GetNbins(),"e"); 
   fYAntiProtons->SetStats(kFALSE);
-  fYAntiProtons->GetYaxis()->SetTitle("dN/dy");
+  fYAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dy)");
   fYAntiProtons->SetTitle("dN/dy antiprotons");
   fYAntiProtons->SetMarkerStyle(kFullCircle);
   fYAntiProtons->SetMarkerColor(4);
+  if(nAnalyzedEvents > 0)
+    fYAntiProtons->Scale(1./nAnalyzedEvents);
 
   return fYAntiProtons;
 }
 
 //____________________________________________________________________//
 TH1D *AliProtonAnalysis::GetProtonPtHistogram() {
+  Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
   TH1D *fPtProtons = (TH1D *)fHistYPtProtons->ProjectionY("fPtProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e"); 
   fPtProtons->SetStats(kFALSE);
-  fPtProtons->GetYaxis()->SetTitle("dN/dP_{T}");
+  fPtProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
   fPtProtons->SetTitle("dN/dPt protons");
   fPtProtons->SetMarkerStyle(kFullCircle);
   fPtProtons->SetMarkerColor(4);
+  if(nAnalyzedEvents > 0)
+    fPtProtons->Scale(1./nAnalyzedEvents);
 
   return fPtProtons;
 }
 
 //____________________________________________________________________//
 TH1D *AliProtonAnalysis::GetAntiProtonPtHistogram() {
+  Int_t nAnalyzedEvents = GetNumberOfAnalyzedEvents();
   TH1D *fPtAntiProtons = (TH1D *)fHistYPtAntiProtons->ProjectionY("fPtAntiProtons",0,fHistYPtProtons->GetXaxis()->GetNbins(),"e"); 
   fPtAntiProtons->SetStats(kFALSE);
-  fPtAntiProtons->GetYaxis()->SetTitle("dN/dP_{T}");
+  fPtAntiProtons->GetYaxis()->SetTitle("(1/N_{events})(dN/dP_{T})");
   fPtAntiProtons->SetTitle("dN/dPt antiprotons");
   fPtAntiProtons->SetMarkerStyle(kFullCircle);
   fPtAntiProtons->SetMarkerColor(4);
+  if(nAnalyzedEvents > 0)
+    fPtAntiProtons->Scale(1./nAnalyzedEvents);
 
   return fPtAntiProtons;
 }
@@ -240,35 +317,150 @@ TH1D *AliProtonAnalysis::GetPtAsymmetryHistogram() {
   return hAsymmetryPt;
 }
 
+//____________________________________________________________________//
+Double_t AliProtonAnalysis::GetParticleFraction(Int_t i, Double_t p) {
+  Double_t partFrac=0;
+  if(fFunctionProbabilityFlag) {
+    if(i == 0) partFrac = fElectronFunction->Eval(p);
+    if(i == 1) partFrac = fMuonFunction->Eval(p);
+    if(i == 2) partFrac = fPionFunction->Eval(p);
+    if(i == 3) partFrac = fKaonFunction->Eval(p);
+    if(i == 4) partFrac = fProtonFunction->Eval(p);
+  }
+  else partFrac = fPartFrac[i];
+
+  return partFrac;
+}
+
 //____________________________________________________________________//
 void AliProtonAnalysis::Analyze(AliESDEvent* fESD) {
-  //Main analysis part
+  //Main analysis part - ESD
+  fHistEvents->Fill(0); //number of analyzed events
+  Double_t Pt = 0.0, P = 0.0;
   Int_t nGoodTracks = fESD->GetNumberOfTracks();
   for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
     AliESDtrack* track = fESD->GetTrack(iTracks);
-    if(IsAccepted(track)) {
-      Double_t Pt = track->Pt();
+    Double_t probability[5];
+
+    if(IsAccepted(track)) {    
+      if(fUseTPCOnly) {
+       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
+       if(!tpcTrack) continue;
+       Pt = tpcTrack->Pt();
+       P = tpcTrack->P();
        
-         //pid
-      Double_t probability[5];
-         track->GetESDpid(probability);
-      Double_t rcc = 0.0;
-      for(Int_t i = 0; i < AliPID::kSPECIES; i++) rcc += probability[i]*fPartFrac[i];
-      if(rcc == 0.0) continue;
-      Double_t w[5];
-      for(Int_t i = 0; i < AliPID::kSPECIES; i++) w[i] = probability[i]*fPartFrac[i]/rcc;
-      Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
-      if(fParticleType == 4) {
-        if(track->Charge() > 0) fHistYPtProtons->Fill(Rapidity(track),Pt);
-        else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(Rapidity(track),Pt);
-      }//proton check
+       //pid
+       track->GetTPCpid(probability);
+       Double_t rcc = 0.0;
+       for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
+         rcc += probability[i]*GetParticleFraction(i,P);
+       if(rcc == 0.0) continue;
+       Double_t w[5];
+       for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
+         w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
+       Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
+       if(fParticleType == 4) {
+         if(tpcTrack->Charge() > 0) 
+           fHistYPtProtons->Fill(Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz()),Pt);
+         else if(tpcTrack->Charge() < 0) 
+           fHistYPtAntiProtons->Fill(Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz()),Pt);
+       }//proton check
+      }//TPC only tracks
+      else if(!fUseTPCOnly) {
+       Pt = track->Pt();
+       P = track->P();
+       
+       //pid
+       track->GetESDpid(probability);
+       Double_t rcc = 0.0;
+       for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
+         rcc += probability[i]*GetParticleFraction(i,P);
+       if(rcc == 0.0) continue;
+       Double_t w[5];
+       for(Int_t i = 0; i < AliPID::kSPECIES; i++) 
+         w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
+       Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIES,w);
+       if(fParticleType == 4) {
+         //cout<<"(Anti)protons found..."<<endl;
+         if(track->Charge() > 0) 
+           fHistYPtProtons->Fill(Rapidity(track->Px(),track->Py(),track->Pz()),Pt);
+         else if(track->Charge() < 0) 
+           fHistYPtAntiProtons->Fill(Rapidity(track->Px(),track->Py(),track->Pz()),Pt);
+       }//proton check 
+      }//combined tracking
     }//cuts
   }//track loop 
 }
 
+//____________________________________________________________________//
+void AliProtonAnalysis::Analyze(AliAODEvent* fAOD) {
+  //Main analysis part - AOD
+  fHistEvents->Fill(0); //number of analyzed events
+  Int_t nGoodTracks = fAOD->GetNumberOfTracks();
+  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) {
+    AliAODTrack* track = fAOD->GetTrack(iTracks);
+    Double_t Pt = track->Pt();
+    Double_t P = track->P();
+    
+    //pid
+    Double_t probability[10];
+    track->GetPID(probability);
+    Double_t rcc = 0.0;
+    for(Int_t i = 0; i < AliPID::kSPECIESN; i++) rcc += probability[i]*GetParticleFraction(i,P);
+    if(rcc == 0.0) continue;
+    Double_t w[10];
+    for(Int_t i = 0; i < AliPID::kSPECIESN; i++) w[i] = probability[i]*GetParticleFraction(i,P)/rcc;
+    Long64_t fParticleType = TMath::LocMax(AliPID::kSPECIESN,w);
+    if(fParticleType == 4) {
+      if(track->Charge() > 0) fHistYPtProtons->Fill(track->Y(fParticleType),Pt);
+      else if(track->Charge() < 0) fHistYPtAntiProtons->Fill(track->Y(fParticleType),Pt);
+    }//proton check
+  }//track loop 
+}
+
+//____________________________________________________________________//
+void AliProtonAnalysis::Analyze(AliStack* stack) {
+  //Main analysis part - MC
+  fHistEvents->Fill(0); //number of analyzed events
+  for(Int_t i = 0; i < stack->GetNprimary(); i++) {
+    TParticle *particle = stack->Particle(i);
+    if(particle->Pt() < 0.1) continue;
+    if(TMath::Abs(particle->Eta()) > 1.0) continue;
+    Int_t pdgcode = particle->GetPdgCode();
+    if(pdgcode == 2212) fHistYPtProtons->Fill(Rapidity(particle->Px(),
+                                                      particle->Py(),
+                                                      particle->Pz()),
+                                             particle->Pt());
+    if(pdgcode == -2212) fHistYPtAntiProtons->Fill(Rapidity(particle->Px(),
+                                                           particle->Py(),
+                                                           particle->Pz()),
+                                                  particle->Pt());
+  }//particle loop                                                                  
+}
+
 //____________________________________________________________________//
 Bool_t AliProtonAnalysis::IsAccepted(AliESDtrack* track) {
   // Checks if the track is excluded from the cuts
+  Double_t Pt = 0.0, Px = 0.0, Py = 0.0, Pz = 0.0;
+  if(fUseTPCOnly) {
+    AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();
+    if(!tpcTrack) {
+      Pt = 0.0; Px = 0.0; Py = 0.0; Pz = 0.0;
+    }
+    else {
+      Pt = tpcTrack->Pt();
+      Px = tpcTrack->Px();
+      Py = tpcTrack->Py();
+      Pz = tpcTrack->Pz();
+    }
+  }
+  else{
+    Pt = track->Pt();
+    Px = track->Px();
+    Py = track->Py();
+    Pz = track->Pz();
+  }
+     
   Int_t  fIdxInt[200];
   Int_t nClustersITS = track->GetITSclusters(fIdxInt);
   Int_t nClustersTPC = track->GetTPCclusters(fIdxInt);
@@ -281,8 +473,6 @@ Bool_t AliProtonAnalysis::IsAccepted(AliESDtrack* track) {
   Double_t extCov[15];
   track->GetExternalCovariance(extCov);
 
-  Double_t Pt = track->Pt();
-
   if(fMinITSClustersFlag)
     if(nClustersITS < fMinITSClusters) return kFALSE;
   if(fMinTPCClustersFlag)
@@ -309,6 +499,7 @@ Bool_t AliProtonAnalysis::IsAccepted(AliESDtrack* track) {
     if ((track->GetStatus() & AliESDtrack::kTPCrefit) == 0) return kFALSE;
 
   if((Pt < fMinPt) || (Pt > fMaxPt)) return kFALSE;
+  if((Rapidity(Px,Py,Pz) < fMinY) || (Rapidity(Px,Py,Pz) > fMaxY)) return kFALSE;
 
   return kTRUE;
 }
@@ -320,7 +511,11 @@ Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) {
   Float_t b[2];
   Float_t bRes[2];
   Float_t bCov[3];
-  esdTrack->GetImpactParameters(b,bCov);
+  if(fUseTPCOnly) 
+    esdTrack->GetImpactParametersTPC(b,bCov);
+  else
+    esdTrack->GetImpactParameters(b,bCov);
+  
   if (bCov[0]<=0 || bCov[2]<=0) {
     //AliDebug(1, "Estimated b resolution lower or equal zero!");
     bCov[0]=0; bCov[2]=0;
@@ -339,16 +534,165 @@ Float_t AliProtonAnalysis::GetSigmaToVertex(AliESDtrack* esdTrack) {
   return d;
 }
 
-Double_t AliProtonAnalysis::Rapidity(AliESDtrack *track) {
+//____________________________________________________________________//
+Double_t AliProtonAnalysis::Rapidity(Double_t Px, Double_t Py, Double_t Pz) {
+  //returns the rapidity of the proton - to be removed
   Double_t fMass = 9.38270000000000048e-01;
   
-  Double_t P = TMath::Sqrt(TMath::Power(track->Px(),2) + 
-                           TMath::Power(track->Py(),2) + 
-                                                  TMath::Power(track->Pz(),2));
+  Double_t P = TMath::Sqrt(TMath::Power(Px,2) + 
+                           TMath::Power(Py,2) + 
+                          TMath::Power(Pz,2));
   Double_t energy = TMath::Sqrt(P*P + fMass*fMass);
   Double_t y = -999;
-  if(energy != track->Pz()
-    y = 0.5*TMath::Log((energy + track->Pz())/(energy - track->Pz()));
+  if(energy != Pz
+    y = 0.5*TMath::Log((energy + Pz)/(energy - Pz));
 
   return y;
 }
+
+//____________________________________________________________________//
+Bool_t AliProtonAnalysis::PrintMean(TH1 *hist, Double_t edge) {
+  //calculates the mean value of the ratio/asymmetry within \pm edge
+  Double_t sum = 0.0;
+  Int_t nentries = 0;
+  //calculate the mean
+  for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
+    Double_t x = hist->GetBinCenter(i+1);
+    Double_t y = hist->GetBinContent(i+1);
+    if(TMath::Abs(x) < edge) {
+      sum += y;
+      nentries += 1;
+    }
+  }
+  Double_t mean = 0.0;
+  if(nentries != 0)
+    mean = sum/nentries;
+
+  //calculate the error
+  for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
+    Double_t x = hist->GetBinCenter(i+1);
+    Double_t y = hist->GetBinContent(i+1);
+    if(TMath::Abs(x) < edge) {
+      sum += TMath::Power((mean - y),2);
+      nentries += 1;
+    }
+  }
+
+  Double_t error = 0.0;
+  if(nentries != 0)
+    error =  TMath::Sqrt(sum)/nentries;
+
+  cout<<"========================================="<<endl;
+  cout<<"Input distribution: "<<hist->GetName()<<endl;
+  cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
+  cout<<"Mean value :"<<mean<<endl;
+  cout<<"Error: "<<error<<endl;
+  cout<<"========================================="<<endl;
+
+  return 0;
+}
+
+//____________________________________________________________________//
+Bool_t AliProtonAnalysis::PrintYields(TH1 *hist, Double_t edge) {
+  //calculates the (anti)proton yields within the \pm edge
+  Double_t sum = 0.0, sumerror = 0.0;
+  Double_t error = 0.0;
+  for(Int_t i = 0; i < hist->GetXaxis()->GetNbins(); i++) {
+    Double_t x = hist->GetBinCenter(i+1);
+    Double_t y = hist->GetBinContent(i+1);
+    if(TMath::Abs(x) < edge) {
+      sum += y;
+      sumerror += TMath::Power(hist->GetBinError(i+1),2); 
+    }
+  }
+
+  error = TMath::Sqrt(sumerror);
+
+  cout<<"========================================="<<endl;
+  cout<<"Input distribution: "<<hist->GetName()<<endl;
+  cout<<"Interval used: -"<<edge<<" -> "<<edge<<endl;
+  cout<<"Yields :"<<sum<<endl;
+  cout<<"Error: "<<error<<endl;
+  cout<<"========================================="<<endl;
+
+  return 0;
+}
+
+//____________________________________________________________________//
+Bool_t AliProtonAnalysis::ReadCorrectionContainer(const char* filename) {
+  // Reads the outout of the correction framework task
+  // Creates the correction maps
+  // Puts the results in the different TList objects
+  Bool_t status = kTRUE;
+
+  TFile *file = TFile::Open(filename);
+  if(!file) {
+    cout<<"Could not find the input CORRFW file "<<filename<<endl;
+    status = kFALSE;
+  }
+
+  AliCFContainer *corrfwContainer = (AliCFContainer*) (file->Get("container"));
+  if(!corrfwContainer) {
+    cout<<"CORRFW container not found!"<<endl;
+    status = kFALSE;
+  }
+  
+  Int_t nSteps = corrfwContainer->GetNStep();
+  TH2D *gYPt[4];
+  //currently the GRID is formed by the y-pT parameters
+  //Add Vz as a next step
+  Int_t iRap = 0, iPt = 1;
+  for(Int_t iStep = 0; iStep < nSteps; iStep++) {
+    gYPt[iStep] = corrfwContainer->ShowProjection(iRap,iPt,iStep);
+    //fCorrectionList2D->Add(gYPt[iStep]);
+  }
+
+  //construct the efficiency grid from the data container                                           
+  TString gTitle = 0;
+  AliCFEffGrid *efficiency[3]; //The efficiency array has nStep-1 entries!!!
+  TH1D *gEfficiency[2][3]; //efficiency as a function of pT and of y (raws - [2]) 
+  TH1D *gCorrection[2][3]; //efficiency as a function of pT and of y (raws - [2]) 
+
+  //Get the 2D efficiency maps
+  for(Int_t iStep = 1; iStep < nSteps; iStep++) { 
+    gTitle = "Efficiency_Step0_Step"; gTitle += iStep; 
+    efficiency[iStep] = new AliCFEffGrid(gTitle.Data(),
+                                        gTitle.Data(),*corrfwContainer);
+    efficiency[iStep]->CalculateEfficiency(iStep,0); //eff= step[i]/step0
+    fCorrectionList2D->Add(efficiency[iStep]);  
+  }
+  //Get the projection of the efficiency maps
+  for(Int_t iParameter = 0; iParameter < 2; iParameter++) { 
+    for(Int_t iStep = 1; iStep < nSteps; iStep++) { 
+      gEfficiency[iParameter][iStep-1] = efficiency[iStep]->Project(iParameter);
+      fEfficiencyList1D->Add(gEfficiency[iParameter][iStep-1]);  
+      gTitle = "Correction_Parameter"; gTitle += iParameter+1;
+      gTitle += "_Step0_Step"; gTitle += iStep; 
+      gCorrection[iParameter][iStep-1] = new TH1D(gTitle.Data(),
+                                                  gTitle.Data(),
+                                                  gEfficiency[iParameter][iStep-1]->GetNbinsX(),
+                                                  gEfficiency[iParameter][iStep-1]->GetXaxis()->GetXmin(),
+                                                  gEfficiency[iParameter][iStep-1]->GetXaxis()->GetXmax());
+      //initialisation of the correction
+      for(Int_t iBin = 1; iBin <= gEfficiency[iParameter][iStep-1]->GetNbinsX(); iBin++)
+       gCorrection[iParameter][iStep-1]->SetBinContent(iBin,1.0);
+    }//step loop
+  }//parameter loop
+  //Calculate the 1D correction parameters as a function of y and pT
+  for(Int_t iParameter = 0; iParameter < 2; iParameter++) { 
+    for(Int_t iStep = 1; iStep < nSteps; iStep++) { 
+      gCorrection[iParameter][iStep-1]->Divide(gEfficiency[iParameter][iStep-1]);
+      fCorrectionList1D->Add(gCorrection[iParameter][iStep-1]);  
+    }
+  }
+}
+
+
+
+
+
+
+
+
+
+