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