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