Added PID task, some fixes for coding conventions
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum.cxx
index 75c432c8be1d1b8332279f91bfdef5bebf78b7df..5a2947d4b7a96762566584bca7b27d48e92ebbb6 100644 (file)
@@ -1,3 +1,9 @@
+// **************************************
+// Task used for the correction of determiantion of reconstructed jet spectra
+// Compares input (gen) and output (rec) jets   
+// *******************************************
+
+
 /**************************************************************************
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  *                                                                        *
@@ -13,7 +19,6 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-
  
 #include <TROOT.h>
 #include <TSystem.h>
@@ -23,6 +28,7 @@
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TH3F.h>
+#include <TProfile.h>
 #include <TList.h>
 #include <TLorentzVector.h>
 #include <TClonesArray.h>
@@ -53,6 +59,8 @@
 
 ClassImp(AliAnalysisTaskJetSpectrum)
 
+const Float_t AliAnalysisTaskJetSpectrum::fgkJetNpartCut[AliAnalysisTaskJetSpectrum::kMaxCorrelation] = {5,10,1E+09};
+
 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
   fJetHeaderRec(0x0),
   fJetHeaderGen(0x0),
@@ -66,19 +74,32 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
   fLimitGenJetEta(kFALSE),
   fAnalysisType(0),
   fExternalWeight(1),    
+  fRecEtaWindow(0.5),
+  fh1Xsec(0x0),
+  fh1Trials(0x0),
   fh1PtHard(0x0),
-  fh1PtHard_NoW(0x0),  
-  fh1PtHard_Trials(0x0),
+  fh1PtHardNoW(0x0),  
+  fh1PtHardTrials(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] = fh2PhiFGen[ij] = fh2EtaFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] =   fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[ij] =  fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] = 0;
+      fh3PtRecGenHard[ij] =  fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
+    }
+    for(int ic = 0;ic < kMaxCorrelation;ic++){
+      fhnCorrelation[ic] = 0;
+    }
   
 }
 
@@ -96,26 +117,90 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
   fLimitGenJetEta(kFALSE),
   fAnalysisType(0),
   fExternalWeight(1),    
+  fRecEtaWindow(0.5),
+  fh1Xsec(0x0),
+  fh1Trials(0x0),
   fh1PtHard(0x0),
-  fh1PtHard_NoW(0x0),  
-  fh1PtHard_Trials(0x0),
+  fh1PtHardNoW(0x0),  
+  fh1PtHardTrials(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;
+    fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] =   fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[ij] =0;
 
-    fh3PtRecGenHard[ij] =  fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
+    fh3PtRecGenHard[ij] =  fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
   }
 
+    for(int ic = 0;ic < kMaxCorrelation;ic++){
+      fhnCorrelation[ic] = 0;
+    }
+
   DefineOutput(1, TList::Class());  
 }
 
 
 
+Bool_t AliAnalysisTaskJetSpectrum::Notify()
+{
+  //
+  // Implemented Notify() to read the cross sections
+  // and number of trials from pyxsec.root
+  // 
+  TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
+  Double_t xsection = 0;
+  UInt_t   ntrials  = 0;
+  if(tree){
+    TFile *curfile = tree->GetCurrentFile();
+    if (!curfile) {
+      Error("Notify","No current file");
+      return kFALSE;
+    }
+    if(!fh1Xsec||!fh1Trials){
+      Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
+      return kFALSE;
+    }
+
+    TString fileName(curfile->GetName());
+    if(fileName.Contains("AliESDs.root")){
+        fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
+    }
+    else if(fileName.Contains("AliAOD.root")){
+        fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
+    }
+    else if(fileName.Contains("galice.root")){
+        // for running with galice and kinematics alone...                      
+        fileName.ReplaceAll("galice.root", "pyxsec.root");
+    }
+    TFile *fxsec = TFile::Open(fileName.Data());
+    if(!fxsec){
+      Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
+      // no a severe condition
+      return kTRUE;
+    }
+    TTree *xtree = (TTree*)fxsec->Get("Xsection");
+    if(!xtree){
+      Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
+    }
+    xtree->SetBranchAddress("xsection",&xsection);
+    xtree->SetBranchAddress("ntrials",&ntrials);
+    xtree->GetEntry(0);
+    fh1Xsec->Fill("<#sigma>",xsection);
+    fh1Trials->Fill("#sum{ntrials}",ntrials);
+  }
+  
+  return kTRUE;
+}
 
 void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
 {
@@ -175,21 +260,20 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
       binLimitsPt[iPt] = 0.0;
     }
     else {// 1.0
-      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2;
-    }
-  }
-  const Int_t nBinEta = 22;
-  Double_t binLimitsEta[nBinEta+1];
-  for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
-    if(iEta==0){
-      binLimitsEta[iEta] = -2.2;
-    }
-    else{
-      binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.2;
+      binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
     }
   }
+  
+  const Int_t nBinEta = 26;
+  Double_t binLimitsEta[nBinEta+1] = {
+    -1.6,-1.4,-1.2,-1.0,
+    -0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,
+    0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 
+    1.0, 1.2, 1.4, 1.6
+  };
 
-  const Int_t nBinPhi = 90;
+
+  const Int_t nBinPhi = 30;
   Double_t binLimitsPhi[nBinPhi+1];
   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
     if(iPhi==0){
@@ -203,11 +287,17 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
   const Int_t nBinFrag = 25;
 
 
+  fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
+  fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
+
+  fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
+  fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
+
   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
 
-  fh1PtHard_NoW = new TH1F("fh1PtHard_NoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
+  fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
 
-  fh1PtHard_Trials = new TH1F("fh1PtHard_Trials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
+  fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
 
   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
 
@@ -222,11 +312,27 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
     fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
 
 
+
     fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
                             nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 
+    fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};#phi_{gen}",
+                            nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi);
 
+    fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};#eta_{gen}",
+                            nBinEta,binLimitsEta,nBinEta,binLimitsEta);
 
+    
+    fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaPhi_j%d",ij),"delta phi vs. P_{T,gen};p_{T,gen} (GeV/c);#phi_{gen}-#phi_{rec}",
+                                   nBinPt,binLimitsPt,100,-1.0,1.0);
+    fh2PtGenDeltaEta[ij] = new TH2F(Form("fh2PtGenDeltaEta_j%d",ij),"delta eta vs. p_{T,gen};p_{T,gen} (GeV/c);#eta_{gen}-#eta_{rec}",
+                                   nBinPt,binLimitsPt,100,-1.0,1.0);
+
+    
+    fh2PtRecDeltaR[ij] = new TH2F(Form("fh2PtRecDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,rec,j};#Delta R", 
+                                 nBinPt,binLimitsPt,60,0,6.0);
+    fh2PtGenDeltaR[ij] = new TH2F(Form("fh2PtGenDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,gen,j};#Delta R", 
+                                 nBinPt,binLimitsPt,60,0,6.0);
 
 
 
@@ -234,7 +340,7 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
 
 
 
-    fh3PtRecGenHard_NoW[ij] = new TH3F(Form("fh3PtRecGenHard_NoW_j%d",ij), "Pt hard vs. pt gen vs. pt rec no weight;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
+    fh3PtRecGenHardNoW[ij] = new TH3F(Form("fh3PtRecGenHardNoW_j%d",ij), "Pt hard vs. pt gen vs. pt rec no weight;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
 
        
     fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
@@ -248,28 +354,60 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
 
 
 
-    fh3RecEtaPhiPt_NoGen[ij] = new TH3F(Form("fh3RecEtaPhiPt_NoGen_j%d",ij),"No generated for found jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
+    fh3RecEtaPhiPtNoGen[ij] = new TH3F(Form("fh3RecEtaPhiPtNoGen_j%d",ij),"No generated for found jet Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
                                        nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 
 
-    fh3GenEtaPhiPt_NoFound[ij] = new TH3F(Form("fh3GenEtaPhiPt_NoFound_j%d",ij),"No found for generated jet eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
+    fh3GenEtaPhiPtNoFound[ij] = new TH3F(Form("fh3GenEtaPhiPtNoFound_j%d",ij),"No found for generated jet eta, phi, pt; #eta; #phi; p_{T,gen} (GeV/c)",
                                        nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 
 
 
-    fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
+    fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,} (GeV/c)",
                                 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
 
   }
+  
+  /////////////////////////////////////////////////////////////////
+  fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
+
+  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.);
 
-  const Int_t saveLevel = 2; // large save level more histos
+  // Response map  
+  //arrays for bin limits
+  const Int_t nbin[4] = {100, 100, 100, 100}; 
+  Double_t vLowEdge[4] = {0.,0.,0.,0.};
+  Double_t vUpEdge[4] = {250., 250., 1., 1.};
 
+  for(int ic = 0;ic < kMaxCorrelation;ic++){
+    fhnCorrelation[ic]  = new THnSparseF(Form("fhnCorrelation_%d",ic),  "Response Map", 4, nbin, vLowEdge, vUpEdge);
+    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 = 3; // large save level more histos
   if(saveLevel>0){
+    fHistList->Add(fh1Xsec);
+    fHistList->Add(fh1Trials);
     fHistList->Add(fh1PtHard);
-    fHistList->Add(fh1PtHard_NoW);
-    fHistList->Add(fh1PtHard_Trials);
+    fHistList->Add(fh1PtHardNoW);
+    fHistList->Add(fh1PtHardTrials);
     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]);
@@ -277,11 +415,17 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
       fHistList->Add(fh1PtGenIn[ij]);
       fHistList->Add(fh1PtGenOut[ij]);
       fHistList->Add(fh2PtFGen[ij]);
+      fHistList->Add(fh2PhiFGen[ij]);
+      fHistList->Add(fh2EtaFGen[ij]);
+      fHistList->Add(fh2PtGenDeltaEta[ij]);
+      fHistList->Add(fh2PtGenDeltaPhi[ij]);
+      fHistList->Add(fh2PtRecDeltaR[ij]);
+      fHistList->Add(fh2PtGenDeltaR[ij]);
+      fHistList->Add(fh3RecEtaPhiPt[ij]);
+      fHistList->Add(fh3GenEtaPhiPt[ij]);      
       if(saveLevel>2){
-       fHistList->Add(fh3RecEtaPhiPt[ij]);
-       fHistList->Add(fh3RecEtaPhiPt_NoGen[ij]);
-       fHistList->Add(fh3GenEtaPhiPt_NoFound[ij]);
-       fHistList->Add(fh3GenEtaPhiPt[ij]);      
+       fHistList->Add(fh3RecEtaPhiPtNoGen[ij]);
+       fHistList->Add(fh3GenEtaPhiPtNoFound[ij]);
       }
     }
   }
@@ -292,7 +436,10 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
     if (h1){
       // Printf("%s ",h1->GetName());
       h1->Sumw2();
+      continue;
     }
+    THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
+    if(hn)hn->Sumw2();
   }
 
   TH1::AddDirectory(oldStatus);
@@ -356,6 +503,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; 
@@ -382,7 +532,23 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
 
     nTrials = pythiaGenHeader->Trials();
     ptHard  = pythiaGenHeader->GetPtHard();
-
+    int iProcessType = pythiaGenHeader->ProcessType();
+    // 11 f+f -> f+f
+    // 12 f+barf -> f+barf
+    // 13 f+barf -> g+g
+    // 28 f+g -> f+g
+    // 53 g+g -> f+barf
+    // 68 g+g -> g+g
+    /*
+    if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
+    //    if(iProcessType != 13 && iProcessType != 68){ // allow only glue
+    if(iProcessType != 11 && iProcessType != 12 && iProcessType != 53){ // allow only quark
+    //    if(iProcessType != 28){ // allow only -> f+g
+      PostData(1, fHistList);
+      return;
+    }
+    */
+    if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
 
     if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
 
