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()),
102 fh1SumPtTrackRec(0x0),
103 fh1SumPtTrackAreaRec(0x0),
106 fh1PtJetsLeadingRecIn(0x0),
107 fh1PtTracksRecIn(0x0),
108 fh1PtTracksLeadingRecIn(0x0),
109 fh1PtTracksGenIn(0x0),
111 fh2NRecTracksPt(0x0),
112 fh2JetsLeadingPhiEta(0x0),
113 fh2JetsLeadingPhiPt(0x0),
114 fh2TracksLeadingPhiEta(0x0),
115 fh2TracksLeadingPhiPt(0x0),
116 fh2TracksLeadingJetPhiPt(0x0),
118 fh2TrackPtTrackPhi(0x0),
120 fh2DijetDeltaPhiPt(0x0),
122 fh2DijetAsymPtCut(0x0),
123 fh2DijetDeltaPhiDeltaEta(0x0),
124 fh2DijetPt2vsPt1(0x0),
125 fh2DijetDifvsSum(0x0),
127 fh1DijetMinvCut(0x0),
141 fh1Ptjethardest(0x0),
142 fh1Ptjetsubhardest1(0x0),
143 fh1Ptjetsubhardest2(0x0),
144 fh1Ptjetsubhardest3(0x0),
145 fh1Rhovspthardest1(0x0),
146 fh1Rhovspthardest2(0x0),
147 fh1Rhovspthardest3(0x0),
148 fh1Errorvspthardest1(0x0),
149 fh1Errorvspthardest2(0x0),
150 fh1Errorvspthardest3(0x0),
153 for(int i = 0;i < kMaxStep*2;++i){
154 fhnJetContainer[i] = 0;
156 for(int i = 0;i < kMaxJets;++i){
157 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
173 AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
174 AliAnalysisTaskSE(name),
179 fhnCorrelationPhiZRec(0x0),
184 fUseAODJetInput(kFALSE),
185 fUseAODTrackInput(kFALSE),
186 fUseAODMCInput(kFALSE),
187 fUseGlobalSelection(kFALSE),
188 fUseExternalWeightOnly(kFALSE),
189 fLimitGenJetEta(kFALSE),
190 fBkgSubtraction(kFALSE),
192 fEventSelectionMask(0),
194 fTrackTypeRec(kTrackUndef),
195 fTrackTypeGen(kTrackUndef),
200 fDeltaPhiWindow(20./180.*TMath::Pi()),
205 fh1PtHardTrials(0x0),
210 fh1SumPtTrackRec(0x0),
211 fh1SumPtTrackAreaRec(0x0),
214 fh1PtJetsLeadingRecIn(0x0),
215 fh1PtTracksRecIn(0x0),
216 fh1PtTracksLeadingRecIn(0x0),
217 fh1PtTracksGenIn(0x0),
219 fh2NRecTracksPt(0x0),
220 fh2JetsLeadingPhiEta(0x0),
221 fh2JetsLeadingPhiPt(0x0),
222 fh2TracksLeadingPhiEta(0x0),
223 fh2TracksLeadingPhiPt(0x0),
224 fh2TracksLeadingJetPhiPt(0x0),
226 fh2TrackPtTrackPhi(0x0),
228 fh2DijetDeltaPhiPt(0x0),
230 fh2DijetAsymPtCut(0x0),
231 fh2DijetDeltaPhiDeltaEta(0x0),
232 fh2DijetPt2vsPt1(0x0),
233 fh2DijetDifvsSum(0x0),
235 fh1DijetMinvCut(0x0),
249 fh1Ptjethardest(0x0),
250 fh1Ptjetsubhardest1(0x0),
251 fh1Ptjetsubhardest2(0x0),
252 fh1Ptjetsubhardest3(0x0),
253 fh1Rhovspthardest1(0x0),
254 fh1Rhovspthardest2(0x0),
255 fh1Rhovspthardest3(0x0),
256 fh1Errorvspthardest1(0x0),
257 fh1Errorvspthardest2(0x0),
258 fh1Errorvspthardest3(0x0),
262 for(int i = 0;i < kMaxStep*2;++i){
263 fhnJetContainer[i] = 0;
265 for(int i = 0;i < kMaxJets;++i){
266 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
279 DefineOutput(1, TList::Class());
284 Bool_t AliAnalysisTaskJetSpectrum2::Notify()
287 // Implemented Notify() to read the cross sections
288 // and number of trials from pyxsec.root
291 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
292 Float_t xsection = 0;
297 TFile *curfile = tree->GetCurrentFile();
299 Error("Notify","No current file");
302 if(!fh1Xsec||!fh1Trials){
303 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
306 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
307 fh1Xsec->Fill("<#sigma>",xsection);
308 // construct a poor man average trials
309 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
310 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
315 void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
319 // Create the output container
329 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
332 if(!fHistList)fHistList = new TList();
333 fHistList->SetOwner(kTRUE);
334 Bool_t oldStatus = TH1::AddDirectoryStatus();
335 TH1::AddDirectory(kFALSE);
340 const Int_t nBinPt = 320;
341 Double_t binLimitsPt[nBinPt+1];
342 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
344 binLimitsPt[iPt] = 0.0;
347 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
351 const Int_t nBinPhi = 90;
352 Double_t binLimitsPhi[nBinPhi+1];
353 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
355 binLimitsPhi[iPhi] = -1.*TMath::Pi();
358 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
363 const Int_t nBinPhi2 = 360;
364 Double_t binLimitsPhi2[nBinPhi2+1];
365 for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
367 binLimitsPhi2[iPhi2] = 0.;
370 binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
376 const Int_t nBinEta = 40;
377 Double_t binLimitsEta[nBinEta+1];
378 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
380 binLimitsEta[iEta] = -2.0;
383 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
387 const Int_t nBinFrag = 25;
390 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
391 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
393 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
394 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
396 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
397 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
398 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
400 fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
402 fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
403 fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
405 fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
406 fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
407 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);
409 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
410 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
411 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
412 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
413 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
415 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
416 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
419 fh2JetsLeadingPhiEta = new TH2F("fh2JetsLeadingPhiEta","delta eta vs delta phi to leading jet;#Delta#phi;#Delta#eta",
420 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
421 fh2JetsLeadingPhiPt = new TH2F("fh2JetsLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
422 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
423 fh2TracksLeadingPhiEta = new TH2F("fh2TracksLeadingPhiEta","delta eta vs delta phi to leading track;#Delta#phi;#Delta#eta",
424 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
425 fh2TracksLeadingPhiPt = new TH2F("fh2TracksLeadingPhiPt","leading p_T vs delta phi to leading jet;#Delta#phi;p_{T} (GeV/c)",
426 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
428 fh2JetPtJetPhi = new TH2F("fh2JetPtJetPhi","Reconstructed jet phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
429 fh2TrackPtTrackPhi = new TH2F("fh2TrackPtTrackPhi","Reconstructed track phi vs. pt",nBinPt,binLimitsPt,nBinPhi2,binLimitsPhi2);
431 fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};#Delta p_{T}/p_{T,Gen}",nBinPt,binLimitsPt,120,-1.2,1.2);
433 for(int ij = 0;ij < kMaxJets;++ij){
434 fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
435 fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
437 fh2PhiPt[ij] = new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
438 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
440 fh2PhiEta[ij] = new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
441 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
443 fh2RhoPtRec[ij] = new TH2F(Form("fh2RhoPtRec_j%d",ij),"jet shape rho for jets;r;p_{T,rec};",
444 40,0.,1.,nBinPt,binLimitsPt);
445 fh2PsiPtRec[ij] = new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
446 40,0.,2.,nBinPt,binLimitsPt);
448 fh2RhoPtGen[ij] = new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
449 40,0.,2.,nBinPt,binLimitsPt);
450 fh2PsiPtGen[ij] = new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
451 40,0.,2.,nBinPt,binLimitsPt);
452 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
455 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",
456 nBinFrag,0.,1.,nBinPt,binLimitsPt);
457 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",
458 nBinFrag,0.,10.,nBinPt,binLimitsPt);
460 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",
461 nBinFrag,0.,1.,nBinPt,binLimitsPt);
462 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",
463 nBinFrag,0.,10.,nBinPt,binLimitsPt);
468 fh2DijetDeltaPhiPt = new TH2F("fh2DijetDeltaPhiPt","Difference in the azimuthal angle;#Delta#phi;p_{T,1};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
469 fh2DijetAsymPt = new TH2F("fh2DijetAsym","Pt asymmetry;#Deltap_{T}/(p_{T,1}+p_{T,2});p_{T,1};Entries",50,0.,1.,nBinPt,binLimitsPt);
470 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);
471 fh2DijetDeltaPhiDeltaEta = new TH2F("fh2DijetDeltaPhiDeltaEta","Difference in the azimuthal angle;#Delta#phi;Entries",180,0.,TMath::Pi(),20,-2.,2.);
472 fh2DijetPt2vsPt1 = new TH2F("fh2DijetPt2vsPt1","Pt2 versus Pt1;p_{T,1} (GeV/c);p_{T,2} (GeV/c)",250,0.,250.,250,0.,250.);
473 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.);
474 fh1DijetMinv = new TH1F("fh1DijetMinv","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
475 fh1DijetMinvCut = new TH1F("fh1DijetMinvCut","Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
476 //background histograms
478 fh1Bkg1 = new TH1F("fh1Bkg1","Background estimate 1",100,0.,10.);
479 fh1Bkg2 = new TH1F("fh1Bkg2","Background estimate 2",100,0.,10.);
480 fh1Bkg3 = new TH1F("fh1Bkg3","Background estimate 3",100,0.,10.);
481 fh1Sigma1 = new TH1F("fh1Sigma1","Background fluctuations 1",100,0.,10.);
482 fh1Sigma2 = new TH1F("fh1Sigma2","Background fluctuations 2",100,0.,10.);
483 fh1Sigma3 = new TH1F("fh1Sigma3","Background fluctuations 3",100,0.,10.);
484 fh1Area1 = new TH1F("fh1Area1","Background mean area 1",50,0.,5.);
485 fh1Area2 = new TH1F("fh1Area2","Background mean area 2",50,0.,5.);
486 fh1Area3 = new TH1F("fh1Area3","Background mean area 3",50,0.,5.);
487 fh1Ptjet = new TH1F("fh1Ptjet","Jet spectrum",100,0.,200.);
488 fh1Ptjetsub1 = new TH1F("fh1Ptjetsub1","Subtracted spectrum 1",50,0.,200.);
489 fh1Ptjetsub2 = new TH1F("fh1Ptjetsub2","Subtracted spectrum 2",50,0.,200.);
490 fh1Ptjetsub3 = new TH1F("fh1Ptjetsub3","Subtracted spectrum 3",50,0.,200.);
491 fh1Ptjethardest = new TH1F("fh1Ptjethardest","Hardest jet spectrum",50,0.,200.);
492 fh1Ptjetsubhardest1 = new TH1F("fh1Pthardestsub1","Subtracted hardest jet spectrum 1",100,0.,200.);
493 fh1Ptjetsubhardest2 = new TH1F("fh1Pthardestsub2","Subtracted hardest jet spectrum 2",100,0.,200.);
494 fh1Ptjetsubhardest3 = new TH1F("fh1Pthardestsub3","Subtracted hardest jet spectrum 3",100,0.,200.);
495 fh1Rhovspthardest1 = new TH2F("fh1Rhovspthardest1","Background vs pTjet 1",100,0.,200.,50,0.,5.);
496 fh1Rhovspthardest2 = new TH2F("fh1Rhovspthardest2","Background vs pTjet 2",100,0.,200.,50,0.,5.);
497 fh1Rhovspthardest3 = new TH2F("fh1Rhovspthardest3","Background vs pTjet 3",100,0.,200.,50,0.,5.);
498 fh1Errorvspthardest1 = new TH2F("fhErorvspthardest1","Relative error 1",100,0.,200.,50,0.,5.);
499 fh1Errorvspthardest2 = new TH2F("fh1Errorvspthardest2","Relative error 2",100,0.,200.,50,0.,5.);
500 fh1Errorvspthardest3 = new TH2F("fh1Errorvspthardest3","Relative error 3",100,0.,200.,50,0.,5.);
506 const Int_t saveLevel = 3; // large save level more histos
508 fHistList->Add(fh1Xsec);
509 fHistList->Add(fh1Trials);
510 fHistList->Add(fh1PtHard);
511 fHistList->Add(fh1PtHardNoW);
512 fHistList->Add(fh1PtHardTrials);
513 fHistList->Add(fh1ZVtx);
514 if(fBranchGen.Length()>0){
515 fHistList->Add(fh1NGenJets);
516 fHistList->Add(fh1PtTracksGenIn);
518 fHistList->Add(fh1PtJetsRecIn);
519 fHistList->Add(fh1PtJetsLeadingRecIn);
520 fHistList->Add(fh1PtTracksRecIn);
521 fHistList->Add(fh1PtTracksLeadingRecIn);
522 fHistList->Add(fh1NRecJets);
523 fHistList->Add(fh1PtTrackRec);
524 fHistList->Add(fh1SumPtTrackRec);
525 fHistList->Add(fh1SumPtTrackAreaRec);
526 fHistList->Add(fh2NRecJetsPt);
527 fHistList->Add(fh2NRecTracksPt);
528 fHistList->Add(fh2JetsLeadingPhiEta );
529 fHistList->Add(fh2JetsLeadingPhiPt );
530 fHistList->Add(fh2TracksLeadingPhiEta);
531 fHistList->Add(fh2TracksLeadingPhiPt);
532 for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
533 for(int ij = 0;ij<kMaxJets;++ij){
534 fHistList->Add( fh1PtRecIn[ij]);
536 if(fBranchGen.Length()>0){
537 fHistList->Add(fh1PtGenIn[ij]);
538 fHistList->Add(fh2FragGen[ij]);
539 fHistList->Add(fh2FragLnGen[ij]);
540 fHistList->Add(fh2RhoPtGen[ij]);
541 fHistList->Add(fh2PsiPtGen[ij]);
543 fHistList->Add( fh2PhiPt[ij]);
544 fHistList->Add( fh2PhiEta[ij]);
545 fHistList->Add( fh2RhoPtRec[ij]);
546 fHistList->Add( fh2PsiPtRec[ij]);
547 fHistList->Add( fh2FragRec[ij]);
548 fHistList->Add( fh2FragLnRec[ij]);
550 fHistList->Add(fhnCorrelation);
551 fHistList->Add(fhnCorrelationPhiZRec);
552 fHistList->Add(fh2JetPtJetPhi);
553 fHistList->Add(fh2TrackPtTrackPhi);
554 fHistList->Add(fh2RelPtFGen);
556 fHistList->Add(fh2DijetDeltaPhiPt);
557 fHistList->Add(fh2DijetAsymPt);
558 fHistList->Add(fh2DijetAsymPtCut);
559 fHistList->Add(fh2DijetDeltaPhiDeltaEta);
560 fHistList->Add(fh2DijetPt2vsPt1);
561 fHistList->Add(fh2DijetDifvsSum);
562 fHistList->Add(fh1DijetMinv);
563 fHistList->Add(fh1DijetMinvCut);
565 fHistList->Add(fh1Bkg1);
566 fHistList->Add(fh1Bkg2);
567 fHistList->Add(fh1Bkg3);
568 fHistList->Add(fh1Sigma1);
569 fHistList->Add(fh1Sigma2);
570 fHistList->Add(fh1Sigma3);
571 fHistList->Add(fh1Area1);
572 fHistList->Add(fh1Area2);
573 fHistList->Add(fh1Area3);
574 fHistList->Add(fh1Ptjet);
575 fHistList->Add(fh1Ptjethardest);
576 fHistList->Add(fh1Ptjetsub1);
577 fHistList->Add(fh1Ptjetsub2);
578 fHistList->Add(fh1Ptjetsub3);
579 fHistList->Add(fh1Ptjetsubhardest1);
580 fHistList->Add(fh1Ptjetsubhardest2);
581 fHistList->Add(fh1Ptjetsubhardest3);
582 fHistList->Add(fh1Rhovspthardest1);
583 fHistList->Add(fh1Rhovspthardest2);
584 fHistList->Add(fh1Rhovspthardest3);
585 fHistList->Add(fh1Errorvspthardest1);
586 fHistList->Add(fh1Errorvspthardest2);
587 fHistList->Add(fh1Errorvspthardest3);
591 // =========== Switch on Sumw2 for all histos ===========
592 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
593 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
598 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
601 TH1::AddDirectory(oldStatus);
604 void AliAnalysisTaskJetSpectrum2::Init()
610 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
614 void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/)
617 Bool_t selected = kTRUE;
619 if(fUseGlobalSelection&&fEventSelectionMask==0){
620 selected = AliAnalysisHelperJetTasks::Selected();
622 else if(fUseGlobalSelection&&fEventSelectionMask>0){
623 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
627 // no selection by the service task, we continue
628 if (fDebug > 1)Printf("Not selected %s:%d",(char*)__FILE__,__LINE__);
629 PostData(1, fHistList);
635 // Execute analysis for current event
637 AliESDEvent *fESD = 0;
639 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
641 Printf("%s:%d AODEvent not found in Input Manager %d",(char*)__FILE__,__LINE__,fUseAODJetInput);
647 // assume that the AOD is in the general output...
650 Printf("%s:%d AODEvent not found in the Output",(char*)__FILE__,__LINE__);
653 fESD = dynamic_cast<AliESDEvent*> (InputEvent());
659 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
662 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
665 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
669 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
670 TClonesArray *aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
672 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
677 AliAODJetEventBackground* evBkg=(AliAODJetEventBackground*)fAOD->FindListObject(fBranchBkg.Data());
681 Printf("%s:%d no reconstructed background array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkg.Data());
686 ///just to start: some very simple plots containing rho, sigma and area of the background.
688 Int_t nJets = aodRecJets->GetEntriesFast();
689 Float_t pthardest=0.;
692 Float_t bkg1=evBkg->GetBackground(0);
693 Float_t bkg2=evBkg->GetBackground(1);
694 Float_t bkg3=evBkg->GetBackground(2);
695 Float_t sigma1=evBkg->GetSigma(0);
696 Float_t sigma2=evBkg->GetSigma(1);
697 Float_t sigma3=evBkg->GetSigma(2);
698 Float_t area1=evBkg->GetMeanarea(0);
699 Float_t area2=evBkg->GetMeanarea(1);
700 Float_t area3=evBkg->GetMeanarea(2);
701 fh1Bkg1->Fill(bkg1); //rho computed with all background jets.
702 fh1Bkg2->Fill(bkg2); //rho computed with all background jets but the 2 hardest.
703 fh1Bkg3->Fill(bkg3); //rho computed with randomized jets
704 fh1Sigma1->Fill(sigma1);
705 fh1Sigma2->Fill(sigma2);
706 fh1Sigma3->Fill(sigma3);
707 fh1Area1->Fill(area1);
708 fh1Area2->Fill(area2);
709 fh1Area3->Fill(area3);
712 for(Int_t k=0;k<nJets;k++){
713 AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets->At(k));
714 fh1Ptjet->Fill(jet->Pt());
715 Float_t ptsub1=jet->Pt()-bkg1*jet->EffectiveAreaCharged();
716 Float_t ptsub2=jet->Pt()-bkg2*jet->EffectiveAreaCharged();
717 Float_t ptsub3=jet->Pt()-bkg3*jet->EffectiveAreaCharged();
718 Float_t err1=sigma1*sqrt(area1);
719 Float_t err2=sigma2*sqrt(area2);
720 Float_t err3=sigma3*sqrt(area3);
721 fh1Ptjetsub1->Fill(ptsub1);
722 fh1Ptjetsub2->Fill(ptsub2);
723 fh1Ptjetsub3->Fill(ptsub3);
726 fh1Ptjethardest->Fill(pthardest);
727 fh1Ptjetsubhardest1->Fill(ptsub1);
728 fh1Ptjetsubhardest2->Fill(ptsub2);
729 fh1Ptjetsubhardest3->Fill(ptsub3);
730 fh1Errorvspthardest1->Fill(ptsub1,err1/ptsub1);
731 fh1Errorvspthardest2->Fill(ptsub2,err2/ptsub2);
732 fh1Errorvspthardest3->Fill(ptsub3,err3/ptsub3);
735 fh1Rhovspthardest1->Fill(pthardest,bkg1);
736 fh1Rhovspthardest2->Fill(pthardest,bkg2);
737 fh1Rhovspthardest3->Fill(pthardest,bkg3);
742 }// background subtraction
745 // ==== General variables needed
748 // We use statice array, not to fragment the memory
749 AliAODJet genJets[kMaxJets];
751 AliAODJet recJets[kMaxJets];
753 ///////////////////////////
758 Double_t nTrials = 1; // Trials for MC trigger
760 if(fUseExternalWeightOnly){
761 eventW = fExternalWeight;
764 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
765 // if(fDebug>0)aodH->SetFillAOD(kFALSE);
766 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
767 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
768 // this is the part we only use when we have MC information
769 AliMCEvent* mcEvent = MCEvent();
770 // AliStack *pStack = 0;
772 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
775 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
778 nTrials = pythiaGenHeader->Trials();
779 ptHard = pythiaGenHeader->GetPtHard();
780 int iProcessType = pythiaGenHeader->ProcessType();
782 // 12 f+barf -> f+barf
787 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
788 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
790 // fetch the pythia generated jets only to be used here
791 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
792 AliAODJet pythiaGenJets[kMaxJets];
793 for(int ip = 0;ip < nPythiaGenJets;++ip){
794 if(iCount>=kMaxJets)continue;
796 pythiaGenHeader->TriggerJet(ip,p);
797 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
799 if(fBranchGen.Length()==0){
802 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
803 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
806 if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
807 // if we have MC particles and we do not read from the aod branch
808 // use the pythia jets
809 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
814 if(fBranchGen.Length()==0)nGenJets = iCount;
815 }// (fAnalysisType&kMCESD)==kMCESD)
818 // we simply fetch the tracks/mc particles as a list of AliVParticles
826 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
827 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
828 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
829 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
832 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
833 fh1PtHard->Fill(ptHard,eventW);
834 fh1PtHardNoW->Fill(ptHard,1);
835 fh1PtHardTrials->Fill(ptHard,nTrials);
836 fh1ZVtx->Fill(fAOD->GetPrimaryVertex()->GetZ());
839 // If we set a second branch for the input jets fetch this
840 if(fBranchGen.Length()>0){
841 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
844 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
845 if(iCount>=kMaxJets)continue;
846 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
850 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
851 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
854 if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
855 genJets[iCount] = *tmp;
861 if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
862 if(fDebug>2)fAOD->Print();
866 fh1NGenJets->Fill(nGenJets);
867 // We do not want to exceed the maximum jet number
868 nGenJets = TMath::Min(nGenJets,kMaxJets);
870 // Fetch the reconstructed jets...
872 nRecJets = aodRecJets->GetEntries();
874 nRecJets = aodRecJets->GetEntries();
875 fh1NRecJets->Fill(nRecJets);
877 // Do something with the all rec jets
878 Int_t nRecOver = nRecJets;
880 // check that the jets are sorted
881 Float_t ptOld = 999999;
882 for(int ir = 0;ir < nRecJets;ir++){
883 AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
884 Float_t tmpPt = tmp->Pt();
886 Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
893 TIterator *recIter = aodRecJets->MakeIterator();
894 AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
895 Float_t pt = tmpRec->Pt();
897 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
898 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
899 while(pt<ptCut&&tmpRec){
901 tmpRec = (AliAODJet*)(recIter->Next());
906 if(nRecOver<=0)break;
907 fh2NRecJetsPt->Fill(ptCut,nRecOver);
911 AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
912 Float_t phi = leading->Phi();
913 if(phi<0)phi+=TMath::Pi()*2.;
914 Float_t eta = leading->Eta();
916 while((tmpRec = (AliAODJet*)(recIter->Next()))){
917 Float_t tmpPt = tmpRec->Pt();
918 fh1PtJetsRecIn->Fill(tmpPt);
920 fh1PtJetsLeadingRecIn->Fill(tmpPt);
924 Float_t tmpPhi = tmpRec->Phi();
926 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
927 Float_t dPhi = phi - tmpRec->Phi();
928 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
929 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
930 Float_t dEta = eta - tmpRec->Eta();
931 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
932 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
937 Int_t nTrackOver = recParticles.GetSize();
938 // do the same for tracks and jets
940 TIterator *recIter = recParticles.MakeIterator();
941 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
942 Float_t pt = tmpRec->Pt();
943 // Printf("Leading track p_t %3.3E",pt);
944 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
945 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
946 while(pt<ptCut&&tmpRec){
948 tmpRec = (AliVParticle*)(recIter->Next());
953 if(nTrackOver<=0)break;
954 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
958 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
959 Float_t phi = leading->Phi();
960 if(phi<0)phi+=TMath::Pi()*2.;
961 Float_t eta = leading->Eta();
963 while((tmpRec = (AliVParticle*)(recIter->Next()))){
964 Float_t tmpPt = tmpRec->Pt();
965 fh1PtTracksRecIn->Fill(tmpPt);
967 fh1PtTracksLeadingRecIn->Fill(tmpPt);
971 Float_t tmpPhi = tmpRec->Phi();
973 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
974 Float_t dPhi = phi - tmpRec->Phi();
975 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
976 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
977 Float_t dEta = eta - tmpRec->Eta();
978 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
979 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
984 if(genParticles.GetSize()){
985 TIterator *genIter = genParticles.MakeIterator();
986 AliVParticle *tmpGen = 0;
987 while((tmpGen = (AliVParticle*)(genIter->Next()))){
988 if(TMath::Abs(tmpGen->Eta())<0.9){
989 Float_t tmpPt = tmpGen->Pt();
990 fh1PtTracksGenIn->Fill(tmpPt);
996 nRecJets = TMath::Min(nRecJets,kMaxJets);
999 for(int ir = 0;ir < nRecJets;++ir){
1000 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
1002 if(tmp->Pt()<fMinJetPt)continue;
1006 nRecJets = iCountRec;
1009 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1011 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
1012 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
1015 for(int i = 0;i<kMaxJets;++i){
1016 iGenIndex[i] = iRecIndex[i] = -1;
1019 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
1020 iGenIndex,iRecIndex,fDebug);
1021 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1024 for(int i = 0;i<kMaxJets;++i){
1025 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
1026 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
1033 Double_t container[6];
1034 Double_t containerPhiZ[6];
1036 // loop over generated jets
1038 // radius; tmp, get from the jet header itself
1039 // at some point, how todeal woht FastJet on MC?
1040 Float_t radiusGen =0.4;
1041 Float_t radiusRec =0.4;
1043 for(int ig = 0;ig < nGenJets;++ig){
1044 Double_t ptGen = genJets[ig].Pt();
1045 Double_t phiGen = genJets[ig].Phi();
1046 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1047 Double_t etaGen = genJets[ig].Eta();
1049 container[3] = ptGen;
1050 container[4] = etaGen;
1051 container[5] = phiGen;
1052 fhnJetContainer[kStep0]->Fill(&container[3],eventW);
1053 Int_t ir = iRecIndex[ig];
1055 if(TMath::Abs(etaGen)<fRecEtaWindow){
1058 fhnJetContainer[kStep1]->Fill(&container[3],eventW);
1059 fh1PtGenIn[ig]->Fill(ptGen,eventW);
1060 // fill the fragmentation function
1061 for(int it = 0;it<genParticles.GetEntries();++it){
1062 AliVParticle *part = (AliVParticle*)genParticles.At(it);
1063 Float_t deltaR = genJets[ig].DeltaR(part);
1064 fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
1065 if(deltaR<radiusGen){
1066 Float_t z = part->Pt()/ptGen;
1067 Float_t lnz = -1.*TMath::Log(z);
1068 fh2FragGen[ig]->Fill(z,ptGen,eventW);
1069 fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
1074 for(int ibx = 1;ibx<=fh2RhoPtGen[ir]->GetNbinsX();ibx++){
1075 Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
1076 Float_t rho = fh1TmpRho->GetBinContent(ibx);
1078 fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
1079 fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
1082 if(ir>=0&&ir<nRecJets){
1083 fhnJetContainer[kStep2]->Fill(&container[3],eventW);
1084 Double_t etaRec = recJets[ir].Eta();
1085 if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
1086 if(TMath::Abs(etaRec)<fRecEtaWindow&&TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep3]->Fill(&container[3],eventW);
1088 }// loop over generated jets
1091 for(int it = 0;it<recParticles.GetEntries();++it){
1092 AliVParticle *part = (AliVParticle*)recParticles.At(it);
1093 // fill sum pt and P_t of all paritles
1094 if(TMath::Abs(part->Eta())<0.9){
1095 Float_t pt = part->Pt();
1096 fh1PtTrackRec->Fill(pt,eventW);
1097 fh2TrackPtTrackPhi->Fill(pt,part->Phi());
1101 fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
1102 fh1SumPtTrackRec->Fill(sumPt,eventW);
1105 // loop over reconstructed jets
1106 for(int ir = 0;ir < nRecJets;++ir){
1107 Double_t etaRec = recJets[ir].Eta();
1108 Double_t ptRec = recJets[ir].Pt();
1109 Double_t phiRec = recJets[ir].Phi();
1110 if(phiRec<0)phiRec+=TMath::Pi()*2.;
1113 // do something with dijets...
1115 Double_t etaRec1 = recJets[0].Eta();
1116 Double_t ptRec1 = recJets[0].Pt();
1117 Double_t phiRec1 = recJets[0].Phi();
1118 if(phiRec1<0)phiRec1+=TMath::Pi()*2.;
1121 if(TMath::Abs(etaRec1)<fRecEtaWindow
1122 &&TMath::Abs(etaRec)<fRecEtaWindow){
1124 Float_t deltaPhi = phiRec1 - phiRec;
1126 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1127 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1128 deltaPhi = TMath::Abs(deltaPhi);
1129 fh2DijetDeltaPhiPt->Fill(deltaPhi,ptRec1);
1130 Float_t asym = (ptRec1-ptRec)/(ptRec1+ptRec);
1131 fh2DijetAsymPt->Fill(asym,ptRec1);
1132 fh2DijetDeltaPhiDeltaEta->Fill(deltaPhi,etaRec1-etaRec);
1133 fh2DijetPt2vsPt1->Fill(ptRec1,ptRec);
1134 fh2DijetDifvsSum->Fill(ptRec1+ptRec,ptRec1-ptRec);
1135 Float_t minv = 2.*(recJets[0].P()*recJets[1].P()-
1136 recJets[0].Px()*recJets[1].Px()-
1137 recJets[0].Py()*recJets[1].Py()-
1138 recJets[0].Pz()*recJets[1].Py());
1139 minv = TMath::Sqrt(minv);
1142 fh1DijetMinv->Fill(minv);
1143 if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
1144 fh1DijetMinvCut->Fill(minv);
1145 fh2DijetAsymPtCut->Fill(asym,ptRec1);
1151 container[0] = ptRec;
1152 container[1] = etaRec;
1153 container[2] = phiRec;
1154 containerPhiZ[0] = ptRec;
1155 containerPhiZ[1] = phiRec;
1156 if(ptRec>30.&&fDebug>0){
1157 // need to cast to int, otherwise the printf overwrites
1158 Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
1159 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
1160 if(fESD)Printf("ESDEvent GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
1161 // aodH->SetFillAOD(kTRUE);
1162 fAOD->GetHeader()->Print();
1163 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
1164 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1165 AliAODTrack *tr = fAOD->GetTrack(it);
1166 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1170 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1178 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
1179 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1181 Float_t zLeading = -1;
1182 if(TMath::Abs(etaRec)<fRecEtaWindow){
1183 fh2JetPtJetPhi->Fill(ptRec,phiRec);
1184 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
1185 fh1PtRecIn[ir]->Fill(ptRec,eventW);
1186 // fill the fragmentation function
1190 for(int it = 0;it<recParticles.GetEntries();++it){
1191 AliVParticle *part = (AliVParticle*)recParticles.At(it);
1192 Float_t eta = part->Eta();
1193 if(TMath::Abs(eta)<0.9){
1194 Float_t phi = part->Phi();
1195 if(phi<0)phi+=TMath::Pi()*2.;
1196 Float_t dPhi = phi - phiRec;
1197 Float_t dEta = eta - etaRec;
1198 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1199 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1200 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1201 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1204 Float_t deltaR = recJets[ir].DeltaR(part);
1205 fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
1208 if(deltaR<radiusRec){
1209 Float_t z = part->Pt()/ptRec;
1210 if(z>zLeading)zLeading=z;
1211 Float_t lnz = -1.*TMath::Log(z);
1212 fh2FragRec[ir]->Fill(z,ptRec,eventW);
1213 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1216 // fill the jet shapes
1218 for(int ibx = 1;ibx<=fh2RhoPtRec[ir]->GetNbinsX();ibx++){
1219 Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1220 Float_t rho = fh1TmpRho->GetBinContent(ibx);
1222 fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1223 fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1228 containerPhiZ[2] = zLeading;
1231 Int_t ig = iGenIndex[ir];
1232 if(ig>=0 && ig<nGenJets){
1233 fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1234 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1235 Double_t ptGen = genJets[ig].Pt();
1236 Double_t phiGen = genJets[ig].Phi();
1237 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1238 Double_t etaGen = genJets[ig].Eta();
1240 container[3] = ptGen;
1241 container[4] = etaGen;
1242 container[5] = phiGen;
1243 containerPhiZ[3] = ptGen;
1245 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1248 if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1249 if(TMath::Abs(etaRec)<fRecEtaWindow){
1250 fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1251 fhnCorrelation->Fill(container);
1253 Float_t delta = (ptRec-ptGen)/ptGen;
1254 fh2RelPtFGen->Fill(ptGen,delta,eventW);
1256 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1258 }// if etarec in window
1262 containerPhiZ[3] = 0;
1263 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1265 }// loop over reconstructed jets
1268 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1269 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1270 PostData(1, fHistList);
1274 void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1276 // Create the particle container for the correction framework manager and
1279 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
1280 const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
1281 const Double_t kEtamin = -3.0, kEtamax = 3.0;
1282 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1283 const Double_t kZmin = 0., kZmax = 1;
1285 // can we neglect migration in eta and phi?
1286 // phi should be no problem since we cover full phi and are phi symmetric
1287 // eta migration is more difficult due to needed acceptance correction
1288 // in limited eta range
1290 //arrays for the number of bins in each dimension
1292 iBin[0] = 320; //bins in pt
1293 iBin[1] = 1; //bins in eta
1294 iBin[2] = 1; // bins in phi
1296 //arrays for lower bounds :
1297 Double_t* binEdges[kNvar];
1298 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1299 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1301 //values for bin lower bounds
1302 // 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);
1303 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1304 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1305 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1308 for(int i = 0;i < kMaxStep*2;++i){
1309 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
1310 for (int k=0; k<kNvar; k++) {
1311 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1314 //create correlation matrix for unfolding
1315 Int_t thnDim[2*kNvar];
1316 for (int k=0; k<kNvar; k++) {
1317 //first half : reconstructed
1319 thnDim[k] = iBin[k];
1320 thnDim[k+kNvar] = iBin[k];
1323 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1324 for (int k=0; k<kNvar; k++) {
1325 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1326 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1328 fhnCorrelation->Sumw2();
1330 // for second correlation histogram
1333 const Int_t kNvarPhiZ = 4;
1334 //arrays for the number of bins in each dimension
1335 Int_t iBinPhiZ[kNvarPhiZ];
1336 iBinPhiZ[0] = 80; //bins in pt
1337 iBinPhiZ[1] = 72; //bins in phi
1338 iBinPhiZ[2] = 20; // bins in Z
1339 iBinPhiZ[3] = 80; //bins in ptgen
1341 //arrays for lower bounds :
1342 Double_t* binEdgesPhiZ[kNvarPhiZ];
1343 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1344 binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1346 for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1347 for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1348 for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1349 for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1351 fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1352 for (int k=0; k<kNvarPhiZ; k++) {
1353 fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1355 fhnCorrelationPhiZRec->Sumw2();
1358 // Add a histogram for Fake jets
1360 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1361 delete [] binEdges[ivar];
1363 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1364 delete [] binEdgesPhiZ[ivar];
1368 void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1370 // Terminate analysis
1372 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1376 Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1378 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1381 if(type==kTrackAOD){
1382 AliAODEvent *aod = 0;
1383 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1384 else aod = AODEvent();
1388 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1389 AliAODTrack *tr = aod->GetTrack(it);
1390 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1391 if(TMath::Abs(tr->Eta())>0.9)continue;
1394 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1395 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1398 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1400 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1410 else if (type == kTrackKineAll||type == kTrackKineCharged){
1411 AliMCEvent* mcEvent = MCEvent();
1412 if(!mcEvent)return iCount;
1413 // we want to have alivpartilces so use get track
1414 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1415 if(!mcEvent->IsPhysicalPrimary(it))continue;
1416 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1417 if(type == kTrackKineAll){
1421 else if(type == kTrackKineCharged){
1422 if(part->Particle()->GetPDG()->Charge()==0)continue;
1428 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1429 AliAODEvent *aod = 0;
1430 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1431 else aod = AODEvent();
1432 if(!aod)return iCount;
1433 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1434 if(!tca)return iCount;
1435 for(int it = 0;it < tca->GetEntriesFast();++it){
1436 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1437 if(!part->IsPhysicalPrimary())continue;
1438 if(type == kTrackAODMCAll){
1442 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1443 if(part->Charge()==0)continue;
1444 if(kTrackAODMCCharged){
1448 if(TMath::Abs(part->Eta())>0.9)continue;