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