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