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