]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
fix to avoid crash when runnign with option 1
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum2.cxx
... / ...
CommitLineData
1// **************************************
2// Task used for the correction of determiantion of reconstructed jet spectra
3// Compares input (gen) and output (rec) jets
4// *******************************************
5
6
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
22
23#include <TROOT.h>
24#include <TRandom.h>
25#include <TSystem.h>
26#include <TInterpreter.h>
27#include <TChain.h>
28#include <TFile.h>
29#include <TKey.h>
30#include <TH1F.h>
31#include <TH2F.h>
32#include <TH3F.h>
33#include <TProfile.h>
34#include <TList.h>
35#include <TLorentzVector.h>
36#include <TClonesArray.h>
37#include <TRefArray.h>
38#include "TDatabasePDG.h"
39
40#include "AliAnalysisTaskJetSpectrum2.h"
41#include "AliAnalysisManager.h"
42#include "AliJetFinder.h"
43#include "AliJetHeader.h"
44#include "AliJetReader.h"
45#include "AliJetReaderHeader.h"
46#include "AliUA1JetHeaderV1.h"
47#include "AliESDEvent.h"
48#include "AliAODEvent.h"
49#include "AliAODHandler.h"
50#include "AliAODTrack.h"
51#include "AliAODJet.h"
52#include "AliAODJetEventBackground.h"
53#include "AliAODMCParticle.h"
54#include "AliMCEventHandler.h"
55#include "AliMCEvent.h"
56#include "AliStack.h"
57#include "AliGenPythiaEventHeader.h"
58#include "AliJetKineReaderHeader.h"
59#include "AliGenCocktailEventHeader.h"
60#include "AliInputEventHandler.h"
61
62
63#include "AliAnalysisHelperJetTasks.h"
64
65ClassImp(AliAnalysisTaskJetSpectrum2)
66
67AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2():
68 AliAnalysisTaskSE(),
69 fJetHeaderRec(0x0),
70 fJetHeaderGen(0x0),
71 fAODIn(0x0),
72 fAODOut(0x0),
73 fAODExtension(0x0),
74 fhnCorrelation(0x0),
75 fhnCorrelationPhiZRec(0x0),
76 f1PtScale(0x0),
77 fBranchRec("jets"),
78 fBranchGen(""),
79 fBranchBkgRec(""),
80 fBranchBkgGen(""),
81 fNonStdFile(""),
82 fUseAODJetInput(kFALSE),
83 fUseAODTrackInput(kFALSE),
84 fUseAODMCInput(kFALSE),
85 fUseGlobalSelection(kFALSE),
86 fUseExternalWeightOnly(kFALSE),
87 fLimitGenJetEta(kFALSE),
88 fNMatchJets(5),
89 fFilterMask(0),
90 fEventSelectionMask(0),
91 fAnalysisType(0),
92 fTrackTypeRec(kTrackUndef),
93 fTrackTypeGen(kTrackUndef),
94 fEventClass(0),
95 fAvgTrials(1),
96 fExternalWeight(1),
97 fJetRecEtaWindow(0.5),
98 fTrackRecEtaWindow(0.5),
99 fMinJetPt(0),
100 fMinTrackPt(0.15),
101 fDeltaPhiWindow(90./180.*TMath::Pi()),
102 fMultRec(0),
103 fMultGen(0),
104 fh1Xsec(0x0),
105 fh1Trials(0x0),
106 fh1PtHard(0x0),
107 fh1PtHardNoW(0x0),
108 fh1PtHardTrials(0x0),
109 fh1ZVtx(0x0),
110 fh1TmpRho(0x0),
111 fh2MultRec(0x0),
112 fh2MultGen(0x0),
113 fh2PtFGen(0x0),
114 fh2RelPtFGen(0x0),
115 fHistList(0x0)
116{
117 for(int i = 0;i < kMaxStep*2;++i){
118 fhnJetContainer[i] = 0;
119 }
120
121 for(int ij = 0;ij <kJetTypes;++ij){
122 fh1NJets[ij] = 0;
123 fh1SumPtTrack[ij] = 0;
124 fh1PtJetsIn[ij] = 0;
125 fh1PtTracksIn[ij] = 0;
126 fh1PtTracksLeadingIn[ij] = 0;
127 fh2MultJetPt[ij] = 0;
128 fh2NJetsPt[ij] = 0;
129 fh2NTracksPt[ij] = 0;
130 fh2LeadingTrackPtTrackPhi[ij] = 0;
131 for(int i = 0;i <= kMaxJets;++i){
132 fh2PhiPt[ij][i] = 0;
133 fh2EtaPt[ij][i] = 0;
134 fh2AreaPt[ij][i] = 0;
135 fh2EtaArea[ij][i] = 0;
136 fh2PhiEta[ij][i] = 0;
137 fh2LTrackPtJetPt[ij][i] = 0;
138 fh1PtIn[ij][i] = 0;
139 fh2RhoPt[ij][i] = 0;
140 fh2PsiPt[ij][i] = 0;
141 }
142
143 fh1DijetMinv[ij] = 0;
144 fh2DijetDeltaPhiPt[ij] = 0;
145 fh2DijetAsymPt[ij] = 0;
146 fh2DijetPt2vsPt1[ij] = 0;
147 fh2DijetDifvsSum[ij] = 0;
148
149 }
150}
151
152AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
153 AliAnalysisTaskSE(name),
154 fJetHeaderRec(0x0),
155 fJetHeaderGen(0x0),
156 fAODIn(0x0),
157 fAODOut(0x0),
158 fAODExtension(0x0),
159 fhnCorrelation(0x0),
160 fhnCorrelationPhiZRec(0x0),
161 f1PtScale(0x0),
162 fBranchRec("jets"),
163 fBranchGen(""),
164 fBranchBkgRec(""),
165 fBranchBkgGen(""),
166 fNonStdFile(""),
167 fUseAODJetInput(kFALSE),
168 fUseAODTrackInput(kFALSE),
169 fUseAODMCInput(kFALSE),
170 fUseGlobalSelection(kFALSE),
171 fUseExternalWeightOnly(kFALSE),
172 fLimitGenJetEta(kFALSE),
173 fNMatchJets(5),
174 fFilterMask(0),
175 fEventSelectionMask(0),
176 fAnalysisType(0),
177 fTrackTypeRec(kTrackUndef),
178 fTrackTypeGen(kTrackUndef),
179 fEventClass(0),
180 fAvgTrials(1),
181 fExternalWeight(1),
182 fJetRecEtaWindow(0.5),
183 fTrackRecEtaWindow(0.5),
184 fMinJetPt(0),
185 fMinTrackPt(0.15),
186 fDeltaPhiWindow(90./180.*TMath::Pi()),
187 fMultRec(0),
188 fMultGen(0),
189 fh1Xsec(0x0),
190 fh1Trials(0x0),
191 fh1PtHard(0x0),
192 fh1PtHardNoW(0x0),
193 fh1PtHardTrials(0x0),
194 fh1ZVtx(0x0),
195 fh1TmpRho(0x0),
196 fh2MultRec(0x0),
197 fh2MultGen(0x0),
198 fh2PtFGen(0x0),
199 fh2RelPtFGen(0x0),
200 fHistList(0x0)
201{
202
203 for(int i = 0;i < kMaxStep*2;++i){
204 fhnJetContainer[i] = 0;
205 }
206
207 for(int ij = 0;ij <kJetTypes;++ij){
208 fh1NJets[ij] = 0;
209 fh1SumPtTrack[ij] = 0;
210 fh1PtJetsIn[ij] = 0;
211 fh1PtTracksIn[ij] = 0;
212 fh1PtTracksLeadingIn[ij] = 0;
213 fh2MultJetPt[ij] = 0;
214 fh2NJetsPt[ij] = 0;
215 fh2NTracksPt[ij] = 0;
216 fh2LeadingTrackPtTrackPhi[ij] = 0;
217 for(int i = 0;i <= kMaxJets;++i){
218 fh2PhiPt[ij][i] = 0;
219 fh2EtaPt[ij][i] = 0;
220 fh2AreaPt[ij][i] = 0;
221 fh2EtaArea[ij][i] = 0;
222 fh2PhiEta[ij][i] = 0;
223 fh2LTrackPtJetPt[ij][i] = 0;
224 fh1PtIn[ij][i] = 0;
225 fh2RhoPt[ij][i] = 0;
226 fh2PsiPt[ij][i] = 0;
227 }
228
229 fh1DijetMinv[ij] = 0;
230 fh2DijetDeltaPhiPt[ij] = 0;
231 fh2DijetAsymPt[ij] = 0;
232 fh2DijetPt2vsPt1[ij] = 0;
233 fh2DijetDifvsSum[ij] = 0;
234 }
235
236 DefineOutput(1, TList::Class());
237}
238
239
240
241Bool_t AliAnalysisTaskJetSpectrum2::Notify()
242{
243
244
245
246 //
247 // Implemented Notify() to read the cross sections
248 // and number of trials from pyxsec.root
249 //
250
251 // Fetch the aod also from the input in,
252 // have todo it in notify
253
254
255 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
256 // assume that the AOD is in the general output...
257 fAODOut = AODEvent();
258
259 if(fNonStdFile.Length()!=0){
260 // case that we have an AOD extension we need can fetch the jets from the extended output
261 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
262 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
263 if(!fAODExtension){
264 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
265 }
266 }
267
268
269 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
270 Float_t xsection = 0;
271 Float_t ftrials = 1;
272
273 fAvgTrials = 1;
274 if(tree){
275 TFile *curfile = tree->GetCurrentFile();
276 if (!curfile) {
277 Error("Notify","No current file");
278 return kFALSE;
279 }
280 if(!fh1Xsec||!fh1Trials){
281 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
282 return kFALSE;
283 }
284 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
285 fh1Xsec->Fill("<#sigma>",xsection);
286 // construct a poor man average trials
287 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
288 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
289 }
290
291 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
292
293 return kTRUE;
294}
295
296void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
297{
298
299
300 // Connect the AOD
301
302 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
303 OpenFile(1);
304 if(!fHistList)fHistList = new TList();
305 fHistList->SetOwner(kTRUE);
306 Bool_t oldStatus = TH1::AddDirectoryStatus();
307 TH1::AddDirectory(kFALSE);
308
309 MakeJetContainer();
310 fHistList->Add(fhnCorrelation);
311 fHistList->Add(fhnCorrelationPhiZRec);
312 for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
313
314 //
315 // Histogram
316
317 const Int_t nBinPt = 320;
318 Double_t binLimitsPt[nBinPt+1];
319 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
320 if(iPt == 0){
321 binLimitsPt[iPt] = 0.0;
322 }
323 else {// 1.0
324 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
325 }
326 }
327 const Int_t nBinPhi = 90;
328 Double_t binLimitsPhi[nBinPhi+1];
329 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
330 if(iPhi==0){
331 binLimitsPhi[iPhi] = 0;
332 }
333 else{
334 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
335 }
336 }
337
338 const Int_t nBinEta = 40;
339 Double_t binLimitsEta[nBinEta+1];
340 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
341 if(iEta==0){
342 binLimitsEta[iEta] = -2.0;
343 }
344 else{
345 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
346 }
347 }
348
349
350 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
351 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
352 fHistList->Add(fh1Xsec);
353 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
354 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
355 fHistList->Add(fh1Trials);
356 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
357 fHistList->Add(fh1PtHard);
358 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
359 fHistList->Add(fh1PtHardNoW);
360 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
361 fHistList->Add(fh1PtHardTrials);
362
363 fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
364 fHistList->Add(fh1ZVtx);
365 fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",400,-0.5,4000,400,0.,4000);
366 fHistList->Add(fh2MultRec);
367 fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,-0.5,4000,400,0.,4000);
368 fHistList->Add(fh2MultGen);
369
370 fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
371 fHistList->Add(fh2PtFGen);
372
373 fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
374 fHistList->Add(fh2RelPtFGen);
375
376 for(int ij = 0;ij <kJetTypes;++ij){
377 TString cAdd = "";
378 TString cJetBranch = "";
379 if(ij==kJetRec){
380 cAdd = "Rec";
381 cJetBranch = fBranchRec.Data();
382 }
383 else if (ij==kJetGen){
384 cAdd = "Gen";
385 cJetBranch = fBranchGen.Data();
386 }
387 else if (ij==kJetRecFull){
388 cAdd = "RecFull";
389 cJetBranch = fBranchRec.Data();
390 }
391 else if (ij==kJetGenFull){
392 cAdd = "GenFull";
393 cJetBranch = fBranchGen.Data();
394 }
395
396 fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
397 fHistList->Add(fh1NJets[ij]);
398
399 fh1PtJetsIn[ij] = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
400 fHistList->Add(fh1PtJetsIn[ij]);
401
402 fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
403 fHistList->Add(fh1PtTracksIn[ij]);
404
405 fh1PtTracksLeadingIn[ij] = new TH1F(Form("fh1PtTracksLeading%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
406 fHistList->Add(fh1PtTracksLeadingIn[ij]);
407
408 fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
409 fHistList->Add(fh1SumPtTrack[ij]);
410
411 fh2MultJetPt[ij] = new TH2F(Form("fh2MultJetPt%s",cAdd.Data()),Form("%s jets p_T;# tracks;;p_{T} (GeV/c)",cAdd.Data()),400,0,4000,nBinPt,binLimitsPt);
412 fHistList->Add(fh2MultJetPt[ij]);
413
414 fh2NJetsPt[ij] = new TH2F(Form("fh2N%sJetsPt",cAdd.Data()),Form("Number of %s jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",cAdd.Data()),nBinPt,binLimitsPt,50,-0.5,49.5);
415 fHistList->Add(fh2NJetsPt[ij]);
416
417 fh2NTracksPt[ij] = new TH2F(Form("fh2N%sTracksPt",cAdd.Data()),Form("Number of %s tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",cAdd.Data()),nBinPt,binLimitsPt,1000,0.,4000);
418 fHistList->Add(fh2NTracksPt[ij]);
419
420 fh2LeadingTrackPtTrackPhi[ij] = new TH2F(Form("fh2Leading%sTrackPtTrackPhi",cAdd.Data()),Form("phi of leading %s track;p_{T};#phi;",cAdd.Data()),
421 nBinPt,binLimitsPt,nBinPhi,binLimitsPhi);
422 fHistList->Add(fh2LeadingTrackPtTrackPhi[ij]);
423
424 for(int i = 0;i <= kMaxJets;++i){
425 fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
426 fHistList->Add(fh1PtIn[ij][i]);
427
428 fh2RhoPt[ij][i] = new TH2F(Form("fh2RhoPt%s_j%d",cAdd.Data(),i),Form("jet shape rho for %s jets;r;p_{T};",cAdd.Data()),
429 40,0.,2.,nBinPt,binLimitsPt);
430 fHistList->Add(fh2RhoPt[ij][i]);
431 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
432 fh2PsiPt[ij][i] = new TH2F(Form("fh2PsiPt%s_j%d",cAdd.Data(),i),Form("jet shape #psi for %s jets;r;p_{T};",cAdd.Data()),
433 40,0.,2.,nBinPt,binLimitsPt);
434 fHistList->Add(fh2PsiPt[ij][i]);
435 fh2PhiPt[ij][i] = new TH2F(Form("fh2PhiPt%s_j%d",cAdd.Data(),i),Form("pt vs phi %s;#phi;p_{T};",cAdd.Data()),
436 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
437 fHistList->Add(fh2PhiPt[ij][i]);
438 fh2EtaPt[ij][i] = new TH2F(Form("fh2EtaPt%s_j%d",cAdd.Data(),i),Form("pt vs eta %s;#eta;p_{T};",cAdd.Data()),
439 50,-1.,1.,nBinPt,binLimitsPt);
440 fHistList->Add(fh2EtaPt[ij][i]);
441 fh2AreaPt[ij][i] = new TH2F(Form("fh2AreaPt%s_j%d",cAdd.Data(),i),
442 Form("pt vs area %s;area;p_{T};",
443 cAdd.Data()),
444 50,0.,1.,nBinPt,binLimitsPt);
445 fHistList->Add(fh2AreaPt[ij][i]);
446 fh2EtaArea[ij][i] = new TH2F(Form("fh2EtaArea%s_j%d",cAdd.Data(),i),
447 Form("area vs eta %s;#eta;area;",
448 cAdd.Data()),
449 50,-1.,1.,50,0,1.);
450 fHistList->Add(fh2EtaArea[ij][i]);
451
452 fh2PhiEta[ij][i] = new TH2F(Form("fh2PhiEta%s_j%d",cAdd.Data(),i),Form("phi vs eta %s ;phi;p_{T};",cAdd.Data()),
453 nBinPhi,binLimitsPhi,20,-1.,1.);
454 fHistList->Add(fh2PhiEta[ij][i]);
455
456 fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i),
457 Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};",
458 cAdd.Data()),
459 200,0.,200.,nBinPt,binLimitsPt);
460 fHistList->Add(fh2LTrackPtJetPt[ij][i]);
461 }
462
463
464 fh1DijetMinv[ij] = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
465 fHistList->Add(fh1DijetMinv[ij]);
466
467 fh2DijetDeltaPhiPt[ij] = new TH2F(Form("fh2Dijet%sDeltaPhiPt",cAdd.Data()),"Difference in the azimuthal angle;#Delta#phi;p_{T,2};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
468 fHistList->Add(fh2DijetDeltaPhiPt[ij]);
469
470 fh2DijetAsymPt[ij] = new TH2F(Form("fh2Dijet%sAsym",cAdd.Data()),"Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
471 fHistList->Add(fh2DijetAsymPt[ij]);
472
473 fh2DijetPt2vsPt1[ij] = new TH2F(Form("fh2Dijet%sPt2vsPt1",cAdd.Data()),"Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
474 fHistList->Add(fh2DijetPt2vsPt1[ij]);
475 fh2DijetDifvsSum[ij] = new TH2F(Form("fh2Dijet%sDifvsSum",cAdd.Data()),"Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.);
476 fHistList->Add( fh2DijetDifvsSum[ij]);
477 }
478 // =========== Switch on Sumw2 for all histos ===========
479 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
480 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
481 if (h1){
482 h1->Sumw2();
483 continue;
484 }
485 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
486 if(hn)hn->Sumw2();
487 }
488 TH1::AddDirectory(oldStatus);
489}
490
491void AliAnalysisTaskJetSpectrum2::Init()
492{
493 //
494 // Initialization
495 //
496
497 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
498
499}
500
501void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
502
503
504 Bool_t selected = kTRUE;
505
506 if(fUseGlobalSelection&&fEventSelectionMask==0){
507 selected = AliAnalysisHelperJetTasks::Selected();
508 }
509 else if(fUseGlobalSelection&&fEventSelectionMask>0){
510 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
511 }
512
513 if(fEventClass>0){
514 selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
515 }
516
517 if(!selected){
518 // no selection by the service task, we continue
519 if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
520 PostData(1, fHistList);
521 return;
522 }
523
524
525 static AliAODEvent* aod = 0;
526
527 // take all other information from the aod we take the tracks from
528 if(!aod){
529 if(fUseAODTrackInput)aod = fAODIn;
530 else aod = fAODOut;
531 }
532
533
534 //
535 // Execute analysis for current event
536 //
537 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
538 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
539 if(!aodH){
540 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
541 return;
542 }
543
544 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
545 TClonesArray *aodRecJets = 0;
546
547 if(fAODOut&&!aodRecJets){
548 aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
549 }
550 if(fAODExtension&&!aodRecJets){
551 aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
552 }
553 if(fAODIn&&!aodRecJets){
554 aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
555 }
556
557
558
559 if(!aodRecJets){
560 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
561 return;
562 }
563
564 TClonesArray *aodGenJets = 0;
565 if(fBranchGen.Length()>0){
566 if(fAODOut&&!aodGenJets){
567 aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
568 }
569 if(fAODExtension&&!aodGenJets){
570 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
571 }
572 if(fAODIn&&!aodGenJets){
573 aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
574 }
575
576 if(!aodGenJets){
577 Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
578 return;
579 }
580 }
581
582 TClonesArray *aodBackRecJets = 0;
583 if(fBranchBkgRec.Length()>0){
584 if(fAODOut&&!aodBackRecJets){
585 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgRec.Data()));
586 }
587 if(fAODExtension&&!aodBackRecJets){
588 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgRec.Data()));
589 }
590 if(fAODIn&&!aodBackRecJets){
591 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgRec.Data()));
592 }
593
594 if(!aodBackRecJets){
595 Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgRec.Data());
596 return;
597 }
598 }
599
600
601 TClonesArray *aodBackGenJets = 0;
602
603 if(fBranchBkgGen.Length()>0){
604 if(fAODOut&&!aodBackGenJets){
605 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgGen.Data()));
606 }
607 if(fAODExtension&&!aodBackGenJets){
608 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgGen.Data()));
609 }
610 if(fAODIn&&!aodBackGenJets){
611 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgGen.Data()));
612 }
613
614 if(!aodBackGenJets){
615 Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgGen.Data());
616 return;
617 }
618 }
619
620
621 // new Scheme
622 // first fill all the pure histograms, i.e. full jets
623 // in case of the correaltion limit to the n-leading jets
624
625 // reconstructed
626
627
628 // generated
629
630
631 // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
632
633
634 TList genJetsList; // full acceptance
635 TList genJetsListCut; // acceptance cut
636 TList recJetsList; // full acceptance
637 TList recJetsListCut; // acceptance cut
638
639
640 GetListOfJets(&recJetsList,aodRecJets,0);
641 GetListOfJets(&recJetsListCut,aodRecJets,1);
642
643 if(fBranchGen.Length()>0){
644 GetListOfJets(&genJetsList,aodGenJets,0);
645 GetListOfJets(&genJetsListCut,aodGenJets,1);
646 }
647
648 Double_t eventW = 1;
649 Double_t ptHard = 0;
650 Double_t nTrials = 1; // Trials for MC trigger
651 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
652
653 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
654 // this is the part we only use when we have MC information
655 AliMCEvent* mcEvent = MCEvent();
656 // AliStack *pStack = 0;
657 if(!mcEvent){
658 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
659 return;
660 }
661 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
662 if(pythiaGenHeader){
663 nTrials = pythiaGenHeader->Trials();
664 ptHard = pythiaGenHeader->GetPtHard();
665 int iProcessType = pythiaGenHeader->ProcessType();
666 // 11 f+f -> f+f
667 // 12 f+barf -> f+barf
668 // 13 f+barf -> g+g
669 // 28 f+g -> f+g
670 // 53 g+g -> f+barf
671 // 68 g+g -> g+g
672 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
673 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
674 }
675 }// (fAnalysisType&kMCESD)==kMCESD)
676
677
678 // we simply fetch the tracks/mc particles as a list of AliVParticles
679
680 TList recParticles;
681 TList genParticles;
682
683 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
684 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
685 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
686 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
687
688 // Event control and counting ...
689 // MC
690 fh1PtHard->Fill(ptHard,eventW);
691 fh1PtHardNoW->Fill(ptHard,1);
692 fh1PtHardTrials->Fill(ptHard,nTrials);
693
694 // Real
695 if(aod->GetPrimaryVertex()){// No vtx for pure MC
696 fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
697 }
698
699
700 Int_t recMult1 = recParticles.GetEntries();
701 Int_t genMult1 = genParticles.GetEntries();
702
703 Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
704 Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
705
706 fh2MultRec->Fill(recMult1,recMult2);
707 fh2MultGen->Fill(genMult1,genMult2);
708 fMultRec = recMult1;
709 if(fMultRec<=0)fMultRec = recMult2;
710 fMultGen = genMult1;
711 if(fMultGen<=0)fMultGen = genMult2;
712
713 // the loops for rec and gen should be indentical... pass it to a separate
714 // function ...
715 // Jet Loop
716 // Track Loop
717 // Jet Jet Loop
718 // Jet Track loop
719
720 FillJetHistos(recJetsListCut,recParticles,kJetRec);
721 FillJetHistos(recJetsList,recParticles,kJetRecFull);
722 FillTrackHistos(recParticles,kJetRec);
723
724 FillJetHistos(genJetsListCut,genParticles,kJetGen);
725 FillJetHistos(genJetsList,genParticles,kJetGenFull);
726 FillTrackHistos(genParticles,kJetGen);
727
728 // Here follows the jet matching part
729 // delegate to a separated method?
730
731 FillMatchHistos(recJetsList,genJetsList);
732
733 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
734 PostData(1, fHistList);
735}
736
737void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
738
739 if(iType>=kJetTypes){
740 return;
741 }
742
743 Int_t refMult = fMultRec;
744 if(iType==kJetGen||iType==kJetGenFull){
745 refMult = fMultGen;
746 }
747
748 Int_t nJets = jetsList.GetEntries();
749 fh1NJets[iType]->Fill(nJets);
750
751 if(nJets<=0)return;
752
753 Float_t ptOld = 1.E+32;
754 Float_t pT0 = 0;
755 Float_t pT1 = 0;
756 Float_t phi0 = 0;
757 Float_t phi1 = 0;
758 Int_t ij0 = -1;
759 Int_t ij1 = -1;
760
761
762 for(int ij = 0;ij < nJets;ij++){
763 AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
764 Float_t ptJet = jet->Pt();
765 fh1PtJetsIn[iType]->Fill(ptJet);
766 fh2MultJetPt[iType]->Fill(refMult,ptJet);
767 if(ptJet>ptOld){
768 Printf("%s:%d Jets Type %d Not Sorted !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,iType,ij,ptJet,ij-1,ptOld);
769 }
770 ptOld = ptJet;
771
772 // find the dijets assume sorting and acceptance cut...
773 if(ij==0){
774 ij0 = ij;
775 pT0 = ptJet;
776 phi0 = jet->Phi();
777 if(phi0<0)phi0+=TMath::Pi()*2.;
778 }
779 else if(ptJet>pT1){
780 // take only the backward hemisphere??
781 phi1 = jet->Phi();
782 if(phi1<0)phi1+=TMath::Pi()*2.;
783 Float_t dPhi = phi1 - phi0;
784 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
785 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
786 if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
787 ij1 = ij;
788 pT1 = ptJet;
789 }
790 }
791 // fill jet histos for kmax jets
792
793 Float_t phiJet = jet->Phi();
794 Float_t etaJet = jet->Eta();
795 if(phiJet<0)phiJet+=TMath::Pi()*2.;
796 fh1TmpRho->Reset();
797 if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
798 fh1PtIn[iType][kMaxJets]->Fill(ptJet);
799 // fill leading jets...
800 AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
801 // AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
802 if(ptJet>10){
803 if(ij<kMaxJets){
804 fh2PhiEta[iType][ij]->Fill(phiJet,etaJet);
805 fh2AreaPt[iType][ij]->Fill(jet->EffectiveAreaCharged(),ptJet);
806 fh2EtaArea[iType][ij]->Fill(etaJet,jet->EffectiveAreaCharged());
807 }
808 fh2PhiEta[iType][kMaxJets]->Fill(phiJet,etaJet);
809 fh2AreaPt[iType][kMaxJets]->Fill(jet->EffectiveAreaCharged(),ptJet);
810 fh2EtaArea[iType][kMaxJets]->Fill(etaJet,jet->EffectiveAreaCharged());
811 }
812 if(ij<kMaxJets){
813 fh2PhiPt[iType][ij]->Fill(phiJet,ptJet);
814 fh2EtaPt[iType][ij]->Fill(etaJet,ptJet);
815 if(leadTrack)fh2LTrackPtJetPt[iType][ij]->Fill(leadTrack->Pt(),ptJet);
816 }
817 fh2PhiPt[iType][kMaxJets]->Fill(phiJet,ptJet);
818 fh2EtaPt[iType][kMaxJets]->Fill(etaJet,ptJet);
819 if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet);
820
821 if(particlesList.GetSize()&&ij<kMaxJets){
822 // Particles... correlated with jets...
823 for(int it = 0;it<particlesList.GetEntries();++it){
824 AliVParticle *part = (AliVParticle*)particlesList.At(it);
825 Float_t deltaR = jet->DeltaR(part);
826 if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
827 }
828 // fill the jet shapes
829 Float_t rhoSum = 0;
830 for(int ibx = 1;ibx<=fh2RhoPt[iType][ij]->GetNbinsX();ibx++){
831 Float_t r = fh2RhoPt[iType][ij]->GetXaxis()->GetBinCenter(ibx);
832 Float_t rho = fh1TmpRho->GetBinContent(ibx);
833 rhoSum += rho;
834 fh2RhoPt[iType][ij]->Fill(r,ptJet,rho);
835 fh2PsiPt[iType][ij]->Fill(r,ptJet,rhoSum);
836 }
837 }// if we have particles
838 }// Jet Loop
839
840
841 // do something with dijets...
842 if(ij0>=0&&ij1>0){
843 AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
844 Double_t ptJet0 = jet0->Pt();
845 Double_t phiJet0 = jet0->Phi();
846 if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
847
848 AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
849 Double_t ptJet1 = jet1->Pt();
850 Double_t phiJet1 = jet1->Phi();
851 if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
852
853 Float_t deltaPhi = phiJet0 - phiJet1;
854 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
855 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
856 deltaPhi = TMath::Abs(deltaPhi);
857 fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
858
859 Float_t asym = 9999;
860 if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
861 fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
862 fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
863 fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
864 Float_t minv = 2.*(jet0->P()*jet1->P()-
865 jet0->Px()*jet1->Px()-
866 jet0->Py()*jet1->Py()-
867 jet0->Pz()*jet1->Pz()); // assume mass == 0;
868 if(minv<0)minv=0; // prevent numerical instabilities
869 minv = TMath::Sqrt(minv);
870 fh1DijetMinv[iType]->Fill(minv);
871 }
872
873
874
875 // count down the jets above thrueshold
876 Int_t nOver = nJets;
877 if(nOver>0){
878 TIterator *jetIter = jetsList.MakeIterator();
879 AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
880 if(tmpJet){
881 Float_t pt = tmpJet->Pt();
882 for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
883 Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
884 while(pt<ptCut&&tmpJet){
885 nOver--;
886 tmpJet = (AliAODJet*)(jetIter->Next());
887 if(tmpJet){
888 pt = tmpJet->Pt();
889 }
890 }
891 if(nOver<=0)break;
892 fh2NJetsPt[iType]->Fill(ptCut,nOver);
893 }
894 }
895 delete jetIter;
896 }
897}
898
899void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
900
901 Int_t nTrackOver = particlesList.GetSize();
902 // do the same for tracks and jets
903 if(nTrackOver>0){
904 TIterator *trackIter = particlesList.MakeIterator();
905 AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
906 Float_t pt = tmpTrack->Pt();
907 for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
908 Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
909 while(pt<ptCut&&tmpTrack){
910 nTrackOver--;
911 tmpTrack = (AliVParticle*)(trackIter->Next());
912 if(tmpTrack){
913 pt = tmpTrack->Pt();
914 }
915 }
916 if(nTrackOver<=0)break;
917 fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
918 }
919
920
921 trackIter->Reset();
922 AliVParticle *leading = (AliVParticle*)particlesList.At(0);
923 Float_t sumPt = 0;
924
925 while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
926 Float_t tmpPt = tmpTrack->Pt();
927 fh1PtTracksIn[iType]->Fill(tmpPt);
928 sumPt += tmpPt;
929 Float_t tmpPhi = tmpTrack->Phi();
930 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
931 if(tmpTrack==leading){
932 fh1PtTracksLeadingIn[iType]->Fill(tmpPt);
933 fh2LeadingTrackPtTrackPhi[iType]->Fill(tmpPt,tmpPhi);
934 continue;
935 }
936 }
937 fh1SumPtTrack[iType]->Fill(sumPt);
938 delete trackIter;
939 }
940
941}
942
943
944void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
945
946
947 // Fill al the matching histos
948 // It is important that the acceptances for the mathing are as large as possible
949 // to avoid false matches on the edge of acceptance
950 // therefore we add some extra matching jets as overhead
951
952 static TArrayI aGenIndex(recJetsList.GetEntries());
953 if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
954
955 static TArrayI aRecIndex(genJetsList.GetEntries());
956 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
957
958 if(fDebug){
959 Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
960 Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
961 }
962 AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
963 &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
964 aGenIndex,aRecIndex,fDebug);
965
966 if(fDebug){
967 for(int i = 0;i< aGenIndex.GetSize();++i){
968 if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
969 }
970 for(int i = 0;i< aRecIndex.GetSize();++i){
971 if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
972 }
973 }
974
975 Double_t container[6];
976 Double_t containerPhiZ[6];
977
978 // loop over generated jets
979 // consider the
980 for(int ig = 0;ig < genJetsList.GetEntries();++ig){
981 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
982 Double_t ptGen = genJet->Pt();
983 Double_t phiGen = genJet->Phi();
984 if(phiGen<0)phiGen+=TMath::Pi()*2.;
985 Double_t etaGen = genJet->Eta();
986 container[3] = ptGen;
987 container[4] = etaGen;
988 container[5] = phiGen;
989 fhnJetContainer[kStep0]->Fill(&container[3]);
990 if(JetSelected(genJet)){
991 fhnJetContainer[kStep1]->Fill(&container[3]);
992 Int_t ir = aRecIndex[ig];
993 if(ir>=0&&ir<recJetsList.GetEntries()){
994 fhnJetContainer[kStep2]->Fill(&container[3]);
995 AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
996 if(JetSelected(recJet))fhnJetContainer[kStep4]->Fill(&container[3]);
997 if(JetSelected(recJet))fhnJetContainer[kStep3]->Fill(&container[3]);
998 }
999 }
1000 }// loop over generated jets used for matching...
1001
1002
1003
1004 // loop over reconstructed jets
1005 for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1006 AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1007 Double_t etaRec = recJet->Eta();
1008 Double_t ptRec = recJet->Pt();
1009 Double_t phiRec = recJet->Phi();
1010 if(phiRec<0)phiRec+=TMath::Pi()*2.;
1011 // do something with dijets...
1012
1013 container[0] = ptRec;
1014 container[1] = etaRec;
1015 container[2] = phiRec;
1016 containerPhiZ[0] = ptRec;
1017 containerPhiZ[1] = phiRec;
1018
1019 fhnJetContainer[kStep0+kMaxStep]->Fill(container);
1020 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1021
1022 if(JetSelected(recJet)){
1023 fhnJetContainer[kStep1+kMaxStep]->Fill(container);
1024 // Fill Correlation
1025 Int_t ig = aGenIndex[ir];
1026 if(ig>=0 && ig<genJetsList.GetEntries()){
1027 fhnJetContainer[kStep2+kMaxStep]->Fill(container);
1028 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1029 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1030 Double_t ptGen = genJet->Pt();
1031 Double_t phiGen = genJet->Phi();
1032 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1033 Double_t etaGen = genJet->Eta();
1034
1035 container[3] = ptGen;
1036 container[4] = etaGen;
1037 container[5] = phiGen;
1038 containerPhiZ[3] = ptGen;
1039 //
1040 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1041 //
1042 if(JetSelected(genJet))fhnJetContainer[kStep4+kMaxStep]->Fill(container);
1043 fhnJetContainer[kStep3+kMaxStep]->Fill(container);
1044 fhnCorrelation->Fill(container);
1045 if(ptGen>0){
1046 Float_t delta = (ptRec-ptGen)/ptGen;
1047 fh2RelPtFGen->Fill(ptGen,delta);
1048 fh2PtFGen->Fill(ptGen,ptRec);
1049 }
1050 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1051 }
1052 else{
1053 containerPhiZ[3] = 0;
1054 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1055 }
1056 }// loop over reconstructed jets
1057 }
1058 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1059}
1060
1061
1062void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1063 //
1064 // Create the particle container for the correction framework manager and
1065 // link it
1066 //
1067 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
1068 const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1069 const Double_t kEtamin = -3.0, kEtamax = 3.0;
1070 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1071 const Double_t kZmin = 0., kZmax = 1;
1072
1073 // can we neglect migration in eta and phi?
1074 // phi should be no problem since we cover full phi and are phi symmetric
1075 // eta migration is more difficult due to needed acceptance correction
1076 // in limited eta range
1077
1078 //arrays for the number of bins in each dimension
1079 Int_t iBin[kNvar];
1080 iBin[0] = 320; //bins in pt
1081 iBin[1] = 1; //bins in eta
1082 iBin[2] = 1; // bins in phi
1083
1084 //arrays for lower bounds :
1085 Double_t* binEdges[kNvar];
1086 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1087 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1088
1089 //values for bin lower bounds
1090 // for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
1091 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1092 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1093 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1094
1095
1096 for(int i = 0;i < kMaxStep*2;++i){
1097 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1098 for (int k=0; k<kNvar; k++) {
1099 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1100 }
1101 }
1102 //create correlation matrix for unfolding
1103 Int_t thnDim[2*kNvar];
1104 for (int k=0; k<kNvar; k++) {
1105 //first half : reconstructed
1106 //second half : MC
1107 thnDim[k] = iBin[k];
1108 thnDim[k+kNvar] = iBin[k];
1109 }
1110
1111 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1112 for (int k=0; k<kNvar; k++) {
1113 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1114 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1115 }
1116 fhnCorrelation->Sumw2();
1117
1118 // for second correlation histogram
1119
1120
1121 const Int_t kNvarPhiZ = 4;
1122 //arrays for the number of bins in each dimension
1123 Int_t iBinPhiZ[kNvarPhiZ];
1124 iBinPhiZ[0] = 80; //bins in pt
1125 iBinPhiZ[1] = 72; //bins in phi
1126 iBinPhiZ[2] = 20; // bins in Z
1127 iBinPhiZ[3] = 80; //bins in ptgen
1128
1129 //arrays for lower bounds :
1130 Double_t* binEdgesPhiZ[kNvarPhiZ];
1131 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1132 binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1133
1134 for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1135 for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1136 for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1137 for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1138
1139 fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1140 for (int k=0; k<kNvarPhiZ; k++) {
1141 fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1142 }
1143 fhnCorrelationPhiZRec->Sumw2();
1144
1145
1146 // Add a histogram for Fake jets
1147
1148 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1149 delete [] binEdges[ivar];
1150
1151 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1152 delete [] binEdgesPhiZ[ivar];
1153
1154}
1155
1156void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1157{
1158// Terminate analysis
1159//
1160 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1161}
1162
1163
1164Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1165
1166 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1167
1168 Int_t iCount = 0;
1169 if(type==kTrackAOD){
1170 AliAODEvent *aod = 0;
1171 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1172 else aod = AODEvent();
1173 if(!aod){
1174 return iCount;
1175 }
1176 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1177 AliAODTrack *tr = aod->GetTrack(it);
1178 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1179 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1180 if(tr->Pt()<fMinTrackPt)continue;
1181 if(fDebug>0){
1182 if(tr->Pt()>20){
1183 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1184 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1185 tr->Print();
1186 // tr->Dump();
1187 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1188 if(fESD){
1189 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1190 esdTr->Print("");
1191 // esdTr->Dump();
1192 }
1193 }
1194 }
1195 list->Add(tr);
1196 iCount++;
1197 }
1198 }
1199 else if (type == kTrackKineAll||type == kTrackKineCharged){
1200 AliMCEvent* mcEvent = MCEvent();
1201 if(!mcEvent)return iCount;
1202 // we want to have alivpartilces so use get track
1203 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1204 if(!mcEvent->IsPhysicalPrimary(it))continue;
1205 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1206 if(part->Pt()<fMinTrackPt)continue;
1207 if(type == kTrackKineAll){
1208 list->Add(part);
1209 iCount++;
1210 }
1211 else if(type == kTrackKineCharged){
1212 if(part->Particle()->GetPDG()->Charge()==0)continue;
1213 list->Add(part);
1214 iCount++;
1215 }
1216 }
1217 }
1218 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1219 AliAODEvent *aod = 0;
1220 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1221 else aod = AODEvent();
1222 if(!aod)return iCount;
1223 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1224 if(!tca)return iCount;
1225 for(int it = 0;it < tca->GetEntriesFast();++it){
1226 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1227 if(!part)continue;
1228 if(part->Pt()<fMinTrackPt)continue;
1229 if(!part->IsPhysicalPrimary())continue;
1230 if(type == kTrackAODMCAll){
1231 list->Add(part);
1232 iCount++;
1233 }
1234 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1235 if(part->Charge()==0)continue;
1236 if(kTrackAODMCCharged){
1237 list->Add(part);
1238 }
1239 else {
1240 if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1241 list->Add(part);
1242 }
1243 iCount++;
1244 }
1245 }
1246 }// AODMCparticle
1247 list->Sort();
1248 return iCount;
1249
1250}
1251
1252
1253
1254Bool_t AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1255 Bool_t selected = false;
1256
1257 if(!jet)return selected;
1258
1259 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1260 selected = kTRUE;
1261 }
1262 return selected;
1263
1264}
1265
1266Int_t AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1267
1268 if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1269 Int_t iCount = 0;
1270
1271 if(!jarray){
1272 Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1273 return iCount;
1274 }
1275
1276
1277 for(int ij=0;ij<jarray->GetEntries();ij++){
1278 AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1279 if(!jet)continue;
1280 if(type==0){
1281 // No cut at all, main purpose here is sorting
1282 list->Add(jet);
1283 iCount++;
1284 }
1285 else if(type == 1){
1286 // eta cut
1287 if(JetSelected(jet)){
1288 list->Add(jet);
1289 iCount++;
1290 }
1291 }
1292 }
1293
1294 list->Sort();
1295 return iCount;
1296
1297}
1298
1299
1300Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1301 if(!jets)return 0;
1302
1303 Int_t refMult = 0;
1304 for(int ij = 0;ij < jets->GetEntries();++ij){
1305 AliAODJet* jet = (AliAODJet*)jets->At(ij);
1306 if(!jet)continue;
1307 TRefArray *refs = jet->GetRefTracks();
1308 if(!refs)continue;
1309 refMult += refs->GetEntries();
1310 }
1311 return refMult;
1312
1313}
1314
1315
1316AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1317 if(!jet)return 0;
1318 TRefArray *refs = jet->GetRefTracks();
1319 if(!refs) return 0;
1320 AliVParticle *leading = 0;
1321 Float_t fMaxPt = 0;
1322 for(int i = 0;i<refs->GetEntries();i++){
1323 AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1324 if(!tmp)continue;
1325 if(tmp->Pt()>fMaxPt){
1326 leading = tmp;
1327 fMaxPt = tmp->Pt();
1328 }
1329 }
1330 return leading;
1331}
1332
1333
1334AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1335 if(!jet)return 0;
1336 if(!list) return 0;
1337 AliVParticle *leading = 0;
1338 Float_t fMaxPt = 0;
1339 for(int i = 0;i<list->GetEntries();i++){
1340 AliVParticle *tmp = (AliVParticle*)(list->At(i));
1341 if(!tmp)continue;
1342 if(jet->DeltaR(tmp)>r)continue;
1343 if(tmp->Pt()>fMaxPt){
1344 leading = tmp;
1345 fMaxPt = tmp->Pt();
1346 }
1347 }
1348 return leading;
1349}