]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGJE/AliAnalysisTaskJetSpectrum2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
A prototype macro for the TDR7 layout
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetSpectrum2.cxx
... / ...
CommitLineData
1// **************************************
2// used for the correction of determiantion of reconstructed jet spectra
3// Compares input (gen) and output (rec) jets
4// *******************************************
5
6
7/**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9 * *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
12 * *
13 * Permission to use, copy, modify and distribute this software and its *
14 * documentation strictly for non-commercial purposes is hereby granted *
15 * without fee, provided that the above copyright notice appears in all *
16 * copies and that both the copyright notice and this permission notice *
17 * appear in the supporting documentation. The authors make no claims *
18 * about the suitability of this software for any purpose. It is *
19 * provided "as is" without express or implied warranty. *
20 **************************************************************************/
21
22
23#include <TROOT.h>
24#include <TRandom.h>
25#include <TRandom3.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 <TProfile2D.h>
36#include <TList.h>
37#include <TLorentzVector.h>
38#include <TClonesArray.h>
39#include <TRefArray.h>
40#include "TDatabasePDG.h"
41
42#include "AliAnalysisTaskJetSpectrum2.h"
43#include "AliAnalysisManager.h"
44#include "AliJetFinder.h"
45#include "AliTHn.h"
46#include "AliJetHeader.h"
47#include "AliJetReader.h"
48#include "AliJetReaderHeader.h"
49#include "AliUA1JetHeaderV1.h"
50#include "AliESDEvent.h"
51#include "AliAODEvent.h"
52#include "AliAODHandler.h"
53#include "AliAODTrack.h"
54#include "AliAODJet.h"
55#include "AliAODJetEventBackground.h"
56#include "AliAODMCParticle.h"
57#include "AliMCEventHandler.h"
58#include "AliMCEvent.h"
59#include "AliStack.h"
60#include "AliGenPythiaEventHeader.h"
61#include "AliJetKineReaderHeader.h"
62#include "AliGenCocktailEventHeader.h"
63#include "AliInputEventHandler.h"
64#include "AliCFContainer.h"
65
66#include "AliAnalysisHelperJetTasks.h"
67
68ClassImp(AliAnalysisTaskJetSpectrum2)
69
70AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2():
71AliAnalysisTaskSE(),
72 fJetHeaderRec(0x0),
73 fJetHeaderGen(0x0),
74 fAODIn(0x0),
75 fAODOut(0x0),
76 fAODExtension(0x0),
77 fhnJetContainer(0x0),
78 fhnCorrelation(0x0),
79 fhnEvent(0x0),
80 f1PtScale(0x0),
81 fBranchRec("jets"),
82 fBranchGen(""),
83 fBranchBkgRec(""),
84 fBranchBkgGen(""),
85 fNonStdFile(""),
86 fRandomizer(0x0),
87 fUseAODJetInput(kFALSE),
88 fUseAODTrackInput(kFALSE),
89 fUseAODMCInput(kFALSE),
90 fUseGlobalSelection(kFALSE),
91 fUseExternalWeightOnly(kFALSE),
92 fLimitGenJetEta(kFALSE),
93 fDoMatching(kFALSE),
94 fNMatchJets(5),
95 fNRPBins(3),
96 fTRP(0),
97 fDebug(0),
98 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
99 fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
100 fFilterMask(0),
101 fEventSelectionMask(0),
102 fNTrigger(0),
103 fTriggerBit(0x0),
104 fNAcceptance(0),
105 fNBinsLeadingTrackPt(10),
106 fNBinsMult(20),
107 fAnalysisType(0),
108 fTrackTypeRec(kTrackUndef),
109 fTrackTypeGen(kTrackUndef),
110 fEventClass(0),
111 fRPMethod(0),
112 fAvgTrials(1),
113 fExternalWeight(1),
114 fJetRecEtaWindow(0.5),
115 fTrackRecEtaWindow(0.5),
116 fMinJetPt(0),
117 fMinTrackPt(0.15),
118 fDeltaPhiWindow(90./180.*TMath::Pi()),
119 fAcceptancePhiMin(0x0),
120 fAcceptancePhiMax(0x0),
121 fAcceptanceEtaMin(0x0),
122 fAcceptanceEtaMax(0x0),
123 fCentrality(100),
124 fRPAngle(0),
125 fMultRec(0),
126 fMultGen(0),
127 fTriggerName(0x0),
128 fh1Xsec(0x0),
129 fh1Trials(0x0),
130 fh1AvgTrials(0x0),
131 fh1PtHard(0x0),
132 fh1PtHardNoW(0x0),
133 fh1PtHardTrials(0x0),
134 fh1ZVtx(0x0),
135 fh1RP(0x0),
136 fh1Centrality(0x0),
137 fh1TmpRho(0x0),
138 fh2MultRec(0x0),
139 fh2MultGen(0x0),
140 fh2RPCentrality(0x0),
141 fh2PtFGen(0x0),
142 fh2deltaPt1Pt2(0x0),
143 fh2RelPtFGen(0x0),
144 fh3RelPtFGenLeadTrkPt(0x0),
145 fHistList(0x0)
146{
147
148 for(int ij = 0;ij <kJetTypes;++ij){
149 fFlagJetType[ij] = 1; // default = on
150 fh1NJets[ij] = 0;
151 fh1SumPtTrack[ij] = 0;
152 fh1PtJetsIn[ij] = 0;
153 fh1PtJetsInRej[ij] = 0;
154 fh1PtJetsInBest[ij] = 0;
155 fh1PtTracksIn[ij] = 0;
156 fh1PtTracksInLow[ij] = 0;
157 fh2NJetsPt[ij] = 0;
158 fh2NTracksPt[ij] = 0;
159 fp2MultRPPhiTrackPt[ij] = 0;
160 fp2CentRPPhiTrackPt[ij] = 0;
161 fhnJetPt[ij] = 0;
162 fhnJetPtBest[ij] = 0;
163 fhnJetPtRej[ij] = 0;
164 fhnJetPtQA[ij] = 0;
165 fhnTrackPt[ij] = 0;
166 fhnTrackPtQA[ij] = 0;
167 for(int i = 0;i <= kMaxJets;++i){
168 fh2LTrackPtJetPt[ij][i] = 0;
169 fh1PtIn[ij][i] = 0;
170 }
171
172 fh1DijetMinv[ij] = 0;
173 fh2DijetDeltaPhiPt[ij] = 0;
174 fh2DijetAsymPt[ij] = 0;
175 fh2DijetPt2vsPt1[ij] = 0;
176 fh2DijetDifvsSum[ij] = 0;
177
178 }
179}
180
181AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
182 AliAnalysisTaskSE(name),
183 fJetHeaderRec(0x0),
184 fJetHeaderGen(0x0),
185 fAODIn(0x0),
186 fAODOut(0x0),
187 fAODExtension(0x0),
188 fhnJetContainer(0x0),
189 fhnCorrelation(0x0),
190 fhnEvent(0x0),
191 f1PtScale(0x0),
192 fBranchRec("jets"),
193 fBranchGen(""),
194 fBranchBkgRec(""),
195 fBranchBkgGen(""),
196 fNonStdFile(""),
197 fRandomizer(0x0),
198 fUseAODJetInput(kFALSE),
199 fUseAODTrackInput(kFALSE),
200 fUseAODMCInput(kFALSE),
201 fUseGlobalSelection(kFALSE),
202 fUseExternalWeightOnly(kFALSE),
203 fLimitGenJetEta(kFALSE),
204 fDoMatching(kFALSE),
205 fNMatchJets(5),
206 fNRPBins(3),
207 fTRP(0),
208 fDebug(0),
209 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
210 fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
211 fFilterMask(0),
212 fEventSelectionMask(0),
213 fNTrigger(0),
214 fTriggerBit(0x0),
215 fNAcceptance(0),
216 fNBinsLeadingTrackPt(10),
217 fNBinsMult(20),
218 fAnalysisType(0),
219 fTrackTypeRec(kTrackUndef),
220 fTrackTypeGen(kTrackUndef),
221 fEventClass(0),
222 fRPMethod(0),
223 fAvgTrials(1),
224 fExternalWeight(1),
225 fJetRecEtaWindow(0.5),
226 fTrackRecEtaWindow(0.5),
227 fMinJetPt(0),
228 fMinTrackPt(0.15),
229 fDeltaPhiWindow(90./180.*TMath::Pi()),
230 fAcceptancePhiMin(0x0),
231 fAcceptancePhiMax(0x0),
232 fAcceptanceEtaMin(0x0),
233 fAcceptanceEtaMax(0x0),
234 fCentrality(100),
235 fRPAngle(0),
236 fMultRec(0),
237 fMultGen(0),
238 fTriggerName(0x0),
239 fh1Xsec(0x0),
240 fh1Trials(0x0),
241 fh1AvgTrials(0x0),
242 fh1PtHard(0x0),
243 fh1PtHardNoW(0x0),
244 fh1PtHardTrials(0x0),
245 fh1ZVtx(0x0),
246 fh1RP(0x0),
247 fh1Centrality(0x0),
248 fh1TmpRho(0x0),
249 fh2MultRec(0x0),
250 fh2MultGen(0x0),
251 fh2RPCentrality(0x0),
252 fh2PtFGen(0x0),
253 fh2deltaPt1Pt2(0x0),
254 fh2RelPtFGen(0x0),
255 fh3RelPtFGenLeadTrkPt(0x0),
256 fHistList(0x0)
257{
258
259 for(int ij = 0;ij <kJetTypes;++ij){
260 fFlagJetType[ij] = 1; // default = on
261 fh1NJets[ij] = 0;
262 fh1SumPtTrack[ij] = 0;
263 fh1PtJetsIn[ij] = 0;
264 fh1PtJetsInRej[ij] = 0;
265 fh1PtJetsInBest[ij] = 0;
266 fh1PtTracksIn[ij] = 0;
267 fh1PtTracksInLow[ij] = 0;
268 fh2NJetsPt[ij] = 0;
269 fh2NTracksPt[ij] = 0;
270 fp2MultRPPhiTrackPt[ij] = 0;
271 fp2CentRPPhiTrackPt[ij] = 0;
272 fhnJetPt[ij] = 0;
273 fhnJetPtBest[ij] = 0;
274 fhnJetPtRej[ij] = 0;
275 fhnJetPtQA[ij] = 0;
276 fhnTrackPt[ij] = 0;
277 fhnTrackPtQA[ij] = 0;
278 for(int i = 0;i <= kMaxJets;++i){
279 fh2LTrackPtJetPt[ij][i] = 0;
280 fh1PtIn[ij][i] = 0;
281 }
282
283 fh1DijetMinv[ij] = 0;
284 fh2DijetDeltaPhiPt[ij] = 0;
285 fh2DijetAsymPt[ij] = 0;
286 fh2DijetPt2vsPt1[ij] = 0;
287 fh2DijetDifvsSum[ij] = 0;
288 }
289
290 DefineOutput(1, TList::Class());
291}
292
293
294
295Bool_t AliAnalysisTaskJetSpectrum2::Notify()
296{
297
298
299
300 //
301 // Implemented Notify() to read the cross sections
302 // and number of trials from pyxsec.root
303 //
304
305 // Fetch the aod also from the input in,
306 // have todo it in notify
307
308
309 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
310 // assume that the AOD is in the general output...
311 fAODOut = AODEvent();
312
313 if(fNonStdFile.Length()!=0){
314 // case that we have an AOD extension we need can fetch the jets from the extended output
315 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
316 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
317 if(!fAODExtension){
318 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
319 }
320 }
321
322
323 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
324 Float_t xsection = 0;
325 Float_t ftrials = 1;
326
327 fAvgTrials = 1;
328 if(tree){
329 TFile *curfile = tree->GetCurrentFile();
330 if (!curfile) {
331 Error("Notify","No current file");
332 return kFALSE;
333 }
334 if(!fh1Xsec||!fh1Trials){
335 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
336 return kFALSE;
337 }
338 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
339 fh1Xsec->Fill("<#sigma>",xsection);
340 // construct a poor man average trials
341 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
342 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
343 fh1Trials->Fill("#sum{ntrials}",ftrials);
344 }
345
346 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
347
348 return kTRUE;
349}
350
351void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
352{
353
354
355 // Connect the AOD
356
357 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
358 OpenFile(1);
359 if(!fHistList)fHistList = new TList();
360 PostData(1, fHistList); // post data in any case once
361
362 if(!fRandomizer)fRandomizer = new TRandom3(0);
363
364 fHistList->SetOwner(kTRUE);
365 Bool_t oldStatus = TH1::AddDirectoryStatus();
366 TH1::AddDirectory(kFALSE);
367
368
369
370 // event npsparse cent, mult
371 const Int_t nBinsSparse0 = 4;
372 const Int_t nBins0[nBinsSparse0] = { 100, 500,fNTrigger,125};
373 const Double_t xmin0[nBinsSparse0] = { 0, 0, -0.5,-2};
374 const Double_t xmax0[nBinsSparse0] = { 100,5000,fNTrigger-0.5,248};
375
376
377 fhnEvent = new THnSparseF("fhnEvent",";cent;mult:trigger;#rho",nBinsSparse0,
378 nBins0,xmin0,xmax0);
379 fHistList->Add(fhnEvent);
380
381 if(fDoMatching){
382 MakeJetContainer();
383 fHistList->Add(fhnCorrelation);
384 fHistList->Add(fhnJetContainer);
385 }
386 //
387 // Histogram
388
389
390
391 const Int_t nBinPt = 150;
392 Double_t binLimitsPt[nBinPt+1];
393 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
394 if(iPt == 0){
395 binLimitsPt[iPt] = -50.0;
396 }
397 else {// 1.0
398 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
399 }
400 }
401 const Int_t nBinPhi = 90;
402 Double_t binLimitsPhi[nBinPhi+1];
403 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
404 if(iPhi==0){
405 binLimitsPhi[iPhi] = 0;
406 }
407 else{
408 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
409 }
410 }
411
412 const Int_t nBinEta = 40;
413 Double_t binLimitsEta[nBinEta+1];
414 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
415 if(iEta==0){
416 binLimitsEta[iEta] = -2.0;
417 }
418 else{
419 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
420 }
421 }
422
423
424 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
425 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
426 fHistList->Add(fh1Xsec);
427 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
428 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
429 fHistList->Add(fh1Trials);
430 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
431 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
432 fHistList->Add(fh1AvgTrials);
433 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
434 fHistList->Add(fh1PtHard);
435 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
436 fHistList->Add(fh1PtHardNoW);
437 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
438 fHistList->Add(fh1PtHardTrials);
439
440 fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
441 fHistList->Add(fh1ZVtx);
442
443
444 fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
445 fHistList->Add(fh1RP);
446
447 fh1Centrality = new TH1F("fh1Centrality","cent;cent (%)",103,-1,102);
448 fHistList->Add(fh1Centrality);
449
450 fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",500,0,5000,500,0.,5000);
451 fHistList->Add(fh2MultRec);
452 fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,0,5000,500,0.,5000);
453 fHistList->Add(fh2MultGen);
454
455
456 fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
457 fHistList->Add(fh2RPCentrality);
458
459 fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
460 fHistList->Add(fh2PtFGen);
461
462 // fh2deltaPt1Pt2(0x0),
463 fh2deltaPt1Pt2 = new TH3F("fh2deltaPt1Pt2",Form("%s vs. %s;delta;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
464 fHistList->Add(fh2deltaPt1Pt2);
465
466 fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
467 fHistList->Add(fh2RelPtFGen);
468
469 const Int_t nBinsLeadingTrackPt = 10;
470 const Double_t binArrayLeadingTrackPt[nBinsLeadingTrackPt+1] = {0.,1.,2.,3.,4.,5.,6.,8.,10.,12.,200.}; //store pT of leading track in jet
471
472 const Int_t nBinDeltaPt = 241;
473 Double_t binLimitsDeltaPt[nBinDeltaPt+1];
474 for(Int_t iDeltaPt = 0;iDeltaPt <= nBinDeltaPt;iDeltaPt++){
475 if(iDeltaPt == 0){
476 binLimitsDeltaPt[iDeltaPt] = -2.41;
477 }
478 else {
479 binLimitsDeltaPt[iDeltaPt] = binLimitsDeltaPt[iDeltaPt-1] + 0.02;
480 }
481 }
482
483 fh3RelPtFGenLeadTrkPt = new TH3F("fh3RelPtFGenLeadTrkPt",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen};p_{T}^{leading track}",nBinPt,binLimitsPt,nBinDeltaPt,binLimitsDeltaPt,nBinsLeadingTrackPt,binArrayLeadingTrackPt);
484 fHistList->Add(fh3RelPtFGenLeadTrkPt);
485
486 for(int ij = 0;ij <kJetTypes;++ij){
487 TString cAdd = "";
488 TString cJetBranch = "";
489 if(ij==kJetRec){
490 cAdd = "Rec";
491 cJetBranch = fBranchRec.Data();
492 }
493 else if (ij==kJetGen){
494 cAdd = "Gen";
495 cJetBranch = fBranchGen.Data();
496 }
497 else if (ij==kJetRecFull){
498 cAdd = "RecFull";
499 cJetBranch = fBranchRec.Data();
500 }
501 else if (ij==kJetGenFull){
502 cAdd = "GenFull";
503 cJetBranch = fBranchGen.Data();
504 }
505
506 if(cJetBranch.Length()==0)fFlagJetType[ij] = 0;
507 if(!fFlagJetType[ij])continue;
508
509 fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
510 fHistList->Add(fh1NJets[ij]);
511
512 fh1PtJetsIn[ij] = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
513 fHistList->Add(fh1PtJetsIn[ij]);
514
515 fh1PtJetsInRej[ij] = new TH1F(Form("fh1PtJets%sInRej",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
516 fHistList->Add(fh1PtJetsInRej[ij]);
517
518 fh1PtJetsInBest[ij] = new TH1F(Form("fh1PtJets%sInBest",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
519 fHistList->Add(fh1PtJetsInBest[ij]);
520
521 fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
522 fHistList->Add(fh1PtTracksIn[ij]);
523
524 fh1PtTracksInLow[ij] = new TH1F(Form("fh1PtTracks%sInLow",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),100,0.,5.);
525 fHistList->Add(fh1PtTracksInLow[ij]);
526
527 fh1SumPtTrack[ij] = new TH1F(Form("fh1SumPtTrack%s",cAdd.Data()),Form("Sum %s track p_T;p_{T} (GeV/c)",cAdd.Data()),1000,0.,3000.);
528 fHistList->Add(fh1SumPtTrack[ij]);
529
530 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);
531 fHistList->Add(fh2NJetsPt[ij]);
532
533 fh2NTracksPt[ij] = new TH2F(Form("fh2N%sTracksPt",cAdd.Data()),Form("Number of %s tracks above threshhold;p_{T,cut} (GeV/c);N_{tracks}",cAdd.Data()),nBinPt,binLimitsPt,1000,0.,5000);
534 fHistList->Add(fh2NTracksPt[ij]);
535
536
537 fp2MultRPPhiTrackPt[ij] = new TProfile2D(Form("fp2MultRPPhiTrackPt%s",cAdd.Data()),"RP phi vs Number of tracks;# tracks;#Delta#phi_{RP}; <p_{T}>",20,0,4000,181,-1./180.*TMath::Pi(),TMath::Pi(),"S");
538 fHistList->Add(fp2MultRPPhiTrackPt[ij]);
539 fp2CentRPPhiTrackPt[ij] = new TProfile2D(Form("fp2CentRPPhiTrackPt%s",cAdd.Data()),"RP phi vs cent;# cent;#Delta#phi_{RP}; <p_{T}>",10,0,100,181,-1./180.*TMath::Pi(),TMath::Pi(),"S");
540 fHistList->Add(fp2CentRPPhiTrackPt[ij]);
541
542 // Bins: Jet number: pTJet, cent, mult, RP, Area, trigger, leading track pT total bins = 4.5M
543 const Int_t nBinsSparse1 = 9;
544 const Int_t nBinsArea = 10;
545 Int_t nbinr5=120;
546 Int_t nbinr5min=-50;
547
548 if(cJetBranch.Contains("KT05")){
549 nbinr5=132;
550 nbinr5min=-80;}
551
552 Int_t nBins1[nBinsSparse1] = { kMaxJets+1,nbinr5, 10, fNBinsMult, fNRPBins, nBinsArea,fNTrigger,fNBinsLeadingTrackPt,fNAcceptance+1};
553 if(cJetBranch.Contains("RandomCone")){
554 nBins1[1] = 600;
555 nBins1[5] = 1;
556 }
557 const Double_t xmin1[nBinsSparse1] = { -0.5,nbinr5min, 0, 0, -0.5, 0., -0.5, 0., -0.5,};
558 const Double_t xmax1[nBinsSparse1] = {kMaxJets+0.5,250,100,4000,fNRPBins-0.5,1.0,fNTrigger-0.5,200.,fNAcceptance+0.5};
559
560 const Double_t binArrayArea[nBinsArea+1] = {xmin1[5],0.07,0.1,0.15,0.2,0.25,0.3,0.4,0.5,0.6,xmax1[5]};
561
562
563 fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leading track p_{T};acceptance bin",nBinsSparse1,nBins1,xmin1,xmax1);
564 fhnJetPt[ij]->SetBinEdges(5,binArrayArea);
565 if(fNBinsLeadingTrackPt>1) fhnJetPt[ij]->SetBinEdges(7,binArrayLeadingTrackPt);
566 fHistList->Add(fhnJetPt[ij]);
567
568
569 // Bins: pTJet, cent, trigger,
570 const Int_t nBinsSparse1b = 3;
571 Int_t nBins1b[nBinsSparse1b] = {120, 10,fNTrigger};
572 const Double_t xmin1b[nBinsSparse1b] = {-50, 0,-0.5};
573 const Double_t xmax1b[nBinsSparse1b] = {250,100,fNTrigger-0.5};
574
575 fhnJetPtBest[ij] = new THnSparseF(Form("fhnJetPtBest%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b);
576 fHistList->Add(fhnJetPtBest[ij]);
577
578 fhnJetPtRej[ij] = new THnSparseF(Form("fhnJetPtRej%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b);
579 fHistList->Add(fhnJetPtRej[ij]);
580
581 // Bins: Jet number: pTJet, cent, eta, phi, Area, trigger, acceptance, signed pT leading
582 const Int_t nBinsSparse2 = 9;
583 Int_t nBins2[nBinsSparse2] = { kMaxJets+1, 60, 8, 6, 90, 100,fNTrigger,fNAcceptance+1,20};
584 if(cJetBranch.Contains("RandomCone")){
585 nBins2[5] = 1;
586 }
587 const Double_t xmin2[nBinsSparse2] = { -0.5, -50, 0,-0.9, 0, 0., -0.5, -0.5, -100};
588 const Double_t xmax2[nBinsSparse2] = {kMaxJets+0.5, 250, 80, 0.9, 2.*TMath::Pi(),1.0,fNTrigger-0.5,fNAcceptance+0.5, 100};
589 fhnJetPtQA[ij] = new THnSparseF(Form("fhnJetPtQA%s",cAdd.Data()),";jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin;signed pt leading",nBinsSparse2,nBins2,xmin2,xmax2);
590 // fhnJetPt[ij]->SetBinEdges(5,binArrayArea);
591 fHistList->Add(fhnJetPtQA[ij]);
592
593 // Bins:track number pTtrack, cent, mult, RP. charge total bins = 224 k
594 const Int_t nBinsSparse3 = 7;
595 const Int_t nRPBinsSparse3 = 3;
596 const Int_t nBins3[nBinsSparse3] = { 2, 100, 10, 1, nRPBinsSparse3,fNTrigger,2};
597 const Double_t xmin3[nBinsSparse3] = { -0.5, 0, 0, 0, -0.5,-0.5,-1.5};
598 const Double_t xmax3[nBinsSparse3] = { 1.5, 200, 100, 4000,nRPBinsSparse3-0.5,fNTrigger-0.5,1.5};
599
600 // change the binning of the pT axis:
601 Double_t *xPt3 = new Double_t[nBins3[1]+1];
602 xPt3[0] = 0.;
603 for(int i = 1; i<=nBins3[1];i++){
604 if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.05; // 1 - 40
605 else if(xPt3[i-1]<4)xPt3[i] = xPt3[i-1] + 0.2; // 41 - 50
606 else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 0.5; // 50 - 62
607 else if(xPt3[i-1]<15)xPt3[i] = xPt3[i-1] + 1.; // 62 - 67
608 else if(xPt3[i-1]<20)xPt3[i] = xPt3[i-1] + 2.; // 67 - 72
609 else if(xPt3[i-1]<60)xPt3[i] = xPt3[i-1] + 5; // 72 - 76
610 else xPt3[i] = xPt3[i-1] + 10; // 76 - 100 = 140
611 }
612
613 fhnTrackPt[ij] = new THnSparseF(Form("fhnTrackPt%s",cAdd.Data()),";track number;p_{T};cent;#tracks;RP;trigger",nBinsSparse3,nBins3,xmin3,xmax3);
614 fhnTrackPt[ij]->SetBinEdges(1,xPt3);
615 fHistList->Add(fhnTrackPt[ij]);
616 delete [] xPt3;
617
618 // Track QA bins track nr, pTrack, cent, eta, phi bins 5.4 M
619 const Int_t nBinsSparse4 = 6;
620 const Int_t nBins4[nBinsSparse4] = { 2, 50, 10, 20, 180,2};
621 const Double_t xmin4[nBinsSparse4] = { -0.5, 0, 0, -1.0, 0.,-1.5};
622 const Double_t xmax4[nBinsSparse4] = { 1.5,150, 100, 1.0,2.*TMath::Pi(),1.5};
623
624 // change the binning ot the pT axis:
625 Double_t *xPt4 = new Double_t[nBins4[1]+1];
626 xPt4[0] = 0.;
627 for(int i = 1; i<=nBins4[1];i++){
628 if(xPt4[i-1]<2)xPt4[i] = xPt4[i-1] + 0.1;
629 else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 0.5;
630 else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] + 1.;
631 else if(xPt4[i-1]<30)xPt4[i] = xPt4[i-1] + 2.5;
632 else xPt4[i] = xPt4[i-1] + 5.;
633 }
634 fhnTrackPtQA[ij] = new THnSparseF(Form("fhnTrackPtQA%s",cAdd.Data()),";track number;p_{T};cent;#eta;#phi;sign",nBinsSparse4,nBins4,xmin4,xmax4);
635 fhnTrackPtQA[ij]->SetBinEdges(1,xPt4);
636 fHistList->Add(fhnTrackPtQA[ij]);
637 delete [] xPt4;
638
639 for(int i = 0;i <= kMaxJets;++i){
640 fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
641 fHistList->Add(fh1PtIn[ij][i]);
642
643
644 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
645 fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i),
646 Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};",
647 cAdd.Data()),
648 200,0.,200.,nBinPt,binLimitsPt);
649 fHistList->Add(fh2LTrackPtJetPt[ij][i]);
650 }
651
652
653 fh1DijetMinv[ij] = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
654 fHistList->Add(fh1DijetMinv[ij]);
655
656 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);
657 fHistList->Add(fh2DijetDeltaPhiPt[ij]);
658
659 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);
660 fHistList->Add(fh2DijetAsymPt[ij]);
661
662 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.);
663 fHistList->Add(fh2DijetPt2vsPt1[ij]);
664 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.);
665 fHistList->Add( fh2DijetDifvsSum[ij]);
666 }
667 // =========== Switch on Sumw2 for all histos ===========
668 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
669 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
670 if (h1){
671 h1->Sumw2();
672 continue;
673 }
674 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
675 if(hn)hn->Sumw2();
676 }
677 TH1::AddDirectory(oldStatus);
678}
679
680void AliAnalysisTaskJetSpectrum2::Init()
681{
682 //
683 // Initialization
684 //
685
686 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
687
688}
689
690void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
691
692
693 Bool_t selected = kTRUE;
694
695 if(fUseGlobalSelection&&fEventSelectionMask==0){
696 selected = AliAnalysisHelperJetTasks::Selected();
697 }
698 else if(fUseGlobalSelection&&fEventSelectionMask>0){
699 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
700 }
701
702 if(fEventClass>0){
703 selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
704 }
705
706 if(!selected){
707 // no selection by the service task, we continue
708 if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
709 PostData(1, fHistList);
710 return;
711 }
712
713
714 static AliAODEvent* aod = 0;
715
716 // take all other information from the aod we take the tracks from
717 if(!aod){
718 if(fUseAODTrackInput)aod = fAODIn;
719 else aod = fAODOut;
720 }
721
722
723 //
724 // Execute analysis for current event
725 //
726 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
727 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
728 if(!aodH){
729 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
730 return;
731 }
732
733 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
734 TClonesArray *aodRecJets = 0;
735
736 if(fAODOut&&!aodRecJets){
737 aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
738 }
739 if(fAODExtension&&!aodRecJets){
740 aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
741 }
742 if(fAODIn&&!aodRecJets){
743 aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
744 }
745
746
747
748 if(!aodRecJets){
749 if(fDebug){
750
751 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
752 if(fAODIn){
753 Printf("Input AOD >>>>");
754 fAODIn->Print();
755 }
756 if(fAODExtension){
757 Printf("AOD Extension >>>>");
758 fAODExtension->Print();
759 }
760 if(fAODOut){
761 Printf("Output AOD >>>>");
762 fAODOut->Print();
763 }
764 }
765 return;
766 }
767
768 TClonesArray *aodGenJets = 0;
769 if(fBranchGen.Length()>0){
770 if(fAODOut&&!aodGenJets){
771 aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
772 }
773 if(fAODExtension&&!aodGenJets){
774 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
775 }
776 if(fAODIn&&!aodGenJets){
777 aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
778 }
779
780 if(!aodGenJets){
781 Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
782 return;
783 }
784 }
785
786 TClonesArray *aodBackRecJets = 0;
787 if(fBranchBkgRec.Length()>0){
788 if(fAODOut&&!aodBackRecJets){
789 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgRec.Data()));
790 }
791 if(fAODExtension&&!aodBackRecJets){
792 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgRec.Data()));
793 }
794 if(fAODIn&&!aodBackRecJets){
795 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgRec.Data()));
796 }
797
798 if(!aodBackRecJets){
799 Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgRec.Data());
800 return;
801 }
802 }
803
804
805 TClonesArray *aodBackGenJets = 0;
806
807 if(fBranchBkgGen.Length()>0){
808 if(fAODOut&&!aodBackGenJets){
809 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgGen.Data()));
810 }
811 if(fAODExtension&&!aodBackGenJets){
812 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgGen.Data()));
813 }
814 if(fAODIn&&!aodBackGenJets){
815 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgGen.Data()));
816 }
817
818 if(!aodBackGenJets){
819 Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgGen.Data());
820 return;
821 }
822 }
823
824
825 // new Scheme
826 // first fill all the pure histograms, i.e. full jets
827 // in case of the correaltion limit to the n-leading jets
828
829 // reconstructed
830
831
832 // generated
833
834
835 // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
836
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 fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
857
858 // Getting some global properties
859 fCentrality = GetCentrality();
860 if(fCentrality<=0)fCentrality = 0;
861 fh1Centrality->Fill(fCentrality);
862
863
864 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
865 // this is the part we only use when we have MC information
866 AliMCEvent* mcEvent = MCEvent();
867 // AliStack *pStack = 0;
868 if(!mcEvent){
869 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
870 return;
871 }
872 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
873 if(pythiaGenHeader){
874 nTrials = pythiaGenHeader->Trials();
875 ptHard = pythiaGenHeader->GetPtHard();
876 int iProcessType = pythiaGenHeader->ProcessType();
877 // 11 f+f -> f+f
878 // 12 f+barf -> f+barf
879 // 13 f+barf -> g+g
880 // 28 f+g -> f+g
881 // 53 g+g -> f+barf
882 // 68 g+g -> g+g
883 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
884 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
885 }
886 }// (fAnalysisType&kMCESD)==kMCESD)
887
888
889 // we simply fetch the tracks/mc particles as a list of AliVParticles
890
891 TList recParticles;
892 TList genParticles;
893
894 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
895
896 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
897 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
898 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
899
900 // CalculateReactionPlaneAngle(&recParticles);
901 fRPAngle = 0;
902
903 if(fRPMethod==0)fRPAngle = aod->GetHeader()->GetEventplane();
904 else if(fRPMethod==1||fRPMethod==2){
905 fRPAngle = aod->GetHeader()->GetQTheta(fRPMethod);
906 }
907 fh1RP->Fill(fRPAngle);
908 fh2RPCentrality->Fill(fCentrality,fRPAngle);
909 // Event control and counting ...
910 // MC
911 fh1PtHard->Fill(ptHard,eventW);
912 fh1PtHardNoW->Fill(ptHard,1);
913 fh1PtHardTrials->Fill(ptHard,nTrials);
914
915 // Real
916 if(aod->GetPrimaryVertex()){// No vtx for pure MC
917 fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
918 }
919
920
921 Int_t recMult1 = recParticles.GetEntries();
922 Int_t genMult1 = genParticles.GetEntries();
923
924 Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
925 Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
926
927 fh2MultRec->Fill(recMult1,recMult2);
928 fh2MultGen->Fill(genMult1,genMult2);
929 fMultRec = recMult1;
930 if(fMultRec<=0)fMultRec = recMult2;
931 fMultGen = genMult1;
932 if(fMultGen<=0)fMultGen = genMult2;
933
934 Double_t var0[4] = {0,};
935 var0[0] = fCentrality;
936 var0[1] = fMultRec;
937 for(int it=0;it<fNTrigger;it++){
938 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
939 var0[2] = it;
940 var0[3] = GetRho(recJetsList);
941 fhnEvent->Fill(var0);
942 }
943 }
944 // the loops for rec and gen should be indentical... pass it to a separate
945 // function ...
946 // Jet Loop
947 // Track Loop
948 // Jet Jet Loop
949 // Jet Track loop
950
951 FillJetHistos(recJetsListCut,recParticles,kJetRec);
952 FillJetHistos(recJetsList,recParticles,kJetRecFull);
953 FillTrackHistos(recParticles,kJetRec);
954
955 FillJetHistos(genJetsListCut,genParticles,kJetGen);
956 FillJetHistos(genJetsList,genParticles,kJetGenFull);
957 FillTrackHistos(genParticles,kJetGen);
958
959 // Here follows the jet matching part
960 // delegate to a separated method?
961
962 if(fDoMatching){
963 FillMatchHistos(recJetsList,genJetsList);
964 }
965
966 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
967 PostData(1, fHistList);
968}
969
970void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
971
972 if(iType>=kJetTypes){
973 return;
974 }
975 if(!fFlagJetType[iType])return;
976
977 Int_t refMult = fMultRec;
978 if(iType==kJetGen||iType==kJetGenFull){
979 refMult = fMultGen;
980 }
981
982 Int_t nJets = jetsList.GetEntries();
983 fh1NJets[iType]->Fill(nJets);
984
985 if(nJets<=0)return;
986
987 Float_t ptOld = 1.E+32;
988 Float_t pT0 = 0;
989 Float_t pT1 = 0;
990 Float_t phi0 = 0;
991 Float_t phi1 = 0;
992 Int_t ij0 = -1;
993 Int_t ij1 = -1;
994
995 Double_t var1[9] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leadingTrackPt
996 Double_t var1b[3] = {0,}; // p_{T,jet};cent;trigger;
997
998 var1[2] = fCentrality;
999 var1b[1] = fCentrality;
1000 var1[3] = refMult;
1001
1002
1003
1004 Double_t var2[9] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area;trigger
1005 var2[2] = fCentrality;
1006
1007 for(int ij = 0;ij < nJets;ij++){
1008 AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
1009 Float_t ptJet = jet->Pt();
1010
1011
1012 if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0);
1013
1014 var1b[0] = ptJet;
1015 if(jet->Trigger()&fJetTriggerBestMask){
1016 fh1PtJetsInBest[iType]->Fill(ptJet);
1017 for(int it = 0;it <fNTrigger;it++){
1018 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1019 var1b[2] = it;
1020 fhnJetPtBest[iType]->Fill(var1b);
1021 }
1022 }
1023 }
1024 if(jet->Trigger()&fJetTriggerExcludeMask){
1025 fh1PtJetsInRej[iType]->Fill(ptJet);
1026 for(int it = 0;it <fNTrigger;it++){
1027 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1028 var1b[2] = it;
1029 fhnJetPtRej[iType]->Fill(var1b);
1030 }
1031 }
1032 continue;
1033 }
1034 fh1PtJetsIn[iType]->Fill(ptJet);
1035 ptOld = ptJet;
1036
1037 // find the dijets assume sorting and acceptance cut...
1038 if(ij==0){
1039 ij0 = ij;
1040 pT0 = ptJet;
1041 phi0 = jet->Phi();
1042 if(phi0<0)phi0+=TMath::Pi()*2.;
1043 }
1044 else if(ptJet>pT1){
1045 // take only the backward hemisphere??
1046 phi1 = jet->Phi();
1047 if(phi1<0)phi1+=TMath::Pi()*2.;
1048 Float_t dPhi = phi1 - phi0;
1049 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1050 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1051 if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
1052 ij1 = ij;
1053 pT1 = ptJet;
1054 }
1055 }
1056 // fill jet histos for kmax jets
1057 Float_t phiJet = jet->Phi();
1058 Float_t etaJet = jet->Eta();
1059 if(phiJet<0)phiJet+=TMath::Pi()*2.;
1060 fh1TmpRho->Reset();
1061 if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
1062
1063 fh1PtIn[iType][kMaxJets]->Fill(ptJet);
1064 // fill leading jets...
1065 AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
1066 //AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
1067 Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
1068 Double_t ptLead = jet->GetPtLeading(); //pT of leading jet
1069 AliVParticle *tt = (AliVParticle*)particlesList.At(0);
1070 Float_t ttphi=tt->Phi();
1071 if(ttphi<0)ttphi+=TMath::Pi()*2.;
1072 Float_t ttpt=tt->Pt();
1073 Int_t phiBintt = GetPhiBin(ttphi-fRPAngle);
1074 Double_t dphitrigjet=RelativePhi(ttphi,phiJet);
1075 if(fTRP==1){
1076 if(TMath::Abs(dphitrigjet)<TMath::Pi()-0.6) continue;
1077 var1[1] = ptJet;
1078 var1[4] = phiBintt;
1079 var1[5] = jet->EffectiveAreaCharged();
1080 var1[7] = ttpt;
1081 var1[8] = CheckAcceptance(phiJet,etaJet);
1082 }
1083
1084 if(fTRP==0){
1085 var1[1] = ptJet;
1086 var1[4] = phiBin;
1087 var1[5] = jet->EffectiveAreaCharged();
1088 var1[7] = ptLead;
1089 var1[8] = CheckAcceptance(phiJet,etaJet);
1090 }
1091
1092 //jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin;signed pt leading
1093 var2[1] = ptJet;
1094 var2[3] = etaJet;
1095 var2[4] = phiJet;
1096 var2[5] = jet->EffectiveAreaCharged();
1097 var2[7] = var1[8];
1098 var2[8] = (leadTrack?leadTrack->Charge()*ptLead:0);//pT of leading jet x charge
1099
1100 if(ij<kMaxJets){
1101 fh2LTrackPtJetPt[iType][ij]->Fill(jet->GetPtLeading(),ptJet);
1102 var1[0] = ij;
1103 var2[0] = ij;
1104 for(int it = 0;it <fNTrigger;it++){
1105 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1106 var1[6] = it;
1107 var2[6] = it;
1108 fhnJetPt[iType]->Fill(var1);
1109 fhnJetPtQA[iType]->Fill(var2);
1110 }
1111 }
1112 }
1113 var1[0] = kMaxJets;// fill for all jets
1114 var2[0] = kMaxJets;// fill for all jets
1115 for(int it = 0;it <fNTrigger;it++){
1116 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1117 var1[6] = it;
1118 var2[6] = it;
1119 fhnJetPt[iType]->Fill(var1);
1120 fhnJetPtQA[iType]->Fill(var2);
1121 }
1122 }
1123
1124 fh2LTrackPtJetPt[iType][kMaxJets]->Fill(jet->GetPtLeading(),ptJet);
1125
1126 if(particlesList.GetSize()&&ij<kMaxJets){
1127 // Particles... correlated with jets...
1128 for(int it = 0;it<particlesList.GetEntries();++it){
1129 AliVParticle *part = (AliVParticle*)particlesList.At(it);
1130 Float_t deltaR = jet->DeltaR(part);
1131 if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
1132 }
1133 // fill the jet shapes
1134 }// if we have particles
1135 }// Jet Loop
1136
1137
1138 // do something with dijets...
1139 if(ij0>=0&&ij1>0){
1140 AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1141 Double_t ptJet0 = jet0->Pt();
1142 Double_t phiJet0 = jet0->Phi();
1143 if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1144
1145 AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1146 Double_t ptJet1 = jet1->Pt();
1147 Double_t phiJet1 = jet1->Phi();
1148 if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1149
1150 Float_t deltaPhi = phiJet0 - phiJet1;
1151 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1152 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1153 deltaPhi = TMath::Abs(deltaPhi);
1154 fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
1155
1156 Float_t asym = 9999;
1157 if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
1158 fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1159 fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1160 fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1161 Float_t minv = 2.*(jet0->P()*jet1->P()-
1162 jet0->Px()*jet1->Px()-
1163 jet0->Py()*jet1->Py()-
1164 jet0->Pz()*jet1->Pz()); // assume mass == 0;
1165 if(minv<0)minv=0; // prevent numerical instabilities
1166 minv = TMath::Sqrt(minv);
1167 fh1DijetMinv[iType]->Fill(minv);
1168 }
1169
1170
1171
1172 // count down the jets above thrueshold
1173 Int_t nOver = nJets;
1174 if(nOver>0){
1175 TIterator *jetIter = jetsList.MakeIterator();
1176 AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
1177 if(tmpJet){
1178 Float_t pt = tmpJet->Pt();
1179 for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1180 Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1181 while(pt<ptCut&&tmpJet){
1182 nOver--;
1183 tmpJet = (AliAODJet*)(jetIter->Next());
1184 if(tmpJet){
1185 pt = tmpJet->Pt();
1186 }
1187 }
1188 if(nOver<=0)break;
1189 fh2NJetsPt[iType]->Fill(ptCut,nOver);
1190 }
1191 }
1192 delete jetIter;
1193 }
1194}
1195
1196void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1197
1198 if(fFlagJetType[iType]<=0)return;
1199 Int_t refMult = fMultRec;
1200 if(iType==kJetGen||iType==kJetGenFull){
1201 refMult = fMultGen;
1202
1203 }
1204
1205 //
1206 Double_t var3[7]; // track number;p_{T};cent;#tracks;RP;cahrge
1207 var3[2] = fCentrality;
1208 var3[3] = refMult;
1209 Double_t var4[6]; // track number;p_{T};cent;#eta;#phi;sign
1210 var4[2] = fCentrality;
1211 Int_t nTrackOver = particlesList.GetSize();
1212 // do the same for tracks and jets
1213 if(nTrackOver>0){
1214 TIterator *trackIter = particlesList.MakeIterator();
1215 AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1216 Float_t pt = tmpTrack->Pt();
1217 for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1218 Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1219 while(pt<ptCut&&tmpTrack){
1220 nTrackOver--;
1221 tmpTrack = (AliVParticle*)(trackIter->Next());
1222 if(tmpTrack){
1223 pt = tmpTrack->Pt();
1224 }
1225 }
1226 if(nTrackOver<=0)break;
1227 fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1228 }
1229
1230
1231 trackIter->Reset();
1232 AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1233 Float_t sumPt = 0;
1234
1235 while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1236 Float_t tmpPt = tmpTrack->Pt();
1237 fh1PtTracksIn[iType]->Fill(tmpPt);
1238 fh1PtTracksInLow[iType]->Fill(tmpPt);
1239
1240 sumPt += tmpPt;
1241 Float_t tmpPhi = tmpTrack->Phi();
1242 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
1243
1244
1245 Float_t phiRP = tmpPhi-fRPAngle;
1246 if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1247 if(phiRP<0)phiRP += TMath::Pi();
1248 if(phiRP<0)phiRP += TMath::Pi();
1249 const float allPhi = -1./180.*TMath::Pi();
1250
1251 if(tmpPt<100){
1252 fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1253 fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1254
1255 fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1256 fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1257 }
1258 Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
1259 var3[0] = 1;
1260 var3[1] = tmpPt;
1261 var3[4] = phiBin;
1262 var3[6] = tmpTrack->Charge();
1263
1264 var4[0] = 1;
1265 var4[1] = tmpPt;
1266 var4[3] = tmpTrack->Eta();
1267 var4[4] = tmpPhi;
1268 var4[5] = tmpTrack->Charge();
1269
1270
1271 for(int it = 0;it <fNTrigger;it++){
1272 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1273 var3[0] = 1;
1274 var3[5] = it;
1275 fhnTrackPt[iType]->Fill(var3);
1276 if(tmpTrack==leading){
1277 var3[0] = 0;
1278 fhnTrackPt[iType]->Fill(var3);
1279 }
1280 }
1281 }
1282
1283
1284
1285
1286 fhnTrackPtQA[iType]->Fill(var4);
1287 if(tmpTrack==leading){
1288 var4[0] = 0;
1289 fhnTrackPtQA[iType]->Fill(var4);
1290 continue;
1291 }
1292
1293
1294
1295
1296
1297 }
1298 fh1SumPtTrack[iType]->Fill(sumPt);
1299 delete trackIter;
1300 }
1301
1302}
1303
1304
1305void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1306
1307
1308 // Fill al the matching histos
1309 // It is important that the acceptances for the mathing are as large as possible
1310 // to avoid false matches on the edge of acceptance
1311 // therefore we add some extra matching jets as overhead
1312
1313 static TArrayI aGenIndex(recJetsList.GetEntries());
1314 if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
1315
1316 static TArrayI aRecIndex(genJetsList.GetEntries());
1317 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
1318
1319 if(fDebug){
1320 Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1321 Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1322 }
1323 AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1324 &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1325 aGenIndex,aRecIndex,fDebug);
1326
1327 if(fDebug){
1328 for(int i = 0;i< aGenIndex.GetSize();++i){
1329 if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1330 }
1331 for(int i = 0;i< aRecIndex.GetSize();++i){
1332 if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1333 }
1334 }
1335
1336 Double_t container[8];
1337 Double_t containerRec[4];
1338 Double_t containerGen[4];
1339
1340 // loop over generated jets
1341 // consider the
1342 for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1343 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1344 Double_t ptGen = genJet->Pt();
1345 Double_t phiGen = genJet->Phi();
1346 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1347 Double_t etaGen = genJet->Eta();
1348 containerGen[0] = ptGen;
1349 containerGen[1] = etaGen;
1350 containerGen[2] = phiGen;
1351 containerGen[3] = genJet->GetPtLeading();
1352
1353 fhnJetContainer->Fill(containerGen,kStep0); //all generated jets
1354
1355 if(JetSelected(genJet))
1356 fhnJetContainer->Fill(containerGen,kStep1); // all generated jets in eta window
1357
1358 Int_t ir = aRecIndex[ig];
1359 if(ir>=0&&ir<recJetsList.GetEntries()){
1360 AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
1361
1362 fhnJetContainer->Fill(containerGen,kStep2); // all generated jets with reconstructed partner
1363
1364 if(JetSelected(genJet)){
1365 fhnJetContainer->Fill(containerGen,kStep3); // all generated jets in eta window with reconstructed partner
1366
1367 if(JetSelected(recJet)) {
1368 fhnJetContainer->Fill(containerGen,kStep4); // all generated jets in eta window with reconstructed partner in eta window
1369
1370 }
1371 containerGen[3] = recJet->GetPtLeading();
1372 fhnJetContainer->Fill(containerGen,kStep5); // all generated jets in eta window with reconstructed partner with leading track on reconstructed level
1373 }
1374 }
1375 }// loop over generated jets used for matching...
1376
1377
1378
1379 // loop over reconstructed jets
1380 for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1381 AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1382 Double_t etaRec = recJet->Eta();
1383 Double_t ptRec = recJet->Pt();
1384 Double_t phiRec = recJet->Phi();
1385 if(phiRec<0)phiRec+=TMath::Pi()*2.;
1386
1387 containerRec[0] = ptRec;
1388 containerRec[1] = etaRec;
1389 containerRec[2] = phiRec;
1390 containerRec[3] = recJet->GetPtLeading();
1391
1392 container[0] = ptRec;
1393 container[1] = etaRec;
1394 container[2] = phiRec;
1395 container[3] = recJet->GetPtLeading();
1396
1397 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1398
1399 if(JetSelected(recJet)){
1400 fhnJetContainer->Fill(containerRec,kStep7); //all rec jets in eta window
1401 // Fill Correlation
1402 Int_t ig = aGenIndex[ir];
1403 if(ig>=0 && ig<genJetsList.GetEntries()){
1404 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1405 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1406 Double_t ptGen = genJet->Pt();
1407 Double_t phiGen = genJet->Phi();
1408 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1409 Double_t etaGen = genJet->Eta();
1410
1411 containerGen[0] = ptGen;
1412 containerGen[1] = etaGen;
1413 containerGen[2] = phiGen;
1414 containerGen[3] = genJet->GetPtLeading();
1415
1416 container[4] = ptGen;
1417 container[5] = etaGen;
1418 container[6] = phiGen;
1419 container[7] = genJet->GetPtLeading();
1420
1421 fhnJetContainer->Fill(containerGen,kStep6); // all rec jets in eta window with generated partner
1422
1423 fhnCorrelation->Fill(container);
1424
1425 Float_t delta = (ptRec-ptGen)/ptGen;
1426 fh2RelPtFGen->Fill(ptGen,delta);
1427 fh3RelPtFGenLeadTrkPt->Fill(ptGen,delta,recJet->GetPtLeading());
1428 fh2PtFGen->Fill(ptGen,ptRec);
1429
1430 fh2deltaPt1Pt2->Fill(ptRec-ptGen,ptRec,ptGen);
1431
1432 }
1433 }// loop over reconstructed jets
1434 }
1435 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1436}
1437
1438
1439void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1440 //
1441 // Create the particle container for the correction framework manager and
1442 // link it
1443 //
1444 const Int_t kNvar = 4 ; //number of variables on the grid:pt,eta, phi
1445 const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
1446 const Double_t kEtamin = -1.0, kEtamax = 1.0;
1447 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1448 const Double_t kPtLeadingTrackPtMin = 0., kPtLeadingTrackPtMax = 200.;
1449
1450 // can we neglect migration in eta and phi?
1451 // phi should be no problem since we cover full phi and are phi symmetric
1452 // eta migration is more difficult due to needed acceptance correction
1453 // in limited eta range
1454
1455 //arrays for the number of bins in each dimension
1456 Int_t iBin[kNvar];
1457 iBin[0] = 125; //bins in pt
1458 iBin[1] = 4; //bins in eta
1459 iBin[2] = 1; // bins in phi
1460 iBin[3] = 10; // bins in leading track Pt
1461
1462 //arrays for lower bounds :
1463 Double_t* binEdges[kNvar];
1464 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1465 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1466
1467 //values for bin lower bounds
1468 // 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);
1469 for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
1470 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1471 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1472 binEdges[3][0]= kPtLeadingTrackPtMin;
1473 binEdges[3][1]= 1.;
1474 binEdges[3][2]= 2.;
1475 binEdges[3][3]= 3.;
1476 binEdges[3][4]= 4.;
1477 binEdges[3][5]= 5.;
1478 binEdges[3][6]= 6.;
1479 binEdges[3][7]= 8.;
1480 binEdges[3][8]= 10.;
1481 binEdges[3][9]= 12.;
1482 binEdges[3][10]= kPtLeadingTrackPtMax;
1483
1484 fhnJetContainer = new AliCFContainer(Form("fhnJetContainer"),Form("AliCFContainer jet info"),kMaxStep,kNvar,iBin);
1485 for (int k=0; k<kNvar; k++) {
1486 fhnJetContainer->SetBinLimits(k,binEdges[k]);
1487 }
1488
1489 //create correlation matrix for unfolding
1490 Int_t thnDim[2*kNvar];
1491 for (int k=0; k<kNvar; k++) {
1492 //first half : reconstructed
1493 //second half : MC
1494 thnDim[k] = iBin[k];
1495 thnDim[k+kNvar] = iBin[k];
1496 }
1497
1498 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparseF with correlations",2*kNvar,thnDim);
1499 for (int k=0; k<kNvar; k++) {
1500 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1501 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
1502 }
1503
1504 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1505 delete [] binEdges[ivar];
1506
1507
1508}
1509
1510void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1511{
1512 // Terminate analysis
1513 //
1514 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1515}
1516
1517
1518Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1519
1520 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
1521
1522 Int_t iCount = 0;
1523 if(type==kTrackAOD){
1524 AliAODEvent *aod = 0;
1525 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1526 else aod = AODEvent();
1527 if(!aod){
1528 return iCount;
1529 }
1530 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1531 AliAODTrack *tr = aod->GetTrack(it);
1532 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
1533 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1534 if(tr->Pt()<fMinTrackPt)continue;
1535 if(fDebug>0){
1536 if(tr->Pt()>20){
1537 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
1538 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
1539 tr->Print();
1540 // tr->Dump();
1541 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1542 if(fESD){
1543 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
1544 esdTr->Print("");
1545 // esdTr->Dump();
1546 }
1547 }
1548 }
1549 list->Add(tr);
1550 iCount++;
1551 }
1552 }
1553 else if (type == kTrackKineAll||type == kTrackKineCharged){
1554 AliMCEvent* mcEvent = MCEvent();
1555 if(!mcEvent)return iCount;
1556 // we want to have alivpartilces so use get track
1557 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1558 if(!mcEvent->IsPhysicalPrimary(it))continue;
1559 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
1560 if(part->Pt()<fMinTrackPt)continue;
1561 if(type == kTrackKineAll){
1562 list->Add(part);
1563 iCount++;
1564 }
1565 else if(type == kTrackKineCharged){
1566 if(part->Particle()->GetPDG()->Charge()==0)continue;
1567 list->Add(part);
1568 iCount++;
1569 }
1570 }
1571 }
1572 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1573 AliAODEvent *aod = 0;
1574 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1575 else aod = AODEvent();
1576 if(!aod)return iCount;
1577 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1578 if(!tca)return iCount;
1579 for(int it = 0;it < tca->GetEntriesFast();++it){
1580 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
1581 if(!part)continue;
1582 if(part->Pt()<fMinTrackPt)continue;
1583 if(!part->IsPhysicalPrimary())continue;
1584 if(type == kTrackAODMCAll){
1585 list->Add(part);
1586 iCount++;
1587 }
1588 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
1589 if(part->Charge()==0)continue;
1590 if(kTrackAODMCCharged){
1591 list->Add(part);
1592 }
1593 else {
1594 if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
1595 list->Add(part);
1596 }
1597 iCount++;
1598 }
1599 }
1600 }// AODMCparticle
1601 list->Sort();
1602 return iCount;
1603
1604}
1605
1606
1607Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1608 AliAODEvent *aod = 0;
1609 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1610 else aod = AODEvent();
1611 if(!aod){
1612 return 101;
1613 }
1614 return aod->GetHeader()->GetCentrality();
1615}
1616
1617
1618
1619Bool_t AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1620 Bool_t selected = false;
1621
1622 if(!jet)return selected;
1623
1624 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1625 selected = kTRUE;
1626 }
1627 return selected;
1628
1629}
1630
1631Int_t AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1632
1633 if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1634 Int_t iCount = 0;
1635
1636 if(!jarray){
1637 Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1638 return iCount;
1639 }
1640
1641
1642 for(int ij=0;ij<jarray->GetEntries();ij++){
1643 AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1644 if(!jet)continue;
1645 if(type==0){
1646 // No cut at all, main purpose here is sorting
1647 list->Add(jet);
1648 iCount++;
1649 }
1650 else if(type == 1){
1651 // eta cut
1652 if(JetSelected(jet)){
1653 list->Add(jet);
1654 iCount++;
1655 }
1656 }
1657 }
1658
1659 list->Sort();
1660 return iCount;
1661
1662}
1663
1664
1665Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1666 if(!jets)return 0;
1667
1668 Int_t refMult = 0;
1669 for(int ij = 0;ij < jets->GetEntries();++ij){
1670 AliAODJet* jet = (AliAODJet*)jets->At(ij);
1671 if(!jet)continue;
1672 TRefArray *refs = jet->GetRefTracks();
1673 if(!refs)continue;
1674 refMult += refs->GetEntries();
1675 }
1676 return refMult;
1677
1678}
1679
1680
1681AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1682 if(!jet)return 0;
1683 TRefArray *refs = jet->GetRefTracks();
1684 if(!refs) return 0;
1685 AliVParticle *leading = 0;
1686 Float_t fMaxPt = 0;
1687 for(int i = 0;i<refs->GetEntries();i++){
1688 AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1689 if(!tmp)continue;
1690 if(tmp->Pt()>fMaxPt){
1691 leading = tmp;
1692 fMaxPt = tmp->Pt();
1693 }
1694 }
1695 return leading;
1696}
1697
1698
1699AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1700 if(!jet)return 0;
1701 if(!list) return 0;
1702 AliVParticle *leading = 0;
1703 Float_t fMaxPt = 0;
1704 for(int i = 0;i<list->GetEntries();i++){
1705 AliVParticle *tmp = (AliVParticle*)(list->At(i));
1706 if(!tmp)continue;
1707 if(jet->DeltaR(tmp)>r)continue;
1708 if(tmp->Pt()>fMaxPt){
1709 leading = tmp;
1710 fMaxPt = tmp->Pt();
1711 }
1712 }
1713 return leading;
1714}
1715
1716
1717Double_t AliAnalysisTaskJetSpectrum2::RelativePhi(Double_t mphi,Double_t vphi){
1718
1719 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1720 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1721 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1722 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1723 double dphi = mphi-vphi;
1724 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1725 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1726 return dphi;//dphi in [-Pi, Pi]
1727}
1728
1729
1730
1731
1732
1733
1734Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1735{
1736 Int_t phibin=-1;
1737 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1738 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1739 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1740 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1741 return phibin;
1742}
1743
1744void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){
1745 //
1746 // set number of trigger bins
1747 //
1748 if(n>0){
1749 fNTrigger = n;
1750 delete [] fTriggerName;
1751 fTriggerName = new TString [fNTrigger];
1752 delete [] fTriggerBit;fTriggerBit = 0;
1753 fTriggerBit = new UInt_t [fNTrigger];
1754 }
1755 else{
1756 fNTrigger = 0;
1757 }
1758}
1759
1760
1761void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){
1762 //
1763 // set trigger bin
1764 //
1765 if(i<fNTrigger){
1766 fTriggerBit[i] = it;
1767 fTriggerName[i] = c;
1768 }
1769}
1770
1771
1772void AliAnalysisTaskJetSpectrum2::SetNAcceptance(Int_t n){
1773 //
1774 // Set number of acceptance bins
1775 //
1776
1777 if(n>0){
1778 fNAcceptance = n;
1779 delete [] fAcceptancePhiMin;
1780 delete [] fAcceptancePhiMax;
1781 delete [] fAcceptanceEtaMin;
1782 delete [] fAcceptanceEtaMax;
1783
1784 fAcceptancePhiMin = new Float_t[fNAcceptance];
1785 fAcceptancePhiMax = new Float_t[fNAcceptance];
1786 fAcceptanceEtaMin = new Float_t[fNAcceptance];
1787 fAcceptanceEtaMax = new Float_t[fNAcceptance];
1788 }
1789 else{
1790 fNTrigger = 0;
1791 }
1792}
1793
1794void AliAnalysisTaskJetSpectrum2::SetAcceptance(Int_t i,Float_t phiMin,Float_t phiMax,Float_t etaMin,Float_t etaMax){
1795 //
1796 // Set acceptance bins
1797 //
1798
1799 if(i<fNAcceptance){
1800 Printf("%s:%d Setting acceptance %d phi %3.3f-%3.3f eta %3.3f-%3.3f",(char*)__FILE__,__LINE__,i,phiMin,phiMax,etaMin,etaMax);
1801
1802 fAcceptancePhiMin[i] = phiMin;
1803 fAcceptancePhiMax[i] = phiMax;
1804 fAcceptanceEtaMin[i] = etaMin;
1805 fAcceptanceEtaMax[i] = etaMax;
1806 }
1807 else{
1808 fNTrigger = 0;
1809 }
1810}
1811
1812
1813Int_t AliAnalysisTaskJetSpectrum2::CheckAcceptance(Float_t phi,Float_t eta){
1814
1815 //
1816 // return aceptnace bin do no use ovelapping bins
1817 //
1818
1819 for(int i = 0;i<fNAcceptance;i++){
1820 if(phi<fAcceptancePhiMin[i])continue;
1821 if(phi>fAcceptancePhiMax[i])continue;
1822 if(eta<fAcceptanceEtaMin[i])continue;
1823 if(eta>fAcceptanceEtaMax[i])continue;
1824 return i;
1825 }
1826 // catch the rest
1827 return fNAcceptance;
1828}
1829
1830Float_t AliAnalysisTaskJetSpectrum2::GetRho(TList &list){
1831
1832 // invert the correction
1833 AliAODJet *jet = (AliAODJet*)list.At(0); // highest pt jet
1834 if(!jet)return -1;
1835 if(jet->EffectiveAreaCharged()<=0)return -1;
1836 Float_t rho = jet->ChargedBgEnergy()/jet->EffectiveAreaCharged();
1837 return rho;
1838
1839}
1840
1841
1842
1843AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){
1844 //
1845 delete [] fTriggerBit;
1846 delete [] fTriggerName;
1847
1848 delete [] fAcceptancePhiMin;
1849 delete [] fAcceptancePhiMax;
1850 delete [] fAcceptanceEtaMin;
1851 delete [] fAcceptanceEtaMax;
1852
1853
1854
1855}