]> 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
limit printouts
[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
3b7ffecf 390 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
391
392}
393
394void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
395{
b5cc0c6d 396
a221560c 397 if(!AliAnalysisHelperJetTasks::Selected()&&fUseGlobalSelection){
b5cc0c6d 398 // no selection by the service task, we continue
a221560c 399 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
b5cc0c6d 400 PostData(1, fHistList);
401 return;
402 }
403
3b7ffecf 404 //
405 // Execute analysis for current event
406 //
7fa8b2da 407 AliESDEvent *fESD = 0;
565584e8 408 if(fUseAODJetInput){
3b7ffecf 409 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
410 if(!fAOD){
565584e8 411 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
3b7ffecf 412 return;
413 }
414 // fethc the header
415 }
416 else{
417 // assume that the AOD is in the general output...
418 fAOD = AODEvent();
419 if(!fAOD){
420 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
421 return;
422 }
7fa8b2da 423 if(fDebug>0){
424 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
425 }
3b7ffecf 426 }
427
428
429
430
431 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
432
433
434 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
435
436 if(!aodH){
437 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
438 return;
439 }
440
441 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
442 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
443 if(!aodRecJets){
444 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
445 return;
446 }
447
448 // ==== General variables needed
449
450
451 // We use statice array, not to fragment the memory
452 AliAODJet genJets[kMaxJets];
453 Int_t nGenJets = 0;
454 AliAODJet recJets[kMaxJets];
455 Int_t nRecJets = 0;
456 ///////////////////////////
599396fa 457
3b7ffecf 458
459 Double_t eventW = 1;
460 Double_t ptHard = 0;
461 Double_t nTrials = 1; // Trials for MC trigger
462
463 if(fUseExternalWeightOnly){
464 eventW = fExternalWeight;
465 }
466
467 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
cf049e94 468 // if(fDebug>0)aodH->SetFillAOD(kFALSE);
3b7ffecf 469 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
470 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
471 // this is the part we only use when we have MC information
472 AliMCEvent* mcEvent = MCEvent();
473 // AliStack *pStack = 0;
474 if(!mcEvent){
475 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
476 return;
477 }
478 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
479 Int_t iCount = 0;
480 if(pythiaGenHeader){
481 nTrials = pythiaGenHeader->Trials();
482 ptHard = pythiaGenHeader->GetPtHard();
483 int iProcessType = pythiaGenHeader->ProcessType();
484 // 11 f+f -> f+f
485 // 12 f+barf -> f+barf
486 // 13 f+barf -> g+g
487 // 28 f+g -> f+g
488 // 53 g+g -> f+barf
489 // 68 g+g -> g+g
490 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
491 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
492
493 // fetch the pythia generated jets only to be used here
494 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
495 AliAODJet pythiaGenJets[kMaxJets];
496 for(int ip = 0;ip < nPythiaGenJets;++ip){
497 if(iCount>=kMaxJets)continue;
498 Float_t p[4];
499 pythiaGenHeader->TriggerJet(ip,p);
500 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
501
5301f754 502 /*
3b7ffecf 503 if(fLimitGenJetEta){
504 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
505 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
506 }
5301f754 507 */
3b7ffecf 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;
5301f754 549 /*
3b7ffecf 550 if(fLimitGenJetEta){
551 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
552 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
553 }
5301f754 554 */
3b7ffecf 555 genJets[iCount] = *tmp;
556 iCount++;
557 }
558 nGenJets = iCount;
559 }
560 else{
5cb0f438 561 if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
562 if(fDebug>2)fAOD->Print();
3b7ffecf 563 }
564 }
565
566 fh1NGenJets->Fill(nGenJets);
567 // We do not want to exceed the maximum jet number
568 nGenJets = TMath::Min(nGenJets,kMaxJets);
569
570 // Fetch the reconstructed jets...
571
572 nRecJets = aodRecJets->GetEntries();
573
cc0649e4 574 nRecJets = aodRecJets->GetEntries();
575 fh1NRecJets->Fill(nRecJets);
576
577 // Do something with the all rec jets
578 Int_t nRecOver = nRecJets;
3b7ffecf 579
cc0649e4 580 // check that the jets are sorted
581 Float_t ptOld = 999999;
582 for(int ir = 0;ir < nRecJets;ir++){
583 AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
584 Float_t tmpPt = tmp->Pt();
585 if(tmpPt>ptOld){
586 Printf("%s:%d Jets Not Sorted!! %d:%.3E %d%.3E",(char*)__FILE__,__LINE__,ir,tmpPt,ir-1,ptOld);
587 }
588 ptOld = tmpPt;
589 }
3b7ffecf 590
cc0649e4 591
592 if(nRecOver>0){
593 TIterator *recIter = aodRecJets->MakeIterator();
594 AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
595 Float_t pt = tmpRec->Pt();
596 if(tmpRec){
597 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
598 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
599 while(pt<ptCut&&tmpRec){
600 nRecOver--;
601 tmpRec = (AliAODJet*)(recIter->Next());
602 if(tmpRec){
603 pt = tmpRec->Pt();
604 }
605 }
606 if(nRecOver<=0)break;
607 fh2NRecJetsPt->Fill(ptCut,nRecOver);
608 }
609 }
610 recIter->Reset();
611 AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
612 Float_t phi = leading->Phi();
613 if(phi<0)phi+=TMath::Pi()*2.;
614 Float_t eta = leading->Eta();
cf049e94 615 pt = leading->Pt();
565584e8 616 while((tmpRec = (AliAODJet*)(recIter->Next()))){
cc0649e4 617 Float_t tmpPt = tmpRec->Pt();
618 fh1PtJetsRecIn->Fill(tmpPt);
565584e8 619 if(tmpRec==leading){
620 fh1PtJetsLeadingRecIn->Fill(tmpPt);
621 continue;
622 }
cc0649e4 623 // correlation
624 Float_t tmpPhi = tmpRec->Phi();
625
626 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
627 Float_t dPhi = phi - tmpRec->Phi();
628 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
629 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
630 Float_t dEta = eta - tmpRec->Eta();
631 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
565584e8 632 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
cc0649e4 633 }
634 delete recIter;
635 }
636
637 Int_t nTrackOver = recParticles.GetSize();
638 // do the same for tracks and jets
639 if(nTrackOver>0){
640 TIterator *recIter = recParticles.MakeIterator();
641 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
642 Float_t pt = tmpRec->Pt();
643 // Printf("Leading track p_t %3.3E",pt);
644 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
645 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
646 while(pt<ptCut&&tmpRec){
647 nTrackOver--;
648 tmpRec = (AliVParticle*)(recIter->Next());
649 if(tmpRec){
650 pt = tmpRec->Pt();
651 }
652 }
653 if(nTrackOver<=0)break;
654 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
655 }
656
657 recIter->Reset();
658 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
659 Float_t phi = leading->Phi();
660 if(phi<0)phi+=TMath::Pi()*2.;
661 Float_t eta = leading->Eta();
cf049e94 662 pt = leading->Pt();
cc0649e4 663 while((tmpRec = (AliVParticle*)(recIter->Next()))){
664 Float_t tmpPt = tmpRec->Pt();
665 fh1PtTracksRecIn->Fill(tmpPt);
565584e8 666 if(tmpRec==leading){
667 fh1PtTracksLeadingRecIn->Fill(tmpPt);
668 continue;
669 }
cc0649e4 670 // correlation
671 Float_t tmpPhi = tmpRec->Phi();
672
673 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
674 Float_t dPhi = phi - tmpRec->Phi();
675 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
676 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
677 Float_t dEta = eta - tmpRec->Eta();
678 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
565584e8 679 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
cc0649e4 680 }
681 delete recIter;
682 }
683
684 if(genParticles.GetSize()){
685 TIterator *genIter = genParticles.MakeIterator();
686 AliVParticle *tmpGen = 0;
687 while((tmpGen = (AliVParticle*)(genIter->Next()))){
688 if(TMath::Abs(tmpGen->Eta())<0.9){
689 Float_t tmpPt = tmpGen->Pt();
690 fh1PtTracksGenIn->Fill(tmpPt);
691 }
692 }
693 delete genIter;
694 }
695
3b7ffecf 696 fh1NRecJets->Fill(nRecJets);
697 nRecJets = TMath::Min(nRecJets,kMaxJets);
698
699 for(int ir = 0;ir < nRecJets;++ir){
700 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
701 if(!tmp)continue;
702 recJets[ir] = *tmp;
703 }
704
705
706 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
707 // Relate the jets
708 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
709 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
710
cc0649e4 711
712 for(int i = 0;i<kMaxJets;++i){
3b7ffecf 713 iGenIndex[i] = iRecIndex[i] = -1;
714 }
715
3b7ffecf 716 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
717 iGenIndex,iRecIndex,fDebug);
718 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
719
720 if(fDebug){
721 for(int i = 0;i<kMaxJets;++i){
722 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
723 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
724 }
725 }
726
727
728
729
730 Double_t container[6];
731
732 // loop over generated jets
733
734 // radius; tmp, get from the jet header itself
735 // at some point, how todeal woht FastJet on MC?
736 Float_t radiusGen =0.4;
737 Float_t radiusRec =0.4;
738
739 for(int ig = 0;ig < nGenJets;++ig){
740 Double_t ptGen = genJets[ig].Pt();
741 Double_t phiGen = genJets[ig].Phi();
742 if(phiGen<0)phiGen+=TMath::Pi()*2.;
743 Double_t etaGen = genJets[ig].Eta();
744
745 container[3] = ptGen;
746 container[4] = etaGen;
747 container[5] = phiGen;
748 fhnJetContainer[kStep0]->Fill(&container[3],eventW);
749 Int_t ir = iRecIndex[ig];
750
751 if(TMath::Abs(etaGen)<fRecEtaWindow){
752 fhnJetContainer[kStep1]->Fill(&container[3],eventW);
753 fh1PtGenIn[ig]->Fill(ptGen,eventW);
754 // fill the fragmentation function
755 for(int it = 0;it<genParticles.GetEntries();++it){
756 AliVParticle *part = (AliVParticle*)genParticles.At(it);
757 if(genJets[ig].DeltaR(part)<radiusGen){
758 Float_t z = part->Pt()/ptGen;
759 Float_t lnz = -1.*TMath::Log(z);
760 fh2FragGen[ig]->Fill(z,ptGen,eventW);
761 fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
762 }
763 }
764 if(ir>=0&&ir<nRecJets){
765 fhnJetContainer[kStep3]->Fill(&container[3],eventW);
766 }
767 }
768
769 if(ir>=0&&ir<nRecJets){
770 fhnJetContainer[kStep2]->Fill(&container[3],eventW);
771 Double_t etaRec = recJets[ir].Eta();
772 if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
773 }
774 }// loop over generated jets
775
7fa8b2da 776
edfbe476 777 Float_t sumPt = 0;
778 for(int it = 0;it<recParticles.GetEntries();++it){
779 AliVParticle *part = (AliVParticle*)recParticles.At(it);
780 // fill sum pt and P_t of all paritles
781 if(TMath::Abs(part->Eta())<0.9){
782 Float_t pt = part->Pt();
783 fh1PtTrackRec->Fill(pt,eventW);
784 sumPt += pt;
785 }
786 }
787 fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
788 fh1SumPtTrackRec->Fill(sumPt,eventW);
789
cf049e94 790
3b7ffecf 791 // loop over reconstructed jets
792 for(int ir = 0;ir < nRecJets;++ir){
793 Double_t etaRec = recJets[ir].Eta();
794 Double_t ptRec = recJets[ir].Pt();
795 Double_t phiRec = recJets[ir].Phi();
796 if(phiRec<0)phiRec+=TMath::Pi()*2.;
797 container[0] = ptRec;
798 container[1] = etaRec;
799 container[2] = phiRec;
800
3e3987f1 801 if(ptRec>10.&&fDebug>0){
7fa8b2da 802 // need to cast to int, otherwise the printf overwrites
803 Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
a923bd34 804 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
cf049e94 805 // aodH->SetFillAOD(kTRUE);
7fa8b2da 806 fAOD->GetHeader()->Print();
a221560c 807 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
7fa8b2da 808 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
809 AliAODTrack *tr = fAOD->GetTrack(it);
810 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
811 tr->Print();
38ecb6a5 812 // tr->Dump();
7fa8b2da 813 if(fESD){
814 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
815 esdTr->Print("");
bb3d13aa 816 // esdTr->Dump();
7fa8b2da 817 }
818 }
819 }
820
821
3b7ffecf 822 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
823 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
824 if(TMath::Abs(etaRec)<fRecEtaWindow){
825 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
826 fh1PtRecIn[ir]->Fill(ptRec,eventW);
827 // fill the fragmentation function
828 for(int it = 0;it<recParticles.GetEntries();++it){
829 AliVParticle *part = (AliVParticle*)recParticles.At(it);
edfbe476 830 Float_t eta = part->Eta();
831 if(TMath::Abs(eta)<0.9){
832 Float_t phi = part->Phi();
833 if(phi<0)phi+=TMath::Pi()*2.;
834 Float_t dPhi = phi - phiRec;
835 Float_t dEta = eta - etaRec;
cc0649e4 836 if(dPhi<(-1.*TMath::Pi()))phiRec+=TMath::Pi()*2.;
edfbe476 837 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
838 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
839 }
3b7ffecf 840 if(recJets[ir].DeltaR(part)<radiusRec){
841 Float_t z = part->Pt()/ptRec;
842 Float_t lnz = -1.*TMath::Log(z);
843 fh2FragRec[ir]->Fill(z,ptRec,eventW);
844 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
845 }
846 }
847 }
848
849
850 // Fill Correlation
851 Int_t ig = iGenIndex[ir];
852 if(ig>=0 && ig<nGenJets){
853 fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
854 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
855 Double_t ptGen = genJets[ig].Pt();
856 Double_t phiGen = genJets[ig].Phi();
857 if(phiGen<0)phiGen+=TMath::Pi()*2.;
858 Double_t etaGen = genJets[ig].Eta();
859
860 container[3] = ptGen;
861 container[4] = etaGen;
862 container[5] = phiGen;
863
864 //
865 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
866 //
867
868 if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
869 if(TMath::Abs(etaRec)<fRecEtaWindow){
870 fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
871 fhnCorrelation->Fill(container);
872 }// if etarec in window
873
874 }
875 ////////////////////////////////////////////////////
876 else{
877
878 }
879 }// loop over reconstructed jets
880
881
882 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
883 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
884 PostData(1, fHistList);
885}
886
887
888void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
889 //
890 // Create the particle container for the correction framework manager and
891 // link it
892 //
893 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
a923bd34 894 const Double_t kPtmin = 5.0, kPtmax = 105.; // we do not want to have empty bins at the beginning...
3b7ffecf 895 const Double_t kEtamin = -3.0, kEtamax = 3.0;
896 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
897
898 // can we neglect migration in eta and phi?
899 // phi should be no problem since we cover full phi and are phi symmetric
900 // eta migration is more difficult due to needed acceptance correction
901 // in limited eta range
902
903
904 //arrays for the number of bins in each dimension
905 Int_t iBin[kNvar];
906 iBin[0] = 100; //bins in pt
907 iBin[1] = 1; //bins in eta
908 iBin[2] = 1; // bins in phi
909
910 //arrays for lower bounds :
911 Double_t* binEdges[kNvar];
912 for(Int_t ivar = 0; ivar < kNvar; ivar++)
913 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
914
915 //values for bin lower bounds
916 // 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 917 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 918 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
919 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
920
921
922 for(int i = 0;i < kMaxStep*2;++i){
923 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d"),kNvar,iBin);
924 for (int k=0; k<kNvar; k++) {
925 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
926 }
927 }
928 //create correlation matrix for unfolding
929 Int_t thnDim[2*kNvar];
930 for (int k=0; k<kNvar; k++) {
931 //first half : reconstructed
932 //second half : MC
933 thnDim[k] = iBin[k];
934 thnDim[k+kNvar] = iBin[k];
935 }
936
937 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
938 for (int k=0; k<kNvar; k++) {
939 fhnCorrelation->SetBinEdges(k,binEdges[k]);
940 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
941 }
942 fhnCorrelation->Sumw2();
943
944 // Add a histogram for Fake jets
945 // thnDim[3] = AliPID::kSPECIES;
946 // fFakeElectrons = new THnSparseF("fakeEkectrons", "Output for Fake Electrons", kNvar + 1, thnDim);
947 // for(Int_t idim = 0; idim < kNvar; idim++)
948 // fFakeElectrons->SetBinEdges(idim, binEdges[idim]);
5cb0f438 949 for(Int_t ivar = 0; ivar < kNvar; ivar++)
950 delete [] binEdges[ivar];
951
3b7ffecf 952}
953
954void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
955{
956// Terminate analysis
957//
958 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
959}
960
961
962Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
963
565584e8 964 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
3b7ffecf 965
966 Int_t iCount = 0;
565584e8 967 if(type==kTrackAOD){
3b7ffecf 968 AliAODEvent *aod = 0;
565584e8 969 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
970 else aod = AODEvent();
3b7ffecf 971 if(!aod){
972 return iCount;
973 }
974 for(int it = 0;it < aod->GetNumberOfTracks();++it){
975 AliAODTrack *tr = aod->GetTrack(it);
976 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
cc0649e4 977 if(TMath::Abs(tr->Eta())>0.9)continue;
38ecb6a5 978 if(fDebug>0){
979 if(tr->Pt()>20){
980 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
981 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),fInputHandler->GetTree()->GetReadEntry());
982 tr->Print();
983 // tr->Dump();
984 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
985 if(fESD){
986 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
987 esdTr->Print("");
bb3d13aa 988 // esdTr->Dump();
38ecb6a5 989 }
990 }
991 }
3b7ffecf 992 list->Add(tr);
993 iCount++;
994 }
995 }
996 else if (type == kTrackKineAll||type == kTrackKineCharged){
997 AliMCEvent* mcEvent = MCEvent();
998 if(!mcEvent)return iCount;
999 // we want to have alivpartilces so use get track
1000 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1001 if(!mcEvent->IsPhysicalPrimary(it))continue;
1002 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1003 if(type == kTrackKineAll){
1004 list->Add(part);
1005 iCount++;
1006 }
1007 else if(type == kTrackKineCharged){
1008 if(part->Particle()->GetPDG()->Charge()==0)continue;
1009 list->Add(part);
1010 iCount++;
1011 }
1012 }
1013 }
565584e8 1014 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1015 AliAODEvent *aod = 0;
5301f754 1016 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
565584e8 1017 else aod = AODEvent();
5010a3f7 1018 if(!aod)return iCount;
3b7ffecf 1019 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1020 if(!tca)return iCount;
1021 for(int it = 0;it < tca->GetEntriesFast();++it){
1022 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1023 if(!part->IsPhysicalPrimary())continue;
c4631cdb 1024 if(type == kTrackAODMCAll){
3b7ffecf 1025 list->Add(part);
1026 iCount++;
1027 }
565584e8 1028 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
3b7ffecf 1029 if(part->Charge()==0)continue;
565584e8 1030 if(kTrackAODMCCharged){
1031 list->Add(part);
1032 }
1033 else {
1034 if(TMath::Abs(part->Eta())>0.9)continue;
1035 list->Add(part);
1036 }
3b7ffecf 1037 iCount++;
1038 }
1039 }
1040 }// AODMCparticle
cc0649e4 1041 list->Sort();
c4631cdb 1042 return iCount;
3b7ffecf 1043
1044}