@@ -419,8 +585,8 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
 
   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
   fh1PtHard->Fill(ptHard,eventW);
-  fh1PtHard_NoW->Fill(ptHard,1);
-  fh1PtHard_Trials->Fill(ptHard,nTrials);
+  fh1PtHardNoW->Fill(ptHard,1);
+  fh1PtHardTrials->Fill(ptHard,nTrials);
 
   // If we set a second branch for the input jets fetch this 
   if(fBranchGen.Length()>0){
@@ -455,6 +621,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));
@@ -490,21 +659,90 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
     Double_t phiRec = recJets[ir].Phi();
     if(phiRec<0)phiRec+=TMath::Pi()*2.;    
     Double_t etaRec = recJets[ir].Eta();
-
+    if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
     fh1E[ir]->Fill(recJets[ir].E(),eventW);
     fh1PtRecIn[ir]->Fill(ptRec,eventW);
     fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
+    for(int irr = ir+1;irr<nRecJets;irr++){
+      fh2PtRecDeltaR[ir]->Fill(recJets[irr].Pt(),recJets[ir].DeltaR(&recJets[irr]));
+    }
     // Fill Correlation
     Int_t ig = iGenIndex[ir];
-    if(ig>=0&&ig<nGenJets){
+    if(ig>=0 && ig<nGenJets){
+      if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+      if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
       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();
+
+      // 
+      // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
+      // 
+
+      if(TMath::Abs(etaRec)<fRecEtaWindow){
+
       fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
+      fh2PhiFGen[ir]->Fill(phiRec,phiGen,eventW);
+      fh2EtaFGen[ir]->Fill(etaRec,etaGen,eventW);
+      fh2PtGenDeltaEta[ir]->Fill(ptGen,etaGen-etaRec,eventW);
+      fh2PtGenDeltaPhi[ir]->Fill(ptGen,phiGen-phiRec,eventW);
       fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
-      fh3PtRecGenHard_NoW[ir]->Fill(ptRec,ptGen,ptHard,1);
-    }
+      fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1);
+      /////////////////////////////////////////////////////
+
+      //      Double_t eRec = recJets[ir].E();
+      //      Double_t eGen = genJets[ig].E();
+      // CKB use p_T not Energy 
+      // TODO recname variabeles and histos
+      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; // CKB. Not allowed! cannot cut on gen properties in real events!
+          // find leading particle
+         //  if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){
+         // TODO implement esd filter flag to be the same as in the jet finder 
+         // allow also for MC particles...
+         Float_t r = recJets[ir].DeltaR(fAOD->GetTrack(it));
+         if (r<0.4 && fAOD->GetTrack(it)->Pt()>ptleading){
+            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) // CKB cannot cut on gen properties 
+          if (fAOD->GetTrack(it)->Pt()>0.03*eRec && fAOD->GetTrack(it)->Pt()<=eRec && r<0.7)
+            nPart++;
+        }
+       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
+
+        // 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]){ // is this corrected for CKB
+           fhnCorrelation[ic]->Fill(var);
+           break;
+         }
+       }
+      }
+
+    }// if etarec in window
+
+    } 
+    ////////////////////////////////////////////////////  
     else{
-      fh3RecEtaPhiPt_NoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
+      fh3RecEtaPhiPtNoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
     }
   }// loop over reconstructed jets
 
