]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fix
[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
de6f8090 263 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
f3050824 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
de6f8090 319 fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};#phi_{gen}",
22ceb537 320 nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi);
f3050824 321
de6f8090 322 fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};#eta_{gen}",
22ceb537 323 nBinEta,binLimitsEta,nBinEta,binLimitsEta);
f3050824 324
22ceb537 325
de6f8090 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}",
77bbf113 327 nBinPt,binLimitsPt,100,-1.0,1.0);
de6f8090 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}",
77bbf113 329 nBinPt,binLimitsPt,100,-1.0,1.0);
22ceb537 330
331
de6f8090 332 fh2PtRecDeltaR[ij] = new TH2F(Form("fh2PtRecDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,rec,j};#Delta R",
d40dc1ba 333 nBinPt,binLimitsPt,60,0,6.0);
de6f8090 334 fh2PtGenDeltaR[ij] = new TH2F(Form("fh2PtGenDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,gen,j};#Delta R",
d40dc1ba 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
de6f8090 366 fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,} (GeV/c)",
f3050824 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();
de6f8090 535 int iProcessType = pythiaGenHeader->ProcessType();
536 // 11 f+f -> f+f
537 // 12 f+barf -> f+barf
538 // 13 f+barf -> g+g
539 // 28 f+g -> f+g
540 // 53 g+g -> f+barf
541 // 68 g+g -> g+g
542 /*
543 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
544 // if(iProcessType != 13 && iProcessType != 68){ // allow only glue
545 if(iProcessType != 11 && iProcessType != 12 && iProcessType != 53){ // allow only quark
546 // if(iProcessType != 28){ // allow only -> f+g
547 PostData(1, fHistList);
548 return;
549 }
550 */
551 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
f3050824 552
553 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
554
555 if(!fUseExternalWeightOnly){
556 // case were we combine more than one p_T hard bin...
557 }
558
559 // fetch the pythia generated jets only to be used here
560 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
561 AliAODJet pythiaGenJets[kMaxJets];
562 Int_t iCount = 0;
563 for(int ip = 0;ip < nPythiaGenJets;++ip){
564 if(iCount>=kMaxJets)continue;
565 Float_t p[4];
566 pythiaGenHeader->TriggerJet(ip,p);
567 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
568
569 if(fLimitGenJetEta){
570 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
571 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
572 }
573
574
575 if(fBranchGen.Length()==0){
576 // if we have MC particles and we do not read from the aod branch
577 // use the pythia jets
578 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
579 }
580 iCount++;
581 }
582 if(fBranchGen.Length()==0)nGenJets = iCount;
583
584 }// (fAnalysisType&kMC)==kMC)
585
586 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
587 fh1PtHard->Fill(ptHard,eventW);
77bbf113 588 fh1PtHardNoW->Fill(ptHard,1);
589 fh1PtHardTrials->Fill(ptHard,nTrials);
f3050824 590
591 // If we set a second branch for the input jets fetch this
592 if(fBranchGen.Length()>0){
593 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
594 if(aodGenJets){
595 Int_t iCount = 0;
596 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
597 if(iCount>=kMaxJets)continue;
598 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
599 if(!tmp)continue;
600 if(fLimitGenJetEta){
601 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
602 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
603 }
604 genJets[iCount] = *tmp;
605 iCount++;
606 }
607 nGenJets = iCount;
608 }
609 else{
610 Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
611 }
612 }
613
614 fh1NGenJets->Fill(nGenJets);
615 // We do not want to exceed the maximum jet number
616 nGenJets = TMath::Min(nGenJets,kMaxJets);
617
618 // Fetch the reconstructed jets...
619
620
621 nRecJets = aodRecJets->GetEntries();
622 fh1NRecJets->Fill(nRecJets);
623 nRecJets = TMath::Min(nRecJets,kMaxJets);
df65bddb 624 //////////////////////////////////////////
625 nTracks = fAOD->GetNumberOfTracks();
626 ///////////////////////////////////////////
f3050824 627
628 for(int ir = 0;ir < nRecJets;++ir){
629 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
630 if(!tmp)continue;
631 recJets[ir] = *tmp;
632 }
633
634
635 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
636 // Relate the jets
637 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
638 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
639
640 for(int i = 0;i<kMaxJets;++i){
641 iGenIndex[i] = iRecIndex[i] = -1;
642 }
643
644
645 GetClosestJets(genJets,nGenJets,recJets,nRecJets,
646 iGenIndex,iRecIndex,fDebug);
647 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
648
649 if(fDebug){
650 for(int i = 0;i<kMaxJets;++i){
651 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
652 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
653 }
654 }
655
656 // loop over reconstructed jets
657 for(int ir = 0;ir < nRecJets;++ir){
658 Double_t ptRec = recJets[ir].Pt();
659 Double_t phiRec = recJets[ir].Phi();
660 if(phiRec<0)phiRec+=TMath::Pi()*2.;
661 Double_t etaRec = recJets[ir].Eta();
77bbf113 662 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
f3050824 663 fh1E[ir]->Fill(recJets[ir].E(),eventW);
664 fh1PtRecIn[ir]->Fill(ptRec,eventW);
665 fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
d40dc1ba 666 for(int irr = ir+1;irr<nRecJets;irr++){
667 fh2PtRecDeltaR[ir]->Fill(recJets[irr].Pt(),recJets[ir].DeltaR(&recJets[irr]));
668 }
f3050824 669 // Fill Correlation
670 Int_t ig = iGenIndex[ir];
df65bddb 671 if(ig>=0 && ig<nGenJets){
77bbf113 672 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
673 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
f3050824 674 fh1PtRecOut[ir]->Fill(ptRec,eventW);
df65bddb 675 Double_t ptGen = genJets[ig].Pt();
676 Double_t phiGen = genJets[ig].Phi();
677 if(phiGen<0)phiGen+=TMath::Pi()*2.;
678 Double_t etaGen = genJets[ig].Eta();
d40dc1ba 679
680 //
681 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
682 //
683
684 if(TMath::Abs(etaRec)<fRecEtaWindow){
685
f3050824 686 fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
22ceb537 687 fh2PhiFGen[ir]->Fill(phiRec,phiGen,eventW);
688 fh2EtaFGen[ir]->Fill(etaRec,etaGen,eventW);
77bbf113 689 fh2PtGenDeltaEta[ir]->Fill(ptGen,etaGen-etaRec,eventW);
77bbf113 690 fh2PtGenDeltaPhi[ir]->Fill(ptGen,phiGen-phiRec,eventW);
f3050824 691 fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
77bbf113 692 fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1);
df65bddb 693 /////////////////////////////////////////////////////
de6f8090 694
695 // Double_t eRec = recJets[ir].E();
696 // Double_t eGen = genJets[ig].E();
697 // CKB use p_T not Energy
698 // TODO recname variabeles and histos
77bbf113 699 Double_t eRec = recJets[ir].E();
700 Double_t eGen = genJets[ig].E();
d40dc1ba 701
77bbf113 702 fh2Efficiency->Fill(eGen, eRec/eGen);
77bbf113 703
704 if (eGen>=0. && eGen<=250.){
705 Double_t eLeading = -1;
df65bddb 706 Double_t ptleading = -1;
707 Int_t nPart=0;
708 // loop over tracks
709 for (Int_t it = 0; it< nTracks; it++){
d40dc1ba 710 // if (fAOD->GetTrack(it)->E() > eGen) continue; // CKB. Not allowed! cannot cut on gen properties in real events!
df65bddb 711 // find leading particle
de6f8090 712 // if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){
713 // TODO implement esd filter flag to be the same as in the jet finder
714 // allow also for MC particles...
715 Float_t r = recJets[ir].DeltaR(fAOD->GetTrack(it));
716 if (r<0.4 && fAOD->GetTrack(it)->Pt()>ptleading){
77bbf113 717 eLeading = fAOD->GetTrack(it)->E();
df65bddb 718 ptleading = fAOD->GetTrack(it)->Pt();
719 }
d40dc1ba 720 // if (fAOD->GetTrack(it)->Pt()>0.03*eGen && fAOD->GetTrack(it)->E()<=eGen && r<0.7) // CKB cannot cut on gen properties
de6f8090 721 if (fAOD->GetTrack(it)->Pt()>0.03*eRec && fAOD->GetTrack(it)->Pt()<=eRec && r<0.7)
df65bddb 722 nPart++;
723 }
77bbf113 724 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
725
df65bddb 726 // fill Response Map (4D histogram) and Energy vs z distributions
77bbf113 727 Double_t var[4] = {eGen, eRec, ptleading/eGen, ptleading/eRec};
df65bddb 728 fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
729 fh2EGenZGen->Fill(var[0],var[2]);
730 fh1JetMultiplicity->Fill(nPart);
77bbf113 731 fh3EGenERecN->Fill(eGen, eRec, nPart);
df65bddb 732 for(int ic = 0;ic <kMaxCorrelation;ic++){
d40dc1ba 733 if (nPart<=fgkJetNpartCut[ic]){ // is this corrected for CKB
77bbf113 734 fhnCorrelation[ic]->Fill(var);
735 break;
736 }
df65bddb 737 }
738 }
d40dc1ba 739
740 }// if etarec in window
741
742 }
df65bddb 743 ////////////////////////////////////////////////////
f3050824 744 else{
77bbf113 745 fh3RecEtaPhiPtNoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
f3050824 746 }
747 }// loop over reconstructed jets
748
749
750 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
751 for(int ig = 0;ig < nGenJets;++ig){
752 Double_t ptGen = genJets[ig].Pt();
753 // Fill Correlation
754 Double_t phiGen = genJets[ig].Phi();
755 if(phiGen<0)phiGen+=TMath::Pi()*2.;
756 Double_t etaGen = genJets[ig].Eta();
757 fh3GenEtaPhiPt[ig]->Fill(etaGen,phiGen,ptGen,eventW);
758 fh1PtGenIn[ig]->Fill(ptGen,eventW);
d40dc1ba 759 for(int igg = ig+1;igg<nGenJets;igg++){
760 fh2PtGenDeltaR[ig]->Fill(genJets[igg].Pt(),genJets[ig].DeltaR(&genJets[igg]));
761 }
f3050824 762 Int_t ir = iRecIndex[ig];
763 if(ir>=0&&ir<nRecJets){
764 fh1PtGenOut[ig]->Fill(ptGen,eventW);
765 }
766 else{
77bbf113 767 fh3GenEtaPhiPtNoFound[ig]->Fill(etaGen,phiGen,ptGen,eventW);
f3050824 768 }
769 }// loop over reconstructed jets
770
771 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
772 PostData(1, fHistList);
773}
774
775void AliAnalysisTaskJetSpectrum::Terminate(Option_t */*option*/)
776{
777// Terminate analysis
778//
779 if (fDebug > 1) printf("AnalysisJetSpectrum: Terminate() \n");
780}
781
782
77bbf113 783void AliAnalysisTaskJetSpectrum::GetClosestJets(AliAODJet *genJets,const Int_t &nGenJets,
784 AliAODJet *recJets,const Int_t &nRecJets,
f3050824 785 Int_t *iGenIndex,Int_t *iRecIndex,Int_t iDebug){
786
787 //
788 // Relate the two input jet Arrays
789 //
790
791 //
792 // The association has to be unique
793 // So check in two directions
794 // find the closest rec to a gen
795 // and check if there is no other rec which is closer
796 // Caveat: Close low energy/split jets may disturb this correlation
797
798 // Idea: search in two directions generated e.g (a--e) and rec (1--3)
799 // Fill a matrix with Flags (1 for closest rec jet, 2 for closest rec jet
800 // in the end we have something like this
801 // 1 2 3
802 // ------------
803 // a| 3 2 0
804 // b| 0 1 0
805 // c| 0 0 3
806 // d| 0 0 1
807 // e| 0 0 1
808 // Topology
809 // 1 2
810 // a b
811 //
812 // d c
813 // 3 e
814 // Only entries with "3" match from both sides
815
816 const int kMode = 3;
817
818 Int_t iFlag[kMaxJets][kMaxJets];
819
820
821
822 for(int i = 0;i < kMaxJets;++i){
823 iRecIndex[i] = -1;
824 iGenIndex[i] = -1;
825 for(int j = 0;j < kMaxJets;++j)iFlag[i][j] = 0;
826 }
827
828 if(nRecJets==0)return;
829 if(nGenJets==0)return;
830
d40dc1ba 831 const Float_t maxDist = 0.5;
f3050824 832 // find the closest distance to the generated
833 for(int ig = 0;ig<nGenJets;++ig){
834 Float_t dist = maxDist;
835 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());
836 for(int ir = 0;ir<nRecJets;++ir){
837 Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
838 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());
839 if(iDebug>1)Printf("Distance (%d)--(%d) %3.3f ",ig,ir,dR);
840 if(dR<dist){
841 iRecIndex[ig] = ir;
842 dist = dR;
843 }
844 }
845 if(iRecIndex[ig]>=0)iFlag[ig][iRecIndex[ig]]+=1;
846 // reset...
847 iRecIndex[ig] = -1;
848 }
849 // other way around
850 for(int ir = 0;ir<nRecJets;++ir){
851 Float_t dist = maxDist;
852 for(int ig = 0;ig<nGenJets;++ig){
853 Double_t dR = genJets[ig].DeltaR(&recJets[ir]);
854 if(dR<dist){
855 iGenIndex[ir] = ig;
856 dist = dR;
857 }
858 }
859 if(iGenIndex[ir]>=0)iFlag[iGenIndex[ir]][ir]+=2;
860 // reset...
861 iGenIndex[ir] = -1;
862 }
863
864 // check for "true" correlations
865
866 if(iDebug>1)Printf(">>>>>> Matrix");
867
868 for(int ig = 0;ig<nGenJets;++ig){
869 for(int ir = 0;ir<nRecJets;++ir){
870 // Print
871 if(iDebug>1)printf("XFL %d ",iFlag[ig][ir]);
872
873 if(kMode==3){
874 // we have a uniqie correlation
875 if(iFlag[ig][ir]==3){
876 iGenIndex[ir] = ig;
877 iRecIndex[ig] = ir;
878 }
879 }
880 else{
881 // we just take the correlation from on side
882 if((iFlag[ig][ir]&2)==2){
883 iGenIndex[ir] = ig;
884 }
885 if((iFlag[ig][ir]&1)==1){
886 iRecIndex[ig] = ir;
887 }
888 }
889 }
890 if(iDebug>1)printf("\n");
891 }
892}
893
894