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