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