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