added 2dim unfolding histograms to jet task, cosmetic changes in steering macro
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 9 Jan 2009 10:51:40 +0000 (10:51 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 9 Jan 2009 10:51:40 +0000 (10:51 +0000)
PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx
PWG4/JetTasks/AliAnalysisTaskJetSpectrum.h
PWG4/JetTasks/AliJetSpectrumUnfolding.cxx
PWG4/macros/AnalysisTrainCAF.C
PWG4/macros/runJetSpectrumUnfolding.C

index 75c432c..26948e5 100644 (file)
@@ -53,6 +53,8 @@
 
 ClassImp(AliAnalysisTaskJetSpectrum)
 
+const Float_t AliAnalysisTaskJetSpectrum::fgkJetNpartCut[AliAnalysisTaskJetSpectrum::kMaxCorrelation] = {5,10,1E+09};
+
 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
   fJetHeaderRec(0x0),
   fJetHeaderGen(0x0),
@@ -71,14 +73,24 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
   fh1PtHard_Trials(0x0),
   fh1NGenJets(0x0),
   fh1NRecJets(0x0),
-  fHistList(0x0)
+  fHistList(0x0) , 
+  ////////////////
+  fh1JetMultiplicity(0x0) ,     
+  fh2ERecZRec(0x0) ,
+  fh2EGenZGen(0x0) ,
+  fh2Efficiency(0x0) ,
+  fh3EGenERecN(0x0) 
+  //////////////// 
 {
   // Default constructor
     for(int ij  = 0;ij<kMaxJets;++ij){
-    fh1E[ij] =  fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
-    fh2PtFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] = 0;
-    fh3PtRecGenHard[ij] =  fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
-  }
+      fh1E[ij] =  fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
+      fh2PtFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] = 0;
+      fh3PtRecGenHard[ij] =  fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
+    }
+    for(int ic = 0;ic < kMaxCorrelation;ic++){
+      fhnCorrelation[ic] = 0;
+    }
   
 }
 
@@ -101,7 +113,14 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
   fh1PtHard_Trials(0x0),
   fh1NGenJets(0x0),
   fh1NRecJets(0x0),
-  fHistList(0x0)
+  fHistList(0x0) ,
+  ////////////////
+  fh1JetMultiplicity(0x0) ,     
+  fh2ERecZRec(0x0) ,
+  fh2EGenZGen(0x0) ,
+  fh2Efficiency(0x0) ,
+  fh3EGenERecN(0x0)
+  //////////////// 
 {
   // Default constructor
   for(int ij  = 0;ij<kMaxJets;++ij){
@@ -111,6 +130,10 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
     fh3PtRecGenHard[ij] =  fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
   }
 
+    for(int ic = 0;ic < kMaxCorrelation;ic++){
+      fhnCorrelation[ic] = 0;
+    }
+
   DefineOutput(1, TList::Class());  
 }
 
@@ -261,15 +284,45 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
                                 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 
   }
+  
+  /////////////////////////////////////////////////////////////////
+  fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
 
-  const Int_t saveLevel = 2; // large save level more histos
+  fh2ERecZRec   = new TH2F("fh2ERecZRec", " ; E^{jet}_{rec} [GeV]; z^{lp}_{rec}", 100, 0., 250., 100, 0., 2.);
+  fh2EGenZGen   = new TH2F("fh2EGenZGen", " ; E^{jet}_{gen} [GeV]; z^{lp}_{gen}", 100, 0., 250., 100, 0., 2.);
+  fh2Efficiency = new TH2F("fh2Efficiency", "ERec/EGen;E^{jet}_{gen} [GeV];E^{jet}_{rec}/E^{jet}_{gen}", 100, 0., 250., 100, 0., 1.5);  
 
