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