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