Jet and Particle identification tasks moved to different directories
[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 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
57   fJetHeaderRec(0x0),
58   fJetHeaderGen(0x0),
59   fAOD(0x0),
60   fBranchRec("jets"),
61   fConfigRec("ConfigJets.C"),
62   fBranchGen(""),
63   fConfigGen(""),
64   fUseAODInput(kFALSE),
65   fUseExternalWeightOnly(kFALSE),
66   fLimitGenJetEta(kFALSE),
67   fAnalysisType(0),
68   fExternalWeight(1),    
69   fh1PtHard(0x0),
70   fh1PtHard_NoW(0x0),  
71   fh1PtHard_Trials(0x0),
72   fh1NGenJets(0x0),
73   fh1NRecJets(0x0),
74   fHistList(0x0)
75 {
76   // Default constructor
77     for(int ij  = 0;ij<kMaxJets;++ij){
78     fh1E[ij] =  fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
79     fh2PtFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] = 0;
80     fh3PtRecGenHard[ij] =  fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
81   }
82   
83 }
84
85 AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
86   AliAnalysisTaskSE(name),
87   fJetHeaderRec(0x0),
88   fJetHeaderGen(0x0),
89   fAOD(0x0),
90   fBranchRec("jets"),
91   fConfigRec("ConfigJets.C"),
92   fBranchGen(""),
93   fConfigGen(""),
94   fUseAODInput(kFALSE),
95   fUseExternalWeightOnly(kFALSE),
96   fLimitGenJetEta(kFALSE),
97   fAnalysisType(0),
98   fExternalWeight(1),    
99   fh1PtHard(0x0),
100   fh1PtHard_NoW(0x0),  
101   fh1PtHard_Trials(0x0),
102   fh1NGenJets(0x0),
103   fh1NRecJets(0x0),
104   fHistList(0x0)
105 {
106   // Default constructor
107   for(int ij  = 0;ij<kMaxJets;++ij){
108     fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
109     fh2PtFGen[ij] =  fh2Frag[ij] = fh2FragLn[ij] = 0;
110
111     fh3PtRecGenHard[ij] =  fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] =  fh3GenEtaPhiPt[ij] = 0;
112   }
113
114   DefineOutput(1, TList::Class());  
115 }
116
117
118
119
120 void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects()
121 {
122
123   //
124   // Create the output container
125   //
126
127   
128   // Connect the AOD
129
130   if(fUseAODInput){
131     fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
132     if(!fAOD){
133       Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODInput);
134       return;
135     }
136     // fethc the header
137     fJetHeaderRec = dynamic_cast<AliJetHeader*>(fInputHandler->GetTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));
138     if(!fJetHeaderRec){
139       Printf("%s:%d Jet Header not found in the Input",(char*)__FILE__,__LINE__);
140     }
141   }
142   else{
143     //  assume that the AOD is in the general output...
144     fAOD  = AODEvent();
145     if(!fAOD){
146       Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
147       return;
148     }
149     fJetHeaderRec = dynamic_cast<AliJetHeader*>(OutputTree()->GetUserInfo()->FindObject(Form("AliJetHeader_%s",fBranchRec.Data())));    
150     if(!fJetHeaderRec){
151       Printf("%s:%d Jet Header not found in the Output",(char*)__FILE__,__LINE__);
152     }
153     else{
154       if(fDebug>10)fJetHeaderRec->Dump();
155     }
156   }
157  
158
159
160   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::UserCreateOutputObjects() \n");
161
162   OpenFile(1);
163   if(!fHistList)fHistList = new TList();
164
165   Bool_t oldStatus = TH1::AddDirectoryStatus();
166   TH1::AddDirectory(kFALSE);
167
168   //
169   //  Histogram
170     
171   const Int_t nBinPt = 100;
172   Double_t binLimitsPt[nBinPt+1];
173   for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
174     if(iPt == 0){
175       binLimitsPt[iPt] = 0.0;
176     }
177     else {// 1.0
178       binLimitsPt[iPt] =  binLimitsPt[iPt-1] + 2;
179     }
180   }
181   const Int_t nBinEta = 22;
182   Double_t binLimitsEta[nBinEta+1];
183   for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
184     if(iEta==0){
185       binLimitsEta[iEta] = -2.2;
186     }
187     else{
188       binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.2;
189     }
190   }
191
192   const Int_t nBinPhi = 90;
193   Double_t binLimitsPhi[nBinPhi+1];
194   for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
195     if(iPhi==0){
196       binLimitsPhi[iPhi] = 0;
197     }
198     else{
199       binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
200     }
201   }
202
203   const Int_t nBinFrag = 25;
204
205
206   fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
207
208   fh1PtHard_NoW = new TH1F("fh1PtHard_NoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
209
210   fh1PtHard_Trials = new TH1F("fh1PtHard_Trials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
211
212   fh1NGenJets  = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
213
214   fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
215
216
217   for(int ij  = 0;ij<kMaxJets;++ij){
218     fh1E[ij] = new TH1F(Form("fh1E_j%d",ij),"Jet Energy;E_{jet} (GeV);N",nBinPt,binLimitsPt);
219     fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
220     fh1PtRecOut[ij] = new TH1F(Form("fh1PtRecOut_j%d",ij),"rec p_T output jets;p_{T,rec}",nBinPt,binLimitsPt);
221     fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
222     fh1PtGenOut[ij] = new TH1F(Form("fh1PtGenOut_j%d",ij),"found p_T output jets;p_{T,gen}",nBinPt,binLimitsPt);
223
224
225     fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)",
226                              nBinPt,binLimitsPt,nBinPt,binLimitsPt);
227
228
229
230
231
232
233     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);
234
235
236
237     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);
238
239         
240     fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
241                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
242
243     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",
244                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
245
246     fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
247                                   nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
248
249
250
251     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)",
252                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
253
254
255     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)",
256                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
257
258
259
260     fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
261                                  nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
262
263   }
264
265   const Int_t saveLevel = 2; // large save level more histos
266
267   if(saveLevel>0){
268     fHistList->Add(fh1PtHard);
269     fHistList->Add(fh1PtHard_NoW);
270     fHistList->Add(fh1PtHard_Trials);
271     fHistList->Add(fh1NGenJets);
272     fHistList->Add(fh1NRecJets);
273     for(int ij  = 0;ij<kMaxJets;++ij){
274       fHistList->Add(fh1E[ij]);
275       fHistList->Add(fh1PtRecIn[ij]);
276       fHistList->Add(fh1PtRecOut[ij]);
277       fHistList->Add(fh1PtGenIn[ij]);
278       fHistList->Add(fh1PtGenOut[ij]);
279       fHistList->Add(fh2PtFGen[ij]);
280       if(saveLevel>2){
281         fHistList->Add(fh3RecEtaPhiPt[ij]);
282         fHistList->Add(fh3RecEtaPhiPt_NoGen[ij]);
283         fHistList->Add(fh3GenEtaPhiPt_NoFound[ij]);
284         fHistList->Add(fh3GenEtaPhiPt[ij]);      
285       }
286     }
287   }
288
289   // =========== Switch on Sumw2 for all histos ===========
290   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
291     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
292     if (h1){
293       // Printf("%s ",h1->GetName());
294       h1->Sumw2();
295     }
296   }
297
298   TH1::AddDirectory(oldStatus);
299
300 }
301
302 void AliAnalysisTaskJetSpectrum::Init()
303 {
304   //
305   // Initialization
306   //
307
308   Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
309   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
310
311 }
312
313 void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
314 {
315   //
316   // Execute analysis for current event
317   //
318
319
320
321   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
322
323   
324   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
325
326   if(!aodH){
327     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
328     return;
329   }
330
331
332   //  aodH->SelectEvent(kTRUE);
333
334   // ========= These pointers need to be valid in any case ======= 
335
336
337   /*
338   AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
339   if(!jhRec){
340     Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
341     return;
342   }
343   */
344   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
345   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
346   if(!aodRecJets){
347     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
348     return;
349   }
350
351   // ==== General variables needed
352
353
354   // We use statice array, not to fragment the memory
355   AliAODJet genJets[kMaxJets];
356   Int_t nGenJets = 0;
357   AliAODJet recJets[kMaxJets];
358   Int_t nRecJets = 0;
359
360   Double_t eventW = 1;
361   Double_t ptHard = 0; 
362   Double_t nTrials = 1; // Trials for MC trigger weigth for real data
363
364   if(fUseExternalWeightOnly){
365     eventW = fExternalWeight;
366   }
367
368
369   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
370   if((fAnalysisType&kAnaMC)==kAnaMC){
371     // this is the part we only use when we have MC information
372     AliMCEvent* mcEvent = MCEvent();
373     //    AliStack *pStack = 0; 
374     if(!mcEvent){
375       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
376       return;
377     }
378     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
379     if(!pythiaGenHeader){
380       return;
381     }
382
383     nTrials = pythiaGenHeader->Trials();
384     ptHard  = pythiaGenHeader->GetPtHard();
385
386
387     if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
388
389     if(!fUseExternalWeightOnly){
390         // case were we combine more than one p_T hard bin...
391     }
392
393     // fetch the pythia generated jets only to be used here
394     Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
395     AliAODJet pythiaGenJets[kMaxJets];
396     Int_t iCount = 0;
397     for(int ip = 0;ip < nPythiaGenJets;++ip){
398       if(iCount>=kMaxJets)continue;
399       Float_t p[4];
400       pythiaGenHeader->TriggerJet(ip,p);
401       pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
402
403       if(fLimitGenJetEta){
404         if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
405            pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
406       }
407
408
409       if(fBranchGen.Length()==0){
410         // if we have MC particles and we do not read from the aod branch
411         // use the pythia jets
412         genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
413       }
414       iCount++;
415     }
416     if(fBranchGen.Length()==0)nGenJets = iCount;    
417
418   }// (fAnalysisType&kMC)==kMC)
419
420   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
421   fh1PtHard->Fill(ptHard,eventW);
422   fh1PtHard_NoW->Fill(ptHard,1);
423   fh1PtHard_Trials->Fill(ptHard,nTrials);
424
425   // If we set a second branch for the input jets fetch this 
426   if(fBranchGen.Length()>0){
427     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
428     if(aodGenJets){
429       Int_t iCount = 0;
430       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
431         if(iCount>=kMaxJets)continue;
432         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
433         if(!tmp)continue;
434         if(fLimitGenJetEta){
435           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
436              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
437         }
438         genJets[iCount] = *tmp;
439         iCount++;
440       }
441       nGenJets = iCount;
442     }
443     else{
444       Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
445     }
446   }
447
448   fh1NGenJets->Fill(nGenJets);
449   // We do not want to exceed the maximum jet number
450   nGenJets = TMath::Min(nGenJets,kMaxJets);
451
452   // Fetch the reconstructed jets...
453
454
455   nRecJets = aodRecJets->GetEntries();
456   fh1NRecJets->Fill(nRecJets);
457   nRecJets = TMath::Min(nRecJets,kMaxJets);
458
459   for(int ir = 0;ir < nRecJets;++ir){
460     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
461     if(!tmp)continue;
462     recJets[ir] = *tmp;
463   }
464
465
466   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
467   // Relate the jets
468   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
469   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
470   
471   for(int i = 0;i<kMaxJets;++i){
472     iGenIndex[i] = iRecIndex[i] = -1;
473   }
474
475
476   GetClosestJets(genJets,nGenJets,recJets,nRecJets,
477                  iGenIndex,iRecIndex,fDebug);
478   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
479
480   if(fDebug){
481     for(int i = 0;i<kMaxJets;++i){
482       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]); 
483       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]); 
484     }
485   }
486
487   // loop over reconstructed jets
488   for(int ir = 0;ir < nRecJets;++ir){
489     Double_t ptRec = recJets[ir].Pt();
490     Double_t phiRec = recJets[ir].Phi();
491     if(phiRec<0)phiRec+=TMath::Pi()*2.;    
492     Double_t etaRec = recJets[ir].Eta();
493
494     fh1E[ir]->Fill(recJets[ir].E(),eventW);
495     fh1PtRecIn[ir]->Fill(ptRec,eventW);
496     fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
497     // Fill Correlation
498     Int_t ig = iGenIndex[ir];
499     if(ig>=0&&ig<nGenJets){
500       fh1PtRecOut[ir]->Fill(ptRec,eventW);
501       Double_t ptGen = genJets[ig].Pt();
502       fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
503       fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
504       fh3PtRecGenHard_NoW[ir]->Fill(ptRec,ptGen,ptHard,1);
505     }
506     else{
507       fh3RecEtaPhiPt_NoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
508     }
509   }// loop over reconstructed jets
510
511
512   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
513   for(int ig = 0;ig < nGenJets;++ig){
514     Double_t ptGen = genJets[ig].Pt();
515     // Fill Correlation
516     Double_t phiGen = genJets[ig].Phi();
517     if(phiGen<0)phiGen+=TMath::Pi()*2.;    
518     Double_t etaGen = genJets[ig].Eta();
519     fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
520     fh1PtGenIn[ig]->Fill(ptGen,eventW);
521     Int_t ir = iRecIndex[ig];
522     if(ir>=0&&ir<nRecJets){
523       fh1PtGenOut[ig]->Fill(ptGen,eventW);
524     }
525     else{
526       fh3GenEtaPhiPt_NoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
527     }
528   }// loop over reconstructed jets
529
530   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
531   PostData(1, fHistList);
532 }
533
534 void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
535 {
536 // Terminate analysis
537 //
538     if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
539 }
540
541
542 void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,Int_t &nGenJets,
543                                                 AliAODJet *recJets,Int_t &nRecJets,
544                                                 Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
545
546   //
547   // Relate the two input jet Arrays
548   //
549
550   //
551   // The association has to be unique
552   // So check in two directions
553   // find the closest rec to a gen
554   // and check if there is no other rec which is closer
555   // Caveat: Close low energy/split jets may disturb this correlation
556
557   // Idea: search in two directions generated e.g (a--e) and rec (1--3)
558   // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
559   // in the end we have something like this
560   //    1   2   3
561   // ------------
562   // a| 3   2   0
563   // b| 0   1   0
564   // c| 0   0   3
565   // d| 0   0   1
566   // e| 0   0   1
567   // Topology
568   //   1     2
569   //     a         b        
570   //
571   //  d      c
572   //        3     e
573   // Only entries with "3" match from both sides
574   
575   const int kMode = 3;
576
577   Int_t iFlag[kMaxJets][kMaxJets];
578
579
580
581   for(int i = 0;i < kMaxJets;++i){
582     iRecIndex[i] = -1;
583     iGenIndex[i] = -1;
584     for(int j = 0;j < kMaxJets;++j)iFlag[i][j] = 0;
585   }
586
587   if(nRecJets==0)return;
588   if(nGenJets==0)return;
589
590   const Float_t maxDist = 1.0;
591   // find the closest distance to the generated
592   for(int ig = 0;ig<nGenJets;++ig){
593     Float_t dist = maxDist;
594     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());
595     for(int ir = 0;ir<nRecJets;++ir){
596       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
597       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());
598       if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
599       if(dR<dist){
600         iRecIndex[ig] = ir;
601         dist = dR;
602       } 
603     }
604     if(iRecIndex[ig]>=0)iFlag[ig][iRecIndex[ig]]+=1;
605     // reset...
606     iRecIndex[ig] = -1;
607   }
608   // other way around
609   for(int ir = 0;ir<nRecJets;++ir){
610     Float_t dist = maxDist;
611     for(int ig = 0;ig<nGenJets;++ig){
612       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
613       if(dR<dist){
614         iGenIndex[ir] = ig;
615         dist = dR;
616       } 
617     }
618     if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]][ir]+=2;
619     // reset...
620     iGenIndex[ir] = -1;
621   }
622
623   // check for "true" correlations
624
625   if(iDebug>1)Printf(">>>>>> Matrix");
626
627   for(int ig = 0;ig<nGenJets;++ig){
628     for(int ir = 0;ir<nRecJets;++ir){
629       // Print
630       if(iDebug>1)printf("XFL %d ",iFlag[ig][ir]);
631
632       if(kMode==3){
633         // we have a uniqie correlation
634         if(iFlag[ig][ir]==3){
635           iGenIndex[ir] = ig;
636           iRecIndex[ig] = ir;
637         }
638       }
639       else{
640         // we just take the correlation from on side
641         if((iFlag[ig][ir]&2)==2){
642           iGenIndex[ir] = ig;
643         }
644         if((iFlag[ig][ir]&1)==1){
645           iRecIndex[ig] = ir;
646         }
647       }
648     }
649     if(iDebug>1)printf("\n");
650   }
651 }
652
653