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