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