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