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