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