+  fh3EGenERecN  = new TH3F("fh3EGenERecN", "Efficiency vs. Jet Multiplicity", 100., 0., 250., 100, 0., 250., 51, 0., 50.);
+
+  // Response map  
+  //arrays for bin limits
+  const Int_t nbin[4] = {100, 100, 100, 100}; 
+  Double_t LowEdge[4] = {0.,0.,0.,0.};
+  Double_t UpEdge[4] = {250., 250., 1., 1.};
+
+  for(int ic = 0;ic < kMaxCorrelation;ic++){
+    fhnCorrelation[ic]  = new THnSparseF(Form("fhnCorrelation_%d",ic),  "Response Map", 4, nbin, LowEdge, UpEdge);
+    if(ic==0) fhnCorrelation[ic]->SetTitle(Form("ResponseMap 0 <= npart <= %.0E",fgkJetNpartCut[ic]));
+    else fhnCorrelation[ic]->SetTitle(Form("ResponseMap %.0E < npart <= %.0E",fgkJetNpartCut[ic-1],fgkJetNpartCut[ic]));
+  }
+  const Int_t saveLevel = 1; // large save level more histos
   if(saveLevel>0){
     fHistList->Add(fh1PtHard);
     fHistList->Add(fh1PtHard_NoW);
     fHistList->Add(fh1PtHard_Trials);
     fHistList->Add(fh1NGenJets);
     fHistList->Add(fh1NRecJets);
+    ////////////////////////
+    fHistList->Add(fh1JetMultiplicity);     
+    fHistList->Add(fh2ERecZRec);
+    fHistList->Add(fh2EGenZGen);
+    fHistList->Add(fh2Efficiency);
+    fHistList->Add(fh3EGenERecN);
+
+    for(int ic = 0;ic < kMaxCorrelation;++ic){
+      fHistList->Add(fhnCorrelation[ic]);
+    }
+    //////////////////////// 
     for(int ij  = 0;ij<kMaxJets;++ij){
       fHistList->Add(fh1E[ij]);
       fHistList->Add(fh1PtRecIn[ij]);
@@ -356,6 +409,9 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
   Int_t nGenJets = 0;
   AliAODJet recJets[kMaxJets];
   Int_t nRecJets = 0;
+  ///////////////////////////
+  Int_t nTracks = 0;
+  //////////////////////////  
 
   Double_t eventW = 1;
   Double_t ptHard = 0; 
@@ -455,6 +511,9 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
   nRecJets = aodRecJets->GetEntries();
   fh1NRecJets->Fill(nRecJets);
   nRecJets = TMath::Min(nRecJets,kMaxJets);
+  //////////////////////////////////////////
+  nTracks  = fAOD->GetNumberOfTracks();
+  ///////////////////////////////////////////  
 
   for(int ir = 0;ir < nRecJets;++ir){
     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
@@ -496,13 +555,55 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
     fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
     // Fill Correlation
     Int_t ig = iGenIndex[ir];
-    if(ig>=0&&ig<nGenJets){
+    if(ig>=0 && ig<nGenJets){
       fh1PtRecOut[ir]->Fill(ptRec,eventW);
-      Double_t ptGen = genJets[ig].Pt();
+      Double_t ptGen  = genJets[ig].Pt();
+      Double_t phiGen = genJets[ig].Phi();
+      if(phiGen<0)phiGen+=TMath::Pi()*2.; 
+      Double_t etaGen = genJets[ig].Eta();
       fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
       fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
       fh3PtRecGenHard_NoW[ir]->Fill(ptRec,ptGen,ptHard,1);
-    }
+      /////////////////////////////////////////////////////
+      Double_t ERec = recJets[ir].E();
+      Double_t EGen = genJets[ig].E();
+      fh2Efficiency->Fill(EGen, ERec/EGen);
+      if (EGen>=0. && EGen<=250.){
+        Double_t Eleading = -1;
+        Double_t ptleading = -1;
+        Int_t nPart=0;
+        // loop over tracks
+        for (Int_t it = 0; it< nTracks; it++){
+          if (fAOD->GetTrack(it)->E() > EGen) continue;
+          Double_t phiTrack = fAOD->GetTrack(it)->Phi();
+          if (phiTrack<0) phiTrack+=TMath::Pi()*2.;            
+          Double_t etaTrack = fAOD->GetTrack(it)->Eta();
+          Float_t deta  = etaRec - etaTrack;
+          Float_t dphi  = TMath::Abs(phiRec - phiTrack);
+          Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);        
+          // find leading particle
+          if (R<0.4 && fAOD->GetTrack(it)->E()>Eleading){
+            Eleading  = fAOD->GetTrack(it)->E();
+            ptleading = fAOD->GetTrack(it)->Pt();            
+          }
+          if (fAOD->GetTrack(it)->Pt()>0.03*EGen && fAOD->GetTrack(it)->E()<=EGen && R<0.7)
+            nPart++;
+        }
+        // fill Response Map (4D histogram) and Energy vs z distributions
+        Double_t var[4] = {EGen, ERec, ptleading/EGen, ptleading/ERec};                       
+        fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
+        fh2EGenZGen->Fill(var[0],var[2]);
+        fh1JetMultiplicity->Fill(nPart);
+        fh3EGenERecN->Fill(EGen, ERec, nPart); 
+       for(int ic = 0;ic <kMaxCorrelation;ic++){
+        if (nPart<=fgkJetNpartCut[ic]){
+         fhnCorrelation[ic]->Fill(var);
+         break;
+       }
+       }
+      }
+    }  
+    ////////////////////////////////////////////////////  
     else{
       fh3RecEtaPhiPt_NoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
     }
index a89345d..47057da 100644 (file)
@@ -5,6 +5,9 @@
  * See cxx source for full Copyright notice                               */
  
 #include "AliAnalysisTaskSE.h"
+////////////////
+#include <THnSparse.h>
+////////////////
 class AliJetHeader;
 class AliESDEvent;
 class AliAODEvent;
@@ -45,14 +48,16 @@ class AliAnalysisTaskJetSpectrum : public AliAnalysisTaskSE
   //
 
     enum {kAnaMC =  0x1};
+    enum {kMaxJets =  5};
+    enum {kMaxCorrelation =  3};
 
  private:
 
     AliAnalysisTaskJetSpectrum(const AliAnalysisTaskJetSpectrum&);
     AliAnalysisTaskJetSpectrum& operator=(const AliAnalysisTaskJetSpectrum&);
     
+    static const Float_t fgkJetNpartCut[kMaxCorrelation];
 
-    enum {kMaxJets =  5};
 
     AliJetHeader *fJetHeaderRec;
     AliJetHeader *fJetHeaderGen;
@@ -102,6 +107,15 @@ class AliAnalysisTaskJetSpectrum : public AliAnalysisTaskSE
     // ============================================                            
   
     TList *fHistList; // Output list
+    
+    ///////// For 2 dimensional unfolding //////////////////
+    TH1F*         fh1JetMultiplicity;   
+    TH2F*         fh2ERecZRec;
+    TH2F*         fh2EGenZGen;
+    TH2F*         fh2Efficiency;
+    TH3F*         fh3EGenERecN;
+    THnSparseF*   fhnCorrelation[kMaxCorrelation];
+    ///////////////////////////////////////////////////////
 
 
     ClassDef(AliAnalysisTaskJetSpectrum, 1) // Analysis task for standard jet analysis
index b1b7fa6..a08c033 100644 (file)
@@ -260,6 +260,8 @@ void AliJetSpectrumUnfolding::SetupCurrentHists(Bool_t createBigBin)
   fCurrentCorrelation = (THnSparseF*)fCorrelation->Clone("fCurrentCorrelation");  
   fCurrentCorrelation->Sumw2();
 
+  Printf("Correlation Matrix has %.0E filled bins", fCurrentCorrelation->GetNbins());
+
   if (createBigBin)
   {
     Int_t maxBinE = 0, maxBinZ = 0;
index 167907e..671c318 100644 (file)
@@ -1,5 +1,5 @@
 //______________________________________________________________________________
-void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
+void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG4/kleinb/LHC08v_jetjet15-50")
 {
 // Example of running analysis train in CAF. To run in debug mode:
 //  - export ROOTSYS=debug  on your local client
@@ -35,12 +35,11 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
     if (iESDfilter) iAODhandler=1;
 
     // Dataset from CAF
-    //    TString dataset = "/PWG4/kleinb/LHC08v_jetjet15-50";
-    TString dataset = "/PWG4/kleinb/LHC08r_jetjet50";
+    TString dataset(ds);
     // CKB quick hack for local analysis
-    gROOT->LoadMacro("CreateESDChain.C");
-    TChain *chain = CreateESDChain("jetjet15-50.txt",1000);
-
+    //    gROOT->LoadMacro("CreateESDChain.C");
+    // TChain *chain = CreateESDChain("jetjet15-50.txt",1000);
+    TChain *chain = 0;
 
     printf("==================================================================\n");
     printf("===========    RUNNING ANALYSIS TRAIN IN CAF MODE    =============\n");
@@ -86,7 +85,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
       TProof::Open(proofNode); 
 
       // Clear packages if changing ROOT version on CAF or local
-      //    gProof->ClearPackages();
+      //      gProof->ClearPackages();
       // Enable proof debugging if needed
       //    gProof->SetLogLevel(5);
       // To debug the train in PROOF mode, type in a root session:
@@ -165,7 +164,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
        AliAODHandler* aodHandler   = new AliAODHandler();
        aodHandler->SetFillAOD(kFALSE);
        mgr->SetOutputEventHandler(aodHandler);       
-       aodHandler->SetOutputFileName(Form("AliAODs_CKB_%07d-%07d.root",nOffset,nOffset+nEvents));
+       aodHandler->SetOutputFileName(Form("AliAODs_pwg4_%07d-%07d.root",nOffset,nOffset+nEvents));
        cout_aod = mgr->CreateContainer("cAOD", TTree::Class(),
                                       AliAnalysisManager::kOutputContainer, "default");
        cout_aod->SetSpecialOutput();
@@ -211,7 +210,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
        mgr->AddTask(jetana);
        // Output histograms list for jet analysis                       
        AliAnalysisDataContainer *cout_jet = mgr->CreateContainer("jethist", TList::Class(),
-                                                             AliAnalysisManager::kOutputContainer,"jethist.root");
+                                                                AliAnalysisManager::kOutputContainer,Form("jethist_%07d-%07d.root",nOffset,nOffset+nEvents));
        // Dummy AOD output container for jet analysis (no client yet)
        c_aodjet = mgr->CreateContainer("cAODjet", TTree::Class(),
                            AliAnalysisManager::kExchangeContainer);
@@ -230,7 +229,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
        mgr->AddTask(jetanaAOD);
        // Output histograms list for jet analysis                       
        AliAnalysisDataContainer *cout_jetAOD = mgr->CreateContainer("jethistAOD", TList::Class(),
-                                                             AliAnalysisManager::kOutputContainer, "jethistAOD.root");
+                                                                   AliAnalysisManager::kOutputContainer, Form("jethistAOD_%07d-%07d.root",nOffset,nOffset+nEvents));
        // Dummy AOD output container for jet analysis (no client yet)
        c_aodjet0 = mgr->CreateContainer("cAODjet0", TTree::Class(),
                            AliAnalysisManager::kExchangeContainer);
@@ -250,7 +249,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
        mgr->AddTask(jetanaMC);
        // Output histograms list for jet analysis                       
        AliAnalysisDataContainer *cout_jetMC = mgr->CreateContainer("jethistMC", TList::Class(),
-                                                             AliAnalysisManager::kOutputContainer, "jethistMC.root");
+                                                                  AliAnalysisManager::kOutputContainer,Form("jethistMC_%07d-%07d.root",nOffset,nOffset+nEvents));
        // Dummy AOD output container for jet analysis (no client yet)
        c_aodjetMC = mgr->CreateContainer("cAODjetMC", TTree::Class(),
                            AliAnalysisManager::kExchangeContainer);
@@ -288,7 +287,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
       //      pwg4spec->SetBranchGen("jetsMC"); 
       mgr->AddTask(pwg4spec);
 
-      AliAnalysisDataContainer *coutput1_Spec = mgr->CreateContainer("pwg4spec", TList::Class(),AliAnalysisManager::kOutputContainer, "histos_pwg4spec.root");
+      AliAnalysisDataContainer *coutput1_Spec = mgr->CreateContainer("pwg4spec", TList::Class(),AliAnalysisManager::kOutputContainer,Form( "pwg4spec_%07d-%07d.root",nOffset,nOffset+nEvents));
       // coutput1_Spec->SetSpecialOutput();
       // Dummy AOD output container for jet analysis (no client yet)
       c_aodSpec = mgr->CreateContainer("cAODjetSpec", TTree::Class(),
@@ -328,7 +327,7 @@ void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
       mgr->AddTask(ueana);
 
 
-      AliAnalysisDataContainer *coutput1_UE = mgr->CreateContainer("histosUE", TList::Class(),AliAnalysisManager::kOutputContainer, "histosUE.root");
+      AliAnalysisDataContainer *coutput1_UE = mgr->CreateContainer("histosUE", TList::Class(),AliAnalysisManager::kOutputContainer, Form("pwg4UE_%07d-%07d.root",nOffset,nOffset+nEvents));
 
       mgr->ConnectInput  (ueana,  0, cinput);    
       //      mgr->ConnectInput  (ueana,  0, c_aodjet);    
index 84d1958..6f57746 100644 (file)
@@ -80,12 +80,12 @@ void unfold(const char* fileNameGen = "gen_pwg4spec.root", const char* folder =
 
   AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
 
-  TFile::Open(fileNameRec);
+  TFile::Open(fileNameGen);
   jetSpec->LoadHistograms();
 
-  TFile::Open(fileNameGen);
-  TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
-  jetSpec->SetGenSpectrum(hist);
+  TFile::Open(fileNameRec);
+  TH2F* hist = (TH2F*) gFile->Get("unfolding/fRecSpectrum");
+  jetSpec->SetRecSpectrum(hist);
 
   jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0);
   // last parameter = calculateErrors  <- this method to calculate the errors takes a lot of time
@@ -301,7 +301,7 @@ void StatisticalUncertainties(const char* fileNameGen = "gen_pwg4spec.root", con
 }
 
 //_______________________________________________________________________________________________________________
-void FillSpecFromFile(const char* fileNameSpec = "histos_pwg4spec.root")
+void FillSpecFromFiles(const char* fileNameReal = "histos_pwg4spec.root",const char* fileNameSim = "histos_pwg4spec.root")
 {
   // This functions avoids problems when the number of bins or the bin limits
   // in the histograms of the AliJetSpectrumUnfolding and AliAnalysisTaskJetSpectrum classes
@@ -312,19 +312,16 @@ void FillSpecFromFile(const char* fileNameSpec = "histos_pwg4spec.root")
   gSystem->Load("libJETAN");
   gSystem->Load("libPWG4JetTasks");
 
-  file = new TFile(fileNameSpec,"read");
+  file = new TFile(fileNameSim,"read");
   tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
-  THnSparseF *fhCorrelation  = (THnSparseF*)(tlist->FindObject("fhCorrelation_less5tracks"));
-  THnSparseF *fhCorrelation2 = (THnSparseF*)(tlist->FindObject("fhCorrelation_5to10tracks"));
-  THnSparseF *fhCorrelation3 = (THnSparseF*)(tlist->FindObject("fhCorrelation_more10tracks"));
-  TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fhEGenZGen"));
-  TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fhERecZRec"));  
 
-  fhCorrelation->Add(fhCorrelation2, 1);
-  fhCorrelation->Add(fhCorrelation3, 1);
-
-  delete fhCorrelation2;
-  delete fhCorrelation3;
+  THnSparseF *fhCorrelation  = 0;
+  for(int ic = 0;ic<3;++ic){
+    THnSparseF *hTmp = (THnSparseF*)(tlist->FindObject(Form("fhnCorrelation_%d",ic)));    
+    if(fhCorrelation==0)fhCorrelation =  (THnSparseF*)hTmp->Clone("fhnCorrelation");
+    else fhCorrelation->Add(hTmp);
+  }
+  TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen"));
 
   AliJetSpectrumUnfolding *jetSpec = new AliJetSpectrumUnfolding("unfolding","unfolding");
 
@@ -341,30 +338,9 @@ void FillSpecFromFile(const char* fileNameSpec = "histos_pwg4spec.root")
        jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
        jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
     }
-  file = TFile::Open("gen_pwg4spec.root", "RECREATE");
-  jetSpec->SaveHistograms();
-  file->Close();
-  jetSpec->GetGenSpectrum()->Reset();
-  printf("true distribution has been set\n");
 
-  // reconstructed jets (measured distribution)
-  for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++)
-    for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++)
-    {
-       Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me);
-       Float_t   zm = fhERecZRec->GetYaxis()->GetBinCenter(mz);
-       Int_t bine   = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec);
-       Int_t binz   = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm);
-       Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz);
-       Float_t err  = jetSpec->GetRecSpectrum()->GetBinError(bine, binz);
-       jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz));
-       jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
-    }
 
