]> 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
fast global merger installed. To use the previous version of the global merger, set...
[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.;
80a790fd 711
c2785065 712 if(nJets!=0){
7f9012c1 713
c2785065 714 Float_t bkg1=evBkg->GetBackground(0);
715 Float_t bkg2=evBkg->GetBackground(1);
7f9012c1 716 Float_t bkg3=evBkg->GetBackground(2);
c2785065 717 Float_t sigma1=evBkg->GetSigma(0);
718 Float_t sigma2=evBkg->GetSigma(1);
7f9012c1 719 Float_t sigma3=evBkg->GetSigma(2);
c2785065 720 Float_t area1=evBkg->GetMeanarea(0);
721 Float_t area2=evBkg->GetMeanarea(1);
7f9012c1 722 Float_t area3=evBkg->GetMeanarea(2);
185fa8e6 723 fh1Bkg1->Fill(bkg1); //rho computed with all background jets.
7f9012c1 724 fh1Bkg2->Fill(bkg2); //rho computed with all background jets but the 2 hardest.
725 fh1Bkg3->Fill(bkg3); //rho computed with randomized jets
c2785065 726 fh1Sigma1->Fill(sigma1);
727 fh1Sigma2->Fill(sigma2);
7f9012c1 728 fh1Sigma3->Fill(sigma3);
185fa8e6 729 fh1Area1->Fill(area1);
730 fh1Area2->Fill(area2);
7f9012c1 731 fh1Area3->Fill(area3);
c2785065 732
80a790fd 733 Int_t iSubJetCounter = 0;
c2785065 734 for(Int_t k=0;k<nJets;k++){
735 AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets->At(k));
736 fh1Ptjet->Fill(jet->Pt());
737 Float_t ptsub1=jet->Pt()-bkg1*jet->EffectiveAreaCharged();
738 Float_t ptsub2=jet->Pt()-bkg2*jet->EffectiveAreaCharged();
7f9012c1 739 Float_t ptsub3=jet->Pt()-bkg3*jet->EffectiveAreaCharged();
6bd3fdae 740 if(ptsub2<0.) ptsub2=0.;
741 if(ptsub1<0.) ptsub1=0.;
742 if(ptsub3<0.) ptsub3=0.;
c2785065 743 Float_t err1=sigma1*sqrt(area1);
744 Float_t err2=sigma2*sqrt(area2);
7f9012c1 745 Float_t err3=sigma3*sqrt(area3);
c2785065 746 fh1Ptjetsub1->Fill(ptsub1);
747 fh1Ptjetsub2->Fill(ptsub2);
7f9012c1 748 fh1Ptjetsub3->Fill(ptsub3);
c2785065 749 if(k==0) {
750 pthardest=jet->Pt();
751 fh1Ptjethardest->Fill(pthardest);
752 fh1Ptjetsubhardest1->Fill(ptsub1);
753 fh1Ptjetsubhardest2->Fill(ptsub2);
7f9012c1 754 fh1Ptjetsubhardest3->Fill(ptsub3);
6bd3fdae 755 fh2Errorvspthardest1->Fill(ptsub1,err1/ptsub1);
756 fh2Errorvspthardest2->Fill(ptsub2,err2/ptsub2);
757 fh2Errorvspthardest3->Fill(ptsub3,err3/ptsub3);
c2785065 758 }
7f9012c1 759
f4132e7d 760 Float_t ptsub=0.;
761 if(fFillCorrBkg==1) ptsub=ptsub1;
762 if(fFillCorrBkg==2) ptsub=ptsub2;
763 if(fFillCorrBkg==3) ptsub=ptsub3;
80a790fd 764 if(ptsub>0){// avoid unphysical jets pT
765 Float_t subphi=jet->Phi();
766 Float_t subtheta=jet->Theta();
767 Float_t subpz = ptsub/TMath::Tan(subtheta);
768 Float_t subpx=ptsub*TMath::Cos(subphi);
769 Float_t subpy=ptsub * TMath::Sin(subphi);
770 Float_t subp = TMath::Sqrt(ptsub*ptsub+subpz*subpz);
771 if(k<kMaxJets){// only store the jets which are assoicated to anohter one
772 genJets[iSubJetCounter].SetPxPyPzE(subpx,subpy,subpz,subp);
773 iSubJetCounter++;
774 }
f4132e7d 775 }
6bd3fdae 776 }
80a790fd 777 nGenJets = iSubJetCounter;
6bd3fdae 778 fh2Rhovspthardest1->Fill(pthardest,bkg1);
779 fh2Rhovspthardest2->Fill(pthardest,bkg2);
780 fh2Rhovspthardest3->Fill(pthardest,bkg3);
c2785065 781 }
782 }// background subtraction
783
784
3b7ffecf 785 Double_t eventW = 1;
786 Double_t ptHard = 0;
787 Double_t nTrials = 1; // Trials for MC trigger
788
789 if(fUseExternalWeightOnly){
790 eventW = fExternalWeight;
791 }
792
793 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
cf049e94 794 // if(fDebug>0)aodH->SetFillAOD(kFALSE);
3b7ffecf 795 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
796 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
797 // this is the part we only use when we have MC information
798 AliMCEvent* mcEvent = MCEvent();
799 // AliStack *pStack = 0;
800 if(!mcEvent){
801 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
802 return;
803 }
804 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
805 Int_t iCount = 0;
806 if(pythiaGenHeader){
807 nTrials = pythiaGenHeader->Trials();
808 ptHard = pythiaGenHeader->GetPtHard();
809 int iProcessType = pythiaGenHeader->ProcessType();
810 // 11 f+f -> f+f
811 // 12 f+barf -> f+barf
812 // 13 f+barf -> g+g
813 // 28 f+g -> f+g
814 // 53 g+g -> f+barf
815 // 68 g+g -> g+g
816 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
817 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
818
819 // fetch the pythia generated jets only to be used here
820 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
821 AliAODJet pythiaGenJets[kMaxJets];
822 for(int ip = 0;ip < nPythiaGenJets;++ip){
823 if(iCount>=kMaxJets)continue;
824 Float_t p[4];
825 pythiaGenHeader->TriggerJet(ip,p);
826 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
827
3b7ffecf 828 if(fBranchGen.Length()==0){
9280dfa6 829 /*
830 if(fLimitGenJetEta){
831 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
832 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
833 }
834 */
835 if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
3b7ffecf 836 // if we have MC particles and we do not read from the aod branch
837 // use the pythia jets
f4132e7d 838 if(!fBkgSubtraction){
839 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
840 iCount++;
841 }
3b7ffecf 842 }
3b7ffecf 843 }
844 }
f4132e7d 845 if(fBranchGen.Length()==0&&!fBkgSubtraction)nGenJets = iCount;
3b7ffecf 846 }// (fAnalysisType&kMCESD)==kMCESD)
847
848
849 // we simply fetch the tracks/mc particles as a list of AliVParticles
850
851 TList recParticles;
852 TList genParticles;
853
565584e8 854
855
856
857 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
858 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
859 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
860 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
3b7ffecf 861
cc0649e4 862
3b7ffecf 863 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
864 fh1PtHard->Fill(ptHard,eventW);
865 fh1PtHardNoW->Fill(ptHard,1);
866 fh1PtHardTrials->Fill(ptHard,nTrials);
57ca1193 867 fh1ZVtx->Fill(fAOD->GetPrimaryVertex()->GetZ());
868
3b7ffecf 869
870 // If we set a second branch for the input jets fetch this
f4132e7d 871 if(fBranchGen.Length()>0&&!fBkgSubtraction){
3b7ffecf 872 TClonesArray *aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
873 if(aodGenJets){
874 Int_t iCount = 0;
875 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
876 if(iCount>=kMaxJets)continue;
877 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
878 if(!tmp)continue;
5301f754 879 /*
3b7ffecf 880 if(fLimitGenJetEta){
881 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
882 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
883 }
5301f754 884 */
9280dfa6 885 if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
3b7ffecf 886 genJets[iCount] = *tmp;
887 iCount++;
888 }
889 nGenJets = iCount;
890 }
891 else{
5cb0f438 892 if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
893 if(fDebug>2)fAOD->Print();
3b7ffecf 894 }
895 }
896
897 fh1NGenJets->Fill(nGenJets);
898 // We do not want to exceed the maximum jet number
899 nGenJets = TMath::Min(nGenJets,kMaxJets);
900
901 // Fetch the reconstructed jets...
902
903 nRecJets = aodRecJets->GetEntries();
904
cc0649e4 905 nRecJets = aodRecJets->GetEntries();
906 fh1NRecJets->Fill(nRecJets);
907
908 // Do something with the all rec jets
909 Int_t nRecOver = nRecJets;
3b7ffecf 910
cc0649e4 911 // check that the jets are sorted
912 Float_t ptOld = 999999;
913 for(int ir = 0;ir < nRecJets;ir++){
914 AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
915 Float_t tmpPt = tmp->Pt();
916 if(tmpPt>ptOld){
3568c3d6 917 Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
cc0649e4 918 }
919 ptOld = tmpPt;
920 }
3b7ffecf 921
cc0649e4 922
923 if(nRecOver>0){
924 TIterator *recIter = aodRecJets->MakeIterator();
925 AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
926 Float_t pt = tmpRec->Pt();
927 if(tmpRec){
928 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
929 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
930 while(pt<ptCut&&tmpRec){
931 nRecOver--;
932 tmpRec = (AliAODJet*)(recIter->Next());
933 if(tmpRec){
934 pt = tmpRec->Pt();
935 }
936 }
937 if(nRecOver<=0)break;
938 fh2NRecJetsPt->Fill(ptCut,nRecOver);
939 }
940 }
941 recIter->Reset();
942 AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
943 Float_t phi = leading->Phi();
944 if(phi<0)phi+=TMath::Pi()*2.;
945 Float_t eta = leading->Eta();
cf049e94 946 pt = leading->Pt();
45a11aef 947 while((tmpRec = (AliAODJet*)(recIter->Next()))){
cc0649e4 948 Float_t tmpPt = tmpRec->Pt();
949 fh1PtJetsRecIn->Fill(tmpPt);
565584e8 950 if(tmpRec==leading){
951 fh1PtJetsLeadingRecIn->Fill(tmpPt);
952 continue;
953 }
cc0649e4 954 // correlation
955 Float_t tmpPhi = tmpRec->Phi();
956
957 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
958 Float_t dPhi = phi - tmpRec->Phi();
959 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
960 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
961 Float_t dEta = eta - tmpRec->Eta();
962 fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
565584e8 963 fh2JetsLeadingPhiPt->Fill(dPhi,pt);
cc0649e4 964 }
965 delete recIter;
966 }
967
968 Int_t nTrackOver = recParticles.GetSize();
969 // do the same for tracks and jets
970 if(nTrackOver>0){
971 TIterator *recIter = recParticles.MakeIterator();
972 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
973 Float_t pt = tmpRec->Pt();
974 // Printf("Leading track p_t %3.3E",pt);
975 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
976 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
977 while(pt<ptCut&&tmpRec){
978 nTrackOver--;
979 tmpRec = (AliVParticle*)(recIter->Next());
980 if(tmpRec){
981 pt = tmpRec->Pt();
982 }
983 }
984 if(nTrackOver<=0)break;
985 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
986 }
987
988 recIter->Reset();
989 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
990 Float_t phi = leading->Phi();
991 if(phi<0)phi+=TMath::Pi()*2.;
992 Float_t eta = leading->Eta();
cf049e94 993 pt = leading->Pt();
cc0649e4 994 while((tmpRec = (AliVParticle*)(recIter->Next()))){
995 Float_t tmpPt = tmpRec->Pt();
996 fh1PtTracksRecIn->Fill(tmpPt);
565584e8 997 if(tmpRec==leading){
998 fh1PtTracksLeadingRecIn->Fill(tmpPt);
999 continue;
1000 }
cc0649e4 1001 // correlation
1002 Float_t tmpPhi = tmpRec->Phi();
1003
1004 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1005 Float_t dPhi = phi - tmpRec->Phi();
1006 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1007 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1008 Float_t dEta = eta - tmpRec->Eta();
1009 fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
565584e8 1010 fh2TracksLeadingPhiPt->Fill(dPhi,pt);
cc0649e4 1011 }
1012 delete recIter;
1013 }
1014
1015 if(genParticles.GetSize()){
1016 TIterator *genIter = genParticles.MakeIterator();
1017 AliVParticle *tmpGen = 0;
1018 while((tmpGen = (AliVParticle*)(genIter->Next()))){
1019 if(TMath::Abs(tmpGen->Eta())<0.9){
1020 Float_t tmpPt = tmpGen->Pt();
1021 fh1PtTracksGenIn->Fill(tmpPt);
1022 }
1023 }
1024 delete genIter;
1025 }
1026
3b7ffecf 1027 nRecJets = TMath::Min(nRecJets,kMaxJets);
6bd3fdae 1028 nJets = TMath::Min(nJets,kMaxJets);
9280dfa6 1029 Int_t iCountRec = 0;
3b7ffecf 1030 for(int ir = 0;ir < nRecJets;++ir){
1031 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
1032 if(!tmp)continue;
9280dfa6 1033 if(tmp->Pt()<fMinJetPt)continue;
3b7ffecf 1034 recJets[ir] = *tmp;
9280dfa6 1035 iCountRec++;
3b7ffecf 1036 }
9280dfa6 1037 nRecJets = iCountRec;
3b7ffecf 1038
1039
1040 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1041 // Relate the jets
1042 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
1043 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
1044
cc0649e4 1045
1046 for(int i = 0;i<kMaxJets;++i){
3b7ffecf 1047 iGenIndex[i] = iRecIndex[i] = -1;
1048 }
1049
6bd3fdae 1050
3b7ffecf 1051 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
1052 iGenIndex,iRecIndex,fDebug);
6bd3fdae 1053
1054
3b7ffecf 1055 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1056
1057 if(fDebug){
1058 for(int i = 0;i<kMaxJets;++i){
1059 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
1060 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
1061 }
1062 }
1063
1064
1065
1066
1067 Double_t container[6];
c8eabe24 1068 Double_t containerPhiZ[6];
3b7ffecf 1069
1070 // loop over generated jets
1071
1072 // radius; tmp, get from the jet header itself
1073 // at some point, how todeal woht FastJet on MC?
1074 Float_t radiusGen =0.4;
1075 Float_t radiusRec =0.4;
1076
1077 for(int ig = 0;ig < nGenJets;++ig){
1078 Double_t ptGen = genJets[ig].Pt();
1079 Double_t phiGen = genJets[ig].Phi();
1080 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1081 Double_t etaGen = genJets[ig].Eta();
1082
1083 container[3] = ptGen;
1084 container[4] = etaGen;
1085 container[5] = phiGen;
1086 fhnJetContainer[kStep0]->Fill(&container[3],eventW);
1087 Int_t ir = iRecIndex[ig];
1088
1089 if(TMath::Abs(etaGen)<fRecEtaWindow){
c8eabe24 1090 fh1TmpRho->Reset();
1091
3b7ffecf 1092 fhnJetContainer[kStep1]->Fill(&container[3],eventW);
1093 fh1PtGenIn[ig]->Fill(ptGen,eventW);
1094 // fill the fragmentation function
1095 for(int it = 0;it<genParticles.GetEntries();++it){
1096 AliVParticle *part = (AliVParticle*)genParticles.At(it);
c8eabe24 1097 Float_t deltaR = genJets[ig].DeltaR(part);
80a790fd 1098 if(ptGen>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
1099 if(deltaR<radiusGen&&ptGen>0){
3b7ffecf 1100 Float_t z = part->Pt()/ptGen;
1101 Float_t lnz = -1.*TMath::Log(z);
1102 fh2FragGen[ig]->Fill(z,ptGen,eventW);
1103 fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
1104 }
c8eabe24 1105
3b7ffecf 1106 }
c8eabe24 1107 Float_t rhoSum = 0;
fb03edbe 1108 for(int ibx = 1;ibx<=fh2RhoPtGen[ir]->GetNbinsX();ibx++){
c8eabe24 1109 Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
1110 Float_t rho = fh1TmpRho->GetBinContent(ibx);
1111 rhoSum += rho;
1112 fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
1113 fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
3b7ffecf 1114 }
c8eabe24 1115 }
3b7ffecf 1116 if(ir>=0&&ir<nRecJets){
1117 fhnJetContainer[kStep2]->Fill(&container[3],eventW);
1118 Double_t etaRec = recJets[ir].Eta();
1119 if(TMath::Abs(etaRec)<fRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
2d1d1b60 1120 if(TMath::Abs(etaRec)<fRecEtaWindow&&TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep3]->Fill(&container[3],eventW);
3b7ffecf 1121 }
1122 }// loop over generated jets
7fa8b2da 1123
edfbe476 1124 Float_t sumPt = 0;
1125 for(int it = 0;it<recParticles.GetEntries();++it){
1126 AliVParticle *part = (AliVParticle*)recParticles.At(it);
1127 // fill sum pt and P_t of all paritles
1128 if(TMath::Abs(part->Eta())<0.9){
1129 Float_t pt = part->Pt();
1130 fh1PtTrackRec->Fill(pt,eventW);
c8eabe24 1131 fh2TrackPtTrackPhi->Fill(pt,part->Phi());
edfbe476 1132 sumPt += pt;
1133 }
1134 }
1135 fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
1136 fh1SumPtTrackRec->Fill(sumPt,eventW);
1137
cf049e94 1138
3b7ffecf 1139 // loop over reconstructed jets
1140 for(int ir = 0;ir < nRecJets;++ir){
1141 Double_t etaRec = recJets[ir].Eta();
1142 Double_t ptRec = recJets[ir].Pt();
1143 Double_t phiRec = recJets[ir].Phi();
1144 if(phiRec<0)phiRec+=TMath::Pi()*2.;
26fa06fb 1145
1146
1147 // do something with dijets...
1148 if(ir==1){
1149 Double_t etaRec1 = recJets[0].Eta();
1150 Double_t ptRec1 = recJets[0].Pt();
1151 Double_t phiRec1 = recJets[0].Phi();
1152 if(phiRec1<0)phiRec1+=TMath::Pi()*2.;
1153
1154
1155 if(TMath::Abs(etaRec1)<fRecEtaWindow
1156 &&TMath::Abs(etaRec)<fRecEtaWindow){
1157
1158 Float_t deltaPhi = phiRec1 - phiRec;
1159
1160 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1161 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1162 deltaPhi = TMath::Abs(deltaPhi);
1163 fh2DijetDeltaPhiPt->Fill(deltaPhi,ptRec1);
80a790fd 1164 Float_t asym = 9999;
1165 if((ptRec1+ptRec)>0)asym = (ptRec1-ptRec)/(ptRec1+ptRec);
26fa06fb 1166 fh2DijetAsymPt->Fill(asym,ptRec1);
1167 fh2DijetDeltaPhiDeltaEta->Fill(deltaPhi,etaRec1-etaRec);
1168 fh2DijetPt2vsPt1->Fill(ptRec1,ptRec);
1169 fh2DijetDifvsSum->Fill(ptRec1+ptRec,ptRec1-ptRec);
1170 Float_t minv = 2.*(recJets[0].P()*recJets[1].P()-
1171 recJets[0].Px()*recJets[1].Px()-
1172 recJets[0].Py()*recJets[1].Py()-
6787fd34 1173 recJets[0].Pz()*recJets[1].Pz());
1174 if(minv<0)minv=0; // prevent numerical instabilities
26fa06fb 1175 minv = TMath::Sqrt(minv);
1176 // with mass == 0;
1177
1178 fh1DijetMinv->Fill(minv);
1179 if((TMath::Pi()-deltaPhi)<fDeltaPhiWindow){
1180 fh1DijetMinvCut->Fill(minv);
1181 fh2DijetAsymPtCut->Fill(asym,ptRec1);
1182 }
1183 }
1184 }
1185
1186
3b7ffecf 1187 container[0] = ptRec;
1188 container[1] = etaRec;
1189 container[2] = phiRec;
c8eabe24 1190 containerPhiZ[0] = ptRec;
1191 containerPhiZ[1] = phiRec;
1277f815 1192 if(ptRec>30.&&fDebug>0){
7fa8b2da 1193 // need to cast to int, otherwise the printf overwrites
1194 Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
1277f815 1195 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
1196 if(fESD)Printf("ESDEvent GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
cf049e94 1197 // aodH->SetFillAOD(kTRUE);
7fa8b2da 1198 fAOD->GetHeader()->Print();
a221560c 1199 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
7fa8b2da 1200 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1201 AliAODTrack *tr = fAOD->GetTrack(it);
1202 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1203 tr->Print();
38ecb6a5 1204 // tr->Dump();
7fa8b2da 1205 if(fESD){
1206 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1207 esdTr->Print("");
bb3d13aa 1208 // esdTr->Dump();
7fa8b2da 1209 }
1210 }
1211 }
1212
1213
3b7ffecf 1214 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
1215 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
45a11aef 1216
c8eabe24 1217 Float_t zLeading = -1;
3b7ffecf 1218 if(TMath::Abs(etaRec)<fRecEtaWindow){
c8eabe24 1219 fh2JetPtJetPhi->Fill(ptRec,phiRec);
3b7ffecf 1220 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
1221 fh1PtRecIn[ir]->Fill(ptRec,eventW);
1222 // fill the fragmentation function
c8eabe24 1223
1224 fh1TmpRho->Reset();
1225
3b7ffecf 1226 for(int it = 0;it<recParticles.GetEntries();++it){
1227 AliVParticle *part = (AliVParticle*)recParticles.At(it);
edfbe476 1228 Float_t eta = part->Eta();
1229 if(TMath::Abs(eta)<0.9){
1230 Float_t phi = part->Phi();
1231 if(phi<0)phi+=TMath::Pi()*2.;
1232 Float_t dPhi = phi - phiRec;
1233 Float_t dEta = eta - etaRec;
26fa06fb 1234 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1235 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
edfbe476 1236 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1237 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1238 }
c8eabe24 1239
1240 Float_t deltaR = recJets[ir].DeltaR(part);
80a790fd 1241 if(ptRec>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
c8eabe24 1242
1243
80a790fd 1244 if(deltaR<radiusRec&&ptRec>0){
3b7ffecf 1245 Float_t z = part->Pt()/ptRec;
c8eabe24 1246 if(z>zLeading)zLeading=z;
3b7ffecf 1247 Float_t lnz = -1.*TMath::Log(z);
1248 fh2FragRec[ir]->Fill(z,ptRec,eventW);
1249 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1250 }
1251 }
c8eabe24 1252 // fill the jet shapes
1253 Float_t rhoSum = 0;
fb03edbe 1254 for(int ibx = 1;ibx<=fh2RhoPtRec[ir]->GetNbinsX();ibx++){
c8eabe24 1255 Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1256 Float_t rho = fh1TmpRho->GetBinContent(ibx);
1257 rhoSum += rho;
1258 fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1259 fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1260 }
3b7ffecf 1261 }
6bd3fdae 1262
3b7ffecf 1263 // Fill Correlation
1264 Int_t ig = iGenIndex[ir];
1265 if(ig>=0 && ig<nGenJets){
1266 fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1267 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1268 Double_t ptGen = genJets[ig].Pt();
1269 Double_t phiGen = genJets[ig].Phi();
1270 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1271 Double_t etaGen = genJets[ig].Eta();
6bd3fdae 1272
3b7ffecf 1273 container[3] = ptGen;
1274 container[4] = etaGen;
1275 container[5] = phiGen;
c8eabe24 1276 containerPhiZ[3] = ptGen;
3b7ffecf 1277 //
1278 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1279 //
6bd3fdae 1280
3b7ffecf 1281 if(TMath::Abs(etaGen)<fRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1282 if(TMath::Abs(etaRec)<fRecEtaWindow){
1283 fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1284 fhnCorrelation->Fill(container);
9a83d4af 1285 if(ptGen>0){
1286 Float_t delta = (ptRec-ptGen)/ptGen;
1287 fh2RelPtFGen->Fill(ptGen,delta,eventW);
1288 }
c8eabe24 1289 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
6bd3fdae 1290 }// if etarec in window
3b7ffecf 1291 }
3b7ffecf 1292 else{
c8eabe24 1293 containerPhiZ[3] = 0;
1294 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
3b7ffecf 1295 }
1296 }// loop over reconstructed jets
1297
1298
1299 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1300 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1301 PostData(1, fHistList);
1302}
1303
1304
1305void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1306 //
1307 // Create the particle container for the correction framework manager and
1308 // link it
1309 //
1310 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
cba109a0 1311 const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
3b7ffecf 1312 const Double_t kEtamin = -3.0, kEtamax = 3.0;
1313 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
c8eabe24 1314 const Double_t kZmin = 0., kZmax = 1;
3b7ffecf 1315
1316 // can we neglect migration in eta and phi?
1317 // phi should be no problem since we cover full phi and are phi symmetric
1318 // eta migration is more difficult due to needed acceptance correction
1319 // in limited eta range
1320
3b7ffecf 1321 //arrays for the number of bins in each dimension
1322 Int_t iBin[kNvar];
cba109a0 1323 iBin[0] = 320; //bins in pt
3b7ffecf 1324 iBin[1] = 1; //bins in eta
1325 iBin[2] = 1; // bins in phi
1326
1327 //arrays for lower bounds :
1328 Double_t* binEdges[kNvar];
1329 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1330 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1331
1332 //values for bin lower bounds
1333 // 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 1334 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 1335 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1336 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1337
1338
1339 for(int i = 0;i < kMaxStep*2;++i){
f51451be 1340 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
3b7ffecf 1341 for (int k=0; k<kNvar; k++) {
1342 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1343 }
1344 }
1345 //create correlation matrix for unfolding
1346 Int_t thnDim[2*kNvar];
1347 for (int k=0; k<kNvar; k++) {
1348 //first half : reconstructed
1349 //second half : MC
1350 thnDim[k] = iBin[k];
1351 thnDim[k+kNvar] = iBin[k];
1352 }
1353
1354 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1355 for (int k=0; k<kNvar; k++) {
1356 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1357 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1358 }
1359 fhnCorrelation->Sumw2();
1360
c8eabe24 1361 // for second correlation histogram
1362
1363
1364 const Int_t kNvarPhiZ = 4;
1365 //arrays for the number of bins in each dimension
1366 Int_t iBinPhiZ[kNvarPhiZ];
1367 iBinPhiZ[0] = 80; //bins in pt
1368 iBinPhiZ[1] = 72; //bins in phi
1369 iBinPhiZ[2] = 20; // bins in Z
1370 iBinPhiZ[3] = 80; //bins in ptgen
1371
1372 //arrays for lower bounds :
1373 Double_t* binEdgesPhiZ[kNvarPhiZ];
1374 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1375 binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1376
1377 for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1378 for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
1379 for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
1380 for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
1381
1382 fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
1383 for (int k=0; k<kNvarPhiZ; k++) {
1384 fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
1385 }
1386 fhnCorrelationPhiZRec->Sumw2();
1387
1388
3b7ffecf 1389 // Add a histogram for Fake jets
c8eabe24 1390
5cb0f438 1391 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1392 delete [] binEdges[ivar];
1393
c8eabe24 1394 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1395 delete [] binEdgesPhiZ[ivar];
1396
3b7ffecf 1397}
1398
1399void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1400{
1401// Terminate analysis
1402//
1403 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1404}
1405
1406
1407Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1408
565584e8 1409 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
3b7ffecf 1410
1411 Int_t iCount = 0;
565584e8 1412 if(type==kTrackAOD){
3b7ffecf 1413 AliAODEvent *aod = 0;
565584e8 1414 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1415 else aod = AODEvent();
3b7ffecf 1416 if(!aod){
1417 return iCount;
1418 }
1419 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1420 AliAODTrack *tr = aod->GetTrack(it);
1421 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
cc0649e4 1422 if(TMath::Abs(tr->Eta())>0.9)continue;
38ecb6a5 1423 if(fDebug>0){
1424 if(tr->Pt()>20){
1425 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
fd5db301 1426 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
38ecb6a5 1427 tr->Print();
1428 // tr->Dump();
1429 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1430 if(fESD){
1431 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1432 esdTr->Print("");
bb3d13aa 1433 // esdTr->Dump();
38ecb6a5 1434 }
1435 }
1436 }
3b7ffecf 1437 list->Add(tr);
1438 iCount++;
1439 }
1440 }
1441 else if (type == kTrackKineAll||type == kTrackKineCharged){
1442 AliMCEvent* mcEvent = MCEvent();
1443 if(!mcEvent)return iCount;
1444 // we want to have alivpartilces so use get track
1445 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1446 if(!mcEvent->IsPhysicalPrimary(it))continue;
1447 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1448 if(type == kTrackKineAll){
1449 list->Add(part);
1450 iCount++;
1451 }
1452 else if(type == kTrackKineCharged){
1453 if(part->Particle()->GetPDG()->Charge()==0)continue;
1454 list->Add(part);
1455 iCount++;
1456 }
1457 }
1458 }
565584e8 1459 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1460 AliAODEvent *aod = 0;
5301f754 1461 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
565584e8 1462 else aod = AODEvent();
5010a3f7 1463 if(!aod)return iCount;
3b7ffecf 1464 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1465 if(!tca)return iCount;
1466 for(int it = 0;it < tca->GetEntriesFast();++it){
1467 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1468 if(!part->IsPhysicalPrimary())continue;
c4631cdb 1469 if(type == kTrackAODMCAll){
3b7ffecf 1470 list->Add(part);
1471 iCount++;
1472 }
565584e8 1473 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
3b7ffecf 1474 if(part->Charge()==0)continue;
565584e8 1475 if(kTrackAODMCCharged){
1476 list->Add(part);
1477 }
1478 else {
1479 if(TMath::Abs(part->Eta())>0.9)continue;
1480 list->Add(part);
1481 }
3b7ffecf 1482 iCount++;
1483 }
1484 }
1485 }// AODMCparticle
cc0649e4 1486 list->Sort();
c4631cdb 1487 return iCount;
3b7ffecf 1488
1489}