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