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