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 "AliAODJetEventBackground.h"
52 #include "AliAODMCParticle.h"
53 #include "AliMCEventHandler.h"
54 #include "AliMCEvent.h"
56 #include "AliGenPythiaEventHeader.h"
57 #include "AliJetKineReaderHeader.h"
58 #include "AliGenCocktailEventHeader.h"
59 #include "AliInputEventHandler.h"
62 #include "AliAnalysisHelperJetTasks.h"
64 ClassImp(AliAnalysisTaskJetSpectrum2)
66 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
71 fhnCorrelationPhiZRec(0x0),
76 fUseAODJetInput(kFALSE),
77 fUseAODTrackInput(kFALSE),
78 fUseAODMCInput(kFALSE),
79 fUseGlobalSelection(kFALSE),
80 fUseExternalWeightOnly(kFALSE),
81 fLimitGenJetEta(kFALSE),
82 fBkgSubtraction(kFALSE),
84 fEventSelectionMask(0),
86 fTrackTypeRec(kTrackUndef),
87 fTrackTypeGen(kTrackUndef),
92 fDeltaPhiWindow(20./180.*TMath::Pi()),
101 fh1SumPtTrackRec(0x0),
102 fh1SumPtTrackAreaRec(0x0),
105 fh1PtJetsLeadingRecIn(0x0),
106 fh1PtTracksRecIn(0x0),
107 fh1PtTracksLeadingRecIn(0x0),
108 fh1PtTracksGenIn(0x0),
110 fh2NRecTracksPt(0x0),
111 fh2JetsLeadingPhiEta(0x0),
112 fh2JetsLeadingPhiPt(0x0),
113 fh2TracksLeadingPhiEta(0x0),
114 fh2TracksLeadingPhiPt(0x0),
115 fh2TracksLeadingJetPhiPt(0x0),
117 fh2TrackPtTrackPhi(0x0),
119 fh2DijetDeltaPhiPt(0x0),
121 fh2DijetAsymPtCut(0x0),
122 fh2DijetDeltaPhiDeltaEta(0x0),
123 fh2DijetPt2vsPt1(0x0),
124 fh2DijetDifvsSum(0x0),
126 fh1DijetMinvCut(0x0),
136 fh1Ptjethardest(0x0),
137 fh1Ptjetsubhardest1(0x0),
138 fh1Ptjetsubhardest2(0x0),
139 fh1Rhovspthardest1(0x0),
140 fh1Rhovspthardest2(0x0),
141 fh1Errorvspthardest1(0x0),
142 fh1Errorvspthardest2(0x0),
145 for(int i = 0;i < kMaxStep*2;++i){
146 fhnJetContainer[i] = 0;
148 for(int i = 0;i < kMaxJets;++i){
149 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
165 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
166 AliAnalysisTaskSE(name),
171 fhnCorrelationPhiZRec(0x0),
176 fUseAODJetInput(kFALSE),
177 fUseAODTrackInput(kFALSE),
178 fUseAODMCInput(kFALSE),
179 fUseGlobalSelection(kFALSE),
180 fUseExternalWeightOnly(kFALSE),
181 fLimitGenJetEta(kFALSE),
182 fBkgSubtraction(kFALSE),
184 fEventSelectionMask(0),
186 fTrackTypeRec(kTrackUndef),
187 fTrackTypeGen(kTrackUndef),
192 fDeltaPhiWindow(20./180.*TMath::Pi()),
197 fh1PtHardTrials(0x0),
201 fh1SumPtTrackRec(0x0),
202 fh1SumPtTrackAreaRec(0x0),
205 fh1PtJetsLeadingRecIn(0x0),
206 fh1PtTracksRecIn(0x0),
207 fh1PtTracksLeadingRecIn(0x0),
208 fh1PtTracksGenIn(0x0),
210 fh2NRecTracksPt(0x0),
211 fh2JetsLeadingPhiEta(0x0),
212 fh2JetsLeadingPhiPt(0x0),
213 fh2TracksLeadingPhiEta(0x0),
214 fh2TracksLeadingPhiPt(0x0),
215 fh2TracksLeadingJetPhiPt(0x0),
217 fh2TrackPtTrackPhi(0x0),
219 fh2DijetDeltaPhiPt(0x0),
221 fh2DijetAsymPtCut(0x0),
222 fh2DijetDeltaPhiDeltaEta(0x0),
223 fh2DijetPt2vsPt1(0x0),
224 fh2DijetDifvsSum(0x0),
226 fh1DijetMinvCut(0x0),
236 fh1Ptjethardest(0x0),
237 fh1Ptjetsubhardest1(0x0),
238 fh1Ptjetsubhardest2(0x0),
239 fh1Rhovspthardest1(0x0),
240 fh1Rhovspthardest2(0x0),
241 fh1Errorvspthardest1(0x0),
242 fh1Errorvspthardest2(0x0),
246 for(int i = 0;i < kMaxStep*2;++i){
247 fhnJetContainer[i] = 0;
249 for(int i = 0;i < kMaxJets;++i){
250 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
263 DefineOutput(1, TList::Class());
268 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
271 // Implemented Notify() to read the cross sections
272 // and number of trials from pyxsec.root
275 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
276 Float_t xsection = 0;
281 TFile *curfile = tree->GetCurrentFile();
283 Error("Notify","No current file");
286 if(!fh1Xsec||!fh1Trials){
287 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
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();
294 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
299 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
303 // Create the output container
313 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
316 if(!fHistList)fHistList = new TList();
317 fHistList->SetOwner(kTRUE);
318 Bool_t oldStatus = TH1::AddDirectoryStatus();
319 TH1::AddDirectory(kFALSE);
324 const Int_t nBinPt = 320;
325 Double_t binLimitsPt[nBinPt+1];
326 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
328 binLimitsPt[iPt] = 0.0;
331 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
335 const Int_t nBinPhi = 90;
336 Double_t binLimitsPhi[nBinPhi+1];
337 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
339 binLimitsPhi[iPhi] = -1.*TMath::Pi();
342 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
347 const Int_t nBinPhi2 = 360;
348 Double_t binLimitsPhi2[nBinPhi2+1];
349 for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
351 binLimitsPhi2[iPhi2] = 0.;
354 binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
360 const Int_t nBinEta = 40;
361 Double_t binLimitsEta[nBinEta+1];
362 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
364 binLimitsEta[iEta] = -2.0;
367 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
371 const Int_t nBinFrag = 25;
374 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
375 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
377 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
378 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
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);
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);
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);
391 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
392 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
393 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
394 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
395 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
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);
401 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
402 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
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);
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);
413 fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};#Delta p_{T}/p_{T,Gen}",nBinPt,binLimitsPt,120,-1.2,1.2);
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);
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);
422 fh2PhiEta[ij] = new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
423 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
425 fh2RhoPtRec[ij] = new TH2F(Form("fh2RhoPtRec_j%d",ij),"jet shape rho for jets;r;p_{T,rec};",
426 40,0.,1.,nBinPt,binLimitsPt);
427 fh2PsiPtRec[ij] = new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
428 40,0.,2.,nBinPt,binLimitsPt);
430 fh2RhoPtGen[ij] = new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
431 40,0.,2.,nBinPt,binLimitsPt);
432 fh2PsiPtGen[ij] = new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
433 40,0.,2.,nBinPt,binLimitsPt);
434 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
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);
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",
443 nBinFrag,0.,1.,nBinPt,binLimitsPt);
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",
445 nBinFrag,0.,10.,nBinPt,binLimitsPt);
450 fh2DijetDeltaPhiPt = new TH2F("fh2DijetDeltaPhiPt","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);
458 //background histograms
460 fh1Bkg1 = new TH1F("fh1Bkg1","Background estimate 1",100,0.,10.);
461 fh1Bkg2 = new TH1F("fh1Bkg2","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 fh1Area1 = new TH1F("fh1Area1","Background mean area 1",50,0.,5.);
465 fh1Area2 = new TH1F("fh1Area2","Background mean area 2",50,0.,5.);
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.);
482 const Int_t saveLevel = 3; // large save level more histos
484 fHistList->Add(fh1Xsec);
485 fHistList->Add(fh1Trials);
486 fHistList->Add(fh1PtHard);
487 fHistList->Add(fh1PtHardNoW);
488 fHistList->Add(fh1PtHardTrials);
489 if(fBranchGen.Length()>0){
490 fHistList->Add(fh1NGenJets);
491 fHistList->Add(fh1PtTracksGenIn);
493 fHistList->Add(fh1PtJetsRecIn);
494 fHistList->Add(fh1PtJetsLeadingRecIn);
495 fHistList->Add(fh1PtTracksRecIn);
496 fHistList->Add(fh1PtTracksLeadingRecIn);
497 fHistList->Add(fh1NRecJets);
498 fHistList->Add(fh1PtTrackRec);
499 fHistList->Add(fh1SumPtTrackRec);
500 fHistList->Add(fh1SumPtTrackAreaRec);
501 fHistList->Add(fh2NRecJetsPt);
502 fHistList->Add(fh2NRecTracksPt);
503 fHistList->Add(fh2JetsLeadingPhiEta );
504 fHistList->Add(fh2JetsLeadingPhiPt );
505 fHistList->Add(fh2TracksLeadingPhiEta);
506 fHistList->Add(fh2TracksLeadingPhiPt);
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]);
511 if(fBranchGen.Length()>0){
512 fHistList->Add(fh1PtGenIn[ij]);
513 fHistList->Add(fh2FragGen[ij]);
514 fHistList->Add(fh2FragLnGen[ij]);
515 fHistList->Add(fh2RhoPtGen[ij]);
516 fHistList->Add(fh2PsiPtGen[ij]);
518 fHistList->Add( fh2PhiPt[ij]);
519 fHistList->Add( fh2PhiEta[ij]);
520 fHistList->Add( fh2RhoPtRec[ij]);
521 fHistList->Add( fh2PsiPtRec[ij]);
522 fHistList->Add( fh2FragRec[ij]);
523 fHistList->Add( fh2FragLnRec[ij]);
525 fHistList->Add(fhnCorrelation);
526 fHistList->Add(fhnCorrelationPhiZRec);
527 fHistList->Add(fh2JetPtJetPhi);
528 fHistList->Add(fh2TrackPtTrackPhi);
529 fHistList->Add(fh2RelPtFGen);
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);
540 fHistList->Add(fh1Bkg1);
541 fHistList->Add(fh1Bkg2);
542 fHistList->Add(fh1Sigma1);
543 fHistList->Add(fh1Sigma2);
544 fHistList->Add(fh1Area1);
545 fHistList->Add(fh1Area2);
546 fHistList->Add(fh1Ptjet);
547 fHistList->Add(fh1Ptjethardest);
548 fHistList->Add(fh1Ptjetsub1);
549 fHistList->Add(fh1Ptjetsub2);
551 fHistList->Add(fh1Ptjetsubhardest1);
552 fHistList->Add(fh1Ptjetsubhardest2);
553 fHistList->Add(fh1Rhovspthardest1);
554 fHistList->Add(fh1Rhovspthardest2);
556 fHistList->Add(fh1Errorvspthardest1);
557 fHistList->Add(fh1Errorvspthardest2);
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));
568 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
571 TH1::AddDirectory(oldStatus);
574 void AliAnalysisTaskJetSpectrum2::Init()
580 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
584 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
587 Bool_t selected = kTRUE;
589 if(fUseGlobalSelection&&fEventSelectionMask==0){
590 selected = AliAnalysisHelperJetTasks::Selected();
592 else if(fUseGlobalSelection&&fEventSelectionMask>0){
593 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
597 // no selection by the service task, we continue
598 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
599 PostData(1, fHistList);
605 // Execute analysis for current event
607 AliESDEvent *fESD = 0;
609 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
611 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
617 // assume that the AOD is in the general output...
620 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
624 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
631 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
634 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
637 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
641 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
642 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
644 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
649 AliAODJetEventBackground* evBkg=(AliAODJetEventBackground*)fAOD->FindListObject(fBranchBkg.Data());
653 Printf("%s:%d no reconstructed background array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkg.Data());
658 ///just to start: some very simple plots containing rho, sigma and area of the background.
660 Int_t nJets = aodRecJets->GetEntriesFast();
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 fh1Bkg1->Fill(bkg1); //rho computed with all background jets.
672 fh1Bkg2->Fill(bkg2); //rho computed with all background jets but the hardest.
673 fh1Sigma1->Fill(sigma1);
674 fh1Sigma2->Fill(sigma2);
675 fh1Area1->Fill(area1);
676 fh1Area2->Fill(area2);
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);
688 fh1Ptjetsub1->Fill(ptsub1);
689 fh1Ptjetsub2->Fill(ptsub2);
693 fh1Ptjethardest->Fill(pthardest);
694 fh1Ptjetsubhardest1->Fill(ptsub1);
695 fh1Ptjetsubhardest2->Fill(ptsub2);
696 fh1Errorvspthardest1->Fill(ptsub1,err1/ptsub1);
697 fh1Errorvspthardest2->Fill(ptsub2,err2/ptsub2);
700 fh1Rhovspthardest1->Fill(pthardest,bkg1);
701 fh1Rhovspthardest2->Fill(pthardest,bkg2);
703 }// background subtraction
706 // ==== General variables needed
709 // We use statice array, not to fragment the memory
710 AliAODJet genJets[kMaxJets];
712 AliAODJet recJets[kMaxJets];
714 ///////////////////////////
719 Double_t nTrials = 1; // Trials for MC trigger
721 if(fUseExternalWeightOnly){
722 eventW = fExternalWeight;
725 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
726 // if(fDebug>0)aodH->SetFillAOD(kFALSE);
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;
733 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
736 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
739 nTrials = pythiaGenHeader->Trials();
740 ptHard = pythiaGenHeader->GetPtHard();
741 int iProcessType = pythiaGenHeader->ProcessType();
743 // 12 f+barf -> f+barf
748 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
749 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
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;
757 pythiaGenHeader->TriggerJet(ip,p);
758 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
760 if(fBranchGen.Length()==0){
763 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
764 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
767 if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
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]);
775 if(fBranchGen.Length()==0)nGenJets = iCount;
776 }// (fAnalysisType&kMCESD)==kMCESD)
779 // we simply fetch the tracks/mc particles as a list of AliVParticles
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());
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);
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()));
803 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
804 if(iCount>=kMaxJets)continue;
805 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
809 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
810 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
813 if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
814 genJets[iCount] = *tmp;
820 if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
821 if(fDebug>2)fAOD->Print();
825 fh1NGenJets->Fill(nGenJets);
826 // We do not want to exceed the maximum jet number
827 nGenJets = TMath::Min(nGenJets,kMaxJets);
829 // Fetch the reconstructed jets...
831 nRecJets = aodRecJets->GetEntries();
833 nRecJets = aodRecJets->GetEntries();
834 fh1NRecJets->Fill(nRecJets);
836 // Do something with the all rec jets
837 Int_t nRecOver = nRecJets;
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();
845 Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
852 TIterator *recIter = aodRecJets->MakeIterator();
853 AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
854 Float_t pt = tmpRec->Pt();
856 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
857 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
858 while(pt<ptCut&&tmpRec){
860 tmpRec = (AliAODJet*)(recIter->Next());
865 if(nRecOver<=0)break;
866 fh2NRecJetsPt->Fill(ptCut,nRecOver);
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();
875 while((tmpRec = (AliAODJet*)(recIter->Next()))){
876 Float_t tmpPt = tmpRec->Pt();
877 fh1PtJetsRecIn->Fill(tmpPt);
879 fh1PtJetsLeadingRecIn->Fill(tmpPt);
883 Float_t tmpPhi = tmpRec->Phi();
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);
891 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
896 Int_t nTrackOver = recParticles.GetSize();
897 // do the same for tracks and jets
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){
907 tmpRec = (AliVParticle*)(recIter->Next());
912 if(nTrackOver<=0)break;
913 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
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();
922 while((tmpRec = (AliVParticle*)(recIter->Next()))){
923 Float_t tmpPt = tmpRec->Pt();
924 fh1PtTracksRecIn->Fill(tmpPt);
926 fh1PtTracksLeadingRecIn->Fill(tmpPt);
930 Float_t tmpPhi = tmpRec->Phi();
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);
938 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
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);
955 nRecJets = TMath::Min(nRecJets,kMaxJets);
958 for(int ir = 0;ir < nRecJets;++ir){
959 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
961 if(tmp->Pt()<fMinJetPt)continue;
965 nRecJets = iCountRec;
968 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
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
974 for(int i = 0;i<kMaxJets;++i){
975 iGenIndex[i] = iRecIndex[i] = -1;
978 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
979 iGenIndex,iRecIndex,fDebug);
980 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
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]);
992 Double_t container[6];
993 Double_t containerPhiZ[6];
995 // loop over generated jets
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;
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();
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];
1014 if(TMath::Abs(etaGen)<fRecEtaWindow){
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);
1022 Float_t deltaR = genJets[ig].DeltaR(part);
1023 fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
1024 if(deltaR<radiusGen){
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);
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);
1037 fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
1038 fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
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);
1046 }// loop over generated jets
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);
1055 fh2TrackPtTrackPhi->Fill(pt,part->Phi());
1059 fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
1060 fh1SumPtTrackRec->Fill(sumPt,eventW);
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.;
1071 // do something with dijets...
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.;
1079 if(TMath::Abs(etaRec1)<fRecEtaWindow
1080 &&TMath::Abs(etaRec)<fRecEtaWindow){
1082 Float_t deltaPhi = phiRec1 - phiRec;
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);
1100 fh1DijetMinv->Fill(minv);
1101 if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
1102 fh1DijetMinvCut->Fill(minv);
1103 fh2DijetAsymPtCut->Fill(asym,ptRec1);
1109 container[0] = ptRec;
1110 container[1] = etaRec;
1111 container[2] = phiRec;
1112 containerPhiZ[0] = ptRec;
1113 containerPhiZ[1] = phiRec;
1114 if(ptRec>30.&&fDebug>0){
1115 // need to cast to int, otherwise the printf overwrites
1116 Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
1117 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
1118 if(fESD)Printf("ESDEvent GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
1119 // aodH->SetFillAOD(kTRUE);
1120 fAOD->GetHeader()->Print();
1121 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
1122 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1123 AliAODTrack *tr = fAOD->GetTrack(it);
1124 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1128 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1136 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
1137 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1139 Float_t zLeading = -1;
1140 if(TMath::Abs(etaRec)<fRecEtaWindow){
1141 fh2JetPtJetPhi->Fill(ptRec,phiRec);
1142 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
1143 fh1PtRecIn[ir]->Fill(ptRec,eventW);
1144 // fill the fragmentation function
1148 for(int it = 0;it<recParticles.GetEntries();++it){
1149 AliVParticle *part = (AliVParticle*)recParticles.At(it);
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;
1156 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1157 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1158 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1159 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1162 Float_t deltaR = recJets[ir].DeltaR(part);
1163 fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
1166 if(deltaR<radiusRec){
1167 Float_t z = part->Pt()/ptRec;
1168 if(z>zLeading)zLeading=z;
1169 Float_t lnz = -1.*TMath::Log(z);
1170 fh2FragRec[ir]->Fill(z,ptRec,eventW);
1171 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1174 // fill the jet shapes
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);
1180 fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1181 fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1186 containerPhiZ[2] = zLeading;
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();
1198 container[3] = ptGen;
1199 container[4] = etaGen;
1200 container[5] = phiGen;
1201 containerPhiZ[3] = ptGen;
1203 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
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);
1211 Float_t delta = (ptRec-ptGen)/ptGen;
1212 fh2RelPtFGen->Fill(ptGen,delta,eventW);
1214 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1216 }// if etarec in window
1220 containerPhiZ[3] = 0;
1221 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1223 }// loop over reconstructed jets
1226 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1227 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1228 PostData(1, fHistList);
1232 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1234 // Create the particle container for the correction framework manager and
1237 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
1238 const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1239 const Double_t kEtamin = -3.0, kEtamax = 3.0;
1240 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1241 const Double_t kZmin = 0., kZmax = 1;
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
1248 //arrays for the number of bins in each dimension
1250 iBin[0] = 320; //bins in pt
1251 iBin[1] = 1; //bins in eta
1252 iBin[2] = 1; // bins in phi
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];
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);
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;
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;
1266 for(int i = 0;i < kMaxStep*2;++i){
1267 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1268 for (int k=0; k<kNvar; k++) {
1269 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1272 //create correlation matrix for unfolding
1273 Int_t thnDim[2*kNvar];
1274 for (int k=0; k<kNvar; k++) {
1275 //first half : reconstructed
1277 thnDim[k] = iBin[k];
1278 thnDim[k+kNvar] = iBin[k];
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]);
1286 fhnCorrelation->Sumw2();
1288 // for second correlation histogram
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
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];
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;
1309 fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1310 for (int k=0; k<kNvarPhiZ; k++) {
1311 fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1313 fhnCorrelationPhiZRec->Sumw2();
1316 // Add a histogram for Fake jets
1318 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1319 delete [] binEdges[ivar];
1321 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1322 delete [] binEdgesPhiZ[ivar];
1326 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1328 // Terminate analysis
1330 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1334 Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1336 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1339 if(type==kTrackAOD){
1340 AliAODEvent *aod = 0;
1341 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1342 else aod = AODEvent();
1346 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1347 AliAODTrack *tr = aod->GetTrack(it);
1348 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1349 if(TMath::Abs(tr->Eta())>0.9)continue;
1352 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1353 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1356 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1358 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
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){
1379 else if(type == kTrackKineCharged){
1380 if(part->Particle()->GetPDG()->Charge()==0)continue;
1386 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1387 AliAODEvent *aod = 0;
1388 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1389 else aod = AODEvent();
1390 if(!aod)return iCount;
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;
1396 if(type == kTrackAODMCAll){
1400 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1401 if(part->Charge()==0)continue;
1402 if(kTrackAODMCCharged){
1406 if(TMath::Abs(part->Eta())>0.9)continue;