@@ -518,12 +756,15 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
     Double_t etaGen = genJets[ig].Eta();
     fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
     fh1PtGenIn[ig]->Fill(ptGen,eventW);
+    for(int igg = ig+1;igg<nGenJets;igg++){
+      fh2PtGenDeltaR[ig]->Fill(genJets[igg].Pt(),genJets[ig].DeltaR(&genJets[igg]));
+    }
     Int_t ir = iRecIndex[ig];
     if(ir>=0&&ir<nRecJets){
       fh1PtGenOut[ig]->Fill(ptGen,eventW);
     }
     else{
-      fh3GenEtaPhiPt_NoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
+      fh3GenEtaPhiPtNoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
     }
   }// loop over reconstructed jets
 
@@ -539,8 +780,8 @@ void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
 }
 
 
-void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,Int_t &nGenJets,
-                                               AliAODJet *recJets,Int_t &nRecJets,
+void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,const Int_t &nGenJets,
+                                               AliAODJet *recJets,const Int_t &nRecJets,
                                                Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
 
   //
@@ -587,7 +828,7 @@ void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,Int_t &nGenJe
   if(nRecJets==0)return;
   if(nGenJets==0)return;
 
-  const Float_t maxDist = 1.0;
+  const Float_t maxDist = 0.5;
   // find the closest distance to the generated
   for(int ig = 0;ig<nGenJets;++ig){
     Float_t dist = maxDist;