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