]> 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
Possibility to use git-svn in the creation of ARVersion.h file (Svein)
[u/mrichter/AliRoot.git] / PWG4 / JetTasks / AliAnalysisTaskJetSpectrum2.cxx
CommitLineData
a31b8a87 1
3b7ffecf 2// **************************************
3// Task used for the correction of determiantion of reconstructed jet spectra
4// Compares input (gen) and output (rec) jets
5// *******************************************
6
7
8/**************************************************************************
9 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
10 * *
11 * Author: The ALICE Off-line Project. *
12 * Contributors are mentioned in the code where appropriate. *
13 * *
14 * Permission to use, copy, modify and distribute this software and its *
15 * documentation strictly for non-commercial purposes is hereby granted *
16 * without fee, provided that the above copyright notice appears in all *
17 * copies and that both the copyright notice and this permission notice *
18 * appear in the supporting documentation. The authors make no claims *
19 * about the suitability of this software for any purpose. It is *
20 * provided "as is" without express or implied warranty. *
21 **************************************************************************/
22
23
24#include <TROOT.h>
25#include <TRandom.h>
26#include <TSystem.h>
27#include <TInterpreter.h>
28#include <TChain.h>
29#include <TFile.h>
30#include <TKey.h>
31#include <TH1F.h>
32#include <TH2F.h>
33#include <TH3F.h>
34#include <TProfile.h>
35#include <TList.h>
36#include <TLorentzVector.h>
37#include <TClonesArray.h>
38#include "TDatabasePDG.h"
39
40#include "AliAnalysisTaskJetSpectrum2.h"
41#include "AliAnalysisManager.h"
42#include "AliJetFinder.h"
43#include "AliJetHeader.h"
44#include "AliJetReader.h"
45#include "AliJetReaderHeader.h"
46#include "AliUA1JetHeaderV1.h"
47#include "AliESDEvent.h"
48#include "AliAODEvent.h"
49#include "AliAODHandler.h"
50#include "AliAODTrack.h"
51#include "AliAODJet.h"
c2785065 52#include "AliAODJetEventBackground.h"
3b7ffecf 53#include "AliAODMCParticle.h"
54#include "AliMCEventHandler.h"
55#include "AliMCEvent.h"
56#include "AliStack.h"
57#include "AliGenPythiaEventHeader.h"
58#include "AliJetKineReaderHeader.h"
59#include "AliGenCocktailEventHeader.h"
60#include "AliInputEventHandler.h"
61
62
63#include "AliAnalysisHelperJetTasks.h"
64
65ClassImp(AliAnalysisTaskJetSpectrum2)
66
67AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(): AliAnalysisTaskSE(),
a31b8a87 68 fJetHeaderRec(0x0),
69 fJetHeaderGen(0x0),
70 fAODIn(0x0),
71 fAODOut(0x0),
72 fAODExtension(0x0),
73 fhnCorrelation(0x0),
c8eabe24 74 fhnCorrelationPhiZRec(0x0),
75 f1PtScale(0x0),
a31b8a87 76 fBranchRec("jets"),
77 fBranchGen(""),
78 fBranchBkg(""),
79 fNonStdFile(""),
80 fUseAODJetInput(kFALSE),
81 fUseAODTrackInput(kFALSE),
82 fUseAODMCInput(kFALSE),
83 fUseGlobalSelection(kFALSE),
84 fUseExternalWeightOnly(kFALSE),
85 fLimitGenJetEta(kFALSE),
86 fBkgSubtraction(kFALSE),
87 fUseOldFill(kFALSE),
88 fNMatchJets(5),
89 fFillCorrBkg(0),
90 fFilterMask(0),
f2dd0695 91 fEventSelectionMask(0),
a31b8a87 92 fAnalysisType(0),
3b7ffecf 93 fTrackTypeRec(kTrackUndef),
94 fTrackTypeGen(kTrackUndef),
f4132e7d 95 fEventClass(0),
3b7ffecf 96 fAvgTrials(1),
97 fExternalWeight(1),
a31b8a87 98 fJetRecEtaWindow(0.5),
99 fTrackRecEtaWindow(0.5),
9280dfa6 100 fMinJetPt(0),
a31b8a87 101 fMinTrackPt(0.15),
102 fDeltaPhiWindow(90./180.*TMath::Pi()),
3b7ffecf 103 fh1Xsec(0x0),
104 fh1Trials(0x0),
105 fh1PtHard(0x0),
106 fh1PtHardNoW(0x0),
107 fh1PtHardTrials(0x0),
57ca1193 108 fh1ZVtx(0x0),
3b7ffecf 109 fh1NGenJets(0x0),
110 fh1NRecJets(0x0),
edfbe476 111 fh1PtTrackRec(0x0),
112 fh1SumPtTrackRec(0x0),
113 fh1SumPtTrackAreaRec(0x0),
c8eabe24 114 fh1TmpRho(0x0),
cc0649e4 115 fh1PtJetsRecIn(0x0),
d8f21f85 116 fh1PtJetsGenIn(0x0),
565584e8 117 fh1PtJetsLeadingRecIn(0x0),
cc0649e4 118 fh1PtTracksRecIn(0x0),
565584e8 119 fh1PtTracksLeadingRecIn(0x0),
cc0649e4 120 fh1PtTracksGenIn(0x0),
121 fh2NRecJetsPt(0x0),
122 fh2NRecTracksPt(0x0),
a31b8a87 123 fh2NGenJetsPt(0x0),
124 fh2NGenTracksPt(0x0),
125 fh2PtFGen(0x0),
9a83d4af 126 fh2RelPtFGen(0x0),
185fa8e6 127 fh1Bkg1(0x0),
128 fh1Bkg2(0x0),
7f9012c1 129 fh1Bkg3(0x0),
c2785065 130 fh1Sigma1(0x0),
131 fh1Sigma2(0x0),
7f9012c1 132 fh1Sigma3(0x0),
185fa8e6 133 fh1Area1(0x0),
134 fh1Area2(0x0),
7f9012c1 135 fh1Area3(0x0),
c2785065 136 fh1Ptjet(0x0),
137 fh1Ptjetsub1(0x0),
138 fh1Ptjetsub2(0x0),
7f9012c1 139 fh1Ptjetsub3(0x0),
c2785065 140 fh1Ptjethardest(0x0),
141 fh1Ptjetsubhardest1(0x0),
142 fh1Ptjetsubhardest2(0x0),
7f9012c1 143 fh1Ptjetsubhardest3(0x0),
6bd3fdae 144 fh2Rhovspthardest1(0x0),
145 fh2Rhovspthardest2(0x0),
146 fh2Rhovspthardest3(0x0),
147 fh2Errorvspthardest1(0x0),
148 fh2Errorvspthardest2(0x0),
149 fh2Errorvspthardest3(0x0),
3b7ffecf 150 fHistList(0x0)
151{
152 for(int i = 0;i < kMaxStep*2;++i){
153 fhnJetContainer[i] = 0;
154 }
155 for(int i = 0;i < kMaxJets;++i){
3b7ffecf 156 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
d931e0ef 157
158 fh2PhiPt[i] = 0;
159 fh2PhiEta[i] = 0;
160 fh2RhoPtRec[i] = 0;
161 fh2RhoPtGen[i] = 0;
162 fh2PsiPtGen[i] = 0;
163 fh2PsiPtRec[i] = 0;
164 fh2FragRec[i] = 0;
165 fh2FragLnRec[i] = 0;
166 fh2FragGen[i] = 0;
167 fh2FragLnGen[i] = 0;
3b7ffecf 168 }
169
a31b8a87 170
171 for(int ij = 0;ij <kJetTypes;++ij){
172 fh1NJets[ij] = 0;
173 fh1SumPtTrack[ij] = 0;
174 fh1PtJetsIn[ij] = 0;
175 fh1PtTracksIn[ij] = 0;
176 fh1PtTracksLeadingIn[ij] = 0;
177 fh2NJetsPt[ij] = 0;
178 fh2NTracksPt[ij] = 0;
179 fh2LeadingJetPtJetPhi[ij] = 0;
180 fh2LeadingTrackPtTrackPhi[ij] = 0;
181 for(int i = 0;i < kMaxJets;++i){
182 fh1PtIn[ij][i] = 0;
183 fh2RhoPt[ij][i] = 0;
184 fh2PsiPt[ij][i] = 0;
185 }
186
187 fh1DijetMinv[ij] = 0;
188 fh2DijetDeltaPhiPt[ij] = 0;
189 fh2DijetAsymPt[ij] = 0;
190 fh2DijetPt2vsPt1[ij] = 0;
191 fh2DijetDifvsSum[ij] = 0;
192
193 }
3b7ffecf 194}
195
196AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
197 AliAnalysisTaskSE(name),
198 fJetHeaderRec(0x0),
199 fJetHeaderGen(0x0),
a31b8a87 200 fAODIn(0x0),
201 fAODOut(0x0),
202 fAODExtension(0x0),
3b7ffecf 203 fhnCorrelation(0x0),
c8eabe24 204 fhnCorrelationPhiZRec(0x0),
205 f1PtScale(0x0),
3b7ffecf 206 fBranchRec("jets"),
207 fBranchGen(""),
c2785065 208 fBranchBkg(""),
a31b8a87 209 fNonStdFile(""),
565584e8 210 fUseAODJetInput(kFALSE),
211 fUseAODTrackInput(kFALSE),
212 fUseAODMCInput(kFALSE),
b5cc0c6d 213 fUseGlobalSelection(kFALSE),
3b7ffecf 214 fUseExternalWeightOnly(kFALSE),
215 fLimitGenJetEta(kFALSE),
c2785065 216 fBkgSubtraction(kFALSE),
a31b8a87 217 fUseOldFill(kFALSE),
218 fNMatchJets(5),
6bd3fdae 219 fFillCorrBkg(0),
3b7ffecf 220 fFilterMask(0),
f2dd0695 221 fEventSelectionMask(0),
3b7ffecf 222 fAnalysisType(0),
223 fTrackTypeRec(kTrackUndef),
224 fTrackTypeGen(kTrackUndef),
f4132e7d 225 fEventClass(0),
3b7ffecf 226 fAvgTrials(1),
227 fExternalWeight(1),
a31b8a87 228 fJetRecEtaWindow(0.5),
229 fTrackRecEtaWindow(0.5),
9280dfa6 230 fMinJetPt(0),
a31b8a87 231 fMinTrackPt(0.15),
232 fDeltaPhiWindow(90./180.*TMath::Pi()),
3b7ffecf 233 fh1Xsec(0x0),
234 fh1Trials(0x0),
235 fh1PtHard(0x0),
236 fh1PtHardNoW(0x0),
237 fh1PtHardTrials(0x0),
57ca1193 238 fh1ZVtx(0x0),
3b7ffecf 239 fh1NGenJets(0x0),
240 fh1NRecJets(0x0),
edfbe476 241 fh1PtTrackRec(0x0),
242 fh1SumPtTrackRec(0x0),
243 fh1SumPtTrackAreaRec(0x0),
c8eabe24 244 fh1TmpRho(0x0),
cc0649e4 245 fh1PtJetsRecIn(0x0),
d8f21f85 246 fh1PtJetsGenIn(0x0),
565584e8 247 fh1PtJetsLeadingRecIn(0x0),
cc0649e4 248 fh1PtTracksRecIn(0x0),
565584e8 249 fh1PtTracksLeadingRecIn(0x0),
cc0649e4 250 fh1PtTracksGenIn(0x0),
251 fh2NRecJetsPt(0x0),
252 fh2NRecTracksPt(0x0),
a31b8a87 253 fh2NGenJetsPt(0x0),
254 fh2NGenTracksPt(0x0),
255 fh2PtFGen(0x0),
9a83d4af 256 fh2RelPtFGen(0x0),
7f9012c1 257 fh1Bkg1(0x0),
185fa8e6 258 fh1Bkg2(0x0),
7f9012c1 259 fh1Bkg3(0x0),
c2785065 260 fh1Sigma1(0x0),
261 fh1Sigma2(0x0),
7f9012c1 262 fh1Sigma3(0x0),
185fa8e6 263 fh1Area1(0x0),
264 fh1Area2(0x0),
7f9012c1 265 fh1Area3(0x0),
c2785065 266 fh1Ptjet(0x0),
7f9012c1 267 fh1Ptjetsub1(0x0),
268 fh1Ptjetsub2(0x0),
269 fh1Ptjetsub3(0x0),
c2785065 270 fh1Ptjethardest(0x0),
271 fh1Ptjetsubhardest1(0x0),
272 fh1Ptjetsubhardest2(0x0),
7f9012c1 273 fh1Ptjetsubhardest3(0x0),
6bd3fdae 274 fh2Rhovspthardest1(0x0),
275 fh2Rhovspthardest2(0x0),
276 fh2Rhovspthardest3(0x0),
277 fh2Errorvspthardest1(0x0),
278 fh2Errorvspthardest2(0x0),
279 fh2Errorvspthardest3(0x0),
3b7ffecf 280 fHistList(0x0)
281{
282
283 for(int i = 0;i < kMaxStep*2;++i){
284 fhnJetContainer[i] = 0;
285 }
286 for(int i = 0;i < kMaxJets;++i){
3b7ffecf 287 fh1PtRecIn[i] = fh1PtGenIn[i] = 0;
d931e0ef 288
289 fh2PhiPt[i] = 0;
290 fh2PhiEta[i] = 0;
291 fh2RhoPtRec[i] = 0;
292 fh2RhoPtGen[i] = 0;
293 fh2PsiPtGen[i] = 0;
294 fh2PsiPtRec[i] = 0;
295 fh2FragRec[i] = 0;
296 fh2FragLnRec[i] = 0;
297 fh2FragGen[i] = 0;
298 fh2FragLnGen[i] = 0;
3b7ffecf 299 }
a31b8a87 300
301 for(int ij = 0;ij <kJetTypes;++ij){
302 fh1NJets[ij] = 0;
303 fh1SumPtTrack[ij] = 0;
304 fh1PtJetsIn[ij] = 0;
305 fh1PtTracksIn[ij] = 0;
306 fh1PtTracksLeadingIn[ij] = 0;
307 fh2NJetsPt[ij] = 0;
308 fh2NTracksPt[ij] = 0;
309 fh2LeadingJetPtJetPhi[ij] = 0;
310 fh2LeadingTrackPtTrackPhi[ij] = 0;
311 for(int i = 0;i < kMaxJets;++i){
312 fh1PtIn[ij][i] = 0;
313 fh2RhoPt[ij][i] = 0;
314 fh2PsiPt[ij][i] = 0;
315 }
316
317 fh1DijetMinv[ij] = 0;
318 fh2DijetDeltaPhiPt[ij] = 0;
319 fh2DijetAsymPt[ij] = 0;
320 fh2DijetPt2vsPt1[ij] = 0;
321 fh2DijetDifvsSum[ij] = 0;
322 }
323
324 DefineOutput(1, TList::Class());
3b7ffecf 325}
326
327
328
329Bool_t AliAnalysisTaskJetSpectrum2::Notify()
330{
a31b8a87 331
332
333
3b7ffecf 334 //
335 // Implemented Notify() to read the cross sections
336 // and number of trials from pyxsec.root
337 //
a31b8a87 338
339 // Fetch the aod also from the input in,
340 // have todo it in notify
341
342
343 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
344 // assume that the AOD is in the general output...
345 fAODOut = AODEvent();
346
347 if(fNonStdFile.Length()!=0){
348 // case that we have an AOD extension we need can fetch the jets from the extended output
349 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
350 fAODExtension = aodH->GetExtension(fNonStdFile.Data());
351 if(!fAODExtension){
352 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
353 }
354 }
355
3b7ffecf 356
357 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
358 Float_t xsection = 0;
359 Float_t ftrials = 1;
360
361 fAvgTrials = 1;
362 if(tree){
363 TFile *curfile = tree->GetCurrentFile();
364 if (!curfile) {
365 Error("Notify","No current file");
366 return kFALSE;
367 }
368 if(!fh1Xsec||!fh1Trials){
369 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
370 return kFALSE;
371 }
372 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
373 fh1Xsec->Fill("<#sigma>",xsection);
374 // construct a poor man average trials
375 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
a221560c 376 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
3b7ffecf 377 }
a31b8a87 378
379
3b7ffecf 380 return kTRUE;
381}
382
383void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
384{
385
3b7ffecf 386
387 // Connect the AOD
388
389
a31b8a87 390
3b7ffecf 391
392
393 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
394
395 OpenFile(1);
396 if(!fHistList)fHistList = new TList();
5cb0f438 397 fHistList->SetOwner(kTRUE);
3b7ffecf 398 Bool_t oldStatus = TH1::AddDirectoryStatus();
399 TH1::AddDirectory(kFALSE);
400
a31b8a87 401
402 MakeJetContainer();
403 fHistList->Add(fhnCorrelation);
404 fHistList->Add(fhnCorrelationPhiZRec);
b5d90baf 405 for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
a31b8a87 406
3b7ffecf 407 //
408 // Histogram
409
cba109a0 410 const Int_t nBinPt = 320;
3b7ffecf 411 Double_t binLimitsPt[nBinPt+1];
412 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
413 if(iPt == 0){
414 binLimitsPt[iPt] = 0.0;
415 }
416 else {// 1.0
edfbe476 417 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
3b7ffecf 418 }
419 }
420
edfbe476 421 const Int_t nBinPhi = 90;
3b7ffecf 422 Double_t binLimitsPhi[nBinPhi+1];
423 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
424 if(iPhi==0){
cc0649e4 425 binLimitsPhi[iPhi] = -1.*TMath::Pi();
3b7ffecf 426 }
427 else{
428 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
429 }
430 }
431
edfbe476 432
c8eabe24 433 const Int_t nBinPhi2 = 360;
434 Double_t binLimitsPhi2[nBinPhi2+1];
435 for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
436 if(iPhi2==0){
437 binLimitsPhi2[iPhi2] = 0.;
438 }
439 else{
440 binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
441 }
442 }
443
edfbe476 444 const Int_t nBinEta = 40;
445 Double_t binLimitsEta[nBinEta+1];
446 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
447 if(iEta==0){
448 binLimitsEta[iEta] = -2.0;
449 }
450 else{
cf049e94 451 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
edfbe476 452 }
453 }
454
3b7ffecf 455 const Int_t nBinFrag = 25;
456
457
a31b8a87 458
459
3b7ffecf 460 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
461 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
a31b8a87 462 fHistList->Add(fh1Xsec);
3b7ffecf 463
464 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
465 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
a31b8a87 466 fHistList->Add(fh1Trials);
3b7ffecf 467
468 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 469 fHistList->Add(fh1PtHard);
470
3b7ffecf 471 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 472 fHistList->Add(fh1PtHardNoW);
473
3b7ffecf 474 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 475 fHistList->Add(fh1PtHardTrials);
3b7ffecf 476
a31b8a87 477
57ca1193 478 fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
a31b8a87 479 fHistList->Add(fh1ZVtx);
480
481 fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
482 fHistList->Add(fh2PtFGen);
57ca1193 483
a31b8a87 484 fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};#Delta p_{T}/p_{T,Gen}",nBinPt,binLimitsPt,120,-2.4,2.4);
485 fHistList->Add(fh2RelPtFGen);
3b7ffecf 486
a31b8a87 487
488 if(fUseOldFill){
489
490 fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
491 fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
492
edfbe476 493 fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
494 fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
495 fh1SumPtTrackAreaRec = new TH1F("fh1SumPtTrackAreaRec","Sum Rec track P_T #eta <0.9 / 1.8 * 2 * 0.4*0.4;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
cc0649e4 496
497 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
d8f21f85 498 fh1PtJetsGenIn = new TH1F("fh1PtJetsGenIn","Gen jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
565584e8 499 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
cc0649e4 500 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
565584e8 501 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
fe77112a 502 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
edfbe476 503
cc0649e4 504 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
505 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
cc0649e4 506
a31b8a87 507 fh2NGenJetsPt = new TH2F("fh2NGenJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
508 fh2NGenTracksPt = new TH2F("fh2NGenTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
509 //
c8eabe24 510
3b7ffecf 511 for(int ij = 0;ij < kMaxJets;++ij){
512 fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
513 fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
514
edfbe476 515 fh2PhiPt[ij] = new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
516 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
517
cc0649e4 518 fh2PhiEta[ij] = new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
edfbe476 519 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
520
c8eabe24 521 fh2RhoPtRec[ij] = new TH2F(Form("fh2RhoPtRec_j%d",ij),"jet shape rho for jets;r;p_{T,rec};",
fb03edbe 522 40,0.,1.,nBinPt,binLimitsPt);
c8eabe24 523 fh2PsiPtRec[ij] = new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
fb03edbe 524 40,0.,2.,nBinPt,binLimitsPt);
c8eabe24 525
526 fh2RhoPtGen[ij] = new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
fb03edbe 527 40,0.,2.,nBinPt,binLimitsPt);
c8eabe24 528 fh2PsiPtGen[ij] = new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
fb03edbe 529 40,0.,2.,nBinPt,binLimitsPt);
530 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
c8eabe24 531
3b7ffecf 532 fh2FragRec[ij] = new TH2F(Form("fh2FragRec_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
533 nBinFrag,0.,1.,nBinPt,binLimitsPt);
534 fh2FragLnRec[ij] = new TH2F(Form("fh2FragLnRec_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
535 nBinFrag,0.,10.,nBinPt,binLimitsPt);
536
537 fh2FragGen[ij] = new TH2F(Form("fh2FragGen_j%d",ij),"Jet Fragmentation;x=p_{T,i}/p_{T,jet};p_{T,jet};1/N_{jet}dN_{ch}/dx",
c8eabe24 538 nBinFrag,0.,1.,nBinPt,binLimitsPt);
3b7ffecf 539 fh2FragLnGen[ij] = new TH2F(Form("fh2FragLnGen_j%d",ij),"Jet Fragmentation Ln;#xi=ln(p_{T,jet}/p_{T,i});p_{T,jet}(GeV);1/N_{jet}dN_{ch}/d#xi",
c8eabe24 540 nBinFrag,0.,10.,nBinPt,binLimitsPt);
3b7ffecf 541 }
542
26fa06fb 543 // Dijet histograms
c2785065 544 //background histograms
545 if(fBkgSubtraction){
185fa8e6 546 fh1Bkg1 = new TH1F("fh1Bkg1","Background estimate 1",100,0.,10.);
547 fh1Bkg2 = new TH1F("fh1Bkg2","Background estimate 2",100,0.,10.);
7f9012c1 548 fh1Bkg3 = new TH1F("fh1Bkg3","Background estimate 3",100,0.,10.);
c2785065 549 fh1Sigma1 = new TH1F("fh1Sigma1","Background fluctuations 1",100,0.,10.);
550 fh1Sigma2 = new TH1F("fh1Sigma2","Background fluctuations 2",100,0.,10.);
7f9012c1 551 fh1Sigma3 = new TH1F("fh1Sigma3","Background fluctuations 3",100,0.,10.);
185fa8e6 552 fh1Area1 = new TH1F("fh1Area1","Background mean area 1",50,0.,5.);
553 fh1Area2 = new TH1F("fh1Area2","Background mean area 2",50,0.,5.);
7f9012c1 554 fh1Area3 = new TH1F("fh1Area3","Background mean area 3",50,0.,5.);
c2785065 555 fh1Ptjet = new TH1F("fh1Ptjet","Jet spectrum",100,0.,200.);
556 fh1Ptjetsub1 = new TH1F("fh1Ptjetsub1","Subtracted spectrum 1",50,0.,200.);
557 fh1Ptjetsub2 = new TH1F("fh1Ptjetsub2","Subtracted spectrum 2",50,0.,200.);
7f9012c1 558 fh1Ptjetsub3 = new TH1F("fh1Ptjetsub3","Subtracted spectrum 3",50,0.,200.);
c2785065 559 fh1Ptjethardest = new TH1F("fh1Ptjethardest","Hardest jet spectrum",50,0.,200.);
560 fh1Ptjetsubhardest1 = new TH1F("fh1Pthardestsub1","Subtracted hardest jet spectrum 1",100,0.,200.);
561 fh1Ptjetsubhardest2 = new TH1F("fh1Pthardestsub2","Subtracted hardest jet spectrum 2",100,0.,200.);
7f9012c1 562 fh1Ptjetsubhardest3 = new TH1F("fh1Pthardestsub3","Subtracted hardest jet spectrum 3",100,0.,200.);
6bd3fdae 563 fh2Rhovspthardest1 = new TH2F("fh2Rhovspthardest1","Background vs pTjet 1",100,0.,200.,50,0.,5.);
564 fh2Rhovspthardest2 = new TH2F("fh2Rhovspthardest2","Background vs pTjet 2",100,0.,200.,50,0.,5.);
565 fh2Rhovspthardest3 = new TH2F("fh2Rhovspthardest3","Background vs pTjet 3",100,0.,200.,50,0.,5.);
566 fh2Errorvspthardest1 = new TH2F("fh2Errorvspthardest1","Relative error 1",100,0.,200.,50,0.,5.);
567 fh2Errorvspthardest2 = new TH2F("fh2Errorvspthardest2","Relative error 2",100,0.,200.,50,0.,5.);
568 fh2Errorvspthardest3 = new TH2F("fh2Errorvspthardest3","Relative error 3",100,0.,200.,50,0.,5.);
c2785065 569 }
570
26fa06fb 571
3b7ffecf 572 const Int_t saveLevel = 3; // large save level more histos
573 if(saveLevel>0){
a31b8a87 574 if(fBranchGen.Length()>0||fBkgSubtraction){
575 fHistList->Add(fh1NGenJets);
576 fHistList->Add(fh1PtTracksGenIn);
cc0649e4 577 }
578 fHistList->Add(fh1PtJetsRecIn);
d8f21f85 579 fHistList->Add(fh1PtJetsGenIn);
565584e8 580 fHistList->Add(fh1PtJetsLeadingRecIn);
cc0649e4 581 fHistList->Add(fh1PtTracksRecIn);
565584e8 582 fHistList->Add(fh1PtTracksLeadingRecIn);
3b7ffecf 583 fHistList->Add(fh1NRecJets);
edfbe476 584 fHistList->Add(fh1PtTrackRec);
585 fHistList->Add(fh1SumPtTrackRec);
586 fHistList->Add(fh1SumPtTrackAreaRec);
cc0649e4 587 fHistList->Add(fh2NRecJetsPt);
588 fHistList->Add(fh2NRecTracksPt);
a31b8a87 589 fHistList->Add(fh2NGenJetsPt);
590 fHistList->Add(fh2NGenTracksPt);
3b7ffecf 591 for(int ij = 0;ij<kMaxJets;++ij){
592 fHistList->Add( fh1PtRecIn[ij]);
c8eabe24 593
f4132e7d 594 if(fBranchGen.Length()>0||fBkgSubtraction){
ffd9d3ff 595 fHistList->Add(fh1PtGenIn[ij]);
596 fHistList->Add(fh2FragGen[ij]);
597 fHistList->Add(fh2FragLnGen[ij]);
c8eabe24 598 fHistList->Add(fh2RhoPtGen[ij]);
599 fHistList->Add(fh2PsiPtGen[ij]);
cc0649e4 600 }
edfbe476 601 fHistList->Add( fh2PhiPt[ij]);
602 fHistList->Add( fh2PhiEta[ij]);
ffd9d3ff 603 fHistList->Add( fh2RhoPtRec[ij]);
604 fHistList->Add( fh2PsiPtRec[ij]);
3b7ffecf 605 fHistList->Add( fh2FragRec[ij]);
606 fHistList->Add( fh2FragLnRec[ij]);
3b7ffecf 607 }
a31b8a87 608 // fHistList->Add(fh2TrackPtTrackPhi);
c2785065 609 if(fBkgSubtraction){
185fa8e6 610 fHistList->Add(fh1Bkg1);
611 fHistList->Add(fh1Bkg2);
7f9012c1 612 fHistList->Add(fh1Bkg3);
c2785065 613 fHistList->Add(fh1Sigma1);
614 fHistList->Add(fh1Sigma2);
7f9012c1 615 fHistList->Add(fh1Sigma3);
185fa8e6 616 fHistList->Add(fh1Area1);
617 fHistList->Add(fh1Area2);
7f9012c1 618 fHistList->Add(fh1Area3);
c2785065 619 fHistList->Add(fh1Ptjet);
620 fHistList->Add(fh1Ptjethardest);
621 fHistList->Add(fh1Ptjetsub1);
622 fHistList->Add(fh1Ptjetsub2);
7f9012c1 623 fHistList->Add(fh1Ptjetsub3);
c2785065 624 fHistList->Add(fh1Ptjetsubhardest1);
625 fHistList->Add(fh1Ptjetsubhardest2);
7f9012c1 626 fHistList->Add(fh1Ptjetsubhardest3);
6bd3fdae 627 fHistList->Add(fh2Rhovspthardest1);
628 fHistList->Add(fh2Rhovspthardest2);
629 fHistList->Add(fh2Rhovspthardest3);
630 fHistList->Add(fh2Errorvspthardest1);
631 fHistList->Add(fh2Errorvspthardest2);
632 fHistList->Add(fh2Errorvspthardest3);
c2785065 633 }
3b7ffecf 634 }
635
a31b8a87 636
637 }// ue old fill
638 else{
639
640 for(int ij = 0;ij <kJetTypes;++ij){
641 TString cAdd = "";
642 TString cJetBranch = "";
643 if(ij==kJetRec){
644 cAdd = "Rec";
645 cJetBranch = fBranchRec.Data();
646 }
647 else if (ij==kJetGen){
648 cAdd = "Gen";
649 cJetBranch = fBranchGen.Data();
650 }
651
652 fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
653 fHistList->Add(fh1NJets[ij]);
654
655 fh1PtJetsIn[ij] = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
656 fHistList->Add(fh1PtJetsIn[ij]);
657
658 fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
659 fHistList->Add(fh1PtTracksIn[ij]);
660
661 fh1PtTracksLeadingIn[ij] = new TH1F(Form("fh1PtTracksLeading%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
662 fHistList->Add(fh1PtTracksLeadingIn[ij]);
663
664 fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
665 fHistList->Add(fh1SumPtTrack[ij]);
666
667 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);
668 fHistList->Add(fh2NJetsPt[ij]);
669
670 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,50,-0.5,49.5);
671 fHistList->Add(fh2NTracksPt[ij]);
672
673 fh2LeadingJetPtJetPhi[ij] = new TH2F(Form("fh2Leading%sJetPtJetPhi",cAdd.Data()),Form("phi of leading %s jet;p_{T};#phi;",cAdd.Data()),
674 nBinPt,binLimitsPt,nBinPhi,binLimitsPhi);
675 fHistList->Add(fh2LeadingJetPtJetPhi[ij]);
676
677 fh2LeadingTrackPtTrackPhi[ij] = new TH2F(Form("fh2Leading%sTrackPtTrackPhi",cAdd.Data()),Form("phi of leading %s track;p_{T};#phi;",cAdd.Data()),
678 nBinPt,binLimitsPt,nBinPhi,binLimitsPhi);
679 fHistList->Add(fh2LeadingTrackPtTrackPhi[ij]);
680
681 for(int i = 0;i < kMaxJets;++i){
682
683 fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
684 fHistList->Add(fh1PtIn[ij][i]);
685
686 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()),
687 40,0.,2.,nBinPt,binLimitsPt);
688 fHistList->Add(fh2RhoPt[ij][i]);
689 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
690
691
692 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()),
693 40,0.,2.,nBinPt,binLimitsPt);
694 fHistList->Add(fh2PsiPt[ij][i]);
695 }
696
697
698 fh1DijetMinv[ij] = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
699 fHistList->Add(fh1DijetMinv[ij]);
700
b5d90baf 701 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 702 fHistList->Add(fh2DijetDeltaPhiPt[ij]);
703
704 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);
705 fHistList->Add(fh2DijetAsymPt[ij]);
706
707 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.);
708 fHistList->Add(fh2DijetPt2vsPt1[ij]);
709
710 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.);
711 fHistList->Add( fh2DijetDifvsSum[ij]);
712 }
713 }// new filling
714
3b7ffecf 715 // =========== Switch on Sumw2 for all histos ===========
716 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
717 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
718 if (h1){
719 h1->Sumw2();
720 continue;
721 }
722 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
723 if(hn)hn->Sumw2();
724 }
725 TH1::AddDirectory(oldStatus);
726}
727
728void AliAnalysisTaskJetSpectrum2::Init()
729{
730 //
731 // Initialization
732 //
733
3b7ffecf 734 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
735
736}
737
a31b8a87 738void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
739
b5cc0c6d 740
f2dd0695 741 Bool_t selected = kTRUE;
742
743 if(fUseGlobalSelection&&fEventSelectionMask==0){
744 selected = AliAnalysisHelperJetTasks::Selected();
745 }
746 else if(fUseGlobalSelection&&fEventSelectionMask>0){
45a11aef 747 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
f2dd0695 748 }
749
f4132e7d 750 if(fEventClass>0){
751 selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
752 }
753
f2dd0695 754 if(!selected){
b5cc0c6d 755 // no selection by the service task, we continue
f4132e7d 756 if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
b5cc0c6d 757 PostData(1, fHistList);
758 return;
759 }
760
f2dd0695 761
a31b8a87 762 static AliAODEvent* aod = 0;
763
764 // take all other information from the aod we take the tracks from
765 if(!aod){
766 if(fUseAODTrackInput)aod = fAODIn;
767 else aod = fAODOut;
768 }
769
770
3b7ffecf 771 //
772 // Execute analysis for current event
773 //
a31b8a87 774 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
775 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
776 if(!aodH){
777 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
778 return;
779 }
780
781 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
782 TClonesArray *aodRecJets = 0;
783
784 if(fAODOut&&!aodRecJets){
785 aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
786 }
787 if(fAODExtension&&!aodRecJets){
788 aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
789 }
790 if(fAODIn&&!aodRecJets){
791 aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
792 }
793
794
795
796 if(!aodRecJets){
797 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
798 return;
799 }
800
801 TClonesArray *aodGenJets = 0;
802 if(fBranchGen.Length()>0){
803 if(fAODOut&&!aodGenJets){
804 aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
805 }
806 if(fAODExtension&&!aodGenJets){
807 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
808 }
809 if(fAODIn&&!aodGenJets){
810 aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
811 }
812
813 if(!aodGenJets){
814 Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
3b7ffecf 815 return;
816 }
3b7ffecf 817 }
a31b8a87 818
819
820 if(false){
821 UserExecOld();
822 return;
823 }
824
825
826 // new Scheme
827 // first fill all the pure histograms, i.e. full jets
828 // in case of the correaltion limit to the n-leading jets
829
830 // reconstructed
831
832
833 // generated
834
835
836 // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
837
838
839 TList genJetsList; // full acceptance
840 TList genJetsListCut; // acceptance cut
841 TList recJetsList; // full acceptance
842 TList recJetsListCut; // acceptance cut
843
844
845 GetListOfJets(&recJetsList,aodRecJets,0);
846 GetListOfJets(&recJetsListCut,aodRecJets,1);
847
848 if(fBranchGen.Length()>0){
849 GetListOfJets(&genJetsList,aodGenJets,0);
850 GetListOfJets(&genJetsListCut,aodGenJets,1);
851 }
852
853 Double_t eventW = 1;
854 Double_t ptHard = 0;
855 Double_t nTrials = 1; // Trials for MC trigger
856 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
857
858 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
859 // this is the part we only use when we have MC information
860 AliMCEvent* mcEvent = MCEvent();
861 // AliStack *pStack = 0;
862 if(!mcEvent){
863 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
3b7ffecf 864 return;
865 }
a31b8a87 866 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
867 if(pythiaGenHeader){
868 nTrials = pythiaGenHeader->Trials();
869 ptHard = pythiaGenHeader->GetPtHard();
870 int iProcessType = pythiaGenHeader->ProcessType();
871 // 11 f+f -> f+f
872 // 12 f+barf -> f+barf
873 // 13 f+barf -> g+g
874 // 28 f+g -> f+g
875 // 53 g+g -> f+barf
876 // 68 g+g -> g+g
877 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
878 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
879 }
880 }// (fAnalysisType&kMCESD)==kMCESD)
881
882
883 // we simply fetch the tracks/mc particles as a list of AliVParticles
884
885 TList recParticles;
886 TList genParticles;
887
888 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
889 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
890 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
891 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
892
893 // Event control and counting ...
894 // MC
895 fh1PtHard->Fill(ptHard,eventW);
896 fh1PtHardNoW->Fill(ptHard,1);
897 fh1PtHardTrials->Fill(ptHard,nTrials);
898
899 // Real
900 fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
901
902 // the loops for rec and gen should be indentical... pass it to a separate
903 // function ...
904 // Jet Loop
905 // Track Loop
906 // Jet Jet Loop
907 // Jet Track loop
908
909 FillJetHistos(recJetsListCut,recParticles,kJetRec);
910 FillTrackHistos(recParticles,kJetRec);
911
912 FillJetHistos(genJetsListCut,genParticles,kJetGen);
913 FillTrackHistos(genParticles,kJetGen);
914
915 // Here follows the jet matching part
916 // delegate to a separated method?
917
918 FillMatchHistos(recJetsList,genJetsList);
919
920 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
921 PostData(1, fHistList);
922}
923
924void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
925
926 if(iType>=kJetTypes){
927 return;
3b7ffecf 928 }
a31b8a87 929
930 Int_t nJets = jetsList.GetEntries();
931 fh1NJets[iType]->Fill(nJets);
932
933 if(nJets<=0)return;
3b7ffecf 934
a31b8a87 935 Float_t ptOld = 1.E+32;
936 Float_t pT0 = 0;
937 Float_t pT1 = 0;
938 Float_t phi0 = 0;
939 Float_t phi1 = 0;
940 Int_t ij0 = -1;
941 Int_t ij1 = -1;
942
943
944 for(int ij = 0;ij < nJets;ij++){
945 AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
946 Float_t ptJet = jet->Pt();
947 fh1PtJetsIn[iType]->Fill(ptJet);
948 if(ptJet>ptOld){
949 Printf("%s:%d Jets Type %d Not Sorted !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,iType,ij,ptJet,ij-1,ptOld);
950 }
951 ptOld = ptJet;
952
b5d90baf 953 // find the dijets assume sorting and acceptance cut...
a31b8a87 954 if(ij==0){
955 ij0 = ij;
956 pT0 = ptJet;
957 phi0 = jet->Phi();
958 if(phi0<0)phi0+=TMath::Pi()*2.;
959 }
960 else if(ptJet>pT1){
961 // take only the backward hemisphere??
962 phi1 = jet->Phi();
b5d90baf 963 if(phi1<0)phi1+=TMath::Pi()*2.;
a31b8a87 964 Float_t dPhi = phi1 - phi0;
965 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
966 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
967 if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
968 ij1 = ij;
969 pT1 = ptJet;
970 }
971 }
972 // fill jet histos for kmax jets
973 if(ij<kMaxJets){
974 Float_t phiJet = jet->Phi();
975 if(phiJet<0)phiJet+=TMath::Pi()*2.;
976 fh1TmpRho->Reset();
977 fh1PtIn[iType][ij]->Fill(ptJet);
978 // fill leading jets...
979 if(ij==0){
980 fh2LeadingJetPtJetPhi[iType]->Fill(ptJet,ptJet);
981 if(ij==0&&iType==0&&fDebug>1){
982 Printf("%d %3.3f %p %s",iType,ptJet,jet,fBranchRec.Data());
983 }
984 }
985 if(particlesList.GetSize()){
986
987 // Particles... correlated with jets...
988 for(int it = 0;it<particlesList.GetEntries();++it){
989 AliVParticle *part = (AliVParticle*)particlesList.At(it);
990 Float_t deltaR = jet->DeltaR(part);
991 if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
992 }
993 // fill the jet shapes
994 Float_t rhoSum = 0;
995 for(int ibx = 1;ibx<=fh2RhoPt[iType][ij]->GetNbinsX();ibx++){
996 Float_t r = fh2RhoPt[iType][ij]->GetXaxis()->GetBinCenter(ibx);
997 Float_t rho = fh1TmpRho->GetBinContent(ibx);
998 rhoSum += rho;
999 fh2RhoPt[iType][ij]->Fill(r,ptJet,rho);
1000 fh2PsiPt[iType][ij]->Fill(r,ptJet,rhoSum);
1001 }
1002 }// if we have particles
1003 }// ij < kMaxJets
1004 }// Jet Loop
1005
1006
1007 // do something with dijets...
1008 if(ij0>=0&&ij1>0){
1009 AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1010 Double_t ptJet0 = jet0->Pt();
1011 Double_t phiJet0 = jet0->Phi();
1012 if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1013
1014 AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1015 Double_t ptJet1 = jet1->Pt();
1016 Double_t phiJet1 = jet1->Phi();
1017 if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1018
1019 Float_t deltaPhi = phiJet0 - phiJet1;
1020 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1021 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1022 deltaPhi = TMath::Abs(deltaPhi);
b5d90baf 1023 fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
a31b8a87 1024
1025 Float_t asym = 9999;
1026 if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1027 fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1028 fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1029 fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1030 Float_t minv = 2.*(jet0->P()*jet1->P()-
1031 jet0->Px()*jet1->Px()-
1032 jet0->Py()*jet1->Py()-
1033 jet0->Pz()*jet1->Pz()); // assume mass == 0;
1034 if(minv<0)minv=0; // prevent numerical instabilities
1035 minv = TMath::Sqrt(minv);
1036 fh1DijetMinv[iType]->Fill(minv);
1037 }
1038
1039
1040
1041 // count down the jets above thrueshold
1042 Int_t nOver = nJets;
1043 if(nOver>0){
1044 TIterator *jetIter = jetsList.MakeIterator();
1045 AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
1046 Float_t pt = tmpJet->Pt();
1047 if(tmpJet){
1048 for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1049 Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1050 while(pt<ptCut&&tmpJet){
1051 nOver--;
1052 tmpJet = (AliAODJet*)(jetIter->Next());
1053 if(tmpJet){
1054 pt = tmpJet->Pt();
1055 }
1056 }
1057 if(nOver<=0)break;
1058 fh2NJetsPt[iType]->Fill(ptCut,nOver);
1059 }
1060 }
1061 delete jetIter;
1062 }
1063}
1064
1065void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1066
1067 Int_t nTrackOver = particlesList.GetSize();
1068 // do the same for tracks and jets
1069 if(nTrackOver>0){
1070 TIterator *trackIter = particlesList.MakeIterator();
1071 AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1072 Float_t pt = tmpTrack->Pt();
1073 for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1074 Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1075 while(pt<ptCut&&tmpTrack){
1076 nTrackOver--;
1077 tmpTrack = (AliVParticle*)(trackIter->Next());
1078 if(tmpTrack){
1079 pt = tmpTrack->Pt();
1080 }
1081 }
1082 if(nTrackOver<=0)break;
1083 fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1084 }
1085
1086
1087 trackIter->Reset();
1088 AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1089 Float_t sumPt = 0;
1090
1091 while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1092 Float_t tmpPt = tmpTrack->Pt();
1093 fh1PtTracksIn[iType]->Fill(tmpPt);
1094 sumPt += tmpPt;
1095 Float_t tmpPhi = tmpTrack->Phi();
1096 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1097 if(tmpTrack==leading){
1098 fh1PtTracksLeadingIn[iType]->Fill(tmpPt);
1099 fh2LeadingTrackPtTrackPhi[iType]->Fill(tmpPt,tmpPhi);
1100 continue;
1101 }
1102 }
1103 fh1SumPtTrack[iType]->Fill(sumPt);
1104 delete trackIter;
1105 }
1106
1107}
1108
1109
1110void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1111
1112
1113 // Fill al the matching histos
1114 // It is important that the acceptances for the mathing are as large as possible
1115 // to avoid false matches on the edge of acceptance
1116 // therefore we add some extra matching jets as overhead
1117
1118 static TArrayI aGenIndex(recJetsList.GetEntries());
1119 if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
3b7ffecf 1120
a31b8a87 1121 static TArrayI aRecIndex(genJetsList.GetEntries());
1122 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
3b7ffecf 1123
a31b8a87 1124 if(fDebug){
1125 Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1126 Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1127 }
1128 AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1129 &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1130 aGenIndex,aRecIndex,fDebug);
1131
1132 if(fDebug){
1133 for(int i = 0;i< aGenIndex.GetSize();++i){
1134 if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1135 }
1136 for(int i = 0;i< aRecIndex.GetSize();++i){
1137 if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1138 }
1139 }
1140
1141 Double_t container[6];
1142 Double_t containerPhiZ[6];
1143
1144 // loop over generated jets
1145 // consider the
1146 for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1147 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1148 Double_t ptGen = genJet->Pt();
1149 Double_t phiGen = genJet->Phi();
1150 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1151 Double_t etaGen = genJet->Eta();
1152 container[3] = ptGen;
1153 container[4] = etaGen;
1154 container[5] = phiGen;
1155 fhnJetContainer[kStep0]->Fill(&container[3]);
1156 if(JetSelected(genJet)){
1157 fhnJetContainer[kStep1]->Fill(&container[3]);
1158 Int_t ir = aRecIndex[ig];
1159 if(ir>=0&&ir<recJetsList.GetEntries()){
1160 fhnJetContainer[kStep2]->Fill(&container[3]);
1161 AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1162 if(JetSelected(recJet))fhnJetContainer[kStep4]->Fill(&container[3]);
1163 if(JetSelected(recJet))fhnJetContainer[kStep3]->Fill(&container[3]);
1164 }
1165 }
1166 }// loop over generated jets used for matching...
1167
1168
1169
1170 // loop over reconstructed jets
1171 for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1172 AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1173 Double_t etaRec = recJet->Eta();
1174 Double_t ptRec = recJet->Pt();
1175 Double_t phiRec = recJet->Phi();
1176 if(phiRec<0)phiRec+=TMath::Pi()*2.;
1177 // do something with dijets...
1178
1179 container[0] = ptRec;
1180 container[1] = etaRec;
1181 container[2] = phiRec;
1182 containerPhiZ[0] = ptRec;
1183 containerPhiZ[1] = phiRec;
1184
1185 fhnJetContainer[kStep0+kMaxStep]->Fill(container);
1186 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1187
1188 if(JetSelected(recJet)){
1189 fhnJetContainer[kStep1+kMaxStep]->Fill(container);
1190 // Fill Correlation
1191 Int_t ig = aGenIndex[ir];
1192 if(ig>=0 && ig<genJetsList.GetEntries()){
1193 fhnJetContainer[kStep2+kMaxStep]->Fill(container);
1194 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1195 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1196 Double_t ptGen = genJet->Pt();
1197 Double_t phiGen = genJet->Phi();
1198 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1199 Double_t etaGen = genJet->Eta();
1200
1201 container[3] = ptGen;
1202 container[4] = etaGen;
1203 container[5] = phiGen;
1204 containerPhiZ[3] = ptGen;
1205 //
1206 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1207 //
1208 if(JetSelected(genJet))fhnJetContainer[kStep4+kMaxStep]->Fill(container);
1209 fhnJetContainer[kStep3+kMaxStep]->Fill(container);
1210 fhnCorrelation->Fill(container);
1211 if(ptGen>0){
1212 Float_t delta = (ptRec-ptGen)/ptGen;
1213 fh2RelPtFGen->Fill(ptGen,delta);
1214 fh2PtFGen->Fill(ptGen,ptRec);
1215 }
1216 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1217 }
1218 else{
1219 containerPhiZ[3] = 0;
1220 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1221 }
1222 }// loop over reconstructed jets
1223 }
1224 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1225}
1226
1227
1228void AliAnalysisTaskJetSpectrum2::UserExecOld(Option_t */*option*/)
1229{
1230
1231
1232 Bool_t selected = kTRUE;
1233
1234
1235
1236
1237 if(fUseGlobalSelection&&fEventSelectionMask==0){
1238 selected = AliAnalysisHelperJetTasks::Selected();
1239 }
1240 else if(fUseGlobalSelection&&fEventSelectionMask>0){
1241 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
1242 }
1243
1244 if(fEventClass>0){
1245 selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
1246 }
1247
1248 if(!selected){
1249 // no selection by the service task, we continue
1250 if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
1251 PostData(1, fHistList);
1252 return;
1253 }
1254
1255
1256 //
1257 // Execute analysis for current event
1258 //
1259 static AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1260
1261
1262 // take all other information from the aod we take the tracks from
1263 static AliAODEvent *aod = 0;
1264 if(!aod){
1265 if(fUseAODTrackInput)aod = fAODIn;
1266 aod = fAODOut;
1267 }
1268 static AliAODEvent *fAOD = 0;
1269 if(!fAOD){
1270 if(fUseAODJetInput)fAOD = fAODIn;
1271 fAOD = fAODOut;
1272 }
3b7ffecf 1273
1274 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
1275
1276
1277 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1278
1279 if(!aodH){
1280 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
1281 return;
1282 }
1283
1284 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
a31b8a87 1285 TClonesArray *aodRecJets = 0;
1286
1287
1288
1289 if(fAOD&&!aodRecJets){
1290 aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
1291 }
1292 if(fAODExtension&&!aodRecJets){
1293 aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
1294 }
1295
3b7ffecf 1296 if(!aodRecJets){
1297 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
1298 return;
1299 }
1300
6bd3fdae 1301 Int_t nJets = aodRecJets->GetEntriesFast();
1302
a31b8a87 1303 TClonesArray *aodGenJets = 0;
1304 if(fBranchGen.Length()>0){
1305 if(fAOD&&!aodGenJets){
1306 aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
1307 }
1308 if(fAODExtension&&!aodGenJets){
1309 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
1310 }
1311
1312 if(!aodGenJets){
1313 Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
1314 return;
1315 }
1316 }
1317
1318
6bd3fdae 1319 // ==== General variables needed
1320 // We use statice array, not to fragment the memory
1321 AliAODJet genJets[kMaxJets];
1322 Int_t nGenJets = 0;
1323 AliAODJet recJets[kMaxJets];
1324 Int_t nRecJets = 0;
a31b8a87 1325
1326 // These will replace the arrays...
1327 TList genJetsList; // full acceptance
1328 TList genJetsListCut; // acceptance cut
1329 TList recJetsList; // full acceptance
1330 TList recJetsListCut; // acceptance cut
6bd3fdae 1331
1332
a31b8a87 1333 GetListOfJets(&genJetsList,aodGenJets,0);
1334 GetListOfJets(&genJetsListCut,aodGenJets,1);
6bd3fdae 1335
a31b8a87 1336 if(fBranchGen.Length()>0){
1337 GetListOfJets(&recJetsList,aodRecJets,0);
1338 GetListOfJets(&recJetsListCut,aodRecJets,1);
1339 }
1340 ///////////////////////////
6bd3fdae 1341
a31b8a87 1342
c2785065 1343
1344
a31b8a87 1345 if(fBkgSubtraction){
1346 AliAODJetEventBackground* evBkg=(AliAODJetEventBackground*)fAOD->FindListObject(fBranchBkg.Data());
c2785065 1347 if(!evBkg){
1348 Printf("%s:%d no reconstructed background array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkg.Data());
a31b8a87 1349 return;
c2785065 1350 }
3b7ffecf 1351
1352
c2785065 1353 ///just to start: some very simple plots containing rho, sigma and area of the background.
1354
6bd3fdae 1355
7f9012c1 1356 Float_t pthardest=0.;
80a790fd 1357
c2785065 1358 if(nJets!=0){
7f9012c1 1359
c2785065 1360 Float_t bkg1=evBkg->GetBackground(0);
1361 Float_t bkg2=evBkg->GetBackground(1);
7f9012c1 1362 Float_t bkg3=evBkg->GetBackground(2);
c2785065 1363 Float_t sigma1=evBkg->GetSigma(0);
1364 Float_t sigma2=evBkg->GetSigma(1);
7f9012c1 1365 Float_t sigma3=evBkg->GetSigma(2);
c2785065 1366 Float_t area1=evBkg->GetMeanarea(0);
1367 Float_t area2=evBkg->GetMeanarea(1);
7f9012c1 1368 Float_t area3=evBkg->GetMeanarea(2);
185fa8e6 1369 fh1Bkg1->Fill(bkg1); //rho computed with all background jets.
7f9012c1 1370 fh1Bkg2->Fill(bkg2); //rho computed with all background jets but the 2 hardest.
1371 fh1Bkg3->Fill(bkg3); //rho computed with randomized jets
c2785065 1372 fh1Sigma1->Fill(sigma1);
1373 fh1Sigma2->Fill(sigma2);
7f9012c1 1374 fh1Sigma3->Fill(sigma3);
185fa8e6 1375 fh1Area1->Fill(area1);
1376 fh1Area2->Fill(area2);
7f9012c1 1377 fh1Area3->Fill(area3);
c2785065 1378
80a790fd 1379 Int_t iSubJetCounter = 0;
c2785065 1380 for(Int_t k=0;k<nJets;k++){
1381 AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets->At(k));
1382 fh1Ptjet->Fill(jet->Pt());
1383 Float_t ptsub1=jet->Pt()-bkg1*jet->EffectiveAreaCharged();
1384 Float_t ptsub2=jet->Pt()-bkg2*jet->EffectiveAreaCharged();
7f9012c1 1385 Float_t ptsub3=jet->Pt()-bkg3*jet->EffectiveAreaCharged();
6bd3fdae 1386 if(ptsub2<0.) ptsub2=0.;
1387 if(ptsub1<0.) ptsub1=0.;
1388 if(ptsub3<0.) ptsub3=0.;
c2785065 1389 Float_t err1=sigma1*sqrt(area1);
1390 Float_t err2=sigma2*sqrt(area2);
7f9012c1 1391 Float_t err3=sigma3*sqrt(area3);
c2785065 1392 fh1Ptjetsub1->Fill(ptsub1);
1393 fh1Ptjetsub2->Fill(ptsub2);
7f9012c1 1394 fh1Ptjetsub3->Fill(ptsub3);
c2785065 1395 if(k==0) {
1396 pthardest=jet->Pt();
1397 fh1Ptjethardest->Fill(pthardest);
1398 fh1Ptjetsubhardest1->Fill(ptsub1);
1399 fh1Ptjetsubhardest2->Fill(ptsub2);
7f9012c1 1400 fh1Ptjetsubhardest3->Fill(ptsub3);
36c36a0c 1401 if (ptsub1 > 0)
1402 fh2Errorvspthardest1->Fill(ptsub1,err1/ptsub1);
1403 if (ptsub2 > 0)
1404 fh2Errorvspthardest2->Fill(ptsub2,err2/ptsub2);
1405 if (ptsub3 > 0)
1406 fh2Errorvspthardest3->Fill(ptsub3,err3/ptsub3);
c2785065 1407 }
7f9012c1 1408
f4132e7d 1409 Float_t ptsub=0.;
1410 if(fFillCorrBkg==1) ptsub=ptsub1;
1411 if(fFillCorrBkg==2) ptsub=ptsub2;
1412 if(fFillCorrBkg==3) ptsub=ptsub3;
80a790fd 1413 if(ptsub>0){// avoid unphysical jets pT
1414 Float_t subphi=jet->Phi();
1415 Float_t subtheta=jet->Theta();
1416 Float_t subpz = ptsub/TMath::Tan(subtheta);
1417 Float_t subpx=ptsub*TMath::Cos(subphi);
1418 Float_t subpy=ptsub * TMath::Sin(subphi);
1419 Float_t subp = TMath::Sqrt(ptsub*ptsub+subpz*subpz);
1420 if(k<kMaxJets){// only store the jets which are assoicated to anohter one
1421 genJets[iSubJetCounter].SetPxPyPzE(subpx,subpy,subpz,subp);
1422 iSubJetCounter++;
1423 }
f4132e7d 1424 }
6bd3fdae 1425 }
80a790fd 1426 nGenJets = iSubJetCounter;
6bd3fdae 1427 fh2Rhovspthardest1->Fill(pthardest,bkg1);
1428 fh2Rhovspthardest2->Fill(pthardest,bkg2);
1429 fh2Rhovspthardest3->Fill(pthardest,bkg3);
c2785065 1430 }
1431 }// background subtraction
1432
1433
3b7ffecf 1434 Double_t eventW = 1;
1435 Double_t ptHard = 0;
1436 Double_t nTrials = 1; // Trials for MC trigger
1437
1438 if(fUseExternalWeightOnly){
1439 eventW = fExternalWeight;
1440 }
1441
1442 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
cf049e94 1443 // if(fDebug>0)aodH->SetFillAOD(kFALSE);
3b7ffecf 1444 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1445 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
1446 // this is the part we only use when we have MC information
1447 AliMCEvent* mcEvent = MCEvent();
1448 // AliStack *pStack = 0;
1449 if(!mcEvent){
1450 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
1451 return;
1452 }
1453 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
1454 Int_t iCount = 0;
1455 if(pythiaGenHeader){
1456 nTrials = pythiaGenHeader->Trials();
1457 ptHard = pythiaGenHeader->GetPtHard();
1458 int iProcessType = pythiaGenHeader->ProcessType();
1459 // 11 f+f -> f+f
1460 // 12 f+barf -> f+barf
1461 // 13 f+barf -> g+g
1462 // 28 f+g -> f+g
1463 // 53 g+g -> f+barf
1464 // 68 g+g -> g+g
1465 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
1466 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
1467
1468 // fetch the pythia generated jets only to be used here
1469 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
1470 AliAODJet pythiaGenJets[kMaxJets];
1471 for(int ip = 0;ip < nPythiaGenJets;++ip){
1472 if(iCount>=kMaxJets)continue;
1473 Float_t p[4];
1474 pythiaGenHeader->TriggerJet(ip,p);
1475 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
1476
3b7ffecf 1477 if(fBranchGen.Length()==0){
9280dfa6 1478 /*
1479 if(fLimitGenJetEta){
1480 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
1481 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
1482 }
a31b8a87 1483 */
9280dfa6 1484 if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
3b7ffecf 1485 // if we have MC particles and we do not read from the aod branch
1486 // use the pythia jets
f4132e7d 1487 if(!fBkgSubtraction){
1488 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
1489 iCount++;
1490 }
3b7ffecf 1491 }
3b7ffecf 1492 }
1493 }
f4132e7d 1494 if(fBranchGen.Length()==0&&!fBkgSubtraction)nGenJets = iCount;
3b7ffecf 1495 }// (fAnalysisType&kMCESD)==kMCESD)
1496
1497
1498 // we simply fetch the tracks/mc particles as a list of AliVParticles
1499
1500 TList recParticles;
1501 TList genParticles;
1502
565584e8 1503
1504
1505
1506 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
1507 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
1508 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
1509 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
3b7ffecf 1510
cc0649e4 1511
3b7ffecf 1512 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1513 fh1PtHard->Fill(ptHard,eventW);
1514 fh1PtHardNoW->Fill(ptHard,1);
1515 fh1PtHardTrials->Fill(ptHard,nTrials);
57ca1193 1516 fh1ZVtx->Fill(fAOD->GetPrimaryVertex()->GetZ());
1517
3b7ffecf 1518
1519 // If we set a second branch for the input jets fetch this
f4132e7d 1520 if(fBranchGen.Length()>0&&!fBkgSubtraction){
3b7ffecf 1521 if(aodGenJets){
1522 Int_t iCount = 0;
1523 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
1524 if(iCount>=kMaxJets)continue;
1525 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
1526 if(!tmp)continue;
5301f754 1527 /*
3b7ffecf 1528 if(fLimitGenJetEta){
1529 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
1530 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
1531 }
5301f754 1532 */
9280dfa6 1533 if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
3b7ffecf 1534 genJets[iCount] = *tmp;
1535 iCount++;
1536 }
1537 nGenJets = iCount;
1538 }
1539 else{
5cb0f438 1540 if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
1541 if(fDebug>2)fAOD->Print();
3b7ffecf 1542 }
1543 }
1544
1545 fh1NGenJets->Fill(nGenJets);
1546 // We do not want to exceed the maximum jet number
1547 nGenJets = TMath::Min(nGenJets,kMaxJets);
1548
1549 // Fetch the reconstructed jets...
1550
1551 nRecJets = aodRecJets->GetEntries();
1552
cc0649e4 1553 nRecJets = aodRecJets->GetEntries();
1554 fh1NRecJets->Fill(nRecJets);
1555
1556 // Do something with the all rec jets
1557 Int_t nRecOver = nRecJets;
3b7ffecf 1558
cc0649e4 1559 // check that the jets are sorted
1560 Float_t ptOld = 999999;
1561 for(int ir = 0;ir < nRecJets;ir++){
1562 AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
1563 Float_t tmpPt = tmp->Pt();
1564 if(tmpPt>ptOld){
3568c3d6 1565 Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
cc0649e4 1566 }
1567 ptOld = tmpPt;
1568 }
3b7ffecf 1569
cc0649e4 1570
1571 if(nRecOver>0){
1572 TIterator *recIter = aodRecJets->MakeIterator();
1573 AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
1574 Float_t pt = tmpRec->Pt();
1575 if(tmpRec){
1576 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1577 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1578 while(pt<ptCut&&tmpRec){
1579 nRecOver--;
1580 tmpRec = (AliAODJet*)(recIter->Next());
1581 if(tmpRec){
1582 pt = tmpRec->Pt();
1583 }
1584 }
1585 if(nRecOver<=0)break;
1586 fh2NRecJetsPt->Fill(ptCut,nRecOver);
1587 }
1588 }
1589 recIter->Reset();
1590 AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
1591 Float_t phi = leading->Phi();
1592 if(phi<0)phi+=TMath::Pi()*2.;
cf049e94 1593 pt = leading->Pt();
45a11aef 1594 while((tmpRec = (AliAODJet*)(recIter->Next()))){
cc0649e4 1595 Float_t tmpPt = tmpRec->Pt();
1596 fh1PtJetsRecIn->Fill(tmpPt);
565584e8 1597 if(tmpRec==leading){
1598 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1599 continue;
1600 }
cc0649e4 1601 // correlation
1602 Float_t tmpPhi = tmpRec->Phi();
1603
1604 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1605 Float_t dPhi = phi - tmpRec->Phi();
1606 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1607 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
a31b8a87 1608 // Float_t dEta = eta - tmpRec->Eta();
1609 // fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1610 // fh2JetsLeadingPhiPt->Fill(dPhi,pt);
cc0649e4 1611 }
1612 delete recIter;
1613 }
1614
1615 Int_t nTrackOver = recParticles.GetSize();
1616 // do the same for tracks and jets
1617 if(nTrackOver>0){
1618 TIterator *recIter = recParticles.MakeIterator();
1619 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1620 Float_t pt = tmpRec->Pt();
1621 // Printf("Leading track p_t %3.3E",pt);
1622 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1623 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1624 while(pt<ptCut&&tmpRec){
1625 nTrackOver--;
1626 tmpRec = (AliVParticle*)(recIter->Next());
1627 if(tmpRec){
1628 pt = tmpRec->Pt();
1629 }
1630 }
1631 if(nTrackOver<=0)break;
1632 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1633 }
1634
1635 recIter->Reset();
1636 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1637 Float_t phi = leading->Phi();
1638 if(phi<0)phi+=TMath::Pi()*2.;
cf049e94 1639 pt = leading->Pt();
cc0649e4 1640 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1641 Float_t tmpPt = tmpRec->Pt();
1642 fh1PtTracksRecIn->Fill(tmpPt);
565584e8 1643 if(tmpRec==leading){
1644 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1645 continue;
1646 }
cc0649e4 1647 // correlation
1648 Float_t tmpPhi = tmpRec->Phi();
1649
1650 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1651 Float_t dPhi = phi - tmpRec->Phi();
1652 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1653 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
a31b8a87 1654 // Float_t dEta = eta - tmpRec->Eta();
1655 // fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1656 // fh2TracksLeadingPhiPt->Fill(dPhi,pt);
cc0649e4 1657 }
1658 delete recIter;
1659 }
1660
1661 if(genParticles.GetSize()){
1662 TIterator *genIter = genParticles.MakeIterator();
1663 AliVParticle *tmpGen = 0;
1664 while((tmpGen = (AliVParticle*)(genIter->Next()))){
1665 if(TMath::Abs(tmpGen->Eta())<0.9){
1666 Float_t tmpPt = tmpGen->Pt();
1667 fh1PtTracksGenIn->Fill(tmpPt);
1668 }
1669 }
1670 delete genIter;
1671 }
1672
3b7ffecf 1673 nRecJets = TMath::Min(nRecJets,kMaxJets);
6bd3fdae 1674 nJets = TMath::Min(nJets,kMaxJets);
9280dfa6 1675 Int_t iCountRec = 0;
3b7ffecf 1676 for(int ir = 0;ir < nRecJets;++ir){
1677 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
1678 if(!tmp)continue;
9280dfa6 1679 if(tmp->Pt()<fMinJetPt)continue;
e6a92928 1680 recJets[iCountRec] = *tmp;
9280dfa6 1681 iCountRec++;
3b7ffecf 1682 }
9280dfa6 1683 nRecJets = iCountRec;
3b7ffecf 1684
1685
1686 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1687 // Relate the jets
1688 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
1689 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
1690
cc0649e4 1691
1692 for(int i = 0;i<kMaxJets;++i){
3b7ffecf 1693 iGenIndex[i] = iRecIndex[i] = -1;
1694 }
1695
6bd3fdae 1696
3b7ffecf 1697 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
1698 iGenIndex,iRecIndex,fDebug);
6bd3fdae 1699
a31b8a87 1700 static TArrayI aGenIndex(recJetsListCut.GetEntries());
1701 if(aGenIndex.GetSize()<recJetsListCut.GetEntries())aGenIndex.Set(recJetsListCut.GetEntries());
1702
1703 static TArrayI aRecIndex(genJetsList.GetEntries());
1704 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1705
1706 if(fDebug){
1707 Printf("New Gens List %d rec index Array %d count %d",genJetsList.GetEntries(),aRecIndex.GetSize(),nGenJets);
1708 Printf("New Rec List %d gen indey Array %d count %d",recJetsListCut.GetEntries(),aGenIndex.GetSize(),nRecJets);
1709 }
1710 AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,nGenJets,&recJetsListCut,nRecJets,
1711 aGenIndex,aRecIndex,fDebug);
1712
6bd3fdae 1713
3b7ffecf 1714 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1715
1716 if(fDebug){
a31b8a87 1717 for(int i = 0;i<kMaxJets;++i){
3b7ffecf 1718 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
1719 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
1720 }
1721 }
1722
1723
1724
1725
1726 Double_t container[6];
c8eabe24 1727 Double_t containerPhiZ[6];
3b7ffecf 1728
1729 // loop over generated jets
1730
1731 // radius; tmp, get from the jet header itself
1732 // at some point, how todeal woht FastJet on MC?
1733 Float_t radiusGen =0.4;
1734 Float_t radiusRec =0.4;
1735
1736 for(int ig = 0;ig < nGenJets;++ig){
1737 Double_t ptGen = genJets[ig].Pt();
1738 Double_t phiGen = genJets[ig].Phi();
1739 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1740 Double_t etaGen = genJets[ig].Eta();
d8f21f85 1741 fh1PtJetsGenIn->Fill(ptGen);
1742
3b7ffecf 1743 container[3] = ptGen;
1744 container[4] = etaGen;
1745 container[5] = phiGen;
1746 fhnJetContainer[kStep0]->Fill(&container[3],eventW);
1747 Int_t ir = iRecIndex[ig];
1748
a31b8a87 1749 if(TMath::Abs(etaGen)<fJetRecEtaWindow){
c8eabe24 1750 fh1TmpRho->Reset();
1751
3b7ffecf 1752 fhnJetContainer[kStep1]->Fill(&container[3],eventW);
1753 fh1PtGenIn[ig]->Fill(ptGen,eventW);
1754 // fill the fragmentation function
1755 for(int it = 0;it<genParticles.GetEntries();++it){
1756 AliVParticle *part = (AliVParticle*)genParticles.At(it);
c8eabe24 1757 Float_t deltaR = genJets[ig].DeltaR(part);
80a790fd 1758 if(ptGen>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
1759 if(deltaR<radiusGen&&ptGen>0){
3b7ffecf 1760 Float_t z = part->Pt()/ptGen;
1761 Float_t lnz = -1.*TMath::Log(z);
1762 fh2FragGen[ig]->Fill(z,ptGen,eventW);
1763 fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
1764 }
c8eabe24 1765
3b7ffecf 1766 }
c8eabe24 1767 Float_t rhoSum = 0;
fb03edbe 1768 for(int ibx = 1;ibx<=fh2RhoPtGen[ir]->GetNbinsX();ibx++){
c8eabe24 1769 Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
1770 Float_t rho = fh1TmpRho->GetBinContent(ibx);
1771 rhoSum += rho;
1772 fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
1773 fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
3b7ffecf 1774 }
c8eabe24 1775 }
3b7ffecf 1776 if(ir>=0&&ir<nRecJets){
1777 fhnJetContainer[kStep2]->Fill(&container[3],eventW);
1778 Double_t etaRec = recJets[ir].Eta();
a31b8a87 1779 if(TMath::Abs(etaRec)<fJetRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
1780 if(TMath::Abs(etaRec)<fJetRecEtaWindow&&TMath::Abs(etaGen)<fJetRecEtaWindow)fhnJetContainer[kStep3]->Fill(&container[3],eventW);
3b7ffecf 1781 }
1782 }// loop over generated jets
7fa8b2da 1783
edfbe476 1784 Float_t sumPt = 0;
1785 for(int it = 0;it<recParticles.GetEntries();++it){
1786 AliVParticle *part = (AliVParticle*)recParticles.At(it);
1787 // fill sum pt and P_t of all paritles
1788 if(TMath::Abs(part->Eta())<0.9){
1789 Float_t pt = part->Pt();
1790 fh1PtTrackRec->Fill(pt,eventW);
a31b8a87 1791 // fh2TrackPtTrackPhi->Fill(pt,part->Phi());
edfbe476 1792 sumPt += pt;
1793 }
1794 }
1795 fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
1796 fh1SumPtTrackRec->Fill(sumPt,eventW);
1797
cf049e94 1798
3b7ffecf 1799 // loop over reconstructed jets
1800 for(int ir = 0;ir < nRecJets;++ir){
1801 Double_t etaRec = recJets[ir].Eta();
1802 Double_t ptRec = recJets[ir].Pt();
1803 Double_t phiRec = recJets[ir].Phi();
1804 if(phiRec<0)phiRec+=TMath::Pi()*2.;
26fa06fb 1805
1806
1807 // do something with dijets...
26fa06fb 1808
3b7ffecf 1809 container[0] = ptRec;
1810 container[1] = etaRec;
1811 container[2] = phiRec;
c8eabe24 1812 containerPhiZ[0] = ptRec;
1813 containerPhiZ[1] = phiRec;
ffab794c 1814 if(ptRec>30.&&fDebug>2){
7fa8b2da 1815 // need to cast to int, otherwise the printf overwrites
1816 Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
1277f815 1817 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
1818 if(fESD)Printf("ESDEvent GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
cf049e94 1819 // aodH->SetFillAOD(kTRUE);
7fa8b2da 1820 fAOD->GetHeader()->Print();
a221560c 1821 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
7fa8b2da 1822 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1823 AliAODTrack *tr = fAOD->GetTrack(it);
1824 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1825 tr->Print();
38ecb6a5 1826 // tr->Dump();
7fa8b2da 1827 if(fESD){
1828 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1829 esdTr->Print("");
bb3d13aa 1830 // esdTr->Dump();
7fa8b2da 1831 }
1832 }
1833 }
1834
1835
3b7ffecf 1836 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
1837 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
45a11aef 1838
c8eabe24 1839 Float_t zLeading = -1;
a31b8a87 1840 if(TMath::Abs(etaRec)<fJetRecEtaWindow){
1841 // fh2JetPtJetPhi->Fill(ptRec,phiRec);
3b7ffecf 1842 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
1843 fh1PtRecIn[ir]->Fill(ptRec,eventW);
1844 // fill the fragmentation function
c8eabe24 1845
1846 fh1TmpRho->Reset();
1847
3b7ffecf 1848 for(int it = 0;it<recParticles.GetEntries();++it){
1849 AliVParticle *part = (AliVParticle*)recParticles.At(it);
edfbe476 1850 Float_t eta = part->Eta();
1851 if(TMath::Abs(eta)<0.9){
1852 Float_t phi = part->Phi();
1853 if(phi<0)phi+=TMath::Pi()*2.;
1854 Float_t dPhi = phi - phiRec;
1855 Float_t dEta = eta - etaRec;
26fa06fb 1856 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1857 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
edfbe476 1858 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1859 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1860 }
c8eabe24 1861
1862 Float_t deltaR = recJets[ir].DeltaR(part);
80a790fd 1863 if(ptRec>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
c8eabe24 1864
1865
80a790fd 1866 if(deltaR<radiusRec&&ptRec>0){
3b7ffecf 1867 Float_t z = part->Pt()/ptRec;
c8eabe24 1868 if(z>zLeading)zLeading=z;
3b7ffecf 1869 Float_t lnz = -1.*TMath::Log(z);
1870 fh2FragRec[ir]->Fill(z,ptRec,eventW);
1871 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1872 }
1873 }
c8eabe24 1874 // fill the jet shapes
1875 Float_t rhoSum = 0;
fb03edbe 1876 for(int ibx = 1;ibx<=fh2RhoPtRec[ir]->GetNbinsX();ibx++){
c8eabe24 1877 Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1878 Float_t rho = fh1TmpRho->GetBinContent(ibx);
1879 rhoSum += rho;
1880 fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1881 fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1882 }
3b7ffecf 1883 }
6bd3fdae 1884
3b7ffecf 1885 // Fill Correlation
1886 Int_t ig = iGenIndex[ir];
1887 if(ig>=0 && ig<nGenJets){
1888 fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1889 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1890 Double_t ptGen = genJets[ig].Pt();
1891 Double_t phiGen = genJets[ig].Phi();
1892 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1893 Double_t etaGen = genJets[ig].Eta();
6bd3fdae 1894
3b7ffecf 1895 container[3] = ptGen;
1896 container[4] = etaGen;
1897 container[5] = phiGen;
c8eabe24 1898 containerPhiZ[3] = ptGen;
3b7ffecf 1899 //
1900 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1901 //
6bd3fdae 1902
a31b8a87 1903 if(TMath::Abs(etaGen)<fJetRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1904 if(TMath::Abs(etaRec)<fJetRecEtaWindow){
3b7ffecf 1905 fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1906 fhnCorrelation->Fill(container);
9a83d4af 1907 if(ptGen>0){
1908 Float_t delta = (ptRec-ptGen)/ptGen;
1909 fh2RelPtFGen->Fill(ptGen,delta,eventW);
1910 }
c8eabe24 1911 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
6bd3fdae 1912 }// if etarec in window
3b7ffecf 1913 }
3b7ffecf 1914 else{
c8eabe24 1915 containerPhiZ[3] = 0;
1916 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
3b7ffecf 1917 }
1918 }// loop over reconstructed jets
1919
1920
1921 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1922 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1923 PostData(1, fHistList);
1924}
1925
1926
1927void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1928 //
1929 // Create the particle container for the correction framework manager and
1930 // link it
1931 //
1932 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
cba109a0 1933 const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
3b7ffecf 1934 const Double_t kEtamin = -3.0, kEtamax = 3.0;
1935 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
c8eabe24 1936 const Double_t kZmin = 0., kZmax = 1;
3b7ffecf 1937
1938 // can we neglect migration in eta and phi?
1939 // phi should be no problem since we cover full phi and are phi symmetric
1940 // eta migration is more difficult due to needed acceptance correction
1941 // in limited eta range
1942
3b7ffecf 1943 //arrays for the number of bins in each dimension
1944 Int_t iBin[kNvar];
cba109a0 1945 iBin[0] = 320; //bins in pt
3b7ffecf 1946 iBin[1] = 1; //bins in eta
1947 iBin[2] = 1; // bins in phi
1948
1949 //arrays for lower bounds :
1950 Double_t* binEdges[kNvar];
1951 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1952 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1953
1954 //values for bin lower bounds
1955 // 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 1956 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 1957 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1958 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1959
1960
1961 for(int i = 0;i < kMaxStep*2;++i){
f51451be 1962 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
3b7ffecf 1963 for (int k=0; k<kNvar; k++) {
1964 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1965 }
1966 }
1967 //create correlation matrix for unfolding
1968 Int_t thnDim[2*kNvar];
1969 for (int k=0; k<kNvar; k++) {
1970 //first half : reconstructed
1971 //second half : MC
1972 thnDim[k] = iBin[k];
1973 thnDim[k+kNvar] = iBin[k];
1974 }
1975
1976 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1977 for (int k=0; k<kNvar; k++) {
1978 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1979 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1980 }
1981 fhnCorrelation->Sumw2();
1982
c8eabe24 1983 // for second correlation histogram
1984
1985
1986 const Int_t kNvarPhiZ = 4;
1987 //arrays for the number of bins in each dimension
1988 Int_t iBinPhiZ[kNvarPhiZ];
1989 iBinPhiZ[0] = 80; //bins in pt
1990 iBinPhiZ[1] = 72; //bins in phi
1991 iBinPhiZ[2] = 20; // bins in Z
1992 iBinPhiZ[3] = 80; //bins in ptgen
1993
1994 //arrays for lower bounds :
1995 Double_t* binEdgesPhiZ[kNvarPhiZ];
1996 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1997 binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1998
1999 for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
2000 for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
2001 for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
2002 for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
2003
2004 fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
2005 for (int k=0; k<kNvarPhiZ; k++) {
2006 fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
2007 }
2008 fhnCorrelationPhiZRec->Sumw2();
2009
2010
3b7ffecf 2011 // Add a histogram for Fake jets
c8eabe24 2012
5cb0f438 2013 for(Int_t ivar = 0; ivar < kNvar; ivar++)
2014 delete [] binEdges[ivar];
2015
c8eabe24 2016 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
2017 delete [] binEdgesPhiZ[ivar];
2018
3b7ffecf 2019}
2020
2021void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
2022{
2023// Terminate analysis
2024//
2025 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
2026}
2027
2028
2029Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
2030
565584e8 2031 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
3b7ffecf 2032
2033 Int_t iCount = 0;
565584e8 2034 if(type==kTrackAOD){
3b7ffecf 2035 AliAODEvent *aod = 0;
565584e8 2036 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
2037 else aod = AODEvent();
3b7ffecf 2038 if(!aod){
2039 return iCount;
2040 }
2041 for(int it = 0;it < aod->GetNumberOfTracks();++it){
2042 AliAODTrack *tr = aod->GetTrack(it);
2043 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
a31b8a87 2044 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
2045 if(tr->Pt()<fMinTrackPt)continue;
38ecb6a5 2046 if(fDebug>0){
2047 if(tr->Pt()>20){
2048 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
fd5db301 2049 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
38ecb6a5 2050 tr->Print();
2051 // tr->Dump();
2052 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
2053 if(fESD){
2054 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
2055 esdTr->Print("");
bb3d13aa 2056 // esdTr->Dump();
38ecb6a5 2057 }
2058 }
2059 }
3b7ffecf 2060 list->Add(tr);
2061 iCount++;
2062 }
2063 }
2064 else if (type == kTrackKineAll||type == kTrackKineCharged){
2065 AliMCEvent* mcEvent = MCEvent();
2066 if(!mcEvent)return iCount;
2067 // we want to have alivpartilces so use get track
2068 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
2069 if(!mcEvent->IsPhysicalPrimary(it))continue;
2070 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
a31b8a87 2071 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 2072 if(type == kTrackKineAll){
2073 list->Add(part);
2074 iCount++;
2075 }
2076 else if(type == kTrackKineCharged){
2077 if(part->Particle()->GetPDG()->Charge()==0)continue;
2078 list->Add(part);
2079 iCount++;
2080 }
2081 }
2082 }
565584e8 2083 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
2084 AliAODEvent *aod = 0;
5301f754 2085 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
565584e8 2086 else aod = AODEvent();
5010a3f7 2087 if(!aod)return iCount;
3b7ffecf 2088 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
2089 if(!tca)return iCount;
2090 for(int it = 0;it < tca->GetEntriesFast();++it){
2091 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
a31b8a87 2092 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 2093 if(!part->IsPhysicalPrimary())continue;
c4631cdb 2094 if(type == kTrackAODMCAll){
3b7ffecf 2095 list->Add(part);
2096 iCount++;
2097 }
565584e8 2098 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
3b7ffecf 2099 if(part->Charge()==0)continue;
565584e8 2100 if(kTrackAODMCCharged){
2101 list->Add(part);
2102 }
2103 else {
a31b8a87 2104 if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
565584e8 2105 list->Add(part);
2106 }
3b7ffecf 2107 iCount++;
2108 }
2109 }
2110 }// AODMCparticle
cc0649e4 2111 list->Sort();
c4631cdb 2112 return iCount;
3b7ffecf 2113
2114}
a31b8a87 2115
2116
2117
2118Bool_t AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
2119 Bool_t selected = false;
2120
2121 if(!jet)return selected;
2122
2123 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
2124 selected = kTRUE;
2125 }
2126 return selected;
2127
2128}
2129
2130Int_t AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
2131
2132 if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
2133 Int_t iCount = 0;
2134
2135 if(!jarray){
2136 Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
2137 return iCount;
2138 }
2139
2140
2141 for(int ij=0;ij<jarray->GetEntries();ij++){
2142 AliAODJet* jet = (AliAODJet*)jarray->At(ij);
2143 if(!jet)continue;
2144 if(type==0){
2145 // No cut at all, main purpose here is sorting
2146 list->Add(jet);
2147 iCount++;
2148 }
2149 else if(type == 1){
2150 // eta cut
2151 if(JetSelected(jet)){
2152 list->Add(jet);
2153 iCount++;
2154 }
2155 }
2156 }
2157
2158 list->Sort();
2159 return iCount;
2160
2161}
2162
2163