added merging of different pt_hard bins in macro and Notify() in task to fill xsection
[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] =  fh2Frag[ij] = fh2FragLn[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] =  fh2Frag[ij] = fh2FragLn[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 = 90;
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
309
310
311
312
313     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);
314
315
316
317     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);
318
319         
320     fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
321                            nBinFrag,0.,1.,nBinPt,binLimitsPt);
322
323     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",
324                              nBinFrag,0.,10.,nBinPt,binLimitsPt);
325
326     fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
327                                   nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
328
329
330
331     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)",
332                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
333
334
335     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)",
336                                         nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
337
338
339
340     fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
341                                  nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
342
343   }
344   
345   /////////////////////////////////////////////////////////////////
346   fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
347
348   fh2ERecZRec   = new TH2F("fh2ERecZRec", " ; E^{jet}_{rec} [GeV]; z^{lp}_{rec}", 100, 0., 250., 100, 0., 2.);
349   fh2EGenZGen   = new TH2F("fh2EGenZGen", " ; E^{jet}_{gen} [GeV]; z^{lp}_{gen}", 100, 0., 250., 100, 0., 2.);
350   fh2Efficiency = new TH2F("fh2Efficiency", "ERec/EGen;E^{jet}_{gen} [GeV];E^{jet}_{rec}/E^{jet}_{gen}", 100, 0., 250., 100, 0., 1.5);  
351
352   fh3EGenERecN  = new TH3F("fh3EGenERecN", "Efficiency vs. Jet Multiplicity", 100., 0., 250., 100, 0., 250., 51, 0., 50.);
353
354   // Response map  
355   //arrays for bin limits
356   const Int_t nbin[4] = {100, 100, 100, 100}; 
357   Double_t LowEdge[4] = {0.,0.,0.,0.};
358   Double_t UpEdge[4] = {250., 250., 1., 1.};
359
360   for(int ic = 0;ic < kMaxCorrelation;ic++){
361     fhnCorrelation[ic]  = new THnSparseF(Form("fhnCorrelation_%d",ic),  "Response Map", 4, nbin, LowEdge, UpEdge);
362     if(ic==0) fhnCorrelation[ic]->SetTitle(Form("ResponseMap 0 <= npart <= %.0E",fgkJetNpartCut[ic]));
363     else fhnCorrelation[ic]->SetTitle(Form("ResponseMap %.0E < npart <= %.0E",fgkJetNpartCut[ic-1],fgkJetNpartCut[ic]));
364   }
365   const Int_t saveLevel = 1; // large save level more histos
366   if(saveLevel>0){
367     fHistList->Add(fh1Xsec);
368     fHistList->Add(fh1Trials);
369     fHistList->Add(fh1PtHard);
370     fHistList->Add(fh1PtHard_NoW);
371     fHistList->Add(fh1PtHard_Trials);
372     fHistList->Add(fh1NGenJets);
373     fHistList->Add(fh1NRecJets);
374     ////////////////////////
375     fHistList->Add(fh1JetMultiplicity);     
376     fHistList->Add(fh2ERecZRec);
377     fHistList->Add(fh2EGenZGen);
378     fHistList->Add(fh2Efficiency);
379     fHistList->Add(fh3EGenERecN);
380
381     for(int ic = 0;ic < kMaxCorrelation;++ic){
382       fHistList->Add(fhnCorrelation[ic]);
383     }
384     //////////////////////// 
385     for(int ij  = 0;ij<kMaxJets;++ij){
386       fHistList->Add(fh1E[ij]);
387       fHistList->Add(fh1PtRecIn[ij]);
388       fHistList->Add(fh1PtRecOut[ij]);
389       fHistList->Add(fh1PtGenIn[ij]);
390       fHistList->Add(fh1PtGenOut[ij]);
391       fHistList->Add(fh2PtFGen[ij]);
392       if(saveLevel>2){
393         fHistList->Add(fh3RecEtaPhiPt[ij]);
394         fHistList->Add(fh3RecEtaPhiPt_NoGen[ij]);
395         fHistList->Add(fh3GenEtaPhiPt_NoFound[ij]);
396         fHistList->Add(fh3GenEtaPhiPt[ij]);      
397       }
398     }
399   }
400
401   // =========== Switch on Sumw2 for all histos ===========
402   for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
403     TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
404     if (h1){
405       // Printf("%s ",h1->GetName());
406       h1->Sumw2();
407       continue;
408     }
409     THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
410     if(hn)hn->Sumw2();
411   }
412
413   TH1::AddDirectory(oldStatus);
414
415 }
416
417 void AliAnalysisTaskJetSpectrum::Init()
418 {
419   //
420   // Initialization
421   //
422
423   Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
424   if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
425
426 }
427
428 void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
429 {
430   //
431   // Execute analysis for current event
432   //
433
434
435
436   if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
437
438   
439   AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
440
441   if(!aodH){
442     Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
443     return;
444   }
445
446
447   //  aodH->SelectEvent(kTRUE);
448
449   // ========= These pointers need to be valid in any case ======= 
450
451
452   /*
453   AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
454   if(!jhRec){
455     Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
456     return;
457   }
458   */
459   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
460   TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
461   if(!aodRecJets){
462     Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
463     return;
464   }
465
466   // ==== General variables needed
467
468
469   // We use statice array, not to fragment the memory
470   AliAODJet genJets[kMaxJets];
471   Int_t nGenJets = 0;
472   AliAODJet recJets[kMaxJets];
473   Int_t nRecJets = 0;
474   ///////////////////////////
475   Int_t nTracks = 0;
476   //////////////////////////  
477
478   Double_t eventW = 1;
479   Double_t ptHard = 0; 
480   Double_t nTrials = 1; // Trials for MC trigger weigth for real data
481
482   if(fUseExternalWeightOnly){
483     eventW = fExternalWeight;
484   }
485
486
487   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
488   if((fAnalysisType&kAnaMC)==kAnaMC){
489     // this is the part we only use when we have MC information
490     AliMCEvent* mcEvent = MCEvent();
491     //    AliStack *pStack = 0; 
492     if(!mcEvent){
493       Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
494       return;
495     }
496     AliGenPythiaEventHeader*  pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
497     if(!pythiaGenHeader){
498       return;
499     }
500
501     nTrials = pythiaGenHeader->Trials();
502     ptHard  = pythiaGenHeader->GetPtHard();
503
504
505     if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
506
507     if(!fUseExternalWeightOnly){
508         // case were we combine more than one p_T hard bin...
509     }
510
511     // fetch the pythia generated jets only to be used here
512     Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
513     AliAODJet pythiaGenJets[kMaxJets];
514     Int_t iCount = 0;
515     for(int ip = 0;ip < nPythiaGenJets;++ip){
516       if(iCount>=kMaxJets)continue;
517       Float_t p[4];
518       pythiaGenHeader->TriggerJet(ip,p);
519       pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
520
521       if(fLimitGenJetEta){
522         if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
523            pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
524       }
525
526
527       if(fBranchGen.Length()==0){
528         // if we have MC particles and we do not read from the aod branch
529         // use the pythia jets
530         genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
531       }
532       iCount++;
533     }
534     if(fBranchGen.Length()==0)nGenJets = iCount;    
535
536   }// (fAnalysisType&kMC)==kMC)
537
538   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
539   fh1PtHard->Fill(ptHard,eventW);
540   fh1PtHard_NoW->Fill(ptHard,1);
541   fh1PtHard_Trials->Fill(ptHard,nTrials);
542
543   // If we set a second branch for the input jets fetch this 
544   if(fBranchGen.Length()>0){
545     TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
546     if(aodGenJets){
547       Int_t iCount = 0;
548       for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
549         if(iCount>=kMaxJets)continue;
550         AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
551         if(!tmp)continue;
552         if(fLimitGenJetEta){
553           if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
554              tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
555         }
556         genJets[iCount] = *tmp;
557         iCount++;
558       }
559       nGenJets = iCount;
560     }
561     else{
562       Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
563     }
564   }
565
566   fh1NGenJets->Fill(nGenJets);
567   // We do not want to exceed the maximum jet number
568   nGenJets = TMath::Min(nGenJets,kMaxJets);
569
570   // Fetch the reconstructed jets...
571
572
573   nRecJets = aodRecJets->GetEntries();
574   fh1NRecJets->Fill(nRecJets);
575   nRecJets = TMath::Min(nRecJets,kMaxJets);
576   //////////////////////////////////////////
577   nTracks  = fAOD->GetNumberOfTracks();
578   ///////////////////////////////////////////  
579
580   for(int ir = 0;ir < nRecJets;++ir){
581     AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
582     if(!tmp)continue;
583     recJets[ir] = *tmp;
584   }
585
586
587   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
588   // Relate the jets
589   Int_t iGenIndex[kMaxJets];    // Index of the generated jet for i-th rec -1 if none
590   Int_t iRecIndex[kMaxJets];    // Index of the rec jet for i-th gen -1 if none
591   
592   for(int i = 0;i<kMaxJets;++i){
593     iGenIndex[i] = iRecIndex[i] = -1;
594   }
595
596
597   GetClosestJets(genJets,nGenJets,recJets,nRecJets,
598                  iGenIndex,iRecIndex,fDebug);
599   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
600
601   if(fDebug){
602     for(int i = 0;i<kMaxJets;++i){
603       if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]); 
604       if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]); 
605     }
606   }
607
608   // loop over reconstructed jets
609   for(int ir = 0;ir < nRecJets;++ir){
610     Double_t ptRec = recJets[ir].Pt();
611     Double_t phiRec = recJets[ir].Phi();
612     if(phiRec<0)phiRec+=TMath::Pi()*2.;    
613     Double_t etaRec = recJets[ir].Eta();
614
615     fh1E[ir]->Fill(recJets[ir].E(),eventW);
616     fh1PtRecIn[ir]->Fill(ptRec,eventW);
617     fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
618     // Fill Correlation
619     Int_t ig = iGenIndex[ir];
620     if(ig>=0 && ig<nGenJets){
621       fh1PtRecOut[ir]->Fill(ptRec,eventW);
622       Double_t ptGen  = genJets[ig].Pt();
623       Double_t phiGen = genJets[ig].Phi();
624       if(phiGen<0)phiGen+=TMath::Pi()*2.; 
625       Double_t etaGen = genJets[ig].Eta();
626       fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
627       fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
628       fh3PtRecGenHard_NoW[ir]->Fill(ptRec,ptGen,ptHard,1);
629       /////////////////////////////////////////////////////
630       Double_t ERec = recJets[ir].E();
631       Double_t EGen = genJets[ig].E();
632       fh2Efficiency->Fill(EGen, ERec/EGen);
633       if (EGen>=0. && EGen<=250.){
634         Double_t Eleading = -1;
635         Double_t ptleading = -1;
636         Int_t nPart=0;
637         // loop over tracks
638         for (Int_t it = 0; it< nTracks; it++){
639           if (fAOD->GetTrack(it)->E() > EGen) continue;
640           Double_t phiTrack = fAOD->GetTrack(it)->Phi();
641           if (phiTrack<0) phiTrack+=TMath::Pi()*2.;            
642           Double_t etaTrack = fAOD->GetTrack(it)->Eta();
643           Float_t deta  = etaRec - etaTrack;
644           Float_t dphi  = TMath::Abs(phiRec - phiTrack);
645           Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);        
646           // find leading particle
647           if (R<0.4 && fAOD->GetTrack(it)->E()>Eleading){
648             Eleading  = fAOD->GetTrack(it)->E();
649             ptleading = fAOD->GetTrack(it)->Pt();            
650           }
651           if (fAOD->GetTrack(it)->Pt()>0.03*EGen && fAOD->GetTrack(it)->E()<=EGen && R<0.7)
652             nPart++;
653         }
654         // fill Response Map (4D histogram) and Energy vs z distributions
655         Double_t var[4] = {EGen, ERec, ptleading/EGen, ptleading/ERec};                       
656         fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
657         fh2EGenZGen->Fill(var[0],var[2]);
658         fh1JetMultiplicity->Fill(nPart);
659         fh3EGenERecN->Fill(EGen, ERec, nPart); 
660         for(int ic = 0;ic <kMaxCorrelation;ic++){
661         if (nPart<=fgkJetNpartCut[ic]){
662           fhnCorrelation[ic]->Fill(var);
663           break;
664         }
665         }
666       }
667     }  
668     ////////////////////////////////////////////////////  
669     else{
670       fh3RecEtaPhiPt_NoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
671     }
672   }// loop over reconstructed jets
673
674
675   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
676   for(int ig = 0;ig < nGenJets;++ig){
677     Double_t ptGen = genJets[ig].Pt();
678     // Fill Correlation
679     Double_t phiGen = genJets[ig].Phi();
680     if(phiGen<0)phiGen+=TMath::Pi()*2.;    
681     Double_t etaGen = genJets[ig].Eta();
682     fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
683     fh1PtGenIn[ig]->Fill(ptGen,eventW);
684     Int_t ir = iRecIndex[ig];
685     if(ir>=0&&ir<nRecJets){
686       fh1PtGenOut[ig]->Fill(ptGen,eventW);
687     }
688     else{
689       fh3GenEtaPhiPt_NoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
690     }
691   }// loop over reconstructed jets
692
693   if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
694   PostData(1, fHistList);
695 }
696
697 void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
698 {
699 // Terminate analysis
700 //
701     if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
702 }
703
704
705 void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,Int_t &nGenJets,
706                                                 AliAODJet *recJets,Int_t &nRecJets,
707                                                 Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
708
709   //
710   // Relate the two input jet Arrays
711   //
712
713   //
714   // The association has to be unique
715   // So check in two directions
716   // find the closest rec to a gen
717   // and check if there is no other rec which is closer
718   // Caveat: Close low energy/split jets may disturb this correlation
719
720   // Idea: search in two directions generated e.g (a--e) and rec (1--3)
721   // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
722   // in the end we have something like this
723   //    1   2   3
724   // ------------
725   // a| 3   2   0
726   // b| 0   1   0
727   // c| 0   0   3
728   // d| 0   0   1
729   // e| 0   0   1
730   // Topology
731   //   1     2
732   //     a         b        
733   //
734   //  d      c
735   //        3     e
736   // Only entries with "3" match from both sides
737   
738   const int kMode = 3;
739
740   Int_t iFlag[kMaxJets][kMaxJets];
741
742
743
744   for(int i = 0;i < kMaxJets;++i){
745     iRecIndex[i] = -1;
746     iGenIndex[i] = -1;
747     for(int j = 0;j < kMaxJets;++j)iFlag[i][j] = 0;
748   }
749
750   if(nRecJets==0)return;
751   if(nGenJets==0)return;
752
753   const Float_t maxDist = 1.0;
754   // find the closest distance to the generated
755   for(int ig = 0;ig<nGenJets;++ig){
756     Float_t dist = maxDist;
757     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());
758     for(int ir = 0;ir<nRecJets;++ir){
759       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
760       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());
761       if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
762       if(dR<dist){
763         iRecIndex[ig] = ir;
764         dist = dR;
765       } 
766     }
767     if(iRecIndex[ig]>=0)iFlag[ig][iRecIndex[ig]]+=1;
768     // reset...
769     iRecIndex[ig] = -1;
770   }
771   // other way around
772   for(int ir = 0;ir<nRecJets;++ir){
773     Float_t dist = maxDist;
774     for(int ig = 0;ig<nGenJets;++ig){
775       Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
776       if(dR<dist){
777         iGenIndex[ir] = ig;
778         dist = dR;
779       } 
780     }
781     if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]][ir]+=2;
782     // reset...
783     iGenIndex[ir] = -1;
784   }
785
786   // check for "true" correlations
787
788   if(iDebug>1)Printf(">>>>>> Matrix");
789
790   for(int ig = 0;ig<nGenJets;++ig){
791     for(int ir = 0;ir<nRecJets;++ir){
792       // Print
793       if(iDebug>1)printf("XFL %d ",iFlag[ig][ir]);
794
795       if(kMode==3){
796         // we have a uniqie correlation
797         if(iFlag[ig][ir]==3){
798           iGenIndex[ir] = ig;
799           iRecIndex[ig] = ir;
800         }
801       }
802       else{
803         // we just take the correlation from on side
804         if((iFlag[ig][ir]&2)==2){
805           iGenIndex[ir] = ig;
806         }
807         if((iFlag[ig][ir]&1)==1){
808           iRecIndex[ig] = ir;
809         }
810       }
811     }
812     if(iDebug>1)printf("\n");
813   }
814 }
815
816