Added PID task, some fixes for coding conventions
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum.cxx
1 // **************************************
2 // Task used for the correction of determiantion of reconstructed jet spectra
3 // Compares input (gen) and output (rec) jets   
4 // *******************************************
5
6
7 /**************************************************************************
8  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9  *                                                                        *
10  * Author: The ALICE Off-line Project.                                    *
11  * Contributors are mentioned in the code where appropriate.              *
12  *                                                                        *
13  * Permission to use, copy, modify and distribute this software and its   *
14  * documentation strictly for non-commercial purposes is hereby granted   *
15  * without fee, provided that the above copyright notice appears in all   *
16  * copies and that both the copyright notice and this permission notice   *
17  * appear in the supporting documentation. The authors make no claims     *
18  * about the suitability of this software for any purpose. It is          *
19  * provided "as is" without express or implied warranty.                  *
20  **************************************************************************/
21
22  
23 #include <TROOT.h>
24 #include <TSystem.h>
25 #include <TInterpreter.h>
26 #include <TChain.h>
27 #include <TFile.h>
28 #include <TH1F.h>
29 #include <TH2F.h>
30 #include <TH3F.h>
31 #include <TProfile.h>
32 #include <TList.h>
33 #include <TLorentzVector.h>
34 #include <TClonesArray.h>
35 #include  "TDatabasePDG.h"
36
37 #include "AliAnalysisTaskJetSpectrum.h"
38 #include "AliAnalysisManager.h"
39 #include "AliJetFinder.h"
40 #include "AliJetHeader.h"
41 #include "AliJetReader.h"
42 #include "AliJetReaderHeader.h"
43 #include "AliUA1JetHeaderV1.h"
44 #include "AliJet.h"
45 #include "AliESDEvent.h"
46 #include "AliAODEvent.h"
47 #include "AliAODHandler.h"
48 #include "AliAODTrack.h"
49 #include "AliAODJet.h"
50 #include "AliMCEventHandler.h"
51 #include "AliMCEvent.h"
52 #include "AliStack.h"
53 #include "AliGenPythiaEventHeader.h"
54 #include "AliJetKineReaderHeader.h"
55 #include "AliGenCocktailEventHeader.h"
56 #include "AliInputEventHandler.h"
57
58 #include "AliAnalysisHelperJetTasks.h"
59
60 ClassImp(AliAnalysisTaskJetSpectrum)
61
62 const Float_t AliAnalysisTaskJetSpectrum::fgkJetNpartCut[AliAnalysisTaskJetSpectrum::kMaxCorrelation] = {5,10,1E+09};
63
64 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
65   fJetHeaderRec(0x0),
66   fJetHeaderGen(0x0),
67   fAOD(0x0),
68   fBranchRec("jets"),
69   fConfigRec("ConfigJets.C"),
70   fBranchGen(""),
71   fConfigGen(""),
72   fUseAODInput(kFALSE),
73   fUseExternalWeightOnly(kFALSE),
74   fLimitGenJetEta(kFALSE),
75   fAnalysisType(0),
76   fExternalWeight(1),    
77   fRecEtaWindow(0.5),
78   fh1Xsec(0x0),
79   fh1Trials(0x0),
80   fh1PtHard(0x0),
81   fh1PtHardNoW(0x0),  
82   fh1PtHardTrials(0x0),
83   fh1NGenJets(0x0),
84   fh1NRecJets(0x0),
85   fHistList(0x0) , 
86   ////////////////
87   fh1JetMultiplicity(0x0) ,     
88   fh2ERecZRec(0x0) ,
89   fh2EGenZGen(0x0) ,
90   fh2Efficiency(0x0) ,
91   fh3EGenERecN(0x0) 
92   //////////////// 
93 {
94   // Default constructor
95     for(int ij  = 0;ij<kMaxJets;++ij){
96       fh1E[ij] =  fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
97       fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] =   fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[ij] =  fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] = 0;
98       fh3PtRecGenHard[ij] =  fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
99     }
100     for(int ic = 0;ic < kMaxCorrelation;ic++){
101       fhnCorrelation[ic] = 0;
102     }
103   
104 }
105
106 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
107   AliAnalysisTaskSE(name),
108   fJetHeaderRec(0x0),
109   fJetHeaderGen(0x0),
110   fAOD(0x0),
111   fBranchRec("jets"),
112   fConfigRec("ConfigJets.C"),
113   fBranchGen(""),
114   fConfigGen(""),
115   fUseAODInput(kFALSE),
116   fUseExternalWeightOnly(kFALSE),
117   fLimitGenJetEta(kFALSE),
118   fAnalysisType(0),
119   fExternalWeight(1),    
120   fRecEtaWindow(0.5),
121   fh1Xsec(0x0),
122   fh1Trials(0x0),
123   fh1PtHard(0x0),
124   fh1PtHardNoW(0x0),  
125   fh1PtHardTrials(0x0),
126   fh1NGenJets(0x0),
127   fh1NRecJets(0x0),
128   fHistList(0x0) ,
129   ////////////////
130   fh1JetMultiplicity(0x0) ,     
131   fh2ERecZRec(0x0) ,
132   fh2EGenZGen(0x0) ,
133   fh2Efficiency(0x0) ,
134   fh3EGenERecN(0x0)
135   //////////////// 
136 {
137   // Default constructor
138   for(int ij  = 0;ij<kMaxJets;++ij){
139     fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
140     fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtGenDeltaPhi[ij] =  fh2PtGenDeltaEta[ij] =   fh2PtRecDeltaR[ij] = fh2PtGenDeltaR[ij] =0;
141
142     fh3PtRecGenHard[ij] =  fh3PtRecGenHardNoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPtNoGen[ij] =fh3GenEtaPhiPtNoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
143   }
144
145     for(int ic = 0;ic < kMaxCorrelation;ic++){
146       fhnCorrelation[ic] = 0;
147     }
148
149   DefineOutput(1, TList::Class());  
150 }
151
152
153
154 Bool_t AliAnalysisTaskJetSpectrum::Notify()
155 {
156   //
157   // Implemented Notify() to read the cross sections
158   // and number of trials from pyxsec.root
159   // 
160   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
161   Double_t xsection = 0;
162   UInt_t   ntrials  = 0;
163   if(tree){
164     TFile *curfile = tree->GetCurrentFile();
165     if (!curfile) {
166       Error("Notify","No current file");
167       return kFALSE;
168     }
169     if(!fh1Xsec||!fh1Trials){
170       Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
171       return kFALSE;
172     }
173
174     TString fileName(curfile->GetName());
175     if(fileName.Contains("AliESDs.root")){
176         fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
177     }
178     else if(fileName.Contains("AliAOD.root")){
179         fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
180     }
181     else if(fileName.Contains("galice.root")){
182         // for running with galice and kinematics alone...                      
183         fileName.ReplaceAll("galice.root", "pyxsec.root");
184     }
185     TFile *fxsec = TFile::Open(fileName.Data());
186     if(!fxsec){
187       Printf("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data());
188       // no a severe condition
189       return kTRUE;
190     }
191     TTree *xtree = (TTree*)fxsec->Get("Xsection");
192     if(!xtree){
193       Printf("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__);
194     }
195     xtree->SetBranchAddress("xsection",&xsection);
196     xtree->SetBranchAddress("ntrials",&ntrials);
197     xtree->GetEntry(0);
198     fh1Xsec->Fill("<#sigma>",xsection);
199     fh1Trials->Fill("#sum{ntrials}",ntrials);
200   }
201   
202   return kTRUE;
203 }
204
205 void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
206 {
207
208   //
209   // Create the output container
210   //
211
212   
213   // Connect the AOD
214
215   if(fUseAODInput){
216     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
217     if(!fAOD){
218       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
219       return;
220     }
221     // fethc the header
222     fJetHeaderRec = dynamic_cast<AliJetHeader*>(fInputHandler->GetTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
223     if(!fJetHeaderRec){
224       Printf("%s:%d Jet Header not found in the Input",(char*)__FILE__,__LINE__);
225     }
226   }
227   else{
228     //  assume that the AOD is in the general output...
229     fAOD  = AODEvent();
230     if(!fAOD){
231       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
232       return;
233     }
234     fJetHeaderRec = dynamic_cast<AliJetHeader*>(OutputTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));    
235     if(!fJetHeaderRec){
236       Printf("%s:%d Jet Header not found in the Output",(char*)__FILE__,__LINE__);
237     }
238     else{
239       if(fDebug>10)fJetHeaderRec->Dump();
240     }
241   }
242  
243
244
245   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
246
247   OpenFile(1);
248   if(!fHistList)fHistList = new TList();
249
250   Bool_t oldStatus = TH1::AddDirectoryStatus();
251   TH1::AddDirectory(kFALSE);
252
253   //
254   //  Histogram
255     
256   const Int_t nBinPt = 100;
257   Double_t binLimitsPt[nBinPt+1];
258   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
259     if(iPt == 0){
260       binLimitsPt[iPt] = 0.0;
261     }
262     else {// 1.0
263       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2.5;
264     }
265   }
266   
267   const Int_t nBinEta = 26;
268   Double_t binLimitsEta[nBinEta+1] = {
269     -1.6,-1.4,-1.2,-1.0,
270     -0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,0.0,
271     0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 
272     1.0, 1.2, 1.4, 1.6
273   };
274
275
276   const Int_t nBinPhi = 30;
277   Double_t binLimitsPhi[nBinPhi+1];
278   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
279     if(iPhi==0){
280       binLimitsPhi[iPhi] = 0;
281     }
282     else{
283       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
284     }
285   }
286
287   const Int_t nBinFrag = 25;
288
289
290   fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
291   fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
292
293   fh1Trials = new TH1F("fh1Trials","trials from pyxsec.root",1,0,1);
294   fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
295
296   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
297
298   fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
299
300   fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
301
302   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
303
304   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
305
306
307   for(int ij  = 0;ij<kMaxJets;++ij){
308     fh1E[ij] = new TH1F(Form("fh1E_j%d",ij),"Jet Energy;E_{jet} (GeV);N",nBinPt,binLimitsPt);
309     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
310     fh1PtRecOut[ij] = new TH1F(Form("fh1PtRecOut_j%d",ij),"rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
311     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
312     fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
313
314
315
316     fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
317                              nBinPt,binLimitsPt,nBinPt,binLimitsPt);
318
319     fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};#phi_{gen}",
320                              nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi);
321
322     fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};#eta_{gen}",
323                              nBinEta,binLimitsEta,nBinEta,binLimitsEta);
324
325     
326     fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaPhi_j%d",ij),"delta phi vs. P_{T,gen};p_{T,gen} (GeV/c);#phi_{gen}-#phi_{rec}",
327                                     nBinPt,binLimitsPt,100,-1.0,1.0);
328     fh2PtGenDeltaEta[ij] = new TH2F(Form("fh2PtGenDeltaEta_j%d",ij),"delta eta vs. p_{T,gen};p_{T,gen} (GeV/c);#eta_{gen}-#eta_{rec}",
329                                     nBinPt,binLimitsPt,100,-1.0,1.0);
330
331     
332     fh2PtRecDeltaR[ij] = new TH2F(Form("fh2PtRecDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,rec,j};#Delta R", 
333                                   nBinPt,binLimitsPt,60,0,6.0);
334     fh2PtGenDeltaR[ij] = new TH2F(Form("fh2PtGenDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,gen,j};#Delta R", 
335                                   nBinPt,binLimitsPt,60,0,6.0);
336
337
338
339     fh3PtRecGenHard[ij] = new TH3F(Form("fh3PtRecGenHard_j%d",ij), "Pt hard vs. pt gen vs. pt rec;p_{T,rec};p_{T,gen} (GeV/c);p_{T,hard} (GeV/c)",nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
340
341
342
343     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);
344
345         
346     fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
347                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
348
349     fh2FragLn[ij] = new TH2F(Form("fh2FragLn_j%d",ij),"Jet Fragmentation Ln;#xi=ln(E_{jet}/E_{i});E_{jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
350                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
351
352     fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
353                                   nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
354
355
356
357     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)",
358                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
359
360
361     fh3GenEtaPhiPtNoFound[ij] = new TH3F(Form("fh3GenEtaPhiPtNoFound_j%d",ij),"No found for generated jet eta, phi, pt; #eta; #phi; p_{T,gen} (GeV/c)",
362                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
363
364
365
366     fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,} (GeV/c)",
367                                  nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
368
369   }
370   
371   /////////////////////////////////////////////////////////////////
372   fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
373
374   fh2ERecZRec   = new TH2F("fh2ERecZRec", " ; E^{jet}_{rec} [GeV]; z^{lp}_{rec}", 100, 0., 250., 100, 0., 2.);
375   fh2EGenZGen   = new TH2F("fh2EGenZGen", " ; E^{jet}_{gen} [GeV]; z^{lp}_{gen}", 100, 0., 250., 100, 0., 2.);
376   fh2Efficiency = new TH2F("fh2Efficiency", "ERec/EGen;E^{jet}_{gen} [GeV];E^{jet}_{rec}/E^{jet}_{gen}", 100, 0., 250., 100, 0., 1.5);  
377
378   fh3EGenERecN  = new TH3F("fh3EGenERecN", "Efficiency vs. Jet Multiplicity", 100, 0., 250., 100, 0., 250., 51, 0., 50.);
379
380   // Response map  
381   //arrays for bin limits
382   const Int_t nbin[4] = {100, 100, 100, 100}; 
383   Double_t vLowEdge[4] = {0.,0.,0.,0.};
384   Double_t vUpEdge[4] = {250., 250., 1., 1.};
385
386   for(int ic = 0;ic < kMaxCorrelation;ic++){
387     fhnCorrelation[ic]  = new THnSparseF(Form("fhnCorrelation_%d",ic),  "Response Map", 4, nbin, vLowEdge, vUpEdge);
388     if(ic==0) fhnCorrelation[ic]->SetTitle(Form("ResponseMap 0 <= npart <= %.0E",fgkJetNpartCut[ic]));
389     else fhnCorrelation[ic]->SetTitle(Form("ResponseMap %.0E < npart <= %.0E",fgkJetNpartCut[ic-1],fgkJetNpartCut[ic]));
390   }
391   const Int_t saveLevel = 3; // large save level more histos
392   if(saveLevel>0){
393     fHistList->Add(fh1Xsec);
394     fHistList->Add(fh1Trials);
395     fHistList->Add(fh1PtHard);
396     fHistList->Add(fh1PtHardNoW);
397     fHistList->Add(fh1PtHardTrials);
398     fHistList->Add(fh1NGenJets);
399     fHistList->Add(fh1NRecJets);
400     ////////////////////////
401     fHistList->Add(fh1JetMultiplicity);     
402     fHistList->Add(fh2ERecZRec);
403     fHistList->Add(fh2EGenZGen);
404     fHistList->Add(fh2Efficiency);
405     fHistList->Add(fh3EGenERecN);
406
407     for(int ic = 0;ic < kMaxCorrelation;++ic){
408       fHistList->Add(fhnCorrelation[ic]);
409     }
410     //////////////////////// 
411     for(int ij  = 0;ij<kMaxJets;++ij){
412       fHistList->Add(fh1E[ij]);
413       fHistList->Add(fh1PtRecIn[ij]);
414       fHistList->Add(fh1PtRecOut[ij]);
415       fHistList->Add(fh1PtGenIn[ij]);
416       fHistList->Add(fh1PtGenOut[ij]);
417       fHistList->Add(fh2PtFGen[ij]);
418       fHistList->Add(fh2PhiFGen[ij]);
419       fHistList->Add(fh2EtaFGen[ij]);
420       fHistList->Add(fh2PtGenDeltaEta[ij]);
421       fHistList->Add(fh2PtGenDeltaPhi[ij]);
422       fHistList->Add(fh2PtRecDeltaR[ij]);
423       fHistList->Add(fh2PtGenDeltaR[ij]);
424       fHistList->Add(fh3RecEtaPhiPt[ij]);
425       fHistList->Add(fh3GenEtaPhiPt[ij]);      
426       if(saveLevel>2){
427         fHistList->Add(fh3RecEtaPhiPtNoGen[ij]);
428         fHistList->Add(fh3GenEtaPhiPtNoFound[ij]);
429       }
430     }
431   }
432
433   // =========== Switch on Sumw2 for all histos ===========
434   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
435     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
436     if (h1){
437       // Printf("%s ",h1->GetName());
438       h1->Sumw2();
439       continue;
440     }
441     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
442     if(hn)hn->Sumw2();
443   }
444
445   TH1::AddDirectory(oldStatus);
446
447 }
448
449 void AliAnalysisTaskJetSpectrum::Init()
450 {
451   //
452   // Initialization
453   //
454
455   Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
456   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
457
458 }
459
460 void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
461 {
462   //
463   // Execute analysis for current event
464   //
465
466
467
468   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
469
470   
471   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
472
473   if(!aodH){
474     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
475     return;
476   }
477
478
479   //  aodH->SelectEvent(kTRUE);
480
481   // ========= These pointers need to be valid in any case ======= 
482
483
484   /*
485   AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
486   if(!jhRec){
487     Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
488     return;
489   }
490   */
491   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
492   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
493   if(!aodRecJets){
494     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
495     return;
496   }
497
498   // ==== General variables needed
499
500
501   // We use statice array, not to fragment the memory
502   AliAODJet genJets[kMaxJets];
503   Int_t nGenJets = 0;
504   AliAODJet recJets[kMaxJets];
505   Int_t nRecJets = 0;
506   ///////////////////////////
507   Int_t nTracks = 0;
508   //////////////////////////  
509
510   Double_t eventW = 1;
511   Double_t ptHard = 0; 
512   Double_t nTrials = 1; // Trials for MC trigger weigth for real data
513
514   if(fUseExternalWeightOnly){
515     eventW = fExternalWeight;
516   }
517
518
519   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
520   if((fAnalysisType&kAnaMC)==kAnaMC){
521     // this is the part we only use when we have MC information
522     AliMCEvent* mcEvent = MCEvent();
523     //    AliStack *pStack = 0; 
524     if(!mcEvent){
525       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
526       return;
527     }
528     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
529     if(!pythiaGenHeader){
530       return;
531     }
532
533     nTrials = pythiaGenHeader->Trials();
534     ptHard  = pythiaGenHeader->GetPtHard();
535     int iProcessType = pythiaGenHeader->ProcessType();
536     // 11 f+f -> f+f
537     // 12 f+barf -> f+barf
538     // 13 f+barf -> g+g
539     // 28 f+g -> f+g
540     // 53 g+g -> f+barf
541     // 68 g+g -> g+g
542     /*
543     if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
544     //    if(iProcessType != 13 && iProcessType != 68){ // allow only glue
545     if(iProcessType != 11 && iProcessType != 12 && iProcessType != 53){ // allow only quark
546     //    if(iProcessType != 28){ // allow only -> f+g
547       PostData(1, fHistList);
548       return;
549     }
550     */
551     if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
552
553     if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
554
555     if(!fUseExternalWeightOnly){
556         // case were we combine more than one p_T hard bin...
557     }
558
559     // fetch the pythia generated jets only to be used here
560     Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
561     AliAODJet pythiaGenJets[kMaxJets];
562     Int_t iCount = 0;
563     for(int ip = 0;ip < nPythiaGenJets;++ip){
564       if(iCount>=kMaxJets)continue;
565       Float_t p[4];
566       pythiaGenHeader->TriggerJet(ip,p);
567       pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
568
569       if(fLimitGenJetEta){
570         if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
571            pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
572       }
573
574
575       if(fBranchGen.Length()==0){
576         // if we have MC particles and we do not read from the aod branch
577         // use the pythia jets
578         genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
579       }
580       iCount++;
581     }
582     if(fBranchGen.Length()==0)nGenJets = iCount;    
583
584   }// (fAnalysisType&kMC)==kMC)
585
586   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
587   fh1PtHard->Fill(ptHard,eventW);
588   fh1PtHardNoW->Fill(ptHard,1);
589   fh1PtHardTrials->Fill(ptHard,nTrials);
590
591   // If we set a second branch for the input jets fetch this 
592   if(fBranchGen.Length()>0){
593     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
594     if(aodGenJets){
595       Int_t iCount = 0;
596       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
597         if(iCount>=kMaxJets)continue;
598         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
599         if(!tmp)continue;
600         if(fLimitGenJetEta){
601           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
602              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
603         }
604         genJets[iCount] = *tmp;
605         iCount++;
606       }
607       nGenJets = iCount;
608     }
609     else{
610       Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
611     }
612   }
613
614   fh1NGenJets->Fill(nGenJets);
615   // We do not want to exceed the maximum jet number
616   nGenJets = TMath::Min(nGenJets,kMaxJets);
617
618   // Fetch the reconstructed jets...
619
620
621   nRecJets = aodRecJets->GetEntries();
622   fh1NRecJets->Fill(nRecJets);
623   nRecJets = TMath::Min(nRecJets,kMaxJets);
624   //////////////////////////////////////////
625   nTracks  = fAOD->GetNumberOfTracks();
626   ///////////////////////////////////////////  
627
628   for(int ir = 0;ir < nRecJets;++ir){
629     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
630     if(!tmp)continue;
631     recJets[ir] = *tmp;
632   }
633
634
635   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
636   // Relate the jets
637   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
638   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
639   
640   for(int i = 0;i<kMaxJets;++i){
641     iGenIndex[i] = iRecIndex[i] = -1;
642   }
643
644
645   GetClosestJets(genJets,nGenJets,recJets,nRecJets,
646                  iGenIndex,iRecIndex,fDebug);
647   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
648
649   if(fDebug){
650     for(int i = 0;i<kMaxJets;++i){
651       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]); 
652       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]); 
653     }
654   }
655
656   // loop over reconstructed jets
657   for(int ir = 0;ir < nRecJets;++ir){
658     Double_t ptRec = recJets[ir].Pt();
659     Double_t phiRec = recJets[ir].Phi();
660     if(phiRec<0)phiRec+=TMath::Pi()*2.;    
661     Double_t etaRec = recJets[ir].Eta();
662     if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
663     fh1E[ir]->Fill(recJets[ir].E(),eventW);
664     fh1PtRecIn[ir]->Fill(ptRec,eventW);
665     fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
666     for(int irr = ir+1;irr<nRecJets;irr++){
667       fh2PtRecDeltaR[ir]->Fill(recJets[irr].Pt(),recJets[ir].DeltaR(&recJets[irr]));
668     }
669     // Fill Correlation
670     Int_t ig = iGenIndex[ir];
671     if(ig>=0 && ig<nGenJets){
672       if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
673       if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
674       fh1PtRecOut[ir]->Fill(ptRec,eventW);
675       Double_t ptGen  = genJets[ig].Pt();
676       Double_t phiGen = genJets[ig].Phi();
677       if(phiGen<0)phiGen+=TMath::Pi()*2.; 
678       Double_t etaGen = genJets[ig].Eta();
679
680       // 
681       // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
682       // 
683
684       if(TMath::Abs(etaRec)<fRecEtaWindow){
685
686       fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
687       fh2PhiFGen[ir]->Fill(phiRec,phiGen,eventW);
688       fh2EtaFGen[ir]->Fill(etaRec,etaGen,eventW);
689       fh2PtGenDeltaEta[ir]->Fill(ptGen,etaGen-etaRec,eventW);
690       fh2PtGenDeltaPhi[ir]->Fill(ptGen,phiGen-phiRec,eventW);
691       fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
692       fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1);
693       /////////////////////////////////////////////////////
694
695       //      Double_t eRec = recJets[ir].E();
696       //      Double_t eGen = genJets[ig].E();
697       // CKB use p_T not Energy 
698       // TODO recname variabeles and histos
699       Double_t eRec = recJets[ir].E();
700       Double_t eGen = genJets[ig].E();
701
702       fh2Efficiency->Fill(eGen, eRec/eGen);
703
704       if (eGen>=0. && eGen<=250.){
705         Double_t eLeading = -1;
706         Double_t ptleading = -1;
707         Int_t nPart=0;
708         // loop over tracks
709         for (Int_t it = 0; it< nTracks; it++){
710           //      if (fAOD->GetTrack(it)->E() > eGen) continue; // CKB. Not allowed! cannot cut on gen properties in real events!
711           // find leading particle
712           //  if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){
713           // TODO implement esd filter flag to be the same as in the jet finder 
714           // allow also for MC particles...
715           Float_t r = recJets[ir].DeltaR(fAOD->GetTrack(it));
716           if (r<0.4 && fAOD->GetTrack(it)->Pt()>ptleading){
717             eLeading  = fAOD->GetTrack(it)->E();
718             ptleading = fAOD->GetTrack(it)->Pt();            
719           }
720           //          if (fAOD->GetTrack(it)->Pt()>0.03*eGen && fAOD->GetTrack(it)->E()<=eGen && r<0.7) // CKB cannot cut on gen properties 
721           if (fAOD->GetTrack(it)->Pt()>0.03*eRec && fAOD->GetTrack(it)->Pt()<=eRec && r<0.7)
722             nPart++;
723         }
724         if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
725
726         // fill Response Map (4D histogram) and Energy vs z distributions
727         Double_t var[4] = {eGen, eRec, ptleading/eGen, ptleading/eRec};                       
728         fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
729         fh2EGenZGen->Fill(var[0],var[2]);
730         fh1JetMultiplicity->Fill(nPart);
731         fh3EGenERecN->Fill(eGen, eRec, nPart); 
732         for(int ic = 0;ic <kMaxCorrelation;ic++){
733           if (nPart<=fgkJetNpartCut[ic]){ // is this corrected for CKB
734             fhnCorrelation[ic]->Fill(var);
735             break;
736           }
737         }
738       }
739
740     }// if etarec in window
741
742     } 
743     ////////////////////////////////////////////////////  
744     else{
745       fh3RecEtaPhiPtNoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
746     }
747   }// loop over reconstructed jets
748
749
750   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
751   for(int ig = 0;ig < nGenJets;++ig){
752     Double_t ptGen = genJets[ig].Pt();
753     // Fill Correlation
754     Double_t phiGen = genJets[ig].Phi();
755     if(phiGen<0)phiGen+=TMath::Pi()*2.;    
756     Double_t etaGen = genJets[ig].Eta();
757     fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
758     fh1PtGenIn[ig]->Fill(ptGen,eventW);
759     for(int igg = ig+1;igg<nGenJets;igg++){
760       fh2PtGenDeltaR[ig]->Fill(genJets[igg].Pt(),genJets[ig].DeltaR(&genJets[igg]));
761     }
762     Int_t ir = iRecIndex[ig];
763     if(ir>=0&&ir<nRecJets){
764       fh1PtGenOut[ig]->Fill(ptGen,eventW);
765     }
766     else{
767       fh3GenEtaPhiPtNoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
768     }
769   }// loop over reconstructed jets
770
771   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
772   PostData(1, fHistList);
773 }
774
775 void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
776 {
777 // Terminate analysis
778 //
779     if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
780 }
781
782
783 void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,const Int_t &nGenJets,
784                                                 AliAODJet *recJets,const Int_t &nRecJets,
785                                                 Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
786
787   //
788   // Relate the two input jet Arrays
789   //
790
791   //
792   // The association has to be unique
793   // So check in two directions
794   // find the closest rec to a gen
795   // and check if there is no other rec which is closer
796   // Caveat: Close low energy/split jets may disturb this correlation
797
798   // Idea: search in two directions generated e.g (a--e) and rec (1--3)
799   // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
800   // in the end we have something like this
801   //    1   2   3
802   // ------------
803   // a| 3   2   0
804   // b| 0   1   0
805   // c| 0   0   3
806   // d| 0   0   1
807   // e| 0   0   1
808   // Topology
809   //   1     2
810   //     a         b        
811   //
812   //  d      c
813   //        3     e
814   // Only entries with "3" match from both sides
815   
816   const int kMode = 3;
817
818   Int_t iFlag[kMaxJets][kMaxJets];
819
820
821
822   for(int i = 0;i < kMaxJets;++i){
823     iRecIndex[i] = -1;
824     iGenIndex[i] = -1;
825     for(int j = 0;j < kMaxJets;++j)iFlag[i][j] = 0;
826   }
827
828   if(nRecJets==0)return;
829   if(nGenJets==0)return;
830
831   const Float_t maxDist = 0.5;
832   // find the closest distance to the generated
833   for(int ig = 0;ig<nGenJets;++ig){
834     Float_t dist = maxDist;
835     if(iDebug>1)Printf("Gen (%d) p_T %3.3f eta %3.3f ph %3.3f ",ig,genJets[ig].Pt(),genJets[ig].Eta(),genJets[ig].Phi());
836     for(int ir = 0;ir<nRecJets;++ir){
837       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
838       if(iDebug>1)Printf("Rec (%d) p_T %3.3f eta %3.3f ph %3.3f ",ir,recJets[ir].Pt(),recJets[ir].Eta(),recJets[ir].Phi());
839       if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
840       if(dR<dist){
841         iRecIndex[ig] = ir;
842         dist = dR;
843       } 
844     }
845     if(iRecIndex[ig]>=0)iFlag[ig][iRecIndex[ig]]+=1;
846     // reset...
847     iRecIndex[ig] = -1;
848   }
849   // other way around
850   for(int ir = 0;ir<nRecJets;++ir){
851     Float_t dist = maxDist;
852     for(int ig = 0;ig<nGenJets;++ig){
853       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
854       if(dR<dist){
855         iGenIndex[ir] = ig;
856         dist = dR;
857       } 
858     }
859     if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]][ir]+=2;
860     // reset...
861     iGenIndex[ir] = -1;
862   }
863
864   // check for "true" correlations
865
866   if(iDebug>1)Printf(">>>>>> Matrix");
867
868   for(int ig = 0;ig<nGenJets;++ig){
869     for(int ir = 0;ir<nRecJets;++ir){
870       // Print
871       if(iDebug>1)printf("XFL %d ",iFlag[ig][ir]);
872
873       if(kMode==3){
874         // we have a uniqie correlation
875         if(iFlag[ig][ir]==3){
876           iGenIndex[ir] = ig;
877           iRecIndex[ig] = ir;
878         }
879       }
880       else{
881         // we just take the correlation from on side
882         if((iFlag[ig][ir]&2)==2){
883           iGenIndex[ir] = ig;
884         }
885         if((iFlag[ig][ir]&1)==1){
886           iRecIndex[ig] = ir;
887         }
888       }
889     }
890     if(iDebug>1)printf("\n");
891   }
892 }
893
894