-  // Response function
-  jetSpec->GetCorrelation()->Reset();
-  jetSpec->GetCorrelation()->Sumw2();
-        
+  Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins());
   for (Int_t idx=1; idx<=fhCorrelation->GetNbins(); idx++)
   {
     //printf("idx: %d\n",idx);
@@ -384,6 +360,45 @@ void FillSpecFromFile(const char* fileNameSpec = "histos_pwg4spec.root")
     jetSpec->GetCorrelation()->SetBinContent(bin, cont + fhCorrelation->GetBinContent(idx));
     jetSpec->GetCorrelation()->SetBinError(bin, err + fhCorrelation->GetBinError(idx));
   }
+  // need number of entries for correct drawing
+  jetSpec->GetCorrelation()->SetEntries(fhCorrelation->GetEntries());
+
+
+  file = TFile::Open("gen_pwg4spec.root", "RECREATE");
+  jetSpec->SaveHistograms();
+  Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins());
+  file->Close();
+  jetSpec->GetGenSpectrum()->Reset();
+  printf("true distribution has been set\n");
+
+  // reconstructed jets (measured distribution)
+
+
+  // Response function
+  jetSpec->GetCorrelation()->Reset();
+  jetSpec->GetCorrelation()->Sumw2();
+  jetSpec->GetGenSpectrum()->Reset();
+  jetSpec->GetRecSpectrum()->Reset();
+
+
+  file = new TFile(fileNameReal,"read");
+  tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
+
+  
+  TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fh2ERecZRec"));  
+
+  for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++)
+    for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++)
+    {
+       Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me);
+       Float_t   zm = fhERecZRec->GetYaxis()->GetBinCenter(mz);
+       Int_t bine   = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec);
+       Int_t binz   = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm);
+       Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz);
+       Float_t err  = jetSpec->GetRecSpectrum()->GetBinError(bine, binz);
+       jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz));
+       jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
+    }
 
   file = TFile::Open("rec_pwg4spec.root", "RECREATE");
   jetSpec->SaveHistograms();
@@ -400,7 +415,7 @@ void FillSpecFromFile(const char* fileNameSpec = "histos_pwg4spec.root")
 void correct(){
   // simple steering to correct a given distribution;
   loadlibs();
-  FillSpecFromFile("/home/ydelgado/pcalice014.cern.ch/macros/jets/CAFdata/histos_pwg4spec.root");
+  FillSpecFromFiles("pwg4spec_0000000-0010000.root","pwg4spec_15-50.root");
 
   char name[100];
   sprintf(name, "unfolded_pwg4spec.root");