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