1 // **************************************
2 // Task used for the correction of determiantion of reconstructed jet spectra
3 // Compares input (gen) and output (rec) jets
4 // *******************************************
7 /**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
26 #include <TInterpreter.h>
35 #include <TLorentzVector.h>
36 #include <TClonesArray.h>
37 #include "TDatabasePDG.h"
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"
55 #include "AliGenPythiaEventHeader.h"
56 #include "AliJetKineReaderHeader.h"
57 #include "AliGenCocktailEventHeader.h"
58 #include "AliInputEventHandler.h"
61 #include "AliAnalysisHelperJetTasks.h"
63 ClassImp(AliAnalysisTaskJetSpectrum2)
65 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
70 fhnCorrelationPhiZRec(0x0),
74 fUseAODJetInput(kFALSE),
75 fUseAODTrackInput(kFALSE),
76 fUseAODMCInput(kFALSE),
77 fUseGlobalSelection(kFALSE),
78 fUseExternalWeightOnly(kFALSE),
79 fLimitGenJetEta(kFALSE),
81 fEventSelectionMask(0),
83 fTrackTypeRec(kTrackUndef),
84 fTrackTypeGen(kTrackUndef),
89 fDeltaPhiWindow(20./180.*TMath::Pi()),
98 fh1SumPtTrackRec(0x0),
99 fh1SumPtTrackAreaRec(0x0),
102 fh1PtJetsLeadingRecIn(0x0),
103 fh1PtTracksRecIn(0x0),
104 fh1PtTracksLeadingRecIn(0x0),
105 fh1PtTracksGenIn(0x0),
107 fh2NRecTracksPt(0x0),
108 fh2JetsLeadingPhiEta(0x0),
109 fh2JetsLeadingPhiPt(0x0),
110 fh2TracksLeadingPhiEta(0x0),
111 fh2TracksLeadingPhiPt(0x0),
112 fh2TracksLeadingJetPhiPt(0x0),
114 fh2TrackPtTrackPhi(0x0),
116 fh2DijetDeltaPhiPt(0x0),
118 fh2DijetAsymPtCut(0x0),
119 fh2DijetDeltaPhiDeltaEta(0x0),
120 fh2DijetPt2vsPt1(0x0),
121 fh2DijetDifvsSum(0x0),
123 fh1DijetMinvCut(0x0),
126 for(int i = 0;i < kMaxStep*2;++i){
127 fhnJetContainer[i] = 0;
129 for(int i = 0;i < kMaxJets;++i){
130 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
146 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
147 AliAnalysisTaskSE(name),
152 fhnCorrelationPhiZRec(0x0),
156 fUseAODJetInput(kFALSE),
157 fUseAODTrackInput(kFALSE),
158 fUseAODMCInput(kFALSE),
159 fUseGlobalSelection(kFALSE),
160 fUseExternalWeightOnly(kFALSE),
161 fLimitGenJetEta(kFALSE),
163 fEventSelectionMask(0),
165 fTrackTypeRec(kTrackUndef),
166 fTrackTypeGen(kTrackUndef),
171 fDeltaPhiWindow(20./180.*TMath::Pi()),
176 fh1PtHardTrials(0x0),
180 fh1SumPtTrackRec(0x0),
181 fh1SumPtTrackAreaRec(0x0),
184 fh1PtJetsLeadingRecIn(0x0),
185 fh1PtTracksRecIn(0x0),
186 fh1PtTracksLeadingRecIn(0x0),
187 fh1PtTracksGenIn(0x0),
189 fh2NRecTracksPt(0x0),
190 fh2JetsLeadingPhiEta(0x0),
191 fh2JetsLeadingPhiPt(0x0),
192 fh2TracksLeadingPhiEta(0x0),
193 fh2TracksLeadingPhiPt(0x0),
194 fh2TracksLeadingJetPhiPt(0x0),
196 fh2TrackPtTrackPhi(0x0),
198 fh2DijetDeltaPhiPt(0x0),
200 fh2DijetAsymPtCut(0x0),
201 fh2DijetDeltaPhiDeltaEta(0x0),
202 fh2DijetPt2vsPt1(0x0),
203 fh2DijetDifvsSum(0x0),
205 fh1DijetMinvCut(0x0),
209 for(int i = 0;i < kMaxStep*2;++i){
210 fhnJetContainer[i] = 0;
212 for(int i = 0;i < kMaxJets;++i){
213 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
226 DefineOutput(1, TList::Class());
231 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
234 // Implemented Notify() to read the cross sections
235 // and number of trials from pyxsec.root
238 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
239 Float_t xsection = 0;
244 TFile *curfile = tree->GetCurrentFile();
246 Error("Notify","No current file");
249 if(!fh1Xsec||!fh1Trials){
250 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
253 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
254 fh1Xsec->Fill("<#sigma>",xsection);
255 // construct a poor man average trials
256 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
257 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
262 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
266 // Create the output container
276 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
279 if(!fHistList)fHistList = new TList();
280 fHistList->SetOwner(kTRUE);
281 Bool_t oldStatus = TH1::AddDirectoryStatus();
282 TH1::AddDirectory(kFALSE);
287 const Int_t nBinPt = 320;
288 Double_t binLimitsPt[nBinPt+1];
289 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
291 binLimitsPt[iPt] = 0.0;
294 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
298 const Int_t nBinPhi = 90;
299 Double_t binLimitsPhi[nBinPhi+1];
300 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
302 binLimitsPhi[iPhi] = -1.*TMath::Pi();
305 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
310 const Int_t nBinPhi2 = 360;
311 Double_t binLimitsPhi2[nBinPhi2+1];
312 for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
314 binLimitsPhi2[iPhi2] = 0.;
317 binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
323 const Int_t nBinEta = 40;
324 Double_t binLimitsEta[nBinEta+1];
325 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
327 binLimitsEta[iEta] = -2.0;
330 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
334 const Int_t nBinFrag = 25;
337 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
338 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
340 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
341 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
343 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
344 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
345 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
347 fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
348 fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
350 fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
351 fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
352 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);
354 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
355 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
356 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
357 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
358 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
360 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
361 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
364 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
365 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
366 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
367 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
368 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
369 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
370 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
371 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
373 fh2JetPtJetPhi = new TH2F("fh2JetPtJetPhi","Reconstructed jet phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
374 fh2TrackPtTrackPhi = new TH2F("fh2TrackPtTrackPhi","Reconstructed track phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
376 fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};#Delta p_{T}/p_{T,Gen}",nBinPt,binLimitsPt,120,-1.2,1.2);
378 for(int ij = 0;ij < kMaxJets;++ij){
379 fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
380 fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
382 fh2PhiPt[ij] = new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
383 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
385 fh2PhiEta[ij] = new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
386 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
388 fh2RhoPtRec[ij] = new TH2F(Form("fh2RhoPtRec_j%d",ij),"jet shape rho for jets;r;p_{T,rec};",
389 20,0.,1.,nBinPt,binLimitsPt);
390 fh2PsiPtRec[ij] = new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
391 20,0.,1.,nBinPt,binLimitsPt);
393 fh2RhoPtGen[ij] = new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
394 20,0.,1.,nBinPt,binLimitsPt);
395 fh2PsiPtGen[ij] = new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
396 20,0.,1.,nBinPt,binLimitsPt);
397 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",20,0.,1);
400 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",
401 nBinFrag,0.,1.,nBinPt,binLimitsPt);
402 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",
403 nBinFrag,0.,10.,nBinPt,binLimitsPt);
405 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",
406 nBinFrag,0.,1.,nBinPt,binLimitsPt);
407 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",
408 nBinFrag,0.,10.,nBinPt,binLimitsPt);
413 fh2DijetDeltaPhiPt = new TH2F("fh2DeltaPhiPt","Difference in the azimuthal angle;#Delta#phi;p_{T,1};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
414 fh2DijetAsymPt = new TH2F("fh2DijetAsym","Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
415 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);
416 fh2DijetDeltaPhiDeltaEta = new TH2F("fh2DijetDeltaPhiDeltaEta","Difference in the azimuthal angle;#Delta#phi;Entries",180,0.,TMath::Pi(),20,-2.,2.);
417 fh2DijetPt2vsPt1 = new TH2F("fh2DijetPt2vsPt1","Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
418 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.);
419 fh1DijetMinv = new TH1F("fh1DijetMinv","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
420 fh1DijetMinvCut = new TH1F("fh1DijetMinvCut","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
425 const Int_t saveLevel = 3; // large save level more histos
427 fHistList->Add(fh1Xsec);
428 fHistList->Add(fh1Trials);
429 fHistList->Add(fh1PtHard);
430 fHistList->Add(fh1PtHardNoW);
431 fHistList->Add(fh1PtHardTrials);
432 if(fBranchGen.Length()>0){
433 fHistList->Add(fh1NGenJets);
434 fHistList->Add(fh1PtTracksGenIn);
436 fHistList->Add(fh1PtJetsRecIn);
437 fHistList->Add(fh1PtJetsLeadingRecIn);
438 fHistList->Add(fh1PtTracksRecIn);
439 fHistList->Add(fh1PtTracksLeadingRecIn);
440 fHistList->Add(fh1NRecJets);
441 fHistList->Add(fh1PtTrackRec);
442 fHistList->Add(fh1SumPtTrackRec);
443 fHistList->Add(fh1SumPtTrackAreaRec);
444 fHistList->Add(fh2NRecJetsPt);
445 fHistList->Add(fh2NRecTracksPt);
446 fHistList->Add(fh2JetsLeadingPhiEta );
447 fHistList->Add(fh2JetsLeadingPhiPt );
448 fHistList->Add(fh2TracksLeadingPhiEta);
449 fHistList->Add(fh2TracksLeadingPhiPt);
450 for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
451 for(int ij = 0;ij<kMaxJets;++ij){
452 fHistList->Add( fh1PtRecIn[ij]);
454 if(fBranchGen.Length()>0){
455 fHistList->Add(fh1PtGenIn[ij]);
456 fHistList->Add(fh2FragGen[ij]);
457 fHistList->Add(fh2FragLnGen[ij]);
458 fHistList->Add(fh2RhoPtGen[ij]);
459 fHistList->Add(fh2PsiPtGen[ij]);
461 fHistList->Add( fh2PhiPt[ij]);
462 fHistList->Add( fh2PhiEta[ij]);
463 fHistList->Add( fh2RhoPtRec[ij]);
464 fHistList->Add( fh2PsiPtRec[ij]);
465 fHistList->Add( fh2FragRec[ij]);
466 fHistList->Add( fh2FragLnRec[ij]);
468 fHistList->Add(fhnCorrelation);
469 fHistList->Add(fhnCorrelationPhiZRec);
470 fHistList->Add(fh2JetPtJetPhi);
471 fHistList->Add(fh2TrackPtTrackPhi);
472 fHistList->Add(fh2RelPtFGen);
474 fHistList->Add(fh2DijetDeltaPhiPt);
475 fHistList->Add(fh2DijetAsymPt);
476 fHistList->Add(fh2DijetAsymPtCut);
477 fHistList->Add(fh2DijetDeltaPhiDeltaEta);
478 fHistList->Add(fh2DijetPt2vsPt1);
479 fHistList->Add(fh2DijetDifvsSum);
480 fHistList->Add(fh1DijetMinv);
481 fHistList->Add(fh1DijetMinvCut);
484 // =========== Switch on Sumw2 for all histos ===========
485 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
486 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
491 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
494 TH1::AddDirectory(oldStatus);
497 void AliAnalysisTaskJetSpectrum2::Init()
503 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
507 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
510 Bool_t selected = kTRUE;
512 if(fUseGlobalSelection&&fEventSelectionMask==0){
513 selected = AliAnalysisHelperJetTasks::Selected();
515 else if(fUseGlobalSelection&&fEventSelectionMask>0){
516 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
520 // no selection by the service task, we continue
521 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
522 PostData(1, fHistList);
528 // Execute analysis for current event
530 AliESDEvent *fESD = 0;
532 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
534 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
540 // assume that the AOD is in the general output...
543 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
547 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
554 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
557 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
560 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
564 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
565 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
567 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
571 // ==== General variables needed
574 // We use statice array, not to fragment the memory
575 AliAODJet genJets[kMaxJets];
577 AliAODJet recJets[kMaxJets];
579 ///////////////////////////
584 Double_t nTrials = 1; // Trials for MC trigger
586 if(fUseExternalWeightOnly){
587 eventW = fExternalWeight;
590 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
591 // if(fDebug>0)aodH->SetFillAOD(kFALSE);
592 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
593 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
594 // this is the part we only use when we have MC information
595 AliMCEvent* mcEvent = MCEvent();
596 // AliStack *pStack = 0;
598 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
601 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
604 nTrials = pythiaGenHeader->Trials();
605 ptHard = pythiaGenHeader->GetPtHard();
606 int iProcessType = pythiaGenHeader->ProcessType();
608 // 12 f+barf -> f+barf
613 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
614 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
616 // fetch the pythia generated jets only to be used here
617 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
618 AliAODJet pythiaGenJets[kMaxJets];
619 for(int ip = 0;ip < nPythiaGenJets;++ip){
620 if(iCount>=kMaxJets)continue;
622 pythiaGenHeader->TriggerJet(ip,p);
623 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
625 if(fBranchGen.Length()==0){
628 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
629 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
632 if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
633 // if we have MC particles and we do not read from the aod branch
634 // use the pythia jets
635 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
640 if(fBranchGen.Length()==0)nGenJets = iCount;
641 }// (fAnalysisType&kMCESD)==kMCESD)
644 // we simply fetch the tracks/mc particles as a list of AliVParticles
652 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
653 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
654 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
655 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
658 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
659 fh1PtHard->Fill(ptHard,eventW);
660 fh1PtHardNoW->Fill(ptHard,1);
661 fh1PtHardTrials->Fill(ptHard,nTrials);
663 // If we set a second branch for the input jets fetch this
664 if(fBranchGen.Length()>0){
665 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
668 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
669 if(iCount>=kMaxJets)continue;
670 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
674 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
675 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
678 if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
679 genJets[iCount] = *tmp;
685 if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
686 if(fDebug>2)fAOD->Print();
690 fh1NGenJets->Fill(nGenJets);
691 // We do not want to exceed the maximum jet number
692 nGenJets = TMath::Min(nGenJets,kMaxJets);
694 // Fetch the reconstructed jets...
696 nRecJets = aodRecJets->GetEntries();
698 nRecJets = aodRecJets->GetEntries();
699 fh1NRecJets->Fill(nRecJets);
701 // Do something with the all rec jets
702 Int_t nRecOver = nRecJets;
704 // check that the jets are sorted
705 Float_t ptOld = 999999;
706 for(int ir = 0;ir < nRecJets;ir++){
707 AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
708 Float_t tmpPt = tmp->Pt();
710 Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
717 TIterator *recIter = aodRecJets->MakeIterator();
718 AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
719 Float_t pt = tmpRec->Pt();
721 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
722 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
723 while(pt<ptCut&&tmpRec){
725 tmpRec = (AliAODJet*)(recIter->Next());
730 if(nRecOver<=0)break;
731 fh2NRecJetsPt->Fill(ptCut,nRecOver);
735 AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
736 Float_t phi = leading->Phi();
737 if(phi<0)phi+=TMath::Pi()*2.;
738 Float_t eta = leading->Eta();
740 while((tmpRec = (AliAODJet*)(recIter->Next()))){
741 Float_t tmpPt = tmpRec->Pt();
742 fh1PtJetsRecIn->Fill(tmpPt);
744 fh1PtJetsLeadingRecIn->Fill(tmpPt);
748 Float_t tmpPhi = tmpRec->Phi();
750 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
751 Float_t dPhi = phi - tmpRec->Phi();
752 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
753 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
754 Float_t dEta = eta - tmpRec->Eta();
755 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
756 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
761 Int_t nTrackOver = recParticles.GetSize();
762 // do the same for tracks and jets
764 TIterator *recIter = recParticles.MakeIterator();
765 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
766 Float_t pt = tmpRec->Pt();
767 // Printf("Leading track p_t %3.3E",pt);
768 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
769 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
770 while(pt<ptCut&&tmpRec){
772 tmpRec = (AliVParticle*)(recIter->Next());
777 if(nTrackOver<=0)break;
778 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
782 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
783 Float_t phi = leading->Phi();
784 if(phi<0)phi+=TMath::Pi()*2.;
785 Float_t eta = leading->Eta();
787 while((tmpRec = (AliVParticle*)(recIter->Next()))){
788 Float_t tmpPt = tmpRec->Pt();
789 fh1PtTracksRecIn->Fill(tmpPt);
791 fh1PtTracksLeadingRecIn->Fill(tmpPt);
795 Float_t tmpPhi = tmpRec->Phi();
797 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
798 Float_t dPhi = phi - tmpRec->Phi();
799 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
800 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
801 Float_t dEta = eta - tmpRec->Eta();
802 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
803 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
808 if(genParticles.GetSize()){
809 TIterator *genIter = genParticles.MakeIterator();
810 AliVParticle *tmpGen = 0;
811 while((tmpGen = (AliVParticle*)(genIter->Next()))){
812 if(TMath::Abs(tmpGen->Eta())<0.9){
813 Float_t tmpPt = tmpGen->Pt();
814 fh1PtTracksGenIn->Fill(tmpPt);
820 nRecJets = TMath::Min(nRecJets,kMaxJets);
823 for(int ir = 0;ir < nRecJets;++ir){
824 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
826 if(tmp->Pt()<fMinJetPt)continue;
830 nRecJets = iCountRec;
833 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
835 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
836 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
839 for(int i = 0;i<kMaxJets;++i){
840 iGenIndex[i] = iRecIndex[i] = -1;
843 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
844 iGenIndex,iRecIndex,fDebug);
845 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
848 for(int i = 0;i<kMaxJets;++i){
849 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
850 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
857 Double_t container[6];
858 Double_t containerPhiZ[6];
860 // loop over generated jets
862 // radius; tmp, get from the jet header itself
863 // at some point, how todeal woht FastJet on MC?
864 Float_t radiusGen =0.4;
865 Float_t radiusRec =0.4;
867 for(int ig = 0;ig < nGenJets;++ig){
868 Double_t ptGen = genJets[ig].Pt();
869 Double_t phiGen = genJets[ig].Phi();
870 if(phiGen<0)phiGen+=TMath::Pi()*2.;
871 Double_t etaGen = genJets[ig].Eta();
873 container[3] = ptGen;
874 container[4] = etaGen;
875 container[5] = phiGen;
876 fhnJetContainer[kStep0]->Fill(&container[3],eventW);
877 Int_t ir = iRecIndex[ig];
879 if(TMath::Abs(etaGen)<fRecEtaWindow){
882 fhnJetContainer[kStep1]->Fill(&container[3],eventW);
883 fh1PtGenIn[ig]->Fill(ptGen,eventW);
884 // fill the fragmentation function
885 for(int it = 0;it<genParticles.GetEntries();++it){
886 AliVParticle *part = (AliVParticle*)genParticles.At(it);
887 Float_t deltaR = genJets[ig].DeltaR(part);
888 fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
889 if(deltaR<radiusGen){
890 Float_t z = part->Pt()/ptGen;
891 Float_t lnz = -1.*TMath::Log(z);
892 fh2FragGen[ig]->Fill(z,ptGen,eventW);
893 fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
898 for(int ibx = 1;ibx<fh2RhoPtGen[ir]->GetNbinsX();ibx++){
899 Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
900 Float_t rho = fh1TmpRho->GetBinContent(ibx);
902 fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
903 fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
906 if(ir>=0&&ir<nRecJets){
907 fhnJetContainer[kStep2]->Fill(&container[3],eventW);
908 Double_t etaRec = recJets[ir].Eta();
909 if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
911 }// loop over generated jets
914 for(int it = 0;it<recParticles.GetEntries();++it){
915 AliVParticle *part = (AliVParticle*)recParticles.At(it);
916 // fill sum pt and P_t of all paritles
917 if(TMath::Abs(part->Eta())<0.9){
918 Float_t pt = part->Pt();
919 fh1PtTrackRec->Fill(pt,eventW);
920 fh2TrackPtTrackPhi->Fill(pt,part->Phi());
924 fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
925 fh1SumPtTrackRec->Fill(sumPt,eventW);
928 // loop over reconstructed jets
929 for(int ir = 0;ir < nRecJets;++ir){
930 Double_t etaRec = recJets[ir].Eta();
931 Double_t ptRec = recJets[ir].Pt();
932 Double_t phiRec = recJets[ir].Phi();
933 if(phiRec<0)phiRec+=TMath::Pi()*2.;
936 // do something with dijets...
938 Double_t etaRec1 = recJets[0].Eta();
939 Double_t ptRec1 = recJets[0].Pt();
940 Double_t phiRec1 = recJets[0].Phi();
941 if(phiRec1<0)phiRec1+=TMath::Pi()*2.;
944 if(TMath::Abs(etaRec1)<fRecEtaWindow
945 &&TMath::Abs(etaRec)<fRecEtaWindow){
947 Float_t deltaPhi = phiRec1 - phiRec;
949 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
950 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
951 deltaPhi = TMath::Abs(deltaPhi);
952 fh2DijetDeltaPhiPt->Fill(deltaPhi,ptRec1);
953 Float_t asym = (ptRec1-ptRec)/(ptRec1+ptRec);
954 fh2DijetAsymPt->Fill(asym,ptRec1);
955 fh2DijetDeltaPhiDeltaEta->Fill(deltaPhi,etaRec1-etaRec);
956 fh2DijetPt2vsPt1->Fill(ptRec1,ptRec);
957 fh2DijetDifvsSum->Fill(ptRec1+ptRec,ptRec1-ptRec);
958 Float_t minv = 2.*(recJets[0].P()*recJets[1].P()-
959 recJets[0].Px()*recJets[1].Px()-
960 recJets[0].Py()*recJets[1].Py()-
961 recJets[0].Pz()*recJets[1].Py());
962 minv = TMath::Sqrt(minv);
965 fh1DijetMinv->Fill(minv);
966 if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
967 fh1DijetMinvCut->Fill(minv);
968 fh2DijetAsymPtCut->Fill(asym,ptRec1);
974 container[0] = ptRec;
975 container[1] = etaRec;
976 container[2] = phiRec;
977 containerPhiZ[0] = ptRec;
978 containerPhiZ[1] = phiRec;
979 if(ptRec>30.&&fDebug>0){
980 // need to cast to int, otherwise the printf overwrites
981 Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
982 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
983 if(fESD)Printf("ESDEvent GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
984 // aodH->SetFillAOD(kTRUE);
985 fAOD->GetHeader()->Print();
986 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
987 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
988 AliAODTrack *tr = fAOD->GetTrack(it);
989 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
993 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1001 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
1002 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1004 Float_t zLeading = -1;
1005 if(TMath::Abs(etaRec)<fRecEtaWindow){
1006 fh2JetPtJetPhi->Fill(ptRec,phiRec);
1007 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
1008 fh1PtRecIn[ir]->Fill(ptRec,eventW);
1009 // fill the fragmentation function
1013 for(int it = 0;it<recParticles.GetEntries();++it){
1014 AliVParticle *part = (AliVParticle*)recParticles.At(it);
1015 Float_t eta = part->Eta();
1016 if(TMath::Abs(eta)<0.9){
1017 Float_t phi = part->Phi();
1018 if(phi<0)phi+=TMath::Pi()*2.;
1019 Float_t dPhi = phi - phiRec;
1020 Float_t dEta = eta - etaRec;
1021 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1022 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1023 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1024 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1027 Float_t deltaR = recJets[ir].DeltaR(part);
1028 fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
1031 if(deltaR<radiusRec){
1032 Float_t z = part->Pt()/ptRec;
1033 if(z>zLeading)zLeading=z;
1034 Float_t lnz = -1.*TMath::Log(z);
1035 fh2FragRec[ir]->Fill(z,ptRec,eventW);
1036 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1039 // fill the jet shapes
1041 for(int ibx = 1;ibx<fh2RhoPtRec[ir]->GetNbinsX();ibx++){
1042 Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1043 Float_t rho = fh1TmpRho->GetBinContent(ibx);
1045 fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1046 fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1051 containerPhiZ[2] = zLeading;
1054 Int_t ig = iGenIndex[ir];
1055 if(ig>=0 && ig<nGenJets){
1056 fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1057 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1058 Double_t ptGen = genJets[ig].Pt();
1059 Double_t phiGen = genJets[ig].Phi();
1060 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1061 Double_t etaGen = genJets[ig].Eta();
1063 container[3] = ptGen;
1064 container[4] = etaGen;
1065 container[5] = phiGen;
1066 containerPhiZ[3] = ptGen;
1068 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1071 if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1072 if(TMath::Abs(etaRec)<fRecEtaWindow){
1073 fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1074 fhnCorrelation->Fill(container);
1076 Float_t delta = (ptRec-ptGen)/ptGen;
1077 fh2RelPtFGen->Fill(ptGen,delta,eventW);
1079 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1081 }// if etarec in window
1085 containerPhiZ[3] = 0;
1086 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1088 }// loop over reconstructed jets
1091 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1092 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1093 PostData(1, fHistList);
1097 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1099 // Create the particle container for the correction framework manager and
1102 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
1103 const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1104 const Double_t kEtamin = -3.0, kEtamax = 3.0;
1105 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1106 const Double_t kZmin = 0., kZmax = 1;
1108 // can we neglect migration in eta and phi?
1109 // phi should be no problem since we cover full phi and are phi symmetric
1110 // eta migration is more difficult due to needed acceptance correction
1111 // in limited eta range
1113 //arrays for the number of bins in each dimension
1115 iBin[0] = 320; //bins in pt
1116 iBin[1] = 1; //bins in eta
1117 iBin[2] = 1; // bins in phi
1119 //arrays for lower bounds :
1120 Double_t* binEdges[kNvar];
1121 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1122 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1124 //values for bin lower bounds
1125 // 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);
1126 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1127 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1128 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1131 for(int i = 0;i < kMaxStep*2;++i){
1132 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1133 for (int k=0; k<kNvar; k++) {
1134 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1137 //create correlation matrix for unfolding
1138 Int_t thnDim[2*kNvar];
1139 for (int k=0; k<kNvar; k++) {
1140 //first half : reconstructed
1142 thnDim[k] = iBin[k];
1143 thnDim[k+kNvar] = iBin[k];
1146 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1147 for (int k=0; k<kNvar; k++) {
1148 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1149 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1151 fhnCorrelation->Sumw2();
1153 // for second correlation histogram
1156 const Int_t kNvarPhiZ = 4;
1157 //arrays for the number of bins in each dimension
1158 Int_t iBinPhiZ[kNvarPhiZ];
1159 iBinPhiZ[0] = 80; //bins in pt
1160 iBinPhiZ[1] = 72; //bins in phi
1161 iBinPhiZ[2] = 20; // bins in Z
1162 iBinPhiZ[3] = 80; //bins in ptgen
1164 //arrays for lower bounds :
1165 Double_t* binEdgesPhiZ[kNvarPhiZ];
1166 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1167 binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1169 for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1170 for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1171 for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1172 for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1174 fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1175 for (int k=0; k<kNvarPhiZ; k++) {
1176 fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1178 fhnCorrelationPhiZRec->Sumw2();
1181 // Add a histogram for Fake jets
1183 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1184 delete [] binEdges[ivar];
1186 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1187 delete [] binEdgesPhiZ[ivar];
1191 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1193 // Terminate analysis
1195 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1199 Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1201 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1204 if(type==kTrackAOD){
1205 AliAODEvent *aod = 0;
1206 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1207 else aod = AODEvent();
1211 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1212 AliAODTrack *tr = aod->GetTrack(it);
1213 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1214 if(TMath::Abs(tr->Eta())>0.9)continue;
1217 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1218 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1221 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1223 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1233 else if (type == kTrackKineAll||type == kTrackKineCharged){
1234 AliMCEvent* mcEvent = MCEvent();
1235 if(!mcEvent)return iCount;
1236 // we want to have alivpartilces so use get track
1237 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1238 if(!mcEvent->IsPhysicalPrimary(it))continue;
1239 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1240 if(type == kTrackKineAll){
1244 else if(type == kTrackKineCharged){
1245 if(part->Particle()->GetPDG()->Charge()==0)continue;
1251 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1252 AliAODEvent *aod = 0;
1253 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1254 else aod = AODEvent();
1255 if(!aod)return iCount;
1256 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1257 if(!tca)return iCount;
1258 for(int it = 0;it < tca->GetEntriesFast();++it){
1259 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1260 if(!part->IsPhysicalPrimary())continue;
1261 if(type == kTrackAODMCAll){
1265 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1266 if(part->Charge()==0)continue;
1267 if(kTrackAODMCCharged){
1271 if(TMath::Abs(part->Eta())>0.9)continue;