Fix minor bug in SDD raw data writing (F. Prino)
[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;
22ceb537 91 fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtGenDeltaPhi[ij] = fh2PtGenDeltaEta[ij] = 0;
df65bddb 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;
22ceb537 133 fh2PtFGen[ij] = fh2PhiFGen[ij] = fh2EtaFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = fh2PtGenDeltaPhi[ij] = fh2PtGenDeltaEta[ij] = 0;
f3050824 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
22ceb537 266 const Int_t nBinPhi = 30;
f3050824 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
22ceb537 308 fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};phi_{gen}",
309 nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi);
f3050824 310
22ceb537 311 fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};eta_{gen}",
312 nBinEta,binLimitsEta,nBinEta,binLimitsEta);
f3050824 313
22ceb537 314
315 fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaPhi_j%d",ij),"delta phi vs. P_tgen;p_{T,gen} (GeV/c);#phi_{gen}-phi_{rec}",
316 100,-1.0,1.0,nBinPt,binLimitsPt);
317 fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaEta_j%d",ij),"delta eta vs. P_tgen;p_{T,gen} (GeV/c);#eta_{gen}-eta_{rec}",
318 100,-1.0,1.0,nBinPt,binLimitsPt);
319
320
f3050824 321
322
323
324 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);
325
326
327
328 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);
329
330
331 fh2Frag[ij] = new TH2F(Form("fh2Frag_j%d",ij),"Jet Fragmentation;x=E_{i}/E_{jet};E_{jet};1/N_{jet}dN_{ch}/dx",
332 nBinFrag,0.,1.,nBinPt,binLimitsPt);
333
334 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",
335 nBinFrag,0.,10.,nBinPt,binLimitsPt);
336
337 fh3RecEtaPhiPt[ij] = new TH3F(Form("fh3RecEtaPhiPt_j%d",ij),"Rec eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
338 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
339
340
341
342 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)",
343 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
344
345
22ceb537 346 fh3GenEtaPhiPt_NoFound[ij] = new TH3F(Form("fh3GenEtaPhiPt_NoFound_j%d",ij),"No found for generated jet eta, phi, pt; #eta; #phi; p_{T,gen} (GeV/c)",
f3050824 347 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
348
349
350
351 fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)",
352 nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
353
354 }
df65bddb 355
356 /////////////////////////////////////////////////////////////////
357 fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
f3050824 358
df65bddb 359 fh2ERecZRec = new TH2F("fh2ERecZRec", " ; E^{jet}_{rec} [GeV]; z^{lp}_{rec}", 100, 0., 250., 100, 0., 2.);
360 fh2EGenZGen = new TH2F("fh2EGenZGen", " ; E^{jet}_{gen} [GeV]; z^{lp}_{gen}", 100, 0., 250., 100, 0., 2.);
361 fh2Efficiency = new TH2F("fh2Efficiency", "ERec/EGen;E^{jet}_{gen} [GeV];E^{jet}_{rec}/E^{jet}_{gen}", 100, 0., 250., 100, 0., 1.5);
f3050824 362
22ceb537 363 fh3EGenERecN = new TH3F("fh3EGenERecN", "Efficiency vs. Jet Multiplicity", 100, 0., 250., 100, 0., 250., 51, 0., 50.);
df65bddb 364
365 // Response map
366 //arrays for bin limits
367 const Int_t nbin[4] = {100, 100, 100, 100};
368 Double_t LowEdge[4] = {0.,0.,0.,0.};
369 Double_t UpEdge[4] = {250., 250., 1., 1.};
370
371 for(int ic = 0;ic < kMaxCorrelation;ic++){
372 fhnCorrelation[ic] = new THnSparseF(Form("fhnCorrelation_%d",ic), "Response Map", 4, nbin, LowEdge, UpEdge);
373 if(ic==0) fhnCorrelation[ic]->SetTitle(Form("ResponseMap 0 <= npart <= %.0E",fgkJetNpartCut[ic]));
374 else fhnCorrelation[ic]->SetTitle(Form("ResponseMap %.0E < npart <= %.0E",fgkJetNpartCut[ic-1],fgkJetNpartCut[ic]));
375 }
22ceb537 376 const Int_t saveLevel = 3; // large save level more histos
f3050824 377 if(saveLevel>0){
4dbfdecc 378 fHistList->Add(fh1Xsec);
379 fHistList->Add(fh1Trials);
f3050824 380 fHistList->Add(fh1PtHard);
381 fHistList->Add(fh1PtHard_NoW);
382 fHistList->Add(fh1PtHard_Trials);
383 fHistList->Add(fh1NGenJets);
384 fHistList->Add(fh1NRecJets);
df65bddb 385 ////////////////////////
386 fHistList->Add(fh1JetMultiplicity);
387 fHistList->Add(fh2ERecZRec);
388 fHistList->Add(fh2EGenZGen);
389 fHistList->Add(fh2Efficiency);
390 fHistList->Add(fh3EGenERecN);
391
392 for(int ic = 0;ic < kMaxCorrelation;++ic){
393 fHistList->Add(fhnCorrelation[ic]);
394 }
395 ////////////////////////
f3050824 396 for(int ij = 0;ij<kMaxJets;++ij){
397 fHistList->Add(fh1E[ij]);
398 fHistList->Add(fh1PtRecIn[ij]);
399 fHistList->Add(fh1PtRecOut[ij]);
400 fHistList->Add(fh1PtGenIn[ij]);
401 fHistList->Add(fh1PtGenOut[ij]);
402 fHistList->Add(fh2PtFGen[ij]);
22ceb537 403 fHistList->Add(fh2PhiFGen[ij]);
404 fHistList->Add(fh2EtaFGen[ij]);
405 fHistList->Add(fh2PtGenDeltaEta[ij]);
406 fHistList->Add(fh2PtGenDeltaPhi[ij]);
407 fHistList->Add(fh3RecEtaPhiPt[ij]);
408 fHistList->Add(fh3GenEtaPhiPt[ij]);
f3050824 409 if(saveLevel>2){
f3050824 410 fHistList->Add(fh3RecEtaPhiPt_NoGen[ij]);
411 fHistList->Add(fh3GenEtaPhiPt_NoFound[ij]);
f3050824 412 }
413 }
414 }
415
416 // =========== Switch on Sumw2 for all histos ===========
417 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
418 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
419 if (h1){
420 // Printf("%s ",h1->GetName());
421 h1->Sumw2();
4dbfdecc 422 continue;
f3050824 423 }
4dbfdecc 424 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
425 if(hn)hn->Sumw2();
f3050824 426 }
427
428 TH1::AddDirectory(oldStatus);
429
430}
431
432void AliAnalysisTaskJetSpectrum::Init()
433{
434 //
435 // Initialization
436 //
437
438 Printf(">>> AnalysisTaskJetSpectrum::Init() debug level %d\n",fDebug);
439 if (fDebug > 1) printf("AnalysisTaskJetSpectrum::Init() \n");
440
441}
442
443void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/)
444{
445 //
446 // Execute analysis for current event
447 //
448
449
450
451 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
452
453
454 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
455
456 if(!aodH){
457 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
458 return;
459 }
460
461
462 // aodH->SelectEvent(kTRUE);
463
464 // ========= These pointers need to be valid in any case =======
465
466
467 /*
468 AliUA1JetHeaderV1 *jhRec = dynamic_cast<AliUA1JetHeaderV1*>(fJetFinderRec->GetHeader());
469 if(!jhRec){
470 Printf("%s:%d No Jet Header found",(char*)__FILE__,__LINE__);
471 return;
472 }
473 */
474 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
475 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
476 if(!aodRecJets){
477 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
478 return;
479 }
480
481 // ==== General variables needed
482
483
484 // We use statice array, not to fragment the memory
485 AliAODJet genJets[kMaxJets];
486 Int_t nGenJets = 0;
487 AliAODJet recJets[kMaxJets];
488 Int_t nRecJets = 0;
df65bddb 489 ///////////////////////////
490 Int_t nTracks = 0;
491 //////////////////////////
f3050824 492
493 Double_t eventW = 1;
494 Double_t ptHard = 0;
495 Double_t nTrials = 1; // Trials for MC trigger weigth for real data
496
497 if(fUseExternalWeightOnly){
498 eventW = fExternalWeight;
499 }
500
501
502 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
503 if((fAnalysisType&kAnaMC)==kAnaMC){
504 // this is the part we only use when we have MC information
505 AliMCEvent* mcEvent = MCEvent();
506 // AliStack *pStack = 0;
507 if(!mcEvent){
508 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
509 return;
510 }
511 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
512 if(!pythiaGenHeader){
513 return;
514 }
515
516 nTrials = pythiaGenHeader->Trials();
517 ptHard = pythiaGenHeader->GetPtHard();
518
519
520 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
521
522 if(!fUseExternalWeightOnly){
523 // case were we combine more than one p_T hard bin...
524 }
525
526 // fetch the pythia generated jets only to be used here
527 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
528 AliAODJet pythiaGenJets[kMaxJets];
529 Int_t iCount = 0;
530 for(int ip = 0;ip < nPythiaGenJets;++ip){
531 if(iCount>=kMaxJets)continue;
532 Float_t p[4];
533 pythiaGenHeader->TriggerJet(ip,p);
534 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
535
536 if(fLimitGenJetEta){
537 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
538 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
539 }
540
541
542 if(fBranchGen.Length()==0){
543 // if we have MC particles and we do not read from the aod branch
544 // use the pythia jets
545 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
546 }
547 iCount++;
548 }
549 if(fBranchGen.Length()==0)nGenJets = iCount;
550
551 }// (fAnalysisType&kMC)==kMC)
552
553 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
554 fh1PtHard->Fill(ptHard,eventW);
555 fh1PtHard_NoW->Fill(ptHard,1);
556 fh1PtHard_Trials->Fill(ptHard,nTrials);
557
558 // If we set a second branch for the input jets fetch this
559 if(fBranchGen.Length()>0){
560 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
561 if(aodGenJets){
562 Int_t iCount = 0;
563 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
564 if(iCount>=kMaxJets)continue;
565 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
566 if(!tmp)continue;
567 if(fLimitGenJetEta){
568 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
569 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
570 }
571 genJets[iCount] = *tmp;
572 iCount++;
573 }
574 nGenJets = iCount;
575 }
576 else{
577 Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
578 }
579 }
580
581 fh1NGenJets->Fill(nGenJets);
582 // We do not want to exceed the maximum jet number
583 nGenJets = TMath::Min(nGenJets,kMaxJets);
584
585 // Fetch the reconstructed jets...
586
587
588 nRecJets = aodRecJets->GetEntries();
589 fh1NRecJets->Fill(nRecJets);
590 nRecJets = TMath::Min(nRecJets,kMaxJets);
df65bddb 591 //////////////////////////////////////////
592 nTracks = fAOD->GetNumberOfTracks();
593 ///////////////////////////////////////////
f3050824 594
595 for(int ir = 0;ir < nRecJets;++ir){
596 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
597 if(!tmp)continue;
598 recJets[ir] = *tmp;
599 }
600
601
602 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
603 // Relate the jets
604 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
605 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
606
607 for(int i = 0;i<kMaxJets;++i){
608 iGenIndex[i] = iRecIndex[i] = -1;
609 }
610
611
612 GetClosestJets(genJets,nGenJets,recJets,nRecJets,
613 iGenIndex,iRecIndex,fDebug);
614 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
615
616 if(fDebug){
617 for(int i = 0;i<kMaxJets;++i){
618 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
619 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
620 }
621 }
622
623 // loop over reconstructed jets
624 for(int ir = 0;ir < nRecJets;++ir){
625 Double_t ptRec = recJets[ir].Pt();
626 Double_t phiRec = recJets[ir].Phi();
627 if(phiRec<0)phiRec+=TMath::Pi()*2.;
628 Double_t etaRec = recJets[ir].Eta();
629
630 fh1E[ir]->Fill(recJets[ir].E(),eventW);
631 fh1PtRecIn[ir]->Fill(ptRec,eventW);
632 fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
633 // Fill Correlation
634 Int_t ig = iGenIndex[ir];
df65bddb 635 if(ig>=0 && ig<nGenJets){
f3050824 636 fh1PtRecOut[ir]->Fill(ptRec,eventW);
df65bddb 637 Double_t ptGen = genJets[ig].Pt();
638 Double_t phiGen = genJets[ig].Phi();
639 if(phiGen<0)phiGen+=TMath::Pi()*2.;
640 Double_t etaGen = genJets[ig].Eta();
f3050824 641 fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
22ceb537 642 fh2PhiFGen[ir]->Fill(phiRec,phiGen,eventW);
643 fh2EtaFGen[ir]->Fill(etaRec,etaGen,eventW);
644 fh2PtGenDeltaEta[ir]->Fill(etaGen-etaRec,eventW);
645 fh2PtGenDeltaEta[ir]->Fill(phiGen-phiRec,eventW);
f3050824 646 fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
647 fh3PtRecGenHard_NoW[ir]->Fill(ptRec,ptGen,ptHard,1);
df65bddb 648 /////////////////////////////////////////////////////
649 Double_t ERec = recJets[ir].E();
650 Double_t EGen = genJets[ig].E();
651 fh2Efficiency->Fill(EGen, ERec/EGen);
652 if (EGen>=0. && EGen<=250.){
653 Double_t Eleading = -1;
654 Double_t ptleading = -1;
655 Int_t nPart=0;
656 // loop over tracks
657 for (Int_t it = 0; it< nTracks; it++){
658 if (fAOD->GetTrack(it)->E() > EGen) continue;
659 Double_t phiTrack = fAOD->GetTrack(it)->Phi();
660 if (phiTrack<0) phiTrack+=TMath::Pi()*2.;
661 Double_t etaTrack = fAOD->GetTrack(it)->Eta();
662 Float_t deta = etaRec - etaTrack;
663 Float_t dphi = TMath::Abs(phiRec - phiTrack);
664 Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
665 // find leading particle
666 if (R<0.4 && fAOD->GetTrack(it)->E()>Eleading){
667 Eleading = fAOD->GetTrack(it)->E();
668 ptleading = fAOD->GetTrack(it)->Pt();
669 }
670 if (fAOD->GetTrack(it)->Pt()>0.03*EGen && fAOD->GetTrack(it)->E()<=EGen && R<0.7)
671 nPart++;
672 }
673 // fill Response Map (4D histogram) and Energy vs z distributions
674 Double_t var[4] = {EGen, ERec, ptleading/EGen, ptleading/ERec};
675 fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
676 fh2EGenZGen->Fill(var[0],var[2]);
677 fh1JetMultiplicity->Fill(nPart);
678 fh3EGenERecN->Fill(EGen, ERec, nPart);
679 for(int ic = 0;ic <kMaxCorrelation;ic++){
680 if (nPart<=fgkJetNpartCut[ic]){
681 fhnCorrelation[ic]->Fill(var);
682 break;
683 }
684 }
685 }
686 }
687 ////////////////////////////////////////////////////
f3050824 688 else{
689 fh3RecEtaPhiPt_NoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
690 }
691 }// loop over reconstructed jets
692
693
694 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
695 for(int ig = 0;ig < nGenJets;++ig){
696 Double_t ptGen = genJets[ig].Pt();
697 // Fill Correlation
698 Double_t phiGen = genJets[ig].Phi();
699 if(phiGen<0)phiGen+=TMath::Pi()*2.;
700 Double_t etaGen = genJets[ig].Eta();
701 fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
702 fh1PtGenIn[ig]->Fill(ptGen,eventW);
703 Int_t ir = iRecIndex[ig];
704 if(ir>=0&&ir<nRecJets){
705 fh1PtGenOut[ig]->Fill(ptGen,eventW);
706 }
707 else{
708 fh3GenEtaPhiPt_NoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
709 }
710 }// loop over reconstructed jets
711
712 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
713 PostData(1, fHistList);
714}
715
716void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
717{
718// Terminate analysis
719//
720 if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
721}
722
723
724void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,Int_t &nGenJets,
725 AliAODJet *recJets,Int_t &nRecJets,
726 Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
727
728 //
729 // Relate the two input jet Arrays
730 //
731
732 //
733 // The association has to be unique
734 // So check in two directions
735 // find the closest rec to a gen
736 // and check if there is no other rec which is closer
737 // Caveat: Close low energy/split jets may disturb this correlation
738
739 // Idea: search in two directions generated e.g (a--e) and rec (1--3)
740 // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
741 // in the end we have something like this
742 // 1 2 3
743 // ------------
744 // a| 3 2 0
745 // b| 0 1 0
746 // c| 0 0 3
747 // d| 0 0 1
748 // e| 0 0 1
749 // Topology
750 // 1 2
751 // a b
752 //
753 // d c
754 // 3 e
755 // Only entries with "3" match from both sides
756
757 const int kMode = 3;
758
759 Int_t iFlag[kMaxJets][kMaxJets];
760
761
762
763 for(int i = 0;i < kMaxJets;++i){
764 iRecIndex[i] = -1;
765 iGenIndex[i] = -1;
766 for(int j = 0;j < kMaxJets;++j)iFlag[i][j] = 0;
767 }
768
769 if(nRecJets==0)return;
770 if(nGenJets==0)return;
771
772 const Float_t maxDist = 1.0;
773 // find the closest distance to the generated
774 for(int ig = 0;ig<nGenJets;++ig){
775 Float_t dist = maxDist;
776 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());
777 for(int ir = 0;ir<nRecJets;++ir){
778 Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
779 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());
780 if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
781 if(dR<dist){
782 iRecIndex[ig] = ir;
783 dist = dR;
784 }
785 }
786 if(iRecIndex[ig]>=0)iFlag[ig][iRecIndex[ig]]+=1;
787 // reset...
788 iRecIndex[ig] = -1;
789 }
790 // other way around
791 for(int ir = 0;ir<nRecJets;++ir){
792 Float_t dist = maxDist;
793 for(int ig = 0;ig<nGenJets;++ig){
794 Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
795 if(dR<dist){
796 iGenIndex[ir] = ig;
797 dist = dR;
798 }
799 }
800 if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]][ir]+=2;
801 // reset...
802 iGenIndex[ir] = -1;
803 }
804
805 // check for "true" correlations
806
807 if(iDebug>1)Printf(">>>>>> Matrix");
808
809 for(int ig = 0;ig<nGenJets;++ig){
810 for(int ir = 0;ir<nRecJets;++ir){
811 // Print
812 if(iDebug>1)printf("XFL %d ",iFlag[ig][ir]);
813
814 if(kMode==3){
815 // we have a uniqie correlation
816 if(iFlag[ig][ir]==3){
817 iGenIndex[ir] = ig;
818 iRecIndex[ig] = ir;
819 }
820 }
821 else{
822 // we just take the correlation from on side
823 if((iFlag[ig][ir]&2)==2){
824 iGenIndex[ir] = ig;
825 }
826 if((iFlag[ig][ir]&1)==1){
827 iRecIndex[ig] = ir;
828 }
829 }
830 }
831 if(iDebug>1)printf("\n");
832 }
833}
834
835