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