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