]> 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
Rewritten spectrum2 task, old running is optional, various minor updates to AddTask...
[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);
405
3b7ffecf 406 //
407 // Histogram
408
cba109a0 409 const Int_t nBinPt = 320;
3b7ffecf 410 Double_t binLimitsPt[nBinPt+1];
411 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
412 if(iPt == 0){
413 binLimitsPt[iPt] = 0.0;
414 }
415 else {// 1.0
edfbe476 416 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 1.0;
3b7ffecf 417 }
418 }
419
edfbe476 420 const Int_t nBinPhi = 90;
3b7ffecf 421 Double_t binLimitsPhi[nBinPhi+1];
422 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
423 if(iPhi==0){
cc0649e4 424 binLimitsPhi[iPhi] = -1.*TMath::Pi();
3b7ffecf 425 }
426 else{
427 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
428 }
429 }
430
edfbe476 431
c8eabe24 432 const Int_t nBinPhi2 = 360;
433 Double_t binLimitsPhi2[nBinPhi2+1];
434 for(Int_t iPhi2 = 0;iPhi2<=nBinPhi2;iPhi2++){
435 if(iPhi2==0){
436 binLimitsPhi2[iPhi2] = 0.;
437 }
438 else{
439 binLimitsPhi2[iPhi2] = binLimitsPhi2[iPhi2-1] + 1/(Float_t)nBinPhi2 * TMath::Pi()*2;
440 }
441 }
442
edfbe476 443 const Int_t nBinEta = 40;
444 Double_t binLimitsEta[nBinEta+1];
445 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
446 if(iEta==0){
447 binLimitsEta[iEta] = -2.0;
448 }
449 else{
cf049e94 450 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
edfbe476 451 }
452 }
453
3b7ffecf 454 const Int_t nBinFrag = 25;
455
456
a31b8a87 457
458
3b7ffecf 459 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
460 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
a31b8a87 461 fHistList->Add(fh1Xsec);
3b7ffecf 462
463 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
464 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
a31b8a87 465 fHistList->Add(fh1Trials);
3b7ffecf 466
467 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 468 fHistList->Add(fh1PtHard);
469
3b7ffecf 470 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 471 fHistList->Add(fh1PtHardNoW);
472
3b7ffecf 473 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 474 fHistList->Add(fh1PtHardTrials);
3b7ffecf 475
a31b8a87 476
57ca1193 477 fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
a31b8a87 478 fHistList->Add(fh1ZVtx);
479
480 fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
481 fHistList->Add(fh2PtFGen);
57ca1193 482
a31b8a87 483 fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};#Delta p_{T}/p_{T,Gen}",nBinPt,binLimitsPt,120,-2.4,2.4);
484 fHistList->Add(fh2RelPtFGen);
3b7ffecf 485
a31b8a87 486
487 if(fUseOldFill){
488
489 fh1NGenJets = new TH1F("fh1NGenJets","N generated jets",20,-0.5,19.5);
490 fh1NRecJets = new TH1F("fh1NRecJets","N reconstructed jets",20,-0.5,19.5);
491
edfbe476 492 fh1PtTrackRec = new TH1F("fh1PtTrackRec","Rec track P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
493 fh1SumPtTrackRec = new TH1F("fh1SumPtTrackRec","Sum Rec track P_T #eta <0.9;p_{T,sum} (GeV/c)",nBinPt,binLimitsPt);
494 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 495
496 fh1PtJetsRecIn = new TH1F("fh1PtJetsRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
d8f21f85 497 fh1PtJetsGenIn = new TH1F("fh1PtJetsGenIn","Gen jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
565584e8 498 fh1PtJetsLeadingRecIn = new TH1F("fh1PtJetsLeadingRecIn","Rec jets P_T;p_{T} (GeV/c)",nBinPt,binLimitsPt);
cc0649e4 499 fh1PtTracksRecIn = new TH1F("fh1PtTracksRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
565584e8 500 fh1PtTracksLeadingRecIn = new TH1F("fh1PtTracksLeadingRecIn","Rec tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
fe77112a 501 fh1PtTracksGenIn = new TH1F("fh1PtTracksGenIn","gen tracks P_T #eta < 0.9;p_{T} (GeV/c)",nBinPt,binLimitsPt);
edfbe476 502
cc0649e4 503 fh2NRecJetsPt = new TH2F("fh2NRecJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
504 fh2NRecTracksPt = new TH2F("fh2NRecTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
cc0649e4 505
a31b8a87 506 fh2NGenJetsPt = new TH2F("fh2NGenJetsPt","Number of jets above threshhold;p_{T,cut} (GeV/c);N_{jets}",nBinPt,binLimitsPt,50,-0.5,49.5);
507 fh2NGenTracksPt = new TH2F("fh2NGenTracksPt","Number of tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",nBinPt,binLimitsPt,50,-0.5,49.5);
508 //
c8eabe24 509
3b7ffecf 510 for(int ij = 0;ij < kMaxJets;++ij){
511 fh1PtRecIn[ij] = new TH1F(Form("fh1PtRecIn_j%d",ij),"rec p_T input ;p_{T,rec}",nBinPt,binLimitsPt);
512 fh1PtGenIn[ij] = new TH1F(Form("fh1PtGenIn_j%d",ij),"found p_T input ;p_{T,gen}",nBinPt,binLimitsPt);
513
edfbe476 514 fh2PhiPt[ij] = new TH2F(Form("fh2PhiPtRec_j%d",ij),"Jet pt vs delta phi;#Delta#phi;p_{T,jet}",
515 nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
516
cc0649e4 517 fh2PhiEta[ij] = new TH2F(Form("fh2PhiEtaRec_j%d",ij),"delta eta vs delta phi for jets;#Delta#phi;#Delta#eta",
edfbe476 518 nBinPhi,binLimitsPhi,nBinEta,binLimitsEta);
519
c8eabe24 520 fh2RhoPtRec[ij] = new TH2F(Form("fh2RhoPtRec_j%d",ij),"jet shape rho for jets;r;p_{T,rec};",
fb03edbe 521 40,0.,1.,nBinPt,binLimitsPt);
c8eabe24 522 fh2PsiPtRec[ij] = new TH2F(Form("fh2PsiPtRec_j%d",ij),"jet shape psi for jets;r;p_{T,rec};",
fb03edbe 523 40,0.,2.,nBinPt,binLimitsPt);
c8eabe24 524
525 fh2RhoPtGen[ij] = new TH2F(Form("fh2RhoPtGen_j%d",ij),"jet shape rho for jets;r;p_{T,gen};",
fb03edbe 526 40,0.,2.,nBinPt,binLimitsPt);
c8eabe24 527 fh2PsiPtGen[ij] = new TH2F(Form("fh2PsiPtGen_j%d",ij),"jet shape psi for jets;r;p_{T,gen};",
fb03edbe 528 40,0.,2.,nBinPt,binLimitsPt);
529 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
c8eabe24 530
3b7ffecf 531 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",
532 nBinFrag,0.,1.,nBinPt,binLimitsPt);
533 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",
534 nBinFrag,0.,10.,nBinPt,binLimitsPt);
535
536 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 537 nBinFrag,0.,1.,nBinPt,binLimitsPt);
3b7ffecf 538 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 539 nBinFrag,0.,10.,nBinPt,binLimitsPt);
3b7ffecf 540 }
541
26fa06fb 542 // Dijet histograms
c2785065 543 //background histograms
544 if(fBkgSubtraction){
185fa8e6 545 fh1Bkg1 = new TH1F("fh1Bkg1","Background estimate 1",100,0.,10.);
546 fh1Bkg2 = new TH1F("fh1Bkg2","Background estimate 2",100,0.,10.);
7f9012c1 547 fh1Bkg3 = new TH1F("fh1Bkg3","Background estimate 3",100,0.,10.);
c2785065 548 fh1Sigma1 = new TH1F("fh1Sigma1","Background fluctuations 1",100,0.,10.);
549 fh1Sigma2 = new TH1F("fh1Sigma2","Background fluctuations 2",100,0.,10.);
7f9012c1 550 fh1Sigma3 = new TH1F("fh1Sigma3","Background fluctuations 3",100,0.,10.);
185fa8e6 551 fh1Area1 = new TH1F("fh1Area1","Background mean area 1",50,0.,5.);
552 fh1Area2 = new TH1F("fh1Area2","Background mean area 2",50,0.,5.);
7f9012c1 553 fh1Area3 = new TH1F("fh1Area3","Background mean area 3",50,0.,5.);
c2785065 554 fh1Ptjet = new TH1F("fh1Ptjet","Jet spectrum",100,0.,200.);
555 fh1Ptjetsub1 = new TH1F("fh1Ptjetsub1","Subtracted spectrum 1",50,0.,200.);
556 fh1Ptjetsub2 = new TH1F("fh1Ptjetsub2","Subtracted spectrum 2",50,0.,200.);
7f9012c1 557 fh1Ptjetsub3 = new TH1F("fh1Ptjetsub3","Subtracted spectrum 3",50,0.,200.);
c2785065 558 fh1Ptjethardest = new TH1F("fh1Ptjethardest","Hardest jet spectrum",50,0.,200.);
559 fh1Ptjetsubhardest1 = new TH1F("fh1Pthardestsub1","Subtracted hardest jet spectrum 1",100,0.,200.);
560 fh1Ptjetsubhardest2 = new TH1F("fh1Pthardestsub2","Subtracted hardest jet spectrum 2",100,0.,200.);
7f9012c1 561 fh1Ptjetsubhardest3 = new TH1F("fh1Pthardestsub3","Subtracted hardest jet spectrum 3",100,0.,200.);
6bd3fdae 562 fh2Rhovspthardest1 = new TH2F("fh2Rhovspthardest1","Background vs pTjet 1",100,0.,200.,50,0.,5.);
563 fh2Rhovspthardest2 = new TH2F("fh2Rhovspthardest2","Background vs pTjet 2",100,0.,200.,50,0.,5.);
564 fh2Rhovspthardest3 = new TH2F("fh2Rhovspthardest3","Background vs pTjet 3",100,0.,200.,50,0.,5.);
565 fh2Errorvspthardest1 = new TH2F("fh2Errorvspthardest1","Relative error 1",100,0.,200.,50,0.,5.);
566 fh2Errorvspthardest2 = new TH2F("fh2Errorvspthardest2","Relative error 2",100,0.,200.,50,0.,5.);
567 fh2Errorvspthardest3 = new TH2F("fh2Errorvspthardest3","Relative error 3",100,0.,200.,50,0.,5.);
c2785065 568 }
569
26fa06fb 570
3b7ffecf 571 const Int_t saveLevel = 3; // large save level more histos
572 if(saveLevel>0){
a31b8a87 573 if(fBranchGen.Length()>0||fBkgSubtraction){
574 fHistList->Add(fh1NGenJets);
575 fHistList->Add(fh1PtTracksGenIn);
cc0649e4 576 }
577 fHistList->Add(fh1PtJetsRecIn);
d8f21f85 578 fHistList->Add(fh1PtJetsGenIn);
565584e8 579 fHistList->Add(fh1PtJetsLeadingRecIn);
cc0649e4 580 fHistList->Add(fh1PtTracksRecIn);
565584e8 581 fHistList->Add(fh1PtTracksLeadingRecIn);
3b7ffecf 582 fHistList->Add(fh1NRecJets);
edfbe476 583 fHistList->Add(fh1PtTrackRec);
584 fHistList->Add(fh1SumPtTrackRec);
585 fHistList->Add(fh1SumPtTrackAreaRec);
cc0649e4 586 fHistList->Add(fh2NRecJetsPt);
587 fHistList->Add(fh2NRecTracksPt);
a31b8a87 588 fHistList->Add(fh2NGenJetsPt);
589 fHistList->Add(fh2NGenTracksPt);
3b7ffecf 590 for(int i = 0;i<kMaxStep*2;++i)fHistList->Add(fhnJetContainer[i]);
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
701 fh2DijetDeltaPhiPt[ij] = new TH2F(Form("fh2Dijet%sDeltaPhiPt",cAdd.Data()),"Difference in the azimuthal angle;#Delta#phi;p_{T,1};Entries",180,0.,TMath::Pi(),nBinPt,binLimitsPt);
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
953 // find the dijets assume soring and acceptance cut...
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();
963 Float_t dPhi = phi1 - phi0;
964 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
965 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
966 if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
967 ij1 = ij;
968 pT1 = ptJet;
969 }
970 }
971 // fill jet histos for kmax jets
972 if(ij<kMaxJets){
973 Float_t phiJet = jet->Phi();
974 if(phiJet<0)phiJet+=TMath::Pi()*2.;
975 fh1TmpRho->Reset();
976 fh1PtIn[iType][ij]->Fill(ptJet);
977 // fill leading jets...
978 if(ij==0){
979 fh2LeadingJetPtJetPhi[iType]->Fill(ptJet,ptJet);
980 if(ij==0&&iType==0&&fDebug>1){
981 Printf("%d %3.3f %p %s",iType,ptJet,jet,fBranchRec.Data());
982 }
983 }
984 if(particlesList.GetSize()){
985
986 // Particles... correlated with jets...
987 for(int it = 0;it<particlesList.GetEntries();++it){
988 AliVParticle *part = (AliVParticle*)particlesList.At(it);
989 Float_t deltaR = jet->DeltaR(part);
990 if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
991 }
992 // fill the jet shapes
993 Float_t rhoSum = 0;
994 for(int ibx = 1;ibx<=fh2RhoPt[iType][ij]->GetNbinsX();ibx++){
995 Float_t r = fh2RhoPt[iType][ij]->GetXaxis()->GetBinCenter(ibx);
996 Float_t rho = fh1TmpRho->GetBinContent(ibx);
997 rhoSum += rho;
998 fh2RhoPt[iType][ij]->Fill(r,ptJet,rho);
999 fh2PsiPt[iType][ij]->Fill(r,ptJet,rhoSum);
1000 }
1001 }// if we have particles
1002 }// ij < kMaxJets
1003 }// Jet Loop
1004
1005
1006 // do something with dijets...
1007 if(ij0>=0&&ij1>0){
1008 AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1009 Double_t ptJet0 = jet0->Pt();
1010 Double_t phiJet0 = jet0->Phi();
1011 if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1012
1013 AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1014 Double_t ptJet1 = jet1->Pt();
1015 Double_t phiJet1 = jet1->Phi();
1016 if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1017
1018 Float_t deltaPhi = phiJet0 - phiJet1;
1019 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1020 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1021 deltaPhi = TMath::Abs(deltaPhi);
1022 fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet0);
1023
1024 Float_t asym = 9999;
1025 if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1026 fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1027 fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1028 fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1029 Float_t minv = 2.*(jet0->P()*jet1->P()-
1030 jet0->Px()*jet1->Px()-
1031 jet0->Py()*jet1->Py()-
1032 jet0->Pz()*jet1->Pz()); // assume mass == 0;
1033 if(minv<0)minv=0; // prevent numerical instabilities
1034 minv = TMath::Sqrt(minv);
1035 fh1DijetMinv[iType]->Fill(minv);
1036 }
1037
1038
1039
1040 // count down the jets above thrueshold
1041 Int_t nOver = nJets;
1042 if(nOver>0){
1043 TIterator *jetIter = jetsList.MakeIterator();
1044 AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
1045 Float_t pt = tmpJet->Pt();
1046 if(tmpJet){
1047 for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1048 Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1049 while(pt<ptCut&&tmpJet){
1050 nOver--;
1051 tmpJet = (AliAODJet*)(jetIter->Next());
1052 if(tmpJet){
1053 pt = tmpJet->Pt();
1054 }
1055 }
1056 if(nOver<=0)break;
1057 fh2NJetsPt[iType]->Fill(ptCut,nOver);
1058 }
1059 }
1060 delete jetIter;
1061 }
1062}
1063
1064void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1065
1066 Int_t nTrackOver = particlesList.GetSize();
1067 // do the same for tracks and jets
1068 if(nTrackOver>0){
1069 TIterator *trackIter = particlesList.MakeIterator();
1070 AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1071 Float_t pt = tmpTrack->Pt();
1072 for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1073 Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1074 while(pt<ptCut&&tmpTrack){
1075 nTrackOver--;
1076 tmpTrack = (AliVParticle*)(trackIter->Next());
1077 if(tmpTrack){
1078 pt = tmpTrack->Pt();
1079 }
1080 }
1081 if(nTrackOver<=0)break;
1082 fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1083 }
1084
1085
1086 trackIter->Reset();
1087 AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1088 Float_t sumPt = 0;
1089
1090 while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1091 Float_t tmpPt = tmpTrack->Pt();
1092 fh1PtTracksIn[iType]->Fill(tmpPt);
1093 sumPt += tmpPt;
1094 Float_t tmpPhi = tmpTrack->Phi();
1095 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1096 if(tmpTrack==leading){
1097 fh1PtTracksLeadingIn[iType]->Fill(tmpPt);
1098 fh2LeadingTrackPtTrackPhi[iType]->Fill(tmpPt,tmpPhi);
1099 continue;
1100 }
1101 }
1102 fh1SumPtTrack[iType]->Fill(sumPt);
1103 delete trackIter;
1104 }
1105
1106}
1107
1108
1109void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1110
1111
1112 // Fill al the matching histos
1113 // It is important that the acceptances for the mathing are as large as possible
1114 // to avoid false matches on the edge of acceptance
1115 // therefore we add some extra matching jets as overhead
1116
1117 static TArrayI aGenIndex(recJetsList.GetEntries());
1118 if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
3b7ffecf 1119
a31b8a87 1120 static TArrayI aRecIndex(genJetsList.GetEntries());
1121 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
3b7ffecf 1122
a31b8a87 1123 if(fDebug){
1124 Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1125 Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1126 }
1127 AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1128 &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1129 aGenIndex,aRecIndex,fDebug);
1130
1131 if(fDebug){
1132 for(int i = 0;i< aGenIndex.GetSize();++i){
1133 if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1134 }
1135 for(int i = 0;i< aRecIndex.GetSize();++i){
1136 if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1137 }
1138 }
1139
1140 Double_t container[6];
1141 Double_t containerPhiZ[6];
1142
1143 // loop over generated jets
1144 // consider the
1145 for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1146 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1147 Double_t ptGen = genJet->Pt();
1148 Double_t phiGen = genJet->Phi();
1149 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1150 Double_t etaGen = genJet->Eta();
1151 container[3] = ptGen;
1152 container[4] = etaGen;
1153 container[5] = phiGen;
1154 fhnJetContainer[kStep0]->Fill(&container[3]);
1155 if(JetSelected(genJet)){
1156 fhnJetContainer[kStep1]->Fill(&container[3]);
1157 Int_t ir = aRecIndex[ig];
1158 if(ir>=0&&ir<recJetsList.GetEntries()){
1159 fhnJetContainer[kStep2]->Fill(&container[3]);
1160 AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1161 if(JetSelected(recJet))fhnJetContainer[kStep4]->Fill(&container[3]);
1162 if(JetSelected(recJet))fhnJetContainer[kStep3]->Fill(&container[3]);
1163 }
1164 }
1165 }// loop over generated jets used for matching...
1166
1167
1168
1169 // loop over reconstructed jets
1170 for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1171 AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1172 Double_t etaRec = recJet->Eta();
1173 Double_t ptRec = recJet->Pt();
1174 Double_t phiRec = recJet->Phi();
1175 if(phiRec<0)phiRec+=TMath::Pi()*2.;
1176 // do something with dijets...
1177
1178 container[0] = ptRec;
1179 container[1] = etaRec;
1180 container[2] = phiRec;
1181 containerPhiZ[0] = ptRec;
1182 containerPhiZ[1] = phiRec;
1183
1184 fhnJetContainer[kStep0+kMaxStep]->Fill(container);
1185 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1186
1187 if(JetSelected(recJet)){
1188 fhnJetContainer[kStep1+kMaxStep]->Fill(container);
1189 // Fill Correlation
1190 Int_t ig = aGenIndex[ir];
1191 if(ig>=0 && ig<genJetsList.GetEntries()){
1192 fhnJetContainer[kStep2+kMaxStep]->Fill(container);
1193 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1194 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1195 Double_t ptGen = genJet->Pt();
1196 Double_t phiGen = genJet->Phi();
1197 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1198 Double_t etaGen = genJet->Eta();
1199
1200 container[3] = ptGen;
1201 container[4] = etaGen;
1202 container[5] = phiGen;
1203 containerPhiZ[3] = ptGen;
1204 //
1205 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1206 //
1207 if(JetSelected(genJet))fhnJetContainer[kStep4+kMaxStep]->Fill(container);
1208 fhnJetContainer[kStep3+kMaxStep]->Fill(container);
1209 fhnCorrelation->Fill(container);
1210 if(ptGen>0){
1211 Float_t delta = (ptRec-ptGen)/ptGen;
1212 fh2RelPtFGen->Fill(ptGen,delta);
1213 fh2PtFGen->Fill(ptGen,ptRec);
1214 }
1215 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1216 }
1217 else{
1218 containerPhiZ[3] = 0;
1219 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
1220 }
1221 }// loop over reconstructed jets
1222 }
1223 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1224}
1225
1226
1227void AliAnalysisTaskJetSpectrum2::UserExecOld(Option_t */*option*/)
1228{
1229
1230
1231 Bool_t selected = kTRUE;
1232
1233
1234
1235
1236 if(fUseGlobalSelection&&fEventSelectionMask==0){
1237 selected = AliAnalysisHelperJetTasks::Selected();
1238 }
1239 else if(fUseGlobalSelection&&fEventSelectionMask>0){
1240 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
1241 }
1242
1243 if(fEventClass>0){
1244 selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
1245 }
1246
1247 if(!selected){
1248 // no selection by the service task, we continue
1249 if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
1250 PostData(1, fHistList);
1251 return;
1252 }
1253
1254
1255 //
1256 // Execute analysis for current event
1257 //
1258 static AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1259
1260
1261 // take all other information from the aod we take the tracks from
1262 static AliAODEvent *aod = 0;
1263 if(!aod){
1264 if(fUseAODTrackInput)aod = fAODIn;
1265 aod = fAODOut;
1266 }
1267 static AliAODEvent *fAOD = 0;
1268 if(!fAOD){
1269 if(fUseAODJetInput)fAOD = fAODIn;
1270 fAOD = fAODOut;
1271 }
3b7ffecf 1272
1273 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
1274
1275
1276 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
1277
1278 if(!aodH){
1279 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
1280 return;
1281 }
1282
1283 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
a31b8a87 1284 TClonesArray *aodRecJets = 0;
1285
1286
1287
1288 if(fAOD&&!aodRecJets){
1289 aodRecJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchRec.Data()));
1290 }
1291 if(fAODExtension&&!aodRecJets){
1292 aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
1293 }
1294
3b7ffecf 1295 if(!aodRecJets){
1296 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
1297 return;
1298 }
1299
6bd3fdae 1300 Int_t nJets = aodRecJets->GetEntriesFast();
1301
a31b8a87 1302 TClonesArray *aodGenJets = 0;
1303 if(fBranchGen.Length()>0){
1304 if(fAOD&&!aodGenJets){
1305 aodGenJets = dynamic_cast<TClonesArray*>(fAOD->FindListObject(fBranchGen.Data()));
1306 }
1307 if(fAODExtension&&!aodGenJets){
1308 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
1309 }
1310
1311 if(!aodGenJets){
1312 Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
1313 return;
1314 }
1315 }
1316
1317
6bd3fdae 1318 // ==== General variables needed
1319 // We use statice array, not to fragment the memory
1320 AliAODJet genJets[kMaxJets];
1321 Int_t nGenJets = 0;
1322 AliAODJet recJets[kMaxJets];
1323 Int_t nRecJets = 0;
a31b8a87 1324
1325 // These will replace the arrays...
1326 TList genJetsList; // full acceptance
1327 TList genJetsListCut; // acceptance cut
1328 TList recJetsList; // full acceptance
1329 TList recJetsListCut; // acceptance cut
6bd3fdae 1330
1331
a31b8a87 1332 GetListOfJets(&genJetsList,aodGenJets,0);
1333 GetListOfJets(&genJetsListCut,aodGenJets,1);
6bd3fdae 1334
a31b8a87 1335 if(fBranchGen.Length()>0){
1336 GetListOfJets(&recJetsList,aodRecJets,0);
1337 GetListOfJets(&recJetsListCut,aodRecJets,1);
1338 }
1339 ///////////////////////////
6bd3fdae 1340
a31b8a87 1341
c2785065 1342
1343
a31b8a87 1344 if(fBkgSubtraction){
1345 AliAODJetEventBackground* evBkg=(AliAODJetEventBackground*)fAOD->FindListObject(fBranchBkg.Data());
c2785065 1346 if(!evBkg){
1347 Printf("%s:%d no reconstructed background array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkg.Data());
a31b8a87 1348 return;
c2785065 1349 }
3b7ffecf 1350
1351
c2785065 1352 ///just to start: some very simple plots containing rho, sigma and area of the background.
1353
6bd3fdae 1354
7f9012c1 1355 Float_t pthardest=0.;
80a790fd 1356
c2785065 1357 if(nJets!=0){
7f9012c1 1358
c2785065 1359 Float_t bkg1=evBkg->GetBackground(0);
1360 Float_t bkg2=evBkg->GetBackground(1);
7f9012c1 1361 Float_t bkg3=evBkg->GetBackground(2);
c2785065 1362 Float_t sigma1=evBkg->GetSigma(0);
1363 Float_t sigma2=evBkg->GetSigma(1);
7f9012c1 1364 Float_t sigma3=evBkg->GetSigma(2);
c2785065 1365 Float_t area1=evBkg->GetMeanarea(0);
1366 Float_t area2=evBkg->GetMeanarea(1);
7f9012c1 1367 Float_t area3=evBkg->GetMeanarea(2);
185fa8e6 1368 fh1Bkg1->Fill(bkg1); //rho computed with all background jets.
7f9012c1 1369 fh1Bkg2->Fill(bkg2); //rho computed with all background jets but the 2 hardest.
1370 fh1Bkg3->Fill(bkg3); //rho computed with randomized jets
c2785065 1371 fh1Sigma1->Fill(sigma1);
1372 fh1Sigma2->Fill(sigma2);
7f9012c1 1373 fh1Sigma3->Fill(sigma3);
185fa8e6 1374 fh1Area1->Fill(area1);
1375 fh1Area2->Fill(area2);
7f9012c1 1376 fh1Area3->Fill(area3);
c2785065 1377
80a790fd 1378 Int_t iSubJetCounter = 0;
c2785065 1379 for(Int_t k=0;k<nJets;k++){
1380 AliAODJet *jet = dynamic_cast<AliAODJet*>(aodRecJets->At(k));
1381 fh1Ptjet->Fill(jet->Pt());
1382 Float_t ptsub1=jet->Pt()-bkg1*jet->EffectiveAreaCharged();
1383 Float_t ptsub2=jet->Pt()-bkg2*jet->EffectiveAreaCharged();
7f9012c1 1384 Float_t ptsub3=jet->Pt()-bkg3*jet->EffectiveAreaCharged();
6bd3fdae 1385 if(ptsub2<0.) ptsub2=0.;
1386 if(ptsub1<0.) ptsub1=0.;
1387 if(ptsub3<0.) ptsub3=0.;
c2785065 1388 Float_t err1=sigma1*sqrt(area1);
1389 Float_t err2=sigma2*sqrt(area2);
7f9012c1 1390 Float_t err3=sigma3*sqrt(area3);
c2785065 1391 fh1Ptjetsub1->Fill(ptsub1);
1392 fh1Ptjetsub2->Fill(ptsub2);
7f9012c1 1393 fh1Ptjetsub3->Fill(ptsub3);
c2785065 1394 if(k==0) {
1395 pthardest=jet->Pt();
1396 fh1Ptjethardest->Fill(pthardest);
1397 fh1Ptjetsubhardest1->Fill(ptsub1);
1398 fh1Ptjetsubhardest2->Fill(ptsub2);
7f9012c1 1399 fh1Ptjetsubhardest3->Fill(ptsub3);
36c36a0c 1400 if (ptsub1 > 0)
1401 fh2Errorvspthardest1->Fill(ptsub1,err1/ptsub1);
1402 if (ptsub2 > 0)
1403 fh2Errorvspthardest2->Fill(ptsub2,err2/ptsub2);
1404 if (ptsub3 > 0)
1405 fh2Errorvspthardest3->Fill(ptsub3,err3/ptsub3);
c2785065 1406 }
7f9012c1 1407
f4132e7d 1408 Float_t ptsub=0.;
1409 if(fFillCorrBkg==1) ptsub=ptsub1;
1410 if(fFillCorrBkg==2) ptsub=ptsub2;
1411 if(fFillCorrBkg==3) ptsub=ptsub3;
80a790fd 1412 if(ptsub>0){// avoid unphysical jets pT
1413 Float_t subphi=jet->Phi();
1414 Float_t subtheta=jet->Theta();
1415 Float_t subpz = ptsub/TMath::Tan(subtheta);
1416 Float_t subpx=ptsub*TMath::Cos(subphi);
1417 Float_t subpy=ptsub * TMath::Sin(subphi);
1418 Float_t subp = TMath::Sqrt(ptsub*ptsub+subpz*subpz);
1419 if(k<kMaxJets){// only store the jets which are assoicated to anohter one
1420 genJets[iSubJetCounter].SetPxPyPzE(subpx,subpy,subpz,subp);
1421 iSubJetCounter++;
1422 }
f4132e7d 1423 }
6bd3fdae 1424 }
80a790fd 1425 nGenJets = iSubJetCounter;
6bd3fdae 1426 fh2Rhovspthardest1->Fill(pthardest,bkg1);
1427 fh2Rhovspthardest2->Fill(pthardest,bkg2);
1428 fh2Rhovspthardest3->Fill(pthardest,bkg3);
c2785065 1429 }
1430 }// background subtraction
1431
1432
3b7ffecf 1433 Double_t eventW = 1;
1434 Double_t ptHard = 0;
1435 Double_t nTrials = 1; // Trials for MC trigger
1436
1437 if(fUseExternalWeightOnly){
1438 eventW = fExternalWeight;
1439 }
1440
1441 fh1Trials->Fill("#sum{ntrials}",fAvgTrials);
cf049e94 1442 // if(fDebug>0)aodH->SetFillAOD(kFALSE);
3b7ffecf 1443 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1444 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
1445 // this is the part we only use when we have MC information
1446 AliMCEvent* mcEvent = MCEvent();
1447 // AliStack *pStack = 0;
1448 if(!mcEvent){
1449 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
1450 return;
1451 }
1452 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
1453 Int_t iCount = 0;
1454 if(pythiaGenHeader){
1455 nTrials = pythiaGenHeader->Trials();
1456 ptHard = pythiaGenHeader->GetPtHard();
1457 int iProcessType = pythiaGenHeader->ProcessType();
1458 // 11 f+f -> f+f
1459 // 12 f+barf -> f+barf
1460 // 13 f+barf -> g+g
1461 // 28 f+g -> f+g
1462 // 53 g+g -> f+barf
1463 // 68 g+g -> g+g
1464 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
1465 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
1466
1467 // fetch the pythia generated jets only to be used here
1468 Int_t nPythiaGenJets = pythiaGenHeader->NTriggerJets();
1469 AliAODJet pythiaGenJets[kMaxJets];
1470 for(int ip = 0;ip < nPythiaGenJets;++ip){
1471 if(iCount>=kMaxJets)continue;
1472 Float_t p[4];
1473 pythiaGenHeader->TriggerJet(ip,p);
1474 pythiaGenJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
1475
3b7ffecf 1476 if(fBranchGen.Length()==0){
9280dfa6 1477 /*
1478 if(fLimitGenJetEta){
1479 if(pythiaGenJets[iCount].Eta()>fJetHeaderRec->GetJetEtaMax()||
1480 pythiaGenJets[iCount].Eta()<fJetHeaderRec->GetJetEtaMin())continue;
1481 }
a31b8a87 1482 */
9280dfa6 1483 if(fMinJetPt>0&&pythiaGenJets[iCount].Pt()<fMinJetPt)continue;
3b7ffecf 1484 // if we have MC particles and we do not read from the aod branch
1485 // use the pythia jets
f4132e7d 1486 if(!fBkgSubtraction){
1487 genJets[iCount].SetPxPyPzE(p[0],p[1],p[2],p[3]);
1488 iCount++;
1489 }
3b7ffecf 1490 }
3b7ffecf 1491 }
1492 }
f4132e7d 1493 if(fBranchGen.Length()==0&&!fBkgSubtraction)nGenJets = iCount;
3b7ffecf 1494 }// (fAnalysisType&kMCESD)==kMCESD)
1495
1496
1497 // we simply fetch the tracks/mc particles as a list of AliVParticles
1498
1499 TList recParticles;
1500 TList genParticles;
1501
565584e8 1502
1503
1504
1505 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
1506 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
1507 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
1508 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
3b7ffecf 1509
cc0649e4 1510
3b7ffecf 1511 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1512 fh1PtHard->Fill(ptHard,eventW);
1513 fh1PtHardNoW->Fill(ptHard,1);
1514 fh1PtHardTrials->Fill(ptHard,nTrials);
57ca1193 1515 fh1ZVtx->Fill(fAOD->GetPrimaryVertex()->GetZ());
1516
3b7ffecf 1517
1518 // If we set a second branch for the input jets fetch this
f4132e7d 1519 if(fBranchGen.Length()>0&&!fBkgSubtraction){
3b7ffecf 1520 if(aodGenJets){
1521 Int_t iCount = 0;
1522 for(int ig = 0;ig < aodGenJets->GetEntries();++ig){
1523 if(iCount>=kMaxJets)continue;
1524 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodGenJets->At(ig));
1525 if(!tmp)continue;
5301f754 1526 /*
3b7ffecf 1527 if(fLimitGenJetEta){
1528 if(tmp->Eta()>fJetHeaderRec->GetJetEtaMax()||
1529 tmp->Eta()<fJetHeaderRec->GetJetEtaMin())continue;
1530 }
5301f754 1531 */
9280dfa6 1532 if(fMinJetPt>0&&tmp->Pt()<fMinJetPt)continue;
3b7ffecf 1533 genJets[iCount] = *tmp;
1534 iCount++;
1535 }
1536 nGenJets = iCount;
1537 }
1538 else{
5cb0f438 1539 if(fDebug>1)Printf("%s:%d Generated jet branch %s not found",(char*)__FILE__,__LINE__,fBranchGen.Data());
1540 if(fDebug>2)fAOD->Print();
3b7ffecf 1541 }
1542 }
1543
1544 fh1NGenJets->Fill(nGenJets);
1545 // We do not want to exceed the maximum jet number
1546 nGenJets = TMath::Min(nGenJets,kMaxJets);
1547
1548 // Fetch the reconstructed jets...
1549
1550 nRecJets = aodRecJets->GetEntries();
1551
cc0649e4 1552 nRecJets = aodRecJets->GetEntries();
1553 fh1NRecJets->Fill(nRecJets);
1554
1555 // Do something with the all rec jets
1556 Int_t nRecOver = nRecJets;
3b7ffecf 1557
cc0649e4 1558 // check that the jets are sorted
1559 Float_t ptOld = 999999;
1560 for(int ir = 0;ir < nRecJets;ir++){
1561 AliAODJet *tmp = (AliAODJet*)(aodRecJets->At(ir));
1562 Float_t tmpPt = tmp->Pt();
1563 if(tmpPt>ptOld){
3568c3d6 1564 Printf("%s:%d Jets Not Sorted %s !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,fBranchRec.Data(),ir,tmpPt,ir-1,ptOld);
cc0649e4 1565 }
1566 ptOld = tmpPt;
1567 }
3b7ffecf 1568
cc0649e4 1569
1570 if(nRecOver>0){
1571 TIterator *recIter = aodRecJets->MakeIterator();
1572 AliAODJet *tmpRec = (AliAODJet*)(recIter->Next());
1573 Float_t pt = tmpRec->Pt();
1574 if(tmpRec){
1575 for(int i = 1;i <= fh2NRecJetsPt->GetNbinsX();i++){
1576 Float_t ptCut = fh2NRecJetsPt->GetXaxis()->GetBinCenter(i);
1577 while(pt<ptCut&&tmpRec){
1578 nRecOver--;
1579 tmpRec = (AliAODJet*)(recIter->Next());
1580 if(tmpRec){
1581 pt = tmpRec->Pt();
1582 }
1583 }
1584 if(nRecOver<=0)break;
1585 fh2NRecJetsPt->Fill(ptCut,nRecOver);
1586 }
1587 }
1588 recIter->Reset();
1589 AliAODJet *leading = (AliAODJet*)aodRecJets->At(0);
1590 Float_t phi = leading->Phi();
1591 if(phi<0)phi+=TMath::Pi()*2.;
cf049e94 1592 pt = leading->Pt();
45a11aef 1593 while((tmpRec = (AliAODJet*)(recIter->Next()))){
cc0649e4 1594 Float_t tmpPt = tmpRec->Pt();
1595 fh1PtJetsRecIn->Fill(tmpPt);
565584e8 1596 if(tmpRec==leading){
1597 fh1PtJetsLeadingRecIn->Fill(tmpPt);
1598 continue;
1599 }
cc0649e4 1600 // correlation
1601 Float_t tmpPhi = tmpRec->Phi();
1602
1603 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1604 Float_t dPhi = phi - tmpRec->Phi();
1605 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1606 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
a31b8a87 1607 // Float_t dEta = eta - tmpRec->Eta();
1608 // fh2JetsLeadingPhiEta->Fill(dPhi,dEta);
1609 // fh2JetsLeadingPhiPt->Fill(dPhi,pt);
cc0649e4 1610 }
1611 delete recIter;
1612 }
1613
1614 Int_t nTrackOver = recParticles.GetSize();
1615 // do the same for tracks and jets
1616 if(nTrackOver>0){
1617 TIterator *recIter = recParticles.MakeIterator();
1618 AliVParticle *tmpRec = (AliVParticle*)(recIter->Next());
1619 Float_t pt = tmpRec->Pt();
1620 // Printf("Leading track p_t %3.3E",pt);
1621 for(int i = 1;i <= fh2NRecTracksPt->GetNbinsX();i++){
1622 Float_t ptCut = fh2NRecTracksPt->GetXaxis()->GetBinCenter(i);
1623 while(pt<ptCut&&tmpRec){
1624 nTrackOver--;
1625 tmpRec = (AliVParticle*)(recIter->Next());
1626 if(tmpRec){
1627 pt = tmpRec->Pt();
1628 }
1629 }
1630 if(nTrackOver<=0)break;
1631 fh2NRecTracksPt->Fill(ptCut,nTrackOver);
1632 }
1633
1634 recIter->Reset();
1635 AliVParticle *leading = (AliVParticle*)recParticles.At(0);
1636 Float_t phi = leading->Phi();
1637 if(phi<0)phi+=TMath::Pi()*2.;
cf049e94 1638 pt = leading->Pt();
cc0649e4 1639 while((tmpRec = (AliVParticle*)(recIter->Next()))){
1640 Float_t tmpPt = tmpRec->Pt();
1641 fh1PtTracksRecIn->Fill(tmpPt);
565584e8 1642 if(tmpRec==leading){
1643 fh1PtTracksLeadingRecIn->Fill(tmpPt);
1644 continue;
1645 }
cc0649e4 1646 // correlation
1647 Float_t tmpPhi = tmpRec->Phi();
1648
1649 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1650 Float_t dPhi = phi - tmpRec->Phi();
1651 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1652 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
a31b8a87 1653 // Float_t dEta = eta - tmpRec->Eta();
1654 // fh2TracksLeadingPhiEta->Fill(dPhi,dEta);
1655 // fh2TracksLeadingPhiPt->Fill(dPhi,pt);
cc0649e4 1656 }
1657 delete recIter;
1658 }
1659
1660 if(genParticles.GetSize()){
1661 TIterator *genIter = genParticles.MakeIterator();
1662 AliVParticle *tmpGen = 0;
1663 while((tmpGen = (AliVParticle*)(genIter->Next()))){
1664 if(TMath::Abs(tmpGen->Eta())<0.9){
1665 Float_t tmpPt = tmpGen->Pt();
1666 fh1PtTracksGenIn->Fill(tmpPt);
1667 }
1668 }
1669 delete genIter;
1670 }
1671
3b7ffecf 1672 nRecJets = TMath::Min(nRecJets,kMaxJets);
6bd3fdae 1673 nJets = TMath::Min(nJets,kMaxJets);
9280dfa6 1674 Int_t iCountRec = 0;
3b7ffecf 1675 for(int ir = 0;ir < nRecJets;++ir){
1676 AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
1677 if(!tmp)continue;
9280dfa6 1678 if(tmp->Pt()<fMinJetPt)continue;
e6a92928 1679 recJets[iCountRec] = *tmp;
9280dfa6 1680 iCountRec++;
3b7ffecf 1681 }
9280dfa6 1682 nRecJets = iCountRec;
3b7ffecf 1683
1684
1685 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1686 // Relate the jets
1687 Int_t iGenIndex[kMaxJets]; // Index of the generated jet for i-th rec -1 if none
1688 Int_t iRecIndex[kMaxJets]; // Index of the rec jet for i-th gen -1 if none
1689
cc0649e4 1690
1691 for(int i = 0;i<kMaxJets;++i){
3b7ffecf 1692 iGenIndex[i] = iRecIndex[i] = -1;
1693 }
1694
6bd3fdae 1695
3b7ffecf 1696 AliAnalysisHelperJetTasks::GetClosestJets(genJets,nGenJets,recJets,nRecJets,
1697 iGenIndex,iRecIndex,fDebug);
6bd3fdae 1698
a31b8a87 1699 static TArrayI aGenIndex(recJetsListCut.GetEntries());
1700 if(aGenIndex.GetSize()<recJetsListCut.GetEntries())aGenIndex.Set(recJetsListCut.GetEntries());
1701
1702 static TArrayI aRecIndex(genJetsList.GetEntries());
1703 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1704
1705 if(fDebug){
1706 Printf("New Gens List %d rec index Array %d count %d",genJetsList.GetEntries(),aRecIndex.GetSize(),nGenJets);
1707 Printf("New Rec List %d gen indey Array %d count %d",recJetsListCut.GetEntries(),aGenIndex.GetSize(),nRecJets);
1708 }
1709 AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,nGenJets,&recJetsListCut,nRecJets,
1710 aGenIndex,aRecIndex,fDebug);
1711
6bd3fdae 1712
3b7ffecf 1713 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1714
1715 if(fDebug){
a31b8a87 1716 for(int i = 0;i<kMaxJets;++i){
3b7ffecf 1717 if(iGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,iGenIndex[i]);
1718 if(iRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,iRecIndex[i]);
1719 }
1720 }
1721
1722
1723
1724
1725 Double_t container[6];
c8eabe24 1726 Double_t containerPhiZ[6];
3b7ffecf 1727
1728 // loop over generated jets
1729
1730 // radius; tmp, get from the jet header itself
1731 // at some point, how todeal woht FastJet on MC?
1732 Float_t radiusGen =0.4;
1733 Float_t radiusRec =0.4;
1734
1735 for(int ig = 0;ig < nGenJets;++ig){
1736 Double_t ptGen = genJets[ig].Pt();
1737 Double_t phiGen = genJets[ig].Phi();
1738 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1739 Double_t etaGen = genJets[ig].Eta();
d8f21f85 1740 fh1PtJetsGenIn->Fill(ptGen);
1741
3b7ffecf 1742 container[3] = ptGen;
1743 container[4] = etaGen;
1744 container[5] = phiGen;
1745 fhnJetContainer[kStep0]->Fill(&container[3],eventW);
1746 Int_t ir = iRecIndex[ig];
1747
a31b8a87 1748 if(TMath::Abs(etaGen)<fJetRecEtaWindow){
c8eabe24 1749 fh1TmpRho->Reset();
1750
3b7ffecf 1751 fhnJetContainer[kStep1]->Fill(&container[3],eventW);
1752 fh1PtGenIn[ig]->Fill(ptGen,eventW);
1753 // fill the fragmentation function
1754 for(int it = 0;it<genParticles.GetEntries();++it){
1755 AliVParticle *part = (AliVParticle*)genParticles.At(it);
c8eabe24 1756 Float_t deltaR = genJets[ig].DeltaR(part);
80a790fd 1757 if(ptGen>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptGen);
1758 if(deltaR<radiusGen&&ptGen>0){
3b7ffecf 1759 Float_t z = part->Pt()/ptGen;
1760 Float_t lnz = -1.*TMath::Log(z);
1761 fh2FragGen[ig]->Fill(z,ptGen,eventW);
1762 fh2FragLnGen[ig]->Fill(lnz,ptGen,eventW);
1763 }
c8eabe24 1764
3b7ffecf 1765 }
c8eabe24 1766 Float_t rhoSum = 0;
fb03edbe 1767 for(int ibx = 1;ibx<=fh2RhoPtGen[ir]->GetNbinsX();ibx++){
c8eabe24 1768 Float_t r = fh2RhoPtGen[ir]->GetXaxis()->GetBinCenter(ibx);
1769 Float_t rho = fh1TmpRho->GetBinContent(ibx);
1770 rhoSum += rho;
1771 fh2RhoPtGen[ig]->Fill(r,ptGen,rho);
1772 fh2PsiPtGen[ig]->Fill(r,ptGen,rhoSum);
3b7ffecf 1773 }
c8eabe24 1774 }
3b7ffecf 1775 if(ir>=0&&ir<nRecJets){
1776 fhnJetContainer[kStep2]->Fill(&container[3],eventW);
1777 Double_t etaRec = recJets[ir].Eta();
a31b8a87 1778 if(TMath::Abs(etaRec)<fJetRecEtaWindow)fhnJetContainer[kStep4]->Fill(&container[3],eventW);
1779 if(TMath::Abs(etaRec)<fJetRecEtaWindow&&TMath::Abs(etaGen)<fJetRecEtaWindow)fhnJetContainer[kStep3]->Fill(&container[3],eventW);
3b7ffecf 1780 }
1781 }// loop over generated jets
7fa8b2da 1782
edfbe476 1783 Float_t sumPt = 0;
1784 for(int it = 0;it<recParticles.GetEntries();++it){
1785 AliVParticle *part = (AliVParticle*)recParticles.At(it);
1786 // fill sum pt and P_t of all paritles
1787 if(TMath::Abs(part->Eta())<0.9){
1788 Float_t pt = part->Pt();
1789 fh1PtTrackRec->Fill(pt,eventW);
a31b8a87 1790 // fh2TrackPtTrackPhi->Fill(pt,part->Phi());
edfbe476 1791 sumPt += pt;
1792 }
1793 }
1794 fh1SumPtTrackAreaRec->Fill(sumPt*0.4*0.4/(2.*1.8),eventW);
1795 fh1SumPtTrackRec->Fill(sumPt,eventW);
1796
cf049e94 1797
3b7ffecf 1798 // loop over reconstructed jets
1799 for(int ir = 0;ir < nRecJets;++ir){
1800 Double_t etaRec = recJets[ir].Eta();
1801 Double_t ptRec = recJets[ir].Pt();
1802 Double_t phiRec = recJets[ir].Phi();
1803 if(phiRec<0)phiRec+=TMath::Pi()*2.;
26fa06fb 1804
1805
1806 // do something with dijets...
26fa06fb 1807
3b7ffecf 1808 container[0] = ptRec;
1809 container[1] = etaRec;
1810 container[2] = phiRec;
c8eabe24 1811 containerPhiZ[0] = ptRec;
1812 containerPhiZ[1] = phiRec;
ffab794c 1813 if(ptRec>30.&&fDebug>2){
7fa8b2da 1814 // need to cast to int, otherwise the printf overwrites
1815 Printf("Jet found in Event %d with p_T, %E",(int)Entry(),ptRec);
1277f815 1816 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(int)fInputHandler->GetTree()->GetTree()->GetReadEntry());
1817 if(fESD)Printf("ESDEvent GetEventNumberInFile(): %d",fESD->GetEventNumberInFile());
cf049e94 1818 // aodH->SetFillAOD(kTRUE);
7fa8b2da 1819 fAOD->GetHeader()->Print();
a221560c 1820 Printf("TriggerClasses: %s",fAOD->GetFiredTriggerClasses().Data());
7fa8b2da 1821 for(int it = 0;it < fAOD->GetNumberOfTracks();++it){
1822 AliAODTrack *tr = fAOD->GetTrack(it);
1823 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1824 tr->Print();
38ecb6a5 1825 // tr->Dump();
7fa8b2da 1826 if(fESD){
1827 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
1828 esdTr->Print("");
bb3d13aa 1829 // esdTr->Dump();
7fa8b2da 1830 }
1831 }
1832 }
1833
1834
3b7ffecf 1835 fhnJetContainer[kStep0+kMaxStep]->Fill(container,eventW);
1836 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
45a11aef 1837
c8eabe24 1838 Float_t zLeading = -1;
a31b8a87 1839 if(TMath::Abs(etaRec)<fJetRecEtaWindow){
1840 // fh2JetPtJetPhi->Fill(ptRec,phiRec);
3b7ffecf 1841 fhnJetContainer[kStep1+kMaxStep]->Fill(container,eventW);
1842 fh1PtRecIn[ir]->Fill(ptRec,eventW);
1843 // fill the fragmentation function
c8eabe24 1844
1845 fh1TmpRho->Reset();
1846
3b7ffecf 1847 for(int it = 0;it<recParticles.GetEntries();++it){
1848 AliVParticle *part = (AliVParticle*)recParticles.At(it);
edfbe476 1849 Float_t eta = part->Eta();
1850 if(TMath::Abs(eta)<0.9){
1851 Float_t phi = part->Phi();
1852 if(phi<0)phi+=TMath::Pi()*2.;
1853 Float_t dPhi = phi - phiRec;
1854 Float_t dEta = eta - etaRec;
26fa06fb 1855 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1856 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
edfbe476 1857 fh2PhiPt[ir]->Fill(dPhi,ptRec,eventW);
1858 fh2PhiEta[ir]->Fill(dPhi,dEta,eventW);
1859 }
c8eabe24 1860
1861 Float_t deltaR = recJets[ir].DeltaR(part);
80a790fd 1862 if(ptRec>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptRec);
c8eabe24 1863
1864
80a790fd 1865 if(deltaR<radiusRec&&ptRec>0){
3b7ffecf 1866 Float_t z = part->Pt()/ptRec;
c8eabe24 1867 if(z>zLeading)zLeading=z;
3b7ffecf 1868 Float_t lnz = -1.*TMath::Log(z);
1869 fh2FragRec[ir]->Fill(z,ptRec,eventW);
1870 fh2FragLnRec[ir]->Fill(lnz,ptRec,eventW);
1871 }
1872 }
c8eabe24 1873 // fill the jet shapes
1874 Float_t rhoSum = 0;
fb03edbe 1875 for(int ibx = 1;ibx<=fh2RhoPtRec[ir]->GetNbinsX();ibx++){
c8eabe24 1876 Float_t r = fh2RhoPtRec[ir]->GetXaxis()->GetBinCenter(ibx);
1877 Float_t rho = fh1TmpRho->GetBinContent(ibx);
1878 rhoSum += rho;
1879 fh2RhoPtRec[ir]->Fill(r,ptRec,rho);
1880 fh2PsiPtRec[ir]->Fill(r,ptRec,rhoSum);
1881 }
3b7ffecf 1882 }
6bd3fdae 1883
3b7ffecf 1884 // Fill Correlation
1885 Int_t ig = iGenIndex[ir];
1886 if(ig>=0 && ig<nGenJets){
1887 fhnJetContainer[kStep2+kMaxStep]->Fill(container,eventW);
1888 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1889 Double_t ptGen = genJets[ig].Pt();
1890 Double_t phiGen = genJets[ig].Phi();
1891 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1892 Double_t etaGen = genJets[ig].Eta();
6bd3fdae 1893
3b7ffecf 1894 container[3] = ptGen;
1895 container[4] = etaGen;
1896 container[5] = phiGen;
c8eabe24 1897 containerPhiZ[3] = ptGen;
3b7ffecf 1898 //
1899 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1900 //
6bd3fdae 1901
a31b8a87 1902 if(TMath::Abs(etaGen)<fJetRecEtaWindow)fhnJetContainer[kStep4+kMaxStep]->Fill(container,eventW);
1903 if(TMath::Abs(etaRec)<fJetRecEtaWindow){
3b7ffecf 1904 fhnJetContainer[kStep3+kMaxStep]->Fill(container,eventW);
1905 fhnCorrelation->Fill(container);
9a83d4af 1906 if(ptGen>0){
1907 Float_t delta = (ptRec-ptGen)/ptGen;
1908 fh2RelPtFGen->Fill(ptGen,delta,eventW);
1909 }
c8eabe24 1910 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
6bd3fdae 1911 }// if etarec in window
3b7ffecf 1912 }
3b7ffecf 1913 else{
c8eabe24 1914 containerPhiZ[3] = 0;
1915 if(fhnCorrelationPhiZRec)fhnCorrelationPhiZRec->Fill(containerPhiZ);
3b7ffecf 1916 }
1917 }// loop over reconstructed jets
1918
1919
1920 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1921 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1922 PostData(1, fHistList);
1923}
1924
1925
1926void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1927 //
1928 // Create the particle container for the correction framework manager and
1929 // link it
1930 //
1931 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
cba109a0 1932 const Double_t kPtmin = 0.0, kPtmax = 320.; // we do not want to have empty bins at the beginning...
3b7ffecf 1933 const Double_t kEtamin = -3.0, kEtamax = 3.0;
1934 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
c8eabe24 1935 const Double_t kZmin = 0., kZmax = 1;
3b7ffecf 1936
1937 // can we neglect migration in eta and phi?
1938 // phi should be no problem since we cover full phi and are phi symmetric
1939 // eta migration is more difficult due to needed acceptance correction
1940 // in limited eta range
1941
3b7ffecf 1942 //arrays for the number of bins in each dimension
1943 Int_t iBin[kNvar];
cba109a0 1944 iBin[0] = 320; //bins in pt
3b7ffecf 1945 iBin[1] = 1; //bins in eta
1946 iBin[2] = 1; // bins in phi
1947
1948 //arrays for lower bounds :
1949 Double_t* binEdges[kNvar];
1950 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1951 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1952
1953 //values for bin lower bounds
1954 // 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 1955 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 1956 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1957 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1958
1959
1960 for(int i = 0;i < kMaxStep*2;++i){
f51451be 1961 fhnJetContainer[i] = new THnSparseF(Form("fhnJetContainer%d",i),Form("THnSparse jet info %d",i),kNvar,iBin);
3b7ffecf 1962 for (int k=0; k<kNvar; k++) {
1963 fhnJetContainer[i]->SetBinEdges(k,binEdges[k]);
1964 }
1965 }
1966 //create correlation matrix for unfolding
1967 Int_t thnDim[2*kNvar];
1968 for (int k=0; k<kNvar; k++) {
1969 //first half : reconstructed
1970 //second half : MC
1971 thnDim[k] = iBin[k];
1972 thnDim[k+kNvar] = iBin[k];
1973 }
1974
1975 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
1976 for (int k=0; k<kNvar; k++) {
1977 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1978 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1979 }
1980 fhnCorrelation->Sumw2();
1981
c8eabe24 1982 // for second correlation histogram
1983
1984
1985 const Int_t kNvarPhiZ = 4;
1986 //arrays for the number of bins in each dimension
1987 Int_t iBinPhiZ[kNvarPhiZ];
1988 iBinPhiZ[0] = 80; //bins in pt
1989 iBinPhiZ[1] = 72; //bins in phi
1990 iBinPhiZ[2] = 20; // bins in Z
1991 iBinPhiZ[3] = 80; //bins in ptgen
1992
1993 //arrays for lower bounds :
1994 Double_t* binEdgesPhiZ[kNvarPhiZ];
1995 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
1996 binEdgesPhiZ[ivar] = new Double_t[iBinPhiZ[ivar] + 1];
1997
1998 for(Int_t i=0; i<=iBinPhiZ[0]; i++) binEdgesPhiZ[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[0]*(Double_t)i;
1999 for(Int_t i=0; i<=iBinPhiZ[1]; i++) binEdgesPhiZ[1][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBinPhiZ[1]*(Double_t)i;
2000 for(Int_t i=0; i<=iBinPhiZ[2]; i++) binEdgesPhiZ[2][i]=(Double_t)kZmin + (kZmax-kZmin)/iBinPhiZ[2]*(Double_t)i;
2001 for(Int_t i=0; i<=iBinPhiZ[3]; i++) binEdgesPhiZ[3][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBinPhiZ[3]*(Double_t)i;
2002
2003 fhnCorrelationPhiZRec = new THnSparseF("fhnCorrelationPhiZRec","THnSparse with correlations",kNvarPhiZ,iBinPhiZ);
2004 for (int k=0; k<kNvarPhiZ; k++) {
2005 fhnCorrelationPhiZRec->SetBinEdges(k,binEdgesPhiZ[k]);
2006 }
2007 fhnCorrelationPhiZRec->Sumw2();
2008
2009
3b7ffecf 2010 // Add a histogram for Fake jets
c8eabe24 2011
5cb0f438 2012 for(Int_t ivar = 0; ivar < kNvar; ivar++)
2013 delete [] binEdges[ivar];
2014
c8eabe24 2015 for(Int_t ivar = 0; ivar < kNvarPhiZ; ivar++)
2016 delete [] binEdgesPhiZ[ivar];
2017
3b7ffecf 2018}
2019
2020void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
2021{
2022// Terminate analysis
2023//
2024 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
2025}
2026
2027
2028Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
2029
565584e8 2030 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
3b7ffecf 2031
2032 Int_t iCount = 0;
565584e8 2033 if(type==kTrackAOD){
3b7ffecf 2034 AliAODEvent *aod = 0;
565584e8 2035 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
2036 else aod = AODEvent();
3b7ffecf 2037 if(!aod){
2038 return iCount;
2039 }
2040 for(int it = 0;it < aod->GetNumberOfTracks();++it){
2041 AliAODTrack *tr = aod->GetTrack(it);
2042 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
a31b8a87 2043 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
2044 if(tr->Pt()<fMinTrackPt)continue;
38ecb6a5 2045 if(fDebug>0){
2046 if(tr->Pt()>20){
2047 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
fd5db301 2048 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
38ecb6a5 2049 tr->Print();
2050 // tr->Dump();
2051 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
2052 if(fESD){
2053 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(tr->GetID());
2054 esdTr->Print("");
bb3d13aa 2055 // esdTr->Dump();
38ecb6a5 2056 }
2057 }
2058 }
3b7ffecf 2059 list->Add(tr);
2060 iCount++;
2061 }
2062 }
2063 else if (type == kTrackKineAll||type == kTrackKineCharged){
2064 AliMCEvent* mcEvent = MCEvent();
2065 if(!mcEvent)return iCount;
2066 // we want to have alivpartilces so use get track
2067 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
2068 if(!mcEvent->IsPhysicalPrimary(it))continue;
2069 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
a31b8a87 2070 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 2071 if(type == kTrackKineAll){
2072 list->Add(part);
2073 iCount++;
2074 }
2075 else if(type == kTrackKineCharged){
2076 if(part->Particle()->GetPDG()->Charge()==0)continue;
2077 list->Add(part);
2078 iCount++;
2079 }
2080 }
2081 }
565584e8 2082 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
2083 AliAODEvent *aod = 0;
5301f754 2084 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
565584e8 2085 else aod = AODEvent();
5010a3f7 2086 if(!aod)return iCount;
3b7ffecf 2087 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
2088 if(!tca)return iCount;
2089 for(int it = 0;it < tca->GetEntriesFast();++it){
2090 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
a31b8a87 2091 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 2092 if(!part->IsPhysicalPrimary())continue;
c4631cdb 2093 if(type == kTrackAODMCAll){
3b7ffecf 2094 list->Add(part);
2095 iCount++;
2096 }
565584e8 2097 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
3b7ffecf 2098 if(part->Charge()==0)continue;
565584e8 2099 if(kTrackAODMCCharged){
2100 list->Add(part);
2101 }
2102 else {
a31b8a87 2103 if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
565584e8 2104 list->Add(part);
2105 }
3b7ffecf 2106 iCount++;
2107 }
2108 }
2109 }// AODMCparticle
cc0649e4 2110 list->Sort();
c4631cdb 2111 return iCount;
3b7ffecf 2112
2113}
a31b8a87 2114
2115
2116
2117Bool_t AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
2118 Bool_t selected = false;
2119
2120 if(!jet)return selected;
2121
2122 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
2123 selected = kTRUE;
2124 }
2125 return selected;
2126
2127}
2128
2129Int_t AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
2130
2131 if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
2132 Int_t iCount = 0;
2133
2134 if(!jarray){
2135 Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
2136 return iCount;
2137 }
2138
2139
2140 for(int ij=0;ij<jarray->GetEntries();ij++){
2141 AliAODJet* jet = (AliAODJet*)jarray->At(ij);
2142 if(!jet)continue;
2143 if(type==0){
2144 // No cut at all, main purpose here is sorting
2145 list->Add(jet);
2146 iCount++;
2147 }
2148 else if(type == 1){
2149 // eta cut
2150 if(JetSelected(jet)){
2151 list->Add(jet);
2152 iCount++;
2153 }
2154 }
2155 }
2156
2157 list->Sort();
2158 return iCount;
2159
2160}
2161
2162