]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/AliAnalysisTaskJetSpectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Correcting the expected string to VZERO, it was V0, for accessing survey objects...
[u/mrichter/AliRoot.git] / PWG4 / AliAnalysisTaskJetSpectrum.cxx
CommitLineData
cfff6259 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"
174cb05f 34#include "AliJetHeader.h"
cfff6259 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"
174cb05f 50#include "AliInputEventHandler.h"
cfff6259 51
52#include "AliAnalysisHelperJetTasks.h"
53
54ClassImp(AliAnalysisTaskJetSpectrum)
55
56AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
174cb05f 57 fJetHeaderRec(0x0),
58 fJetHeaderGen(0x0),
cfff6259 59 fAOD(0x0),
60 fBranchRec("jets"),
61 fConfigRec("ConfigJets.C"),
62 fBranchGen(""),
63 fConfigGen(""),
64 fUseAODInput(kFALSE),
65 fUseExternalWeightOnly(kFALSE),
174cb05f 66 fLimitGenJetEta(kFALSE),
cfff6259 67 fAnalysisType(0),
68 fExternalWeight(1),
69 fh1PtHard(0x0),
70 fh1PtHard_NoW(0x0),
71 fh1PtHard_Trials(0x0),
cfff6259 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;
174cb05f 80 fh3PtRecGenHard[ij] = fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] = fh3GenEtaPhiPt[ij] = 0;
cfff6259 81 }
82
83}
84
85AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(const char* name):
86 AliAnalysisTaskSE(name),
174cb05f 87 fJetHeaderRec(0x0),
88 fJetHeaderGen(0x0),
cfff6259 89 fAOD(0x0),
90 fBranchRec("jets"),
91 fConfigRec("ConfigJets.C"),
92 fBranchGen(""),
93 fConfigGen(""),
94 fUseAODInput(kFALSE),
95 fUseExternalWeightOnly(kFALSE),
174cb05f 96 fLimitGenJetEta(kFALSE),
cfff6259 97 fAnalysisType(0),
98 fExternalWeight(1),
99 fh1PtHard(0x0),
100 fh1PtHard_NoW(0x0),
101 fh1PtHard_Trials(0x0),
cfff6259 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
174cb05f 111 fh3PtRecGenHard[ij] = fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] = fh3GenEtaPhiPt[ij] = 0;
cfff6259 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;
174cb05f 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 }
cfff6259 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;
174cb05f 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 }
cfff6259 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){
174cb05f 185 binLimitsEta[iEta] = -2.2;
cfff6259 186 }
187 else{
174cb05f 188 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.2;
cfff6259 189 }
190 }
191
174cb05f 192 const Int_t nBinPhi = 90;
cfff6259 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
cfff6259 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
174cb05f 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)",
cfff6259 256 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
257
258
259
174cb05f 260 fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
cfff6259 261 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
262
263 }
264
174cb05f 265 const Int_t saveLevel = 2; // large save level more histos
cfff6259 266
267 if(saveLevel>0){
268 fHistList->Add(fh1PtHard);
269 fHistList->Add(fh1PtHard_NoW);
270 fHistList->Add(fh1PtHard_Trials);
cfff6259 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]);
174cb05f 283 fHistList->Add(fh3GenEtaPhiPt_NoFound[ij]);
284 fHistList->Add(fh3GenEtaPhiPt[ij]);
cfff6259 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
174cb05f 319
320
cfff6259 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
174cb05f 386
387 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
388
cfff6259 389 if(!fUseExternalWeightOnly){
390 // case were we combine more than one p_T hard bin...
cfff6259 391 }
392
393 // fetch the pythia generated jets only to be used here
394 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
395 AliAODJet pythiaGenJets[kMaxJets];
174cb05f 396 Int_t iCount = 0;
cfff6259 397 for(int ip = 0;ip < nPythiaGenJets;++ip){
174cb05f 398 if(iCount>=kMaxJets)continue;
cfff6259 399 Float_t p[4];
400 pythiaGenHeader->TriggerJet(ip,p);
174cb05f 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
cfff6259 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
174cb05f 412 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
cfff6259 413 }
174cb05f 414 iCount++;
cfff6259 415 }
174cb05f 416 if(fBranchGen.Length()==0)nGenJets = iCount;
cfff6259 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
cfff6259 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){
174cb05f 429 Int_t iCount = 0;
430 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
431 if(iCount>=kMaxJets)continue;
cfff6259 432 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
433 if(!tmp)continue;
174cb05f 434 if(fLimitGenJetEta){
435 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
436 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
437 }
438 genJets[iCount] = *tmp;
439 iCount++;
cfff6259 440 }
174cb05f 441 nGenJets = iCount;
cfff6259 442 }
443 else{
174cb05f 444 Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
cfff6259 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();
174cb05f 456 fh1NRecJets->Fill(nRecJets);
457 nRecJets = TMath::Min(nRecJets,kMaxJets);
458
cfff6259 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
174cb05f 465
cfff6259 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
174cb05f 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
cfff6259 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 }
174cb05f 506 else{
507 fh3RecEtaPhiPt_NoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
508 }
cfff6259 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
174cb05f 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);
cfff6259 520 fh1PtGenIn[ig]->Fill(ptGen,eventW);
521 Int_t ir = iRecIndex[ig];
174cb05f 522 if(ir>=0&&ir<nRecJets){
cfff6259 523 fh1PtGenOut[ig]->Fill(ptGen,eventW);
524 }
174cb05f 525 else{
526 fh3GenEtaPhiPt_NoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
527 }
cfff6259 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
174cb05f 574
575 const int kMode = 3;
cfff6259 576
577 Int_t iFlag[kMaxJets][kMaxJets];
578
174cb05f 579
580
cfff6259 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
174cb05f 590 const Float_t maxDist = 1.0;
591 // find the closest distance to the generated
cfff6259 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;
174cb05f 605 // reset...
606 iRecIndex[ig] = -1;
cfff6259 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;
174cb05f 619 // reset...
620 iGenIndex[ir] = -1;
cfff6259 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
174cb05f 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 }
cfff6259 638 }
639 else{
174cb05f 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 }
cfff6259 647 }
648 }
649 if(iDebug>1)printf("\n");
650 }
651}
652
653