]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Update of the GPU tracker
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum2.cxx
CommitLineData
3b7ffecf 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 "TDatabasePDG.h"
38
39#include "AliAnalysisTaskJetSpectrum2.h"
40#include "AliAnalysisManager.h"
41#include "AliJetFinder.h"
42#include "AliJetHeader.h"
43#include "AliJetReader.h"
44#include "AliJetReaderHeader.h"
45#include "AliUA1JetHeaderV1.h"
46#include "AliESDEvent.h"
47#include "AliAODEvent.h"
48#include "AliAODHandler.h"
49#include "AliAODTrack.h"
50#include "AliAODJet.h"
51#include "AliAODMCParticle.h"
52#include "AliMCEventHandler.h"
53#include "AliMCEvent.h"
54#include "AliStack.h"
55#include "AliGenPythiaEventHeader.h"
56#include "AliJetKineReaderHeader.h"
57#include "AliGenCocktailEventHeader.h"
58#include "AliInputEventHandler.h"
59
60
61#include "AliAnalysisHelperJetTasks.h"
62
63ClassImp(AliAnalysisTaskJetSpectrum2)
64
65AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
26fa06fb 66 fJetHeaderRec(0x0),
67 fJetHeaderGen(0x0),
68 fAOD(0x0),
69 fhnCorrelation(0x0),
70 fBranchRec("jets"),
71 fBranchGen(""),
72 fUseAODJetInput(kFALSE),
73 fUseAODTrackInput(kFALSE),
74 fUseAODMCInput(kFALSE),
75 fUseGlobalSelection(kFALSE),
76 fUseExternalWeightOnly(kFALSE),
77 fLimitGenJetEta(kFALSE),
78 fFilterMask(0),
79 fAnalysisType(0),
3b7ffecf 80 fTrackTypeRec(kTrackUndef),
81 fTrackTypeGen(kTrackUndef),
82 fAvgTrials(1),
83 fExternalWeight(1),
26fa06fb 84 fRecEtaWindow(0.5),
9280dfa6 85 fMinJetPt(0),
26fa06fb 86 fDeltaPhiWindow(20./180.*TMath::Pi()),
3b7ffecf 87 fh1Xsec(0x0),
88 fh1Trials(0x0),
89 fh1PtHard(0x0),
90 fh1PtHardNoW(0x0),
91 fh1PtHardTrials(0x0),
92 fh1NGenJets(0x0),
93 fh1NRecJets(0x0),
edfbe476 94 fh1PtTrackRec(0x0),
95 fh1SumPtTrackRec(0x0),
96 fh1SumPtTrackAreaRec(0x0),
cc0649e4 97 fh1PtJetsRecIn(0x0),
565584e8 98 fh1PtJetsLeadingRecIn(0x0),
cc0649e4 99 fh1PtTracksRecIn(0x0),
565584e8 100 fh1PtTracksLeadingRecIn(0x0),
cc0649e4 101 fh1PtTracksGenIn(0x0),
102 fh2NRecJetsPt(0x0),
103 fh2NRecTracksPt(0x0),
104 fh2JetsLeadingPhiEta(0x0),
565584e8 105 fh2JetsLeadingPhiPt(0x0),
cc0649e4 106 fh2TracksLeadingPhiEta(0x0),
565584e8 107 fh2TracksLeadingPhiPt(0x0),
cf049e94 108 fh2TracksLeadingJetPhiPt(0x0),
26fa06fb 109 fh2DijetDeltaPhiPt(0x0),
110 fh2DijetAsymPt(0x0),
111 fh2DijetAsymPtCut(0x0),
112 fh2DijetDeltaPhiDeltaEta(0x0),
113 fh2DijetPt2vsPt1(0x0),
114 fh2DijetDifvsSum(0x0),
115 fh1DijetMinv(0x0),
116 fh1DijetMinvCut(0x0),
3b7ffecf 117 fHistList(0x0)
118{
119 for(int i = 0;i < kMaxStep*2;++i){
120 fhnJetContainer[i] = 0;
121 }
122 for(int i = 0;i < kMaxJets;++i){
edfbe476 123 fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
3b7ffecf 124 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
125 }
126
127}
128
129AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
130 AliAnalysisTaskSE(name),
131 fJetHeaderRec(0x0),
132 fJetHeaderGen(0x0),
133 fAOD(0x0),
134 fhnCorrelation(0x0),
135 fBranchRec("jets"),
136 fBranchGen(""),
565584e8 137 fUseAODJetInput(kFALSE),
138 fUseAODTrackInput(kFALSE),
139 fUseAODMCInput(kFALSE),
b5cc0c6d 140 fUseGlobalSelection(kFALSE),
3b7ffecf 141 fUseExternalWeightOnly(kFALSE),
142 fLimitGenJetEta(kFALSE),
143 fFilterMask(0),
144 fAnalysisType(0),
145 fTrackTypeRec(kTrackUndef),
146 fTrackTypeGen(kTrackUndef),
147 fAvgTrials(1),
148 fExternalWeight(1),
149 fRecEtaWindow(0.5),
9280dfa6 150 fMinJetPt(0),
151 fDeltaPhiWindow(20./180.*TMath::Pi()),
3b7ffecf 152 fh1Xsec(0x0),
153 fh1Trials(0x0),
154 fh1PtHard(0x0),
155 fh1PtHardNoW(0x0),
156 fh1PtHardTrials(0x0),
157 fh1NGenJets(0x0),
158 fh1NRecJets(0x0),
edfbe476 159 fh1PtTrackRec(0x0),
160 fh1SumPtTrackRec(0x0),
161 fh1SumPtTrackAreaRec(0x0),
cc0649e4 162 fh1PtJetsRecIn(0x0),
565584e8 163 fh1PtJetsLeadingRecIn(0x0),
cc0649e4 164 fh1PtTracksRecIn(0x0),
565584e8 165 fh1PtTracksLeadingRecIn(0x0),
cc0649e4 166 fh1PtTracksGenIn(0x0),
167 fh2NRecJetsPt(0x0),
168 fh2NRecTracksPt(0x0),
169 fh2JetsLeadingPhiEta(0x0),
565584e8 170 fh2JetsLeadingPhiPt(0x0),
cc0649e4 171 fh2TracksLeadingPhiEta(0x0),
565584e8 172 fh2TracksLeadingPhiPt(0x0),
cf049e94 173 fh2TracksLeadingJetPhiPt(0x0),
26fa06fb 174 fh2DijetDeltaPhiPt(0x0),
175 fh2DijetAsymPt(0x0),
176 fh2DijetAsymPtCut(0x0),
177 fh2DijetDeltaPhiDeltaEta(0x0),
178 fh2DijetPt2vsPt1(0x0),
179 fh2DijetDifvsSum(0x0),
180 fh1DijetMinv(0x0),
181 fh1DijetMinvCut(0x0),
3b7ffecf 182 fHistList(0x0)
183{
184
185 for(int i = 0;i < kMaxStep*2;++i){
186 fhnJetContainer[i] = 0;
187 }
188 for(int i = 0;i < kMaxJets;++i){
edfbe476 189 fh2PhiPt[i] = fh2PhiEta[i] = fh2FragRec[i] = fh2FragLnRec[i] = fh2FragGen[i] = fh2FragLnGen[i] = 0;
3b7ffecf 190 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
191 }
192 DefineOutput(1, TList::Class());
193}
194
195
196
197Bool_t AliAnalysisTaskJetSpectrum2::Notify()
198{
199 //
200 // Implemented Notify() to read the cross sections
201 // and number of trials from pyxsec.root
202 //
203
204 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
205 Float_t xsection = 0;
206 Float_t ftrials = 1;
207
208 fAvgTrials = 1;
209 if(tree){
210 TFile *curfile = tree->GetCurrentFile();
211 if (!curfile) {
212 Error("Notify","No current file");
213 return kFALSE;
214 }
215 if(!fh1Xsec||!fh1Trials){
216 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
217 return kFALSE;
218 }
219 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
220 fh1Xsec->Fill("<#sigma>",xsection);
221 // construct a poor man average trials
222 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
a221560c 223 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
3b7ffecf 224 }
225 return kTRUE;
226}
227
228void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
229{
230
231 //
232 // Create the output container
233 //
234
235
236 // Connect the AOD
237
238
239 MakeJetContainer();
240
241
242 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
243
244 OpenFile(1);
245 if(!fHistList)fHistList = new TList();
5cb0f438 246 fHistList->SetOwner(kTRUE);
3b7ffecf 247 Bool_t oldStatus = TH1::AddDirectoryStatus();
248 TH1::AddDirectory(kFALSE);
249
250 //
251 // Histogram
252
253 const Int_t nBinPt = 100;
254 Double_t binLimitsPt[nBinPt+1];
255 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
256 if(iPt == 0){
257 binLimitsPt[iPt] = 0.0;
258 }
259 else {// 1.0
edfbe476 260 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
3b7ffecf 261 }
262 }
263
edfbe476 264 const Int_t nBinPhi = 90;
3b7ffecf 265 Double_t binLimitsPhi[nBinPhi+1];
266 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
267 if(iPhi==0){
cc0649e4 268 binLimitsPhi[iPhi] = -1.*TMath::Pi();
3b7ffecf 269 }
270 else{
271 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
272 }
273 }
274
edfbe476 275
276
277 const Int_t nBinEta = 40;
278 Double_t binLimitsEta[nBinEta+1];
279 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
280 if(iEta==0){
281 binLimitsEta[iEta] = -2.0;
282 }
283 else{
cf049e94 284 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
edfbe476 285 }
286 }
287
3b7ffecf 288 const Int_t nBinFrag = 25;
289
290
291 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
292 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
293
294 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
295 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
296
297 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
298 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
299 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
300
301 fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
302 fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
303
edfbe476 304 fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
305 fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
306 fh1SumPtTrackAreaRec = new TH1F("fh1SumPtTrackAreaRec","Sum Rec track P_T #eta <0.9 / 1.8 * 2 * 0.4*0.4;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
cc0649e4 307
308 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
565584e8 309 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
cc0649e4 310 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
565584e8 311 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
fe77112a 312 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
edfbe476 313
cc0649e4 314 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
315 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
3b7ffecf 316 //
cc0649e4 317
565584e8 318 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
cc0649e4 319 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
565584e8 320 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
321 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
322 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
323 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
324 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
325 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
cc0649e4 326
3b7ffecf 327 for(int ij = 0;ij < kMaxJets;++ij){
328 fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
329 fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
330
edfbe476 331 fh2PhiPt[ij] = new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
332 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
333
cc0649e4 334 fh2PhiEta[ij] = new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
edfbe476 335 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
336
337
3b7ffecf 338 fh2FragRec[ij] = new TH2F(Form("fh2FragRec_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
339 nBinFrag,0.,1.,nBinPt,binLimitsPt);
340 fh2FragLnRec[ij] = new TH2F(Form("fh2FragLnRec_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
341 nBinFrag,0.,10.,nBinPt,binLimitsPt);
342
343 fh2FragGen[ij] = new TH2F(Form("fh2FragGen_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
344 nBinFrag,0.,1.,nBinPt,binLimitsPt);
345 fh2FragLnGen[ij] = new TH2F(Form("fh2FragLnGen_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
346 nBinFrag,0.,10.,nBinPt,binLimitsPt);
347 }
348
26fa06fb 349 // Dijet histograms
350
351 fh2DijetDeltaPhiPt = new TH2F("fh2DeltaPhiPt","Difference in the azimuthal angle;#Delta#phi;p_{T,1};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
352 fh2DijetAsymPt = new TH2F("fh2DijetAsym","Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
353 fh2DijetAsymPtCut = new TH2F("fh2DijetAsymCut","Pt asymmetry after delta phi cut;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
354 fh2DijetDeltaPhiDeltaEta = new TH2F("fh2DijetDeltaPhiDeltaEta","Difference in the azimuthal angle;#Delta#phi;Entries",180,0.,TMath::Pi(),20,-2.,2.);
355 fh2DijetPt2vsPt1 = new TH2F("fh2DijetPt2vsPt1","Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
356 fh2DijetDifvsSum = new TH2F("fh2DijetDifvsSum","Pt difference vs Pt sum;p_{T,1}+p_{T,2} (GeV/c);#Deltap_{T} (GeV/c)",400,0.,400.,150,0.,150.);
357 fh1DijetMinv = new TH1F("fh1DijetMinv","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
358 fh1DijetMinvCut = new TH1F("fh1DijetMinvCut","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
359
360
361
3b7ffecf 362
363 const Int_t saveLevel = 3; // large save level more histos
364 if(saveLevel>0){
365 fHistList->Add(fh1Xsec);
366 fHistList->Add(fh1Trials);
367 fHistList->Add(fh1PtHard);
368 fHistList->Add(fh1PtHardNoW);
369 fHistList->Add(fh1PtHardTrials);
cc0649e4 370 if(fBranchGen.Length()>0){
371 fHistList->Add(fh1NGenJets);
372 fHistList->Add(fh1PtTracksGenIn);
373 }
374 fHistList->Add(fh1PtJetsRecIn);
565584e8 375 fHistList->Add(fh1PtJetsLeadingRecIn);
cc0649e4 376 fHistList->Add(fh1PtTracksRecIn);
565584e8 377 fHistList->Add(fh1PtTracksLeadingRecIn);
3b7ffecf 378 fHistList->Add(fh1NRecJets);
edfbe476 379 fHistList->Add(fh1PtTrackRec);
380 fHistList->Add(fh1SumPtTrackRec);
381 fHistList->Add(fh1SumPtTrackAreaRec);
cc0649e4 382 fHistList->Add(fh2NRecJetsPt);
383 fHistList->Add(fh2NRecTracksPt);
384 fHistList->Add(fh2JetsLeadingPhiEta );
565584e8 385 fHistList->Add(fh2JetsLeadingPhiPt );
cc0649e4 386 fHistList->Add(fh2TracksLeadingPhiEta);
565584e8 387 fHistList->Add(fh2TracksLeadingPhiPt);
3b7ffecf 388 for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
389 for(int ij = 0;ij<kMaxJets;++ij){
390 fHistList->Add( fh1PtRecIn[ij]);
cc0649e4 391 if(fBranchGen.Length()>0){
392 fHistList->Add( fh1PtGenIn[ij]);
393 fHistList->Add( fh2FragGen[ij]);
394 fHistList->Add( fh2FragLnGen[ij]);
395 }
edfbe476 396 fHistList->Add( fh2PhiPt[ij]);
397 fHistList->Add( fh2PhiEta[ij]);
3b7ffecf 398 fHistList->Add( fh2FragRec[ij]);
399 fHistList->Add( fh2FragLnRec[ij]);
3b7ffecf 400 }
401 fHistList->Add(fhnCorrelation);
26fa06fb 402
403
404 fHistList->Add(fh2DijetDeltaPhiPt);
405 fHistList->Add(fh2DijetAsymPt);
406 fHistList->Add(fh2DijetAsymPtCut);
407 fHistList->Add(fh2DijetDeltaPhiDeltaEta);
408 fHistList->Add(fh2DijetPt2vsPt1);
409 fHistList->Add(fh2DijetDifvsSum);
410 fHistList->Add(fh1DijetMinv);
411 fHistList->Add(fh1DijetMinvCut);
3b7ffecf 412 }
413
414 // =========== Switch on Sumw2 for all histos ===========
415 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
416 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
417 if (h1){
418 h1->Sumw2();
419 continue;
420 }
421 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
422 if(hn)hn->Sumw2();
423 }
424 TH1::AddDirectory(oldStatus);
425}
426
427void AliAnalysisTaskJetSpectrum2::Init()
428{
429 //
430 // Initialization
431 //
432
3b7ffecf 433 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
434
435}
436
437void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
438{
b5cc0c6d 439
a221560c 440 if(!AliAnalysisHelperJetTasks::Selected()&&fUseGlobalSelection){
b5cc0c6d 441 // no selection by the service task, we continue
a221560c 442 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
b5cc0c6d 443 PostData(1, fHistList);
444 return;
445 }
446
3b7ffecf 447 //
448 // Execute analysis for current event
449 //
7fa8b2da 450 AliESDEvent *fESD = 0;
565584e8 451 if(fUseAODJetInput){
3b7ffecf 452 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
453 if(!fAOD){
565584e8 454 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
3b7ffecf 455 return;
456 }
457 // fethc the header
458 }
459 else{
460 // assume that the AOD is in the general output...
461 fAOD = AODEvent();
462 if(!fAOD){
463 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
464 return;
465 }
7fa8b2da 466 if(fDebug>0){
467 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
468 }
3b7ffecf 469 }
470
471
472
473
474 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
475
476
477 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
478
479 if(!aodH){
480 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
481 return;
482 }
483
484 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
485 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
486 if(!aodRecJets){
487 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
488 return;
489 }
490
491 // ==== General variables needed
492
493
494 // We use statice array, not to fragment the memory
495 AliAODJet genJets[kMaxJets];
496 Int_t nGenJets = 0;
497 AliAODJet recJets[kMaxJets];
498 Int_t nRecJets = 0;
499 ///////////////////////////
599396fa 500
3b7ffecf 501
502 Double_t eventW = 1;
503 Double_t ptHard = 0;
504 Double_t nTrials = 1; // Trials for MC trigger
505
506 if(fUseExternalWeightOnly){
507 eventW = fExternalWeight;
508 }
509
510 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
cf049e94 511 // if(fDebug>0)aodH->SetFillAOD(kFALSE);
3b7ffecf 512 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
513 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
514 // this is the part we only use when we have MC information
515 AliMCEvent* mcEvent = MCEvent();
516 // AliStack *pStack = 0;
517 if(!mcEvent){
518 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
519 return;
520 }
521 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
522 Int_t iCount = 0;
523 if(pythiaGenHeader){
524 nTrials = pythiaGenHeader->Trials();
525 ptHard = pythiaGenHeader->GetPtHard();
526 int iProcessType = pythiaGenHeader->ProcessType();
527 // 11 f+f -> f+f
528 // 12 f+barf -> f+barf
529 // 13 f+barf -> g+g
530 // 28 f+g -> f+g
531 // 53 g+g -> f+barf
532 // 68 g+g -> g+g
533 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
534 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
535
536 // fetch the pythia generated jets only to be used here
537 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
538 AliAODJet pythiaGenJets[kMaxJets];
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
3b7ffecf 545 if(fBranchGen.Length()==0){
9280dfa6 546 /*
547 if(fLimitGenJetEta){
548 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
549 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
550 }
551 */
552 if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
3b7ffecf 553 // if we have MC particles and we do not read from the aod branch
554 // use the pythia jets
555 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
9280dfa6 556 iCount++;
3b7ffecf 557 }
3b7ffecf 558 }
559 }
560 if(fBranchGen.Length()==0)nGenJets = iCount;
561 }// (fAnalysisType&kMCESD)==kMCESD)
562
563
564 // we simply fetch the tracks/mc particles as a list of AliVParticles
565
566 TList recParticles;
567 TList genParticles;
568
565584e8 569
570
571
572 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
573 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
574 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
575 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
3b7ffecf 576
cc0649e4 577
3b7ffecf 578 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
579 fh1PtHard->Fill(ptHard,eventW);
580 fh1PtHardNoW->Fill(ptHard,1);
581 fh1PtHardTrials->Fill(ptHard,nTrials);
582
583 // If we set a second branch for the input jets fetch this
584 if(fBranchGen.Length()>0){
585 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
586 if(aodGenJets){
587 Int_t iCount = 0;
588 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
589 if(iCount>=kMaxJets)continue;
590 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
591 if(!tmp)continue;
5301f754 592 /*
3b7ffecf 593 if(fLimitGenJetEta){
594 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
595 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
596 }
5301f754 597 */
9280dfa6 598 if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
3b7ffecf 599 genJets[iCount] = *tmp;
600 iCount++;
601 }
602 nGenJets = iCount;
603 }
604 else{
5cb0f438 605 if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
606 if(fDebug>2)fAOD->Print();
3b7ffecf 607 }
608 }
609
610 fh1NGenJets->Fill(nGenJets);
611 // We do not want to exceed the maximum jet number
612 nGenJets = TMath::Min(nGenJets,kMaxJets);
613
614 // Fetch the reconstructed jets...
615
616 nRecJets = aodRecJets->GetEntries();
617
cc0649e4 618 nRecJets = aodRecJets->GetEntries();
619 fh1NRecJets->Fill(nRecJets);
620
621 // Do something with the all rec jets
622 Int_t nRecOver = nRecJets;
3b7ffecf 623
cc0649e4 624 // check that the jets are sorted
625 Float_t ptOld = 999999;
626 for(int ir = 0;ir < nRecJets;ir++){
627 AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
628 Float_t tmpPt = tmp->Pt();
629 if(tmpPt>ptOld){
3568c3d6 630 Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
cc0649e4 631 }
632 ptOld = tmpPt;
633 }
3b7ffecf 634
cc0649e4 635
636 if(nRecOver>0){
637 TIterator *recIter = aodRecJets->MakeIterator();
638 AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
639 Float_t pt = tmpRec->Pt();
640 if(tmpRec){
641 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
642 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
643 while(pt<ptCut&&tmpRec){
644 nRecOver--;
645 tmpRec = (AliAODJet*)(recIter->Next());
646 if(tmpRec){
647 pt = tmpRec->Pt();
648 }
649 }
650 if(nRecOver<=0)break;
651 fh2NRecJetsPt->Fill(ptCut,nRecOver);
652 }
653 }
654 recIter->Reset();
655 AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
656 Float_t phi = leading->Phi();
657 if(phi<0)phi+=TMath::Pi()*2.;
658 Float_t eta = leading->Eta();
cf049e94 659 pt = leading->Pt();
565584e8 660 while((tmpRec = (AliAODJet*)(recIter->Next()))){
cc0649e4 661 Float_t tmpPt = tmpRec->Pt();
662 fh1PtJetsRecIn->Fill(tmpPt);
565584e8 663 if(tmpRec==leading){
664 fh1PtJetsLeadingRecIn->Fill(tmpPt);
665 continue;
666 }
cc0649e4 667 // correlation
668 Float_t tmpPhi = tmpRec->Phi();
669
670 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
671 Float_t dPhi = phi - tmpRec->Phi();
672 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
673 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
674 Float_t dEta = eta - tmpRec->Eta();
675 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
565584e8 676 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
cc0649e4 677 }
678 delete recIter;
679 }
680
681 Int_t nTrackOver = recParticles.GetSize();
682 // do the same for tracks and jets
683 if(nTrackOver>0){
684 TIterator *recIter = recParticles.MakeIterator();
685 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
686 Float_t pt = tmpRec->Pt();
687 // Printf("Leading track p_t %3.3E",pt);
688 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
689 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
690 while(pt<ptCut&&tmpRec){
691 nTrackOver--;
692 tmpRec = (AliVParticle*)(recIter->Next());
693 if(tmpRec){
694 pt = tmpRec->Pt();
695 }
696 }
697 if(nTrackOver<=0)break;
698 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
699 }
700
701 recIter->Reset();
702 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
703 Float_t phi = leading->Phi();
704 if(phi<0)phi+=TMath::Pi()*2.;
705 Float_t eta = leading->Eta();
cf049e94 706 pt = leading->Pt();
cc0649e4 707 while((tmpRec = (AliVParticle*)(recIter->Next()))){
708 Float_t tmpPt = tmpRec->Pt();
709 fh1PtTracksRecIn->Fill(tmpPt);
565584e8 710 if(tmpRec==leading){
711 fh1PtTracksLeadingRecIn->Fill(tmpPt);
712 continue;
713 }
cc0649e4 714 // correlation
715 Float_t tmpPhi = tmpRec->Phi();
716
717 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
718 Float_t dPhi = phi - tmpRec->Phi();
719 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
720 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
721 Float_t dEta = eta - tmpRec->Eta();
722 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
565584e8 723 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
cc0649e4 724 }
725 delete recIter;
726 }
727
728 if(genParticles.GetSize()){
729 TIterator *genIter = genParticles.MakeIterator();
730 AliVParticle *tmpGen = 0;
731 while((tmpGen = (AliVParticle*)(genIter->Next()))){
732 if(TMath::Abs(tmpGen->Eta())<0.9){
733 Float_t tmpPt = tmpGen->Pt();
734 fh1PtTracksGenIn->Fill(tmpPt);
735 }
736 }
737 delete genIter;
738 }
739
3b7ffecf 740 nRecJets = TMath::Min(nRecJets,kMaxJets);
9280dfa6 741
742 Int_t iCountRec = 0;
3b7ffecf 743 for(int ir = 0;ir < nRecJets;++ir){
744 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
745 if(!tmp)continue;
9280dfa6 746 if(tmp->Pt()<fMinJetPt)continue;
3b7ffecf 747 recJets[ir] = *tmp;
9280dfa6 748 iCountRec++;
3b7ffecf 749 }
9280dfa6 750 nRecJets = iCountRec;
3b7ffecf 751
752
753 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
754 // Relate the jets
755 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
756 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
757
cc0649e4 758
759 for(int i = 0;i<kMaxJets;++i){
3b7ffecf 760 iGenIndex[i] = iRecIndex[i] = -1;
761 }
762
3b7ffecf 763 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
764 iGenIndex,iRecIndex,fDebug);
765 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
766
767 if(fDebug){
768 for(int i = 0;i<kMaxJets;++i){
769 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
770 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
771 }
772 }
773
774
775
776
777 Double_t container[6];
778
779 // loop over generated jets
780
781 // radius; tmp, get from the jet header itself
782 // at some point, how todeal woht FastJet on MC?
783 Float_t radiusGen =0.4;
784 Float_t radiusRec =0.4;
785
786 for(int ig = 0;ig < nGenJets;++ig){
787 Double_t ptGen = genJets[ig].Pt();
788 Double_t phiGen = genJets[ig].Phi();
789 if(phiGen<0)phiGen+=TMath::Pi()*2.;
790 Double_t etaGen = genJets[ig].Eta();
791
792 container[3] = ptGen;
793 container[4] = etaGen;
794 container[5] = phiGen;
795 fhnJetContainer[kStep0]->Fill(&container[3],eventW);
796 Int_t ir = iRecIndex[ig];
797
798 if(TMath::Abs(etaGen)<fRecEtaWindow){
799 fhnJetContainer[kStep1]->Fill(&container[3],eventW);
800 fh1PtGenIn[ig]->Fill(ptGen,eventW);
801 // fill the fragmentation function
802 for(int it = 0;it<genParticles.GetEntries();++it){
803 AliVParticle *part = (AliVParticle*)genParticles.At(it);
804 if(genJets[ig].DeltaR(part)<radiusGen){
805 Float_t z = part->Pt()/ptGen;
806 Float_t lnz = -1.*TMath::Log(z);
807 fh2FragGen[ig]->Fill(z,ptGen,eventW);
808 fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
809 }
810 }
811 if(ir>=0&&ir<nRecJets){
812 fhnJetContainer[kStep3]->Fill(&container[3],eventW);
813 }
814 }
815
816 if(ir>=0&&ir<nRecJets){
817 fhnJetContainer[kStep2]->Fill(&container[3],eventW);
818 Double_t etaRec = recJets[ir].Eta();
819 if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
820 }
821 }// loop over generated jets
822
7fa8b2da 823
edfbe476 824 Float_t sumPt = 0;
825 for(int it = 0;it<recParticles.GetEntries();++it){
826 AliVParticle *part = (AliVParticle*)recParticles.At(it);
827 // fill sum pt and P_t of all paritles
828 if(TMath::Abs(part->Eta())<0.9){
829 Float_t pt = part->Pt();
830 fh1PtTrackRec->Fill(pt,eventW);
831 sumPt += pt;
832 }
833 }
834 fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
835 fh1SumPtTrackRec->Fill(sumPt,eventW);
836
cf049e94 837
3b7ffecf 838 // loop over reconstructed jets
839 for(int ir = 0;ir < nRecJets;++ir){
840 Double_t etaRec = recJets[ir].Eta();
841 Double_t ptRec = recJets[ir].Pt();
842 Double_t phiRec = recJets[ir].Phi();
843 if(phiRec<0)phiRec+=TMath::Pi()*2.;
26fa06fb 844
845
846 // do something with dijets...
847 if(ir==1){
848 Double_t etaRec1 = recJets[0].Eta();
849 Double_t ptRec1 = recJets[0].Pt();
850 Double_t phiRec1 = recJets[0].Phi();
851 if(phiRec1<0)phiRec1+=TMath::Pi()*2.;
852
853
854 if(TMath::Abs(etaRec1)<fRecEtaWindow
855 &&TMath::Abs(etaRec)<fRecEtaWindow){
856
857 Float_t deltaPhi = phiRec1 - phiRec;
858
859 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
860 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
861 deltaPhi = TMath::Abs(deltaPhi);
862 fh2DijetDeltaPhiPt->Fill(deltaPhi,ptRec1);
863 Float_t asym = (ptRec1-ptRec)/(ptRec1+ptRec);
864 fh2DijetAsymPt->Fill(asym,ptRec1);
865 fh2DijetDeltaPhiDeltaEta->Fill(deltaPhi,etaRec1-etaRec);
866 fh2DijetPt2vsPt1->Fill(ptRec1,ptRec);
867 fh2DijetDifvsSum->Fill(ptRec1+ptRec,ptRec1-ptRec);
868 Float_t minv = 2.*(recJets[0].P()*recJets[1].P()-
869 recJets[0].Px()*recJets[1].Px()-
870 recJets[0].Py()*recJets[1].Py()-
871 recJets[0].Pz()*recJets[1].Py());
872 minv = TMath::Sqrt(minv);
873 // with mass == 0;
874
875 fh1DijetMinv->Fill(minv);
876 if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
877 fh1DijetMinvCut->Fill(minv);
878 fh2DijetAsymPtCut->Fill(asym,ptRec1);
879 }
880 }
881 }
882
883
3b7ffecf 884 container[0] = ptRec;
885 container[1] = etaRec;
886 container[2] = phiRec;
887
1277f815 888 if(ptRec>30.&&fDebug>0){
7fa8b2da 889 // need to cast to int, otherwise the printf overwrites
890 Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
1277f815 891 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
892 if(fESD)Printf("ESDEvent GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
cf049e94 893 // aodH->SetFillAOD(kTRUE);
7fa8b2da 894 fAOD->GetHeader()->Print();
a221560c 895 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
7fa8b2da 896 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
897 AliAODTrack *tr = fAOD->GetTrack(it);
898 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
899 tr->Print();
38ecb6a5 900 // tr->Dump();
7fa8b2da 901 if(fESD){
902 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
903 esdTr->Print("");
bb3d13aa 904 // esdTr->Dump();
7fa8b2da 905 }
906 }
907 }
908
909
3b7ffecf 910 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
911 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
912 if(TMath::Abs(etaRec)<fRecEtaWindow){
913 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
914 fh1PtRecIn[ir]->Fill(ptRec,eventW);
915 // fill the fragmentation function
916 for(int it = 0;it<recParticles.GetEntries();++it){
917 AliVParticle *part = (AliVParticle*)recParticles.At(it);
edfbe476 918 Float_t eta = part->Eta();
919 if(TMath::Abs(eta)<0.9){
920 Float_t phi = part->Phi();
921 if(phi<0)phi+=TMath::Pi()*2.;
922 Float_t dPhi = phi - phiRec;
923 Float_t dEta = eta - etaRec;
26fa06fb 924 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
925 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
edfbe476 926 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
927 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
928 }
3b7ffecf 929 if(recJets[ir].DeltaR(part)<radiusRec){
930 Float_t z = part->Pt()/ptRec;
931 Float_t lnz = -1.*TMath::Log(z);
932 fh2FragRec[ir]->Fill(z,ptRec,eventW);
933 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
934 }
935 }
936 }
937
938
939 // Fill Correlation
940 Int_t ig = iGenIndex[ir];
941 if(ig>=0 && ig<nGenJets){
942 fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
943 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
944 Double_t ptGen = genJets[ig].Pt();
945 Double_t phiGen = genJets[ig].Phi();
946 if(phiGen<0)phiGen+=TMath::Pi()*2.;
947 Double_t etaGen = genJets[ig].Eta();
948
949 container[3] = ptGen;
950 container[4] = etaGen;
951 container[5] = phiGen;
952
953 //
954 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
955 //
956
957 if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
958 if(TMath::Abs(etaRec)<fRecEtaWindow){
959 fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
960 fhnCorrelation->Fill(container);
961 }// if etarec in window
962
963 }
964 ////////////////////////////////////////////////////
965 else{
966
967 }
968 }// loop over reconstructed jets
969
970
971 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
972 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
973 PostData(1, fHistList);
974}
975
976
977void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
978 //
979 // Create the particle container for the correction framework manager and
980 // link it
981 //
982 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
383d44fd 983 const Double_t kPtmin = 0.0, kPtmax = 160.; // we do not want to have empty bins at the beginning...
3b7ffecf 984 const Double_t kEtamin = -3.0, kEtamax = 3.0;
985 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
986
987 // can we neglect migration in eta and phi?
988 // phi should be no problem since we cover full phi and are phi symmetric
989 // eta migration is more difficult due to needed acceptance correction
990 // in limited eta range
991
992
993 //arrays for the number of bins in each dimension
994 Int_t iBin[kNvar];
383d44fd 995 iBin[0] = 160; //bins in pt
3b7ffecf 996 iBin[1] = 1; //bins in eta
997 iBin[2] = 1; // bins in phi
998
999 //arrays for lower bounds :
1000 Double_t* binEdges[kNvar];
1001 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1002 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1003
1004 //values for bin lower bounds
1005 // 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);
a923bd34 1006 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
3b7ffecf 1007 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1008 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1009
1010
1011 for(int i = 0;i < kMaxStep*2;++i){
1012 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
1013 for (int k=0; k<kNvar; k++) {
1014 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1015 }
1016 }
1017 //create correlation matrix for unfolding
1018 Int_t thnDim[2*kNvar];
1019 for (int k=0; k<kNvar; k++) {
1020 //first half : reconstructed
1021 //second half : MC
1022 thnDim[k] = iBin[k];
1023 thnDim[k+kNvar] = iBin[k];
1024 }
1025
1026 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1027 for (int k=0; k<kNvar; k++) {
1028 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1029 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1030 }
1031 fhnCorrelation->Sumw2();
1032
1033 // Add a histogram for Fake jets
1034 // thnDim[3] = AliPID::kSPECIES;
1035 // fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
1036 // for(Int_t idim = 0; idim < kNvar; idim++)
1037 // fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
5cb0f438 1038 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1039 delete [] binEdges[ivar];
1040
3b7ffecf 1041}
1042
1043void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1044{
1045// Terminate analysis
1046//
1047 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1048}
1049
1050
1051Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1052
565584e8 1053 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
3b7ffecf 1054
1055 Int_t iCount = 0;
565584e8 1056 if(type==kTrackAOD){
3b7ffecf 1057 AliAODEvent *aod = 0;
565584e8 1058 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1059 else aod = AODEvent();
3b7ffecf 1060 if(!aod){
1061 return iCount;
1062 }
1063 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1064 AliAODTrack *tr = aod->GetTrack(it);
1065 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
cc0649e4 1066 if(TMath::Abs(tr->Eta())>0.9)continue;
38ecb6a5 1067 if(fDebug>0){
1068 if(tr->Pt()>20){
1069 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1070 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
1071 tr->Print();
1072 // tr->Dump();
1073 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1074 if(fESD){
1075 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1076 esdTr->Print("");
bb3d13aa 1077 // esdTr->Dump();
38ecb6a5 1078 }
1079 }
1080 }
3b7ffecf 1081 list->Add(tr);
1082 iCount++;
1083 }
1084 }
1085 else if (type == kTrackKineAll||type == kTrackKineCharged){
1086 AliMCEvent* mcEvent = MCEvent();
1087 if(!mcEvent)return iCount;
1088 // we want to have alivpartilces so use get track
1089 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1090 if(!mcEvent->IsPhysicalPrimary(it))continue;
1091 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1092 if(type == kTrackKineAll){
1093 list->Add(part);
1094 iCount++;
1095 }
1096 else if(type == kTrackKineCharged){
1097 if(part->Particle()->GetPDG()->Charge()==0)continue;
1098 list->Add(part);
1099 iCount++;
1100 }
1101 }
1102 }
565584e8 1103 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1104 AliAODEvent *aod = 0;
5301f754 1105 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
565584e8 1106 else aod = AODEvent();
5010a3f7 1107 if(!aod)return iCount;
3b7ffecf 1108 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1109 if(!tca)return iCount;
1110 for(int it = 0;it < tca->GetEntriesFast();++it){
1111 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1112 if(!part->IsPhysicalPrimary())continue;
c4631cdb 1113 if(type == kTrackAODMCAll){
3b7ffecf 1114 list->Add(part);
1115 iCount++;
1116 }
565584e8 1117 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
3b7ffecf 1118 if(part->Charge()==0)continue;
565584e8 1119 if(kTrackAODMCCharged){
1120 list->Add(part);
1121 }
1122 else {
1123 if(TMath::Abs(part->Eta())>0.9)continue;
1124 list->Add(part);
1125 }
3b7ffecf 1126 iCount++;
1127 }
1128 }
1129 }// AODMCparticle
cc0649e4 1130 list->Sort();
c4631cdb 1131 return iCount;
3b7ffecf 1132
1133}