]> 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
Some clean up in addTask macros, added centrality selection and event classes to...
[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 "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 "AliAODJetEventBackground.h"
52#include "AliAODMCParticle.h"
53#include "AliMCEventHandler.h"
54#include "AliMCEvent.h"
55#include "AliStack.h"
56#include "AliGenPythiaEventHeader.h"
57#include "AliJetKineReaderHeader.h"
58#include "AliGenCocktailEventHeader.h"
59#include "AliInputEventHandler.h"
60
61
62#include "AliAnalysisHelperJetTasks.h"
63
64ClassImp(AliAnalysisTaskJetSpectrum2)
65
66AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
67 fJetHeaderRec(0x0),
68 fJetHeaderGen(0x0),
69 fAOD(0x0),
70 fhnCorrelation(0x0),
71 fhnCorrelationPhiZRec(0x0),
72 f1PtScale(0x0),
73 fBranchRec("jets"),
74 fBranchGen(""),
75 fBranchBkg(""),
76 fUseAODJetInput(kFALSE),
77 fUseAODTrackInput(kFALSE),
78 fUseAODMCInput(kFALSE),
79 fUseGlobalSelection(kFALSE),
80 fUseExternalWeightOnly(kFALSE),
81 fLimitGenJetEta(kFALSE),
82 fBkgSubtraction(kFALSE),
83 fFillCorrBkg(0),
84 fFilterMask(0),
85 fEventSelectionMask(0),
86 fAnalysisType(0),
87 fTrackTypeRec(kTrackUndef),
88 fTrackTypeGen(kTrackUndef),
89 fEventClass(0),
90 fAvgTrials(1),
91 fExternalWeight(1),
92 fRecEtaWindow(0.5),
93 fMinJetPt(0),
94 fDeltaPhiWindow(20./180.*TMath::Pi()),
95 fh1Xsec(0x0),
96 fh1Trials(0x0),
97 fh1PtHard(0x0),
98 fh1PtHardNoW(0x0),
99 fh1PtHardTrials(0x0),
100 fh1ZVtx(0x0),
101 fh1NGenJets(0x0),
102 fh1NRecJets(0x0),
103 fh1PtTrackRec(0x0),
104 fh1SumPtTrackRec(0x0),
105 fh1SumPtTrackAreaRec(0x0),
106 fh1TmpRho(0x0),
107 fh1PtJetsRecIn(0x0),
108 fh1PtJetsLeadingRecIn(0x0),
109 fh1PtTracksRecIn(0x0),
110 fh1PtTracksLeadingRecIn(0x0),
111 fh1PtTracksGenIn(0x0),
112 fh2NRecJetsPt(0x0),
113 fh2NRecTracksPt(0x0),
114 fh2JetsLeadingPhiEta(0x0),
115 fh2JetsLeadingPhiPt(0x0),
116 fh2TracksLeadingPhiEta(0x0),
117 fh2TracksLeadingPhiPt(0x0),
118 fh2TracksLeadingJetPhiPt(0x0),
119 fh2JetPtJetPhi(0x0),
120 fh2TrackPtTrackPhi(0x0),
121 fh2RelPtFGen(0x0),
122 fh2DijetDeltaPhiPt(0x0),
123 fh2DijetAsymPt(0x0),
124 fh2DijetAsymPtCut(0x0),
125 fh2DijetDeltaPhiDeltaEta(0x0),
126 fh2DijetPt2vsPt1(0x0),
127 fh2DijetDifvsSum(0x0),
128 fh1DijetMinv(0x0),
129 fh1DijetMinvCut(0x0),
130 fh1Bkg1(0x0),
131 fh1Bkg2(0x0),
132 fh1Bkg3(0x0),
133 fh1Sigma1(0x0),
134 fh1Sigma2(0x0),
135 fh1Sigma3(0x0),
136 fh1Area1(0x0),
137 fh1Area2(0x0),
138 fh1Area3(0x0),
139 fh1Ptjet(0x0),
140 fh1Ptjetsub1(0x0),
141 fh1Ptjetsub2(0x0),
142 fh1Ptjetsub3(0x0),
143 fh1Ptjethardest(0x0),
144 fh1Ptjetsubhardest1(0x0),
145 fh1Ptjetsubhardest2(0x0),
146 fh1Ptjetsubhardest3(0x0),
147 fh2Rhovspthardest1(0x0),
148 fh2Rhovspthardest2(0x0),
149 fh2Rhovspthardest3(0x0),
150 fh2Errorvspthardest1(0x0),
151 fh2Errorvspthardest2(0x0),
152 fh2Errorvspthardest3(0x0),
153 fHistList(0x0)
154{
155 for(int i = 0;i < kMaxStep*2;++i){
156 fhnJetContainer[i] = 0;
157 }
158 for(int i = 0;i < kMaxJets;++i){
159 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
160
161 fh2PhiPt[i] = 0;
162 fh2PhiEta[i] = 0;
163 fh2RhoPtRec[i] = 0;
164 fh2RhoPtGen[i] = 0;
165 fh2PsiPtGen[i] = 0;
166 fh2PsiPtRec[i] = 0;
167 fh2FragRec[i] = 0;
168 fh2FragLnRec[i] = 0;
169 fh2FragGen[i] = 0;
170 fh2FragLnGen[i] = 0;
171 }
172
173}
174
175AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
176 AliAnalysisTaskSE(name),
177 fJetHeaderRec(0x0),
178 fJetHeaderGen(0x0),
179 fAOD(0x0),
180 fhnCorrelation(0x0),
181 fhnCorrelationPhiZRec(0x0),
182 f1PtScale(0x0),
183 fBranchRec("jets"),
184 fBranchGen(""),
185 fBranchBkg(""),
186 fUseAODJetInput(kFALSE),
187 fUseAODTrackInput(kFALSE),
188 fUseAODMCInput(kFALSE),
189 fUseGlobalSelection(kFALSE),
190 fUseExternalWeightOnly(kFALSE),
191 fLimitGenJetEta(kFALSE),
192 fBkgSubtraction(kFALSE),
193 fFillCorrBkg(0),
194 fFilterMask(0),
195 fEventSelectionMask(0),
196 fAnalysisType(0),
197 fTrackTypeRec(kTrackUndef),
198 fTrackTypeGen(kTrackUndef),
199 fEventClass(0),
200 fAvgTrials(1),
201 fExternalWeight(1),
202 fRecEtaWindow(0.5),
203 fMinJetPt(0),
204 fDeltaPhiWindow(20./180.*TMath::Pi()),
205 fh1Xsec(0x0),
206 fh1Trials(0x0),
207 fh1PtHard(0x0),
208 fh1PtHardNoW(0x0),
209 fh1PtHardTrials(0x0),
210 fh1ZVtx(0x0),
211 fh1NGenJets(0x0),
212 fh1NRecJets(0x0),
213 fh1PtTrackRec(0x0),
214 fh1SumPtTrackRec(0x0),
215 fh1SumPtTrackAreaRec(0x0),
216 fh1TmpRho(0x0),
217 fh1PtJetsRecIn(0x0),
218 fh1PtJetsLeadingRecIn(0x0),
219 fh1PtTracksRecIn(0x0),
220 fh1PtTracksLeadingRecIn(0x0),
221 fh1PtTracksGenIn(0x0),
222 fh2NRecJetsPt(0x0),
223 fh2NRecTracksPt(0x0),
224 fh2JetsLeadingPhiEta(0x0),
225 fh2JetsLeadingPhiPt(0x0),
226 fh2TracksLeadingPhiEta(0x0),
227 fh2TracksLeadingPhiPt(0x0),
228 fh2TracksLeadingJetPhiPt(0x0),
229 fh2JetPtJetPhi(0x0),
230 fh2TrackPtTrackPhi(0x0),
231 fh2RelPtFGen(0x0),
232 fh2DijetDeltaPhiPt(0x0),
233 fh2DijetAsymPt(0x0),
234 fh2DijetAsymPtCut(0x0),
235 fh2DijetDeltaPhiDeltaEta(0x0),
236 fh2DijetPt2vsPt1(0x0),
237 fh2DijetDifvsSum(0x0),
238 fh1DijetMinv(0x0),
239 fh1DijetMinvCut(0x0),
240 fh1Bkg1(0x0),
241 fh1Bkg2(0x0),
242 fh1Bkg3(0x0),
243 fh1Sigma1(0x0),
244 fh1Sigma2(0x0),
245 fh1Sigma3(0x0),
246 fh1Area1(0x0),
247 fh1Area2(0x0),
248 fh1Area3(0x0),
249 fh1Ptjet(0x0),
250 fh1Ptjetsub1(0x0),
251 fh1Ptjetsub2(0x0),
252 fh1Ptjetsub3(0x0),
253 fh1Ptjethardest(0x0),
254 fh1Ptjetsubhardest1(0x0),
255 fh1Ptjetsubhardest2(0x0),
256 fh1Ptjetsubhardest3(0x0),
257 fh2Rhovspthardest1(0x0),
258 fh2Rhovspthardest2(0x0),
259 fh2Rhovspthardest3(0x0),
260 fh2Errorvspthardest1(0x0),
261 fh2Errorvspthardest2(0x0),
262 fh2Errorvspthardest3(0x0),
263 fHistList(0x0)
264{
265
266 for(int i = 0;i < kMaxStep*2;++i){
267 fhnJetContainer[i] = 0;
268 }
269 for(int i = 0;i < kMaxJets;++i){
270 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
271
272 fh2PhiPt[i] = 0;
273 fh2PhiEta[i] = 0;
274 fh2RhoPtRec[i] = 0;
275 fh2RhoPtGen[i] = 0;
276 fh2PsiPtGen[i] = 0;
277 fh2PsiPtRec[i] = 0;
278 fh2FragRec[i] = 0;
279 fh2FragLnRec[i] = 0;
280 fh2FragGen[i] = 0;
281 fh2FragLnGen[i] = 0;
282 }
283 DefineOutput(1, TList::Class());
284}
285
286
287
288Bool_t AliAnalysisTaskJetSpectrum2::Notify()
289{
290 //
291 // Implemented Notify() to read the cross sections
292 // and number of trials from pyxsec.root
293 //
294
295 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
296 Float_t xsection = 0;
297 Float_t ftrials = 1;
298
299 fAvgTrials = 1;
300 if(tree){
301 TFile *curfile = tree->GetCurrentFile();
302 if (!curfile) {
303 Error("Notify","No current file");
304 return kFALSE;
305 }
306 if(!fh1Xsec||!fh1Trials){
307 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
308 return kFALSE;
309 }
310 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
311 fh1Xsec->Fill("<#sigma>",xsection);
312 // construct a poor man average trials
313 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
314 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
315 }
316 return kTRUE;
317}
318
319void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
320{
321
322 //
323 // Create the output container
324 //
325
326
327 // Connect the AOD
328
329
330 MakeJetContainer();
331
332
333 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
334
335 OpenFile(1);
336 if(!fHistList)fHistList = new TList();
337 fHistList->SetOwner(kTRUE);
338 Bool_t oldStatus = TH1::AddDirectoryStatus();
339 TH1::AddDirectory(kFALSE);
340
341 //
342 // Histogram
343
344 const Int_t nBinPt = 320;
345 Double_t binLimitsPt[nBinPt+1];
346 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
347 if(iPt == 0){
348 binLimitsPt[iPt] = 0.0;
349 }
350 else {// 1.0
351 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
352 }
353 }
354
355 const Int_t nBinPhi = 90;
356 Double_t binLimitsPhi[nBinPhi+1];
357 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
358 if(iPhi==0){
359 binLimitsPhi[iPhi] = -1.*TMath::Pi();
360 }
361 else{
362 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
363 }
364 }
365
366
367 const Int_t nBinPhi2 = 360;
368 Double_t binLimitsPhi2[nBinPhi2+1];
369 for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
370 if(iPhi2==0){
371 binLimitsPhi2[iPhi2] = 0.;
372 }
373 else{
374 binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
375 }
376 }
377
378
379
380 const Int_t nBinEta = 40;
381 Double_t binLimitsEta[nBinEta+1];
382 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
383 if(iEta==0){
384 binLimitsEta[iEta] = -2.0;
385 }
386 else{
387 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
388 }
389 }
390
391 const Int_t nBinFrag = 25;
392
393
394 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
395 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
396
397 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
398 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
399
400 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
401 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
402 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
403
404 fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
405
406 fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
407 fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
408
409 fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
410 fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
411 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);
412
413 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
414 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
415 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
416 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
417 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
418
419 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
420 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
421 //
422
423 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
424 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
425 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
426 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
427 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
428 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
429 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
430 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
431
432 fh2JetPtJetPhi = new TH2F("fh2JetPtJetPhi","Reconstructed jet phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
433 fh2TrackPtTrackPhi = new TH2F("fh2TrackPtTrackPhi","Reconstructed track phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
434
435 fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};#Delta p_{T}/p_{T,Gen}",nBinPt,binLimitsPt,120,-1.2,1.2);
436
437 for(int ij = 0;ij < kMaxJets;++ij){
438 fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
439 fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
440
441 fh2PhiPt[ij] = new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
442 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
443
444 fh2PhiEta[ij] = new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
445 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
446
447 fh2RhoPtRec[ij] = new TH2F(Form("fh2RhoPtRec_j%d",ij),"jet shape rho for jets;r;p_{T,rec};",
448 40,0.,1.,nBinPt,binLimitsPt);
449 fh2PsiPtRec[ij] = new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
450 40,0.,2.,nBinPt,binLimitsPt);
451
452 fh2RhoPtGen[ij] = new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
453 40,0.,2.,nBinPt,binLimitsPt);
454 fh2PsiPtGen[ij] = new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
455 40,0.,2.,nBinPt,binLimitsPt);
456 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
457
458
459 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",
460 nBinFrag,0.,1.,nBinPt,binLimitsPt);
461 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",
462 nBinFrag,0.,10.,nBinPt,binLimitsPt);
463
464 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",
465 nBinFrag,0.,1.,nBinPt,binLimitsPt);
466 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",
467 nBinFrag,0.,10.,nBinPt,binLimitsPt);
468 }
469
470 // Dijet histograms
471
472 fh2DijetDeltaPhiPt = new TH2F("fh2DijetDeltaPhiPt","Difference in the azimuthal angle;#Delta#phi;p_{T,1};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
473 fh2DijetAsymPt = new TH2F("fh2DijetAsym","Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
474 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);
475 fh2DijetDeltaPhiDeltaEta = new TH2F("fh2DijetDeltaPhiDeltaEta","Difference in the azimuthal angle;#Delta#phi;Entries",180,0.,TMath::Pi(),20,-2.,2.);
476 fh2DijetPt2vsPt1 = new TH2F("fh2DijetPt2vsPt1","Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
477 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.);
478 fh1DijetMinv = new TH1F("fh1DijetMinv","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
479 fh1DijetMinvCut = new TH1F("fh1DijetMinvCut","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
480 //background histograms
481 if(fBkgSubtraction){
482 fh1Bkg1 = new TH1F("fh1Bkg1","Background estimate 1",100,0.,10.);
483 fh1Bkg2 = new TH1F("fh1Bkg2","Background estimate 2",100,0.,10.);
484 fh1Bkg3 = new TH1F("fh1Bkg3","Background estimate 3",100,0.,10.);
485 fh1Sigma1 = new TH1F("fh1Sigma1","Background fluctuations 1",100,0.,10.);
486 fh1Sigma2 = new TH1F("fh1Sigma2","Background fluctuations 2",100,0.,10.);
487 fh1Sigma3 = new TH1F("fh1Sigma3","Background fluctuations 3",100,0.,10.);
488 fh1Area1 = new TH1F("fh1Area1","Background mean area 1",50,0.,5.);
489 fh1Area2 = new TH1F("fh1Area2","Background mean area 2",50,0.,5.);
490 fh1Area3 = new TH1F("fh1Area3","Background mean area 3",50,0.,5.);
491 fh1Ptjet = new TH1F("fh1Ptjet","Jet spectrum",100,0.,200.);
492 fh1Ptjetsub1 = new TH1F("fh1Ptjetsub1","Subtracted spectrum 1",50,0.,200.);
493 fh1Ptjetsub2 = new TH1F("fh1Ptjetsub2","Subtracted spectrum 2",50,0.,200.);
494 fh1Ptjetsub3 = new TH1F("fh1Ptjetsub3","Subtracted spectrum 3",50,0.,200.);
495 fh1Ptjethardest = new TH1F("fh1Ptjethardest","Hardest jet spectrum",50,0.,200.);
496 fh1Ptjetsubhardest1 = new TH1F("fh1Pthardestsub1","Subtracted hardest jet spectrum 1",100,0.,200.);
497 fh1Ptjetsubhardest2 = new TH1F("fh1Pthardestsub2","Subtracted hardest jet spectrum 2",100,0.,200.);
498 fh1Ptjetsubhardest3 = new TH1F("fh1Pthardestsub3","Subtracted hardest jet spectrum 3",100,0.,200.);
499 fh2Rhovspthardest1 = new TH2F("fh2Rhovspthardest1","Background vs pTjet 1",100,0.,200.,50,0.,5.);
500 fh2Rhovspthardest2 = new TH2F("fh2Rhovspthardest2","Background vs pTjet 2",100,0.,200.,50,0.,5.);
501 fh2Rhovspthardest3 = new TH2F("fh2Rhovspthardest3","Background vs pTjet 3",100,0.,200.,50,0.,5.);
502 fh2Errorvspthardest1 = new TH2F("fh2Errorvspthardest1","Relative error 1",100,0.,200.,50,0.,5.);
503 fh2Errorvspthardest2 = new TH2F("fh2Errorvspthardest2","Relative error 2",100,0.,200.,50,0.,5.);
504 fh2Errorvspthardest3 = new TH2F("fh2Errorvspthardest3","Relative error 3",100,0.,200.,50,0.,5.);
505 }
506
507
508
509
510 const Int_t saveLevel = 3; // large save level more histos
511 if(saveLevel>0){
512 fHistList->Add(fh1Xsec);
513 fHistList->Add(fh1Trials);
514 fHistList->Add(fh1PtHard);
515 fHistList->Add(fh1PtHardNoW);
516 fHistList->Add(fh1PtHardTrials);
517 fHistList->Add(fh1ZVtx);
518 if(fBranchGen.Length()>0||fBkgSubtraction){
519 fHistList->Add(fh1NGenJets);
520 fHistList->Add(fh1PtTracksGenIn);
521 }
522 fHistList->Add(fh1PtJetsRecIn);
523 fHistList->Add(fh1PtJetsLeadingRecIn);
524 fHistList->Add(fh1PtTracksRecIn);
525 fHistList->Add(fh1PtTracksLeadingRecIn);
526 fHistList->Add(fh1NRecJets);
527 fHistList->Add(fh1PtTrackRec);
528 fHistList->Add(fh1SumPtTrackRec);
529 fHistList->Add(fh1SumPtTrackAreaRec);
530 fHistList->Add(fh2NRecJetsPt);
531 fHistList->Add(fh2NRecTracksPt);
532 fHistList->Add(fh2JetsLeadingPhiEta );
533 fHistList->Add(fh2JetsLeadingPhiPt );
534 fHistList->Add(fh2TracksLeadingPhiEta);
535 fHistList->Add(fh2TracksLeadingPhiPt);
536 for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
537 for(int ij = 0;ij<kMaxJets;++ij){
538 fHistList->Add( fh1PtRecIn[ij]);
539
540 if(fBranchGen.Length()>0||fBkgSubtraction){
541 fHistList->Add(fh1PtGenIn[ij]);
542 fHistList->Add(fh2FragGen[ij]);
543 fHistList->Add(fh2FragLnGen[ij]);
544 fHistList->Add(fh2RhoPtGen[ij]);
545 fHistList->Add(fh2PsiPtGen[ij]);
546 }
547 fHistList->Add( fh2PhiPt[ij]);
548 fHistList->Add( fh2PhiEta[ij]);
549 fHistList->Add( fh2RhoPtRec[ij]);
550 fHistList->Add( fh2PsiPtRec[ij]);
551 fHistList->Add( fh2FragRec[ij]);
552 fHistList->Add( fh2FragLnRec[ij]);
553 }
554 fHistList->Add(fhnCorrelation);
555 fHistList->Add(fhnCorrelationPhiZRec);
556 fHistList->Add(fh2JetPtJetPhi);
557 fHistList->Add(fh2TrackPtTrackPhi);
558 fHistList->Add(fh2RelPtFGen);
559
560 fHistList->Add(fh2DijetDeltaPhiPt);
561 fHistList->Add(fh2DijetAsymPt);
562 fHistList->Add(fh2DijetAsymPtCut);
563 fHistList->Add(fh2DijetDeltaPhiDeltaEta);
564 fHistList->Add(fh2DijetPt2vsPt1);
565 fHistList->Add(fh2DijetDifvsSum);
566 fHistList->Add(fh1DijetMinv);
567 fHistList->Add(fh1DijetMinvCut);
568 if(fBkgSubtraction){
569 fHistList->Add(fh1Bkg1);
570 fHistList->Add(fh1Bkg2);
571 fHistList->Add(fh1Bkg3);
572 fHistList->Add(fh1Sigma1);
573 fHistList->Add(fh1Sigma2);
574 fHistList->Add(fh1Sigma3);
575 fHistList->Add(fh1Area1);
576 fHistList->Add(fh1Area2);
577 fHistList->Add(fh1Area3);
578 fHistList->Add(fh1Ptjet);
579 fHistList->Add(fh1Ptjethardest);
580 fHistList->Add(fh1Ptjetsub1);
581 fHistList->Add(fh1Ptjetsub2);
582 fHistList->Add(fh1Ptjetsub3);
583 fHistList->Add(fh1Ptjetsubhardest1);
584 fHistList->Add(fh1Ptjetsubhardest2);
585 fHistList->Add(fh1Ptjetsubhardest3);
586 fHistList->Add(fh2Rhovspthardest1);
587 fHistList->Add(fh2Rhovspthardest2);
588 fHistList->Add(fh2Rhovspthardest3);
589 fHistList->Add(fh2Errorvspthardest1);
590 fHistList->Add(fh2Errorvspthardest2);
591 fHistList->Add(fh2Errorvspthardest3);
592 }
593 }
594
595 // =========== Switch on Sumw2 for all histos ===========
596 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
597 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
598 if (h1){
599 h1->Sumw2();
600 continue;
601 }
602 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
603 if(hn)hn->Sumw2();
604 }
605 TH1::AddDirectory(oldStatus);
606}
607
608void AliAnalysisTaskJetSpectrum2::Init()
609{
610 //
611 // Initialization
612 //
613
614 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
615
616}
617
618void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
619{
620
621 Bool_t selected = kTRUE;
622
623 if(fUseGlobalSelection&&fEventSelectionMask==0){
624 selected = AliAnalysisHelperJetTasks::Selected();
625 }
626 else if(fUseGlobalSelection&&fEventSelectionMask>0){
627 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
628 }
629
630 if(fEventClass>0){
631 selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
632 }
633
634 if(!selected){
635 // no selection by the service task, we continue
636 if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
637 PostData(1, fHistList);
638 return;
639 }
640
641
642 //
643 // Execute analysis for current event
644 //
645 AliESDEvent *fESD = 0;
646 if(fUseAODJetInput){
647 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
648 if(!fAOD){
649 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
650 return;
651 }
652 // fethc the header
653 }
654 else{
655 // assume that the AOD is in the general output...
656 fAOD = AODEvent();
657 if(!fAOD){
658 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
659 return;
660 }
661 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
662 }
663
664
665
666
667 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
668
669
670 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
671
672 if(!aodH){
673 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
674 return;
675 }
676
677 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
678 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
679 if(!aodRecJets){
680 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
681 return;
682 }
683
684 Int_t nJets = aodRecJets->GetEntriesFast();
685
686 // ==== General variables needed
687 // We use statice array, not to fragment the memory
688 AliAODJet genJets[kMaxJets];
689 Int_t nGenJets = 0;
690 AliAODJet recJets[kMaxJets];
691 Int_t nRecJets = 0;
692 ///////////////////////////
693
694
695
696
697 if(fBkgSubtraction){
698 AliAODJetEventBackground* evBkg=(AliAODJetEventBackground*)fAOD->FindListObject(fBranchBkg.Data());
699
700
701 if(!evBkg){
702 Printf("%s:%d no reconstructed background array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkg.Data());
703 return;
704 }
705
706
707 ///just to start: some very simple plots containing rho, sigma and area of the background.
708
709
710 Float_t pthardest=0.;
711 if(nJets!=0){
712
713 Float_t bkg1=evBkg->GetBackground(0);
714 Float_t bkg2=evBkg->GetBackground(1);
715 Float_t bkg3=evBkg->GetBackground(2);
716 Float_t sigma1=evBkg->GetSigma(0);
717 Float_t sigma2=evBkg->GetSigma(1);
718 Float_t sigma3=evBkg->GetSigma(2);
719 Float_t area1=evBkg->GetMeanarea(0);
720 Float_t area2=evBkg->GetMeanarea(1);
721 Float_t area3=evBkg->GetMeanarea(2);
722 fh1Bkg1->Fill(bkg1); //rho computed with all background jets.
723 fh1Bkg2->Fill(bkg2); //rho computed with all background jets but the 2 hardest.
724 fh1Bkg3->Fill(bkg3); //rho computed with randomized jets
725 fh1Sigma1->Fill(sigma1);
726 fh1Sigma2->Fill(sigma2);
727 fh1Sigma3->Fill(sigma3);
728 fh1Area1->Fill(area1);
729 fh1Area2->Fill(area2);
730 fh1Area3->Fill(area3);
731
732
733 for(Int_t k=0;k<nJets;k++){
734 AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets->At(k));
735 fh1Ptjet->Fill(jet->Pt());
736 Float_t ptsub1=jet->Pt()-bkg1*jet->EffectiveAreaCharged();
737 Float_t ptsub2=jet->Pt()-bkg2*jet->EffectiveAreaCharged();
738 Float_t ptsub3=jet->Pt()-bkg3*jet->EffectiveAreaCharged();
739 if(ptsub2<0.) ptsub2=0.;
740 if(ptsub1<0.) ptsub1=0.;
741 if(ptsub3<0.) ptsub3=0.;
742 Float_t err1=sigma1*sqrt(area1);
743 Float_t err2=sigma2*sqrt(area2);
744 Float_t err3=sigma3*sqrt(area3);
745 fh1Ptjetsub1->Fill(ptsub1);
746 fh1Ptjetsub2->Fill(ptsub2);
747 fh1Ptjetsub3->Fill(ptsub3);
748 if(k==0) {
749 pthardest=jet->Pt();
750 fh1Ptjethardest->Fill(pthardest);
751 fh1Ptjetsubhardest1->Fill(ptsub1);
752 fh1Ptjetsubhardest2->Fill(ptsub2);
753 fh1Ptjetsubhardest3->Fill(ptsub3);
754 fh2Errorvspthardest1->Fill(ptsub1,err1/ptsub1);
755 fh2Errorvspthardest2->Fill(ptsub2,err2/ptsub2);
756 fh2Errorvspthardest3->Fill(ptsub3,err3/ptsub3);
757 }
758
759 Float_t ptsub=0.;
760 if(fFillCorrBkg==1) ptsub=ptsub1;
761 if(fFillCorrBkg==2) ptsub=ptsub2;
762 if(fFillCorrBkg==3) ptsub=ptsub3;
763 Float_t subphi=jet->Phi();
764 Float_t subtheta=jet->Theta();
765 Float_t subpz = ptsub/TMath::Tan(subtheta);
766 Float_t subpx=ptsub*TMath::Cos(subphi);
767 Float_t subpy=ptsub * TMath::Sin(subphi);
768 Float_t subp = TMath::Sqrt(ptsub*ptsub+subpz*subpz);
769 if(k<kMaxJets){
770 genJets[k].SetPxPyPzE(subpx,subpy,subpz,subp);
771 nGenJets = k+1;
772 }
773 }
774 fh2Rhovspthardest1->Fill(pthardest,bkg1);
775 fh2Rhovspthardest2->Fill(pthardest,bkg2);
776 fh2Rhovspthardest3->Fill(pthardest,bkg3);
777 }
778 }// background subtraction
779
780
781 Double_t eventW = 1;
782 Double_t ptHard = 0;
783 Double_t nTrials = 1; // Trials for MC trigger
784
785 if(fUseExternalWeightOnly){
786 eventW = fExternalWeight;
787 }
788
789 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
790 // if(fDebug>0)aodH->SetFillAOD(kFALSE);
791 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
792 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
793 // this is the part we only use when we have MC information
794 AliMCEvent* mcEvent = MCEvent();
795 // AliStack *pStack = 0;
796 if(!mcEvent){
797 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
798 return;
799 }
800 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
801 Int_t iCount = 0;
802 if(pythiaGenHeader){
803 nTrials = pythiaGenHeader->Trials();
804 ptHard = pythiaGenHeader->GetPtHard();
805 int iProcessType = pythiaGenHeader->ProcessType();
806 // 11 f+f -> f+f
807 // 12 f+barf -> f+barf
808 // 13 f+barf -> g+g
809 // 28 f+g -> f+g
810 // 53 g+g -> f+barf
811 // 68 g+g -> g+g
812 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
813 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
814
815 // fetch the pythia generated jets only to be used here
816 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
817 AliAODJet pythiaGenJets[kMaxJets];
818 for(int ip = 0;ip < nPythiaGenJets;++ip){
819 if(iCount>=kMaxJets)continue;
820 Float_t p[4];
821 pythiaGenHeader->TriggerJet(ip,p);
822 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
823
824 if(fBranchGen.Length()==0){
825 /*
826 if(fLimitGenJetEta){
827 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
828 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
829 }
830 */
831 if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
832 // if we have MC particles and we do not read from the aod branch
833 // use the pythia jets
834 if(!fBkgSubtraction){
835 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
836 iCount++;
837 }
838 }
839 }
840 }
841 if(fBranchGen.Length()==0&&!fBkgSubtraction)nGenJets = iCount;
842 }// (fAnalysisType&kMCESD)==kMCESD)
843
844
845 // we simply fetch the tracks/mc particles as a list of AliVParticles
846
847 TList recParticles;
848 TList genParticles;
849
850
851
852
853 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
854 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
855 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
856 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
857
858
859 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
860 fh1PtHard->Fill(ptHard,eventW);
861 fh1PtHardNoW->Fill(ptHard,1);
862 fh1PtHardTrials->Fill(ptHard,nTrials);
863 fh1ZVtx->Fill(fAOD->GetPrimaryVertex()->GetZ());
864
865
866 // If we set a second branch for the input jets fetch this
867 if(fBranchGen.Length()>0&&!fBkgSubtraction){
868 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
869 if(aodGenJets){
870 Int_t iCount = 0;
871 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
872 if(iCount>=kMaxJets)continue;
873 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
874 if(!tmp)continue;
875 /*
876 if(fLimitGenJetEta){
877 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
878 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
879 }
880 */
881 if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
882 genJets[iCount] = *tmp;
883 iCount++;
884 }
885 nGenJets = iCount;
886 }
887 else{
888 if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
889 if(fDebug>2)fAOD->Print();
890 }
891 }
892
893 fh1NGenJets->Fill(nGenJets);
894 // We do not want to exceed the maximum jet number
895 nGenJets = TMath::Min(nGenJets,kMaxJets);
896
897 // Fetch the reconstructed jets...
898
899 nRecJets = aodRecJets->GetEntries();
900
901 nRecJets = aodRecJets->GetEntries();
902 fh1NRecJets->Fill(nRecJets);
903
904 // Do something with the all rec jets
905 Int_t nRecOver = nRecJets;
906
907 // check that the jets are sorted
908 Float_t ptOld = 999999;
909 for(int ir = 0;ir < nRecJets;ir++){
910 AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
911 Float_t tmpPt = tmp->Pt();
912 if(tmpPt>ptOld){
913 Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
914 }
915 ptOld = tmpPt;
916 }
917
918
919 if(nRecOver>0){
920 TIterator *recIter = aodRecJets->MakeIterator();
921 AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
922 Float_t pt = tmpRec->Pt();
923 if(tmpRec){
924 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
925 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
926 while(pt<ptCut&&tmpRec){
927 nRecOver--;
928 tmpRec = (AliAODJet*)(recIter->Next());
929 if(tmpRec){
930 pt = tmpRec->Pt();
931 }
932 }
933 if(nRecOver<=0)break;
934 fh2NRecJetsPt->Fill(ptCut,nRecOver);
935 }
936 }
937 recIter->Reset();
938 AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
939 Float_t phi = leading->Phi();
940 if(phi<0)phi+=TMath::Pi()*2.;
941 Float_t eta = leading->Eta();
942 pt = leading->Pt();
943 while((tmpRec = (AliAODJet*)(recIter->Next()))){
944 Float_t tmpPt = tmpRec->Pt();
945 fh1PtJetsRecIn->Fill(tmpPt);
946 if(tmpRec==leading){
947 fh1PtJetsLeadingRecIn->Fill(tmpPt);
948 continue;
949 }
950 // correlation
951 Float_t tmpPhi = tmpRec->Phi();
952
953 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
954 Float_t dPhi = phi - tmpRec->Phi();
955 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
956 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
957 Float_t dEta = eta - tmpRec->Eta();
958 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
959 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
960 }
961 delete recIter;
962 }
963
964 Int_t nTrackOver = recParticles.GetSize();
965 // do the same for tracks and jets
966 if(nTrackOver>0){
967 TIterator *recIter = recParticles.MakeIterator();
968 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
969 Float_t pt = tmpRec->Pt();
970 // Printf("Leading track p_t %3.3E",pt);
971 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
972 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
973 while(pt<ptCut&&tmpRec){
974 nTrackOver--;
975 tmpRec = (AliVParticle*)(recIter->Next());
976 if(tmpRec){
977 pt = tmpRec->Pt();
978 }
979 }
980 if(nTrackOver<=0)break;
981 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
982 }
983
984 recIter->Reset();
985 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
986 Float_t phi = leading->Phi();
987 if(phi<0)phi+=TMath::Pi()*2.;
988 Float_t eta = leading->Eta();
989 pt = leading->Pt();
990 while((tmpRec = (AliVParticle*)(recIter->Next()))){
991 Float_t tmpPt = tmpRec->Pt();
992 fh1PtTracksRecIn->Fill(tmpPt);
993 if(tmpRec==leading){
994 fh1PtTracksLeadingRecIn->Fill(tmpPt);
995 continue;
996 }
997 // correlation
998 Float_t tmpPhi = tmpRec->Phi();
999
1000 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1001 Float_t dPhi = phi - tmpRec->Phi();
1002 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1003 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1004 Float_t dEta = eta - tmpRec->Eta();
1005 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1006 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
1007 }
1008 delete recIter;
1009 }
1010
1011 if(genParticles.GetSize()){
1012 TIterator *genIter = genParticles.MakeIterator();
1013 AliVParticle *tmpGen = 0;
1014 while((tmpGen = (AliVParticle*)(genIter->Next()))){
1015 if(TMath::Abs(tmpGen->Eta())<0.9){
1016 Float_t tmpPt = tmpGen->Pt();
1017 fh1PtTracksGenIn->Fill(tmpPt);
1018 }
1019 }
1020 delete genIter;
1021 }
1022
1023 nRecJets = TMath::Min(nRecJets,kMaxJets);
1024 nJets = TMath::Min(nJets,kMaxJets);
1025 Int_t iCountRec = 0;
1026 for(int ir = 0;ir < nRecJets;++ir){
1027 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
1028 if(!tmp)continue;
1029 if(tmp->Pt()<fMinJetPt)continue;
1030 recJets[ir] = *tmp;
1031 iCountRec++;
1032 }
1033 nRecJets = iCountRec;
1034
1035
1036 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1037 // Relate the jets
1038 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
1039 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
1040
1041
1042 for(int i = 0;i<kMaxJets;++i){
1043 iGenIndex[i] = iRecIndex[i] = -1;
1044 }
1045
1046
1047 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
1048 iGenIndex,iRecIndex,fDebug);
1049
1050
1051 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1052
1053 if(fDebug){
1054 for(int i = 0;i<kMaxJets;++i){
1055 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
1056 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
1057 }
1058 }
1059
1060
1061
1062
1063 Double_t container[6];
1064 Double_t containerPhiZ[6];
1065
1066 // loop over generated jets
1067
1068 // radius; tmp, get from the jet header itself
1069 // at some point, how todeal woht FastJet on MC?
1070 Float_t radiusGen =0.4;
1071 Float_t radiusRec =0.4;
1072
1073 for(int ig = 0;ig < nGenJets;++ig){
1074 Double_t ptGen = genJets[ig].Pt();
1075 Double_t phiGen = genJets[ig].Phi();
1076 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1077 Double_t etaGen = genJets[ig].Eta();
1078
1079 container[3] = ptGen;
1080 container[4] = etaGen;
1081 container[5] = phiGen;
1082 fhnJetContainer[kStep0]->Fill(&container[3],eventW);
1083 Int_t ir = iRecIndex[ig];
1084
1085 if(TMath::Abs(etaGen)<fRecEtaWindow){
1086 fh1TmpRho->Reset();
1087
1088 fhnJetContainer[kStep1]->Fill(&container[3],eventW);
1089 fh1PtGenIn[ig]->Fill(ptGen,eventW);
1090 // fill the fragmentation function
1091 for(int it = 0;it<genParticles.GetEntries();++it){
1092 AliVParticle *part = (AliVParticle*)genParticles.At(it);
1093 Float_t deltaR = genJets[ig].DeltaR(part);
1094 fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
1095 if(deltaR<radiusGen){
1096 Float_t z = part->Pt()/ptGen;
1097 Float_t lnz = -1.*TMath::Log(z);
1098 fh2FragGen[ig]->Fill(z,ptGen,eventW);
1099 fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
1100 }
1101
1102 }
1103 Float_t rhoSum = 0;
1104 for(int ibx = 1;ibx<=fh2RhoPtGen[ir]->GetNbinsX();ibx++){
1105 Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
1106 Float_t rho = fh1TmpRho->GetBinContent(ibx);
1107 rhoSum += rho;
1108 fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
1109 fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
1110 }
1111 }
1112 if(ir>=0&&ir<nRecJets){
1113 fhnJetContainer[kStep2]->Fill(&container[3],eventW);
1114 Double_t etaRec = recJets[ir].Eta();
1115 if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
1116 if(TMath::Abs(etaRec)<fRecEtaWindow&&TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep3]->Fill(&container[3],eventW);
1117 }
1118 }// loop over generated jets
1119
1120 Float_t sumPt = 0;
1121 for(int it = 0;it<recParticles.GetEntries();++it){
1122 AliVParticle *part = (AliVParticle*)recParticles.At(it);
1123 // fill sum pt and P_t of all paritles
1124 if(TMath::Abs(part->Eta())<0.9){
1125 Float_t pt = part->Pt();
1126 fh1PtTrackRec->Fill(pt,eventW);
1127 fh2TrackPtTrackPhi->Fill(pt,part->Phi());
1128 sumPt += pt;
1129 }
1130 }
1131 fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
1132 fh1SumPtTrackRec->Fill(sumPt,eventW);
1133
1134
1135 // loop over reconstructed jets
1136 for(int ir = 0;ir < nRecJets;++ir){
1137 Double_t etaRec = recJets[ir].Eta();
1138 Double_t ptRec = recJets[ir].Pt();
1139 Double_t phiRec = recJets[ir].Phi();
1140 if(phiRec<0)phiRec+=TMath::Pi()*2.;
1141
1142
1143 // do something with dijets...
1144 if(ir==1){
1145 Double_t etaRec1 = recJets[0].Eta();
1146 Double_t ptRec1 = recJets[0].Pt();
1147 Double_t phiRec1 = recJets[0].Phi();
1148 if(phiRec1<0)phiRec1+=TMath::Pi()*2.;
1149
1150
1151 if(TMath::Abs(etaRec1)<fRecEtaWindow
1152 &&TMath::Abs(etaRec)<fRecEtaWindow){
1153
1154 Float_t deltaPhi = phiRec1 - phiRec;
1155
1156 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1157 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1158 deltaPhi = TMath::Abs(deltaPhi);
1159 fh2DijetDeltaPhiPt->Fill(deltaPhi,ptRec1);
1160 Float_t asym = (ptRec1-ptRec)/(ptRec1+ptRec);
1161 fh2DijetAsymPt->Fill(asym,ptRec1);
1162 fh2DijetDeltaPhiDeltaEta->Fill(deltaPhi,etaRec1-etaRec);
1163 fh2DijetPt2vsPt1->Fill(ptRec1,ptRec);
1164 fh2DijetDifvsSum->Fill(ptRec1+ptRec,ptRec1-ptRec);
1165 Float_t minv = 2.*(recJets[0].P()*recJets[1].P()-
1166 recJets[0].Px()*recJets[1].Px()-
1167 recJets[0].Py()*recJets[1].Py()-
1168 recJets[0].Pz()*recJets[1].Py());
1169 minv = TMath::Sqrt(minv);
1170 // with mass == 0;
1171
1172 fh1DijetMinv->Fill(minv);
1173 if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
1174 fh1DijetMinvCut->Fill(minv);
1175 fh2DijetAsymPtCut->Fill(asym,ptRec1);
1176 }
1177 }
1178 }
1179
1180
1181 container[0] = ptRec;
1182 container[1] = etaRec;
1183 container[2] = phiRec;
1184 containerPhiZ[0] = ptRec;
1185 containerPhiZ[1] = phiRec;
1186 if(ptRec>30.&&fDebug>0){
1187 // need to cast to int, otherwise the printf overwrites
1188 Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
1189 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
1190 if(fESD)Printf("ESDEvent GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
1191 // aodH->SetFillAOD(kTRUE);
1192 fAOD->GetHeader()->Print();
1193 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
1194 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1195 AliAODTrack *tr = fAOD->GetTrack(it);
1196 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1197 tr->Print();
1198 // tr->Dump();
1199 if(fESD){
1200 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1201 esdTr->Print("");
1202 // esdTr->Dump();
1203 }
1204 }
1205 }
1206
1207
1208 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
1209 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1210
1211 Float_t zLeading = -1;
1212 if(TMath::Abs(etaRec)<fRecEtaWindow){
1213 fh2JetPtJetPhi->Fill(ptRec,phiRec);
1214 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
1215 fh1PtRecIn[ir]->Fill(ptRec,eventW);
1216 // fill the fragmentation function
1217
1218 fh1TmpRho->Reset();
1219
1220 for(int it = 0;it<recParticles.GetEntries();++it){
1221 AliVParticle *part = (AliVParticle*)recParticles.At(it);
1222 Float_t eta = part->Eta();
1223 if(TMath::Abs(eta)<0.9){
1224 Float_t phi = part->Phi();
1225 if(phi<0)phi+=TMath::Pi()*2.;
1226 Float_t dPhi = phi - phiRec;
1227 Float_t dEta = eta - etaRec;
1228 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1229 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1230 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1231 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1232 }
1233
1234 Float_t deltaR = recJets[ir].DeltaR(part);
1235 fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
1236
1237
1238 if(deltaR<radiusRec){
1239 Float_t z = part->Pt()/ptRec;
1240 if(z>zLeading)zLeading=z;
1241 Float_t lnz = -1.*TMath::Log(z);
1242 fh2FragRec[ir]->Fill(z,ptRec,eventW);
1243 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1244 }
1245 }
1246 // fill the jet shapes
1247 Float_t rhoSum = 0;
1248 for(int ibx = 1;ibx<=fh2RhoPtRec[ir]->GetNbinsX();ibx++){
1249 Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1250 Float_t rho = fh1TmpRho->GetBinContent(ibx);
1251 rhoSum += rho;
1252 fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1253 fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1254 }
1255 }
1256
1257 // Fill Correlation
1258 Int_t ig = iGenIndex[ir];
1259 if(ig>=0 && ig<nGenJets){
1260 fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1261 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1262 Double_t ptGen = genJets[ig].Pt();
1263 Double_t phiGen = genJets[ig].Phi();
1264 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1265 Double_t etaGen = genJets[ig].Eta();
1266
1267 container[3] = ptGen;
1268 container[4] = etaGen;
1269 container[5] = phiGen;
1270 containerPhiZ[3] = ptGen;
1271 //
1272 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1273 //
1274
1275 if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1276 if(TMath::Abs(etaRec)<fRecEtaWindow){
1277 fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1278 fhnCorrelation->Fill(container);
1279 if(ptGen>0){
1280 Float_t delta = (ptRec-ptGen)/ptGen;
1281 fh2RelPtFGen->Fill(ptGen,delta,eventW);
1282 }
1283 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1284 }// if etarec in window
1285 }
1286 else{
1287 containerPhiZ[3] = 0;
1288 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1289 }
1290 }// loop over reconstructed jets
1291
1292
1293 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1294 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1295 PostData(1, fHistList);
1296}
1297
1298
1299void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1300 //
1301 // Create the particle container for the correction framework manager and
1302 // link it
1303 //
1304 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
1305 const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1306 const Double_t kEtamin = -3.0, kEtamax = 3.0;
1307 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1308 const Double_t kZmin = 0., kZmax = 1;
1309
1310 // can we neglect migration in eta and phi?
1311 // phi should be no problem since we cover full phi and are phi symmetric
1312 // eta migration is more difficult due to needed acceptance correction
1313 // in limited eta range
1314
1315 //arrays for the number of bins in each dimension
1316 Int_t iBin[kNvar];
1317 iBin[0] = 320; //bins in pt
1318 iBin[1] = 1; //bins in eta
1319 iBin[2] = 1; // bins in phi
1320
1321 //arrays for lower bounds :
1322 Double_t* binEdges[kNvar];
1323 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1324 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1325
1326 //values for bin lower bounds
1327 // 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);
1328 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1329 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1330 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1331
1332
1333 for(int i = 0;i < kMaxStep*2;++i){
1334 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1335 for (int k=0; k<kNvar; k++) {
1336 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1337 }
1338 }
1339 //create correlation matrix for unfolding
1340 Int_t thnDim[2*kNvar];
1341 for (int k=0; k<kNvar; k++) {
1342 //first half : reconstructed
1343 //second half : MC
1344 thnDim[k] = iBin[k];
1345 thnDim[k+kNvar] = iBin[k];
1346 }
1347
1348 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1349 for (int k=0; k<kNvar; k++) {
1350 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1351 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1352 }
1353 fhnCorrelation->Sumw2();
1354
1355 // for second correlation histogram
1356
1357
1358 const Int_t kNvarPhiZ = 4;
1359 //arrays for the number of bins in each dimension
1360 Int_t iBinPhiZ[kNvarPhiZ];
1361 iBinPhiZ[0] = 80; //bins in pt
1362 iBinPhiZ[1] = 72; //bins in phi
1363 iBinPhiZ[2] = 20; // bins in Z
1364 iBinPhiZ[3] = 80; //bins in ptgen
1365
1366 //arrays for lower bounds :
1367 Double_t* binEdgesPhiZ[kNvarPhiZ];
1368 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1369 binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1370
1371 for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1372 for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1373 for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1374 for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1375
1376 fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1377 for (int k=0; k<kNvarPhiZ; k++) {
1378 fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1379 }
1380 fhnCorrelationPhiZRec->Sumw2();
1381
1382
1383 // Add a histogram for Fake jets
1384
1385 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1386 delete [] binEdges[ivar];
1387
1388 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1389 delete [] binEdgesPhiZ[ivar];
1390
1391}
1392
1393void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1394{
1395// Terminate analysis
1396//
1397 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1398}
1399
1400
1401Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1402
1403 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1404
1405 Int_t iCount = 0;
1406 if(type==kTrackAOD){
1407 AliAODEvent *aod = 0;
1408 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1409 else aod = AODEvent();
1410 if(!aod){
1411 return iCount;
1412 }
1413 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1414 AliAODTrack *tr = aod->GetTrack(it);
1415 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1416 if(TMath::Abs(tr->Eta())>0.9)continue;
1417 if(fDebug>0){
1418 if(tr->Pt()>20){
1419 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1420 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1421 tr->Print();
1422 // tr->Dump();
1423 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1424 if(fESD){
1425 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1426 esdTr->Print("");
1427 // esdTr->Dump();
1428 }
1429 }
1430 }
1431 list->Add(tr);
1432 iCount++;
1433 }
1434 }
1435 else if (type == kTrackKineAll||type == kTrackKineCharged){
1436 AliMCEvent* mcEvent = MCEvent();
1437 if(!mcEvent)return iCount;
1438 // we want to have alivpartilces so use get track
1439 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1440 if(!mcEvent->IsPhysicalPrimary(it))continue;
1441 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1442 if(type == kTrackKineAll){
1443 list->Add(part);
1444 iCount++;
1445 }
1446 else if(type == kTrackKineCharged){
1447 if(part->Particle()->GetPDG()->Charge()==0)continue;
1448 list->Add(part);
1449 iCount++;
1450 }
1451 }
1452 }
1453 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1454 AliAODEvent *aod = 0;
1455 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1456 else aod = AODEvent();
1457 if(!aod)return iCount;
1458 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1459 if(!tca)return iCount;
1460 for(int it = 0;it < tca->GetEntriesFast();++it){
1461 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1462 if(!part->IsPhysicalPrimary())continue;
1463 if(type == kTrackAODMCAll){
1464 list->Add(part);
1465 iCount++;
1466 }
1467 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1468 if(part->Charge()==0)continue;
1469 if(kTrackAODMCCharged){
1470 list->Add(part);
1471 }
1472 else {
1473 if(TMath::Abs(part->Eta())>0.9)continue;
1474 list->Add(part);
1475 }
1476 iCount++;
1477 }
1478 }
1479 }// AODMCparticle
1480 list->Sort();
1481 return iCount;
1482
1483}