Added PID task, some fixes for coding conventions
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum.cxx
index 26948e55ae21b654e755a80a1d326a365c7c5fe9..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>
@@ -68,9 +74,12 @@ 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) , 
@@ -85,8 +94,8 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
   // 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;
+      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;
@@ -108,9 +117,12 @@ 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) ,
@@ -125,9 +137,9 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
   // 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++){
@@ -139,6 +151,56 @@ AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
 
 
 
+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()
 {
@@ -198,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){
@@ -226,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);
 
@@ -245,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);
 
 
 
@@ -257,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",
@@ -271,16 +354,16 @@ 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);
 
   }
@@ -292,24 +375,26 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
   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.);
+  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.};
+  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, LowEdge, UpEdge);
+    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 = 1; // large save level more histos
+  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);
     ////////////////////////
@@ -330,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]);
       }
     }
   }
@@ -345,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);
@@ -438,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);
 
@@ -475,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){
@@ -549,63 +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 (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 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();
-      fh2Efficiency->Fill(EGen, ERec/EGen);
-      if (EGen>=0. && EGen<=250.){
-        Double_t Eleading = -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;
-          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);        
+         //      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){
-            Eleading  = fAOD->GetTrack(it)->E();
+         //  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)
+         //          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};                       
+        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); 
+        fh3EGenERecN->Fill(eGen, eRec, nPart); 
        for(int ic = 0;ic <kMaxCorrelation;ic++){
-        if (nPart<=fgkJetNpartCut[ic]){
-         fhnCorrelation[ic]->Fill(var);
-         break;
-       }
+         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
 
@@ -619,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
 
@@ -640,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){
 
   //
@@ -688,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;