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