]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGJE/AliAnalysisTaskJetSpectrum2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coverity 18434
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetSpectrum2.cxx
CommitLineData
3b7ffecf 1// **************************************
3170a3f8 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>
742ee86c 25#include <TRandom3.h>
3b7ffecf 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>
2eded560 35#include <TProfile2D.h>
3b7ffecf 36#include <TList.h>
37#include <TLorentzVector.h>
38#include <TClonesArray.h>
37eb26ea 39#include <TRefArray.h>
3b7ffecf 40#include "TDatabasePDG.h"
41
42#include "AliAnalysisTaskJetSpectrum2.h"
43#include "AliAnalysisManager.h"
44#include "AliJetFinder.h"
b99018c6 45#include "AliTHn.h"
3b7ffecf 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"
c2785065 55#include "AliAODJetEventBackground.h"
3b7ffecf 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
37eb26ea 70AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2():
71 AliAnalysisTaskSE(),
a31b8a87 72 fJetHeaderRec(0x0),
73 fJetHeaderGen(0x0),
74 fAODIn(0x0),
75 fAODOut(0x0),
76 fAODExtension(0x0),
b99018c6 77 fhnJetContainer(0x0),
a31b8a87 78 fhnCorrelation(0x0),
90a006c2 79 fhnEvent(0x0),
c8eabe24 80 f1PtScale(0x0),
a31b8a87 81 fBranchRec("jets"),
82 fBranchGen(""),
37eb26ea 83 fBranchBkgRec(""),
84 fBranchBkgGen(""),
a31b8a87 85 fNonStdFile(""),
742ee86c 86 fRandomizer(0x0),
a31b8a87 87 fUseAODJetInput(kFALSE),
88 fUseAODTrackInput(kFALSE),
89 fUseAODMCInput(kFALSE),
90 fUseGlobalSelection(kFALSE),
91 fUseExternalWeightOnly(kFALSE),
92 fLimitGenJetEta(kFALSE),
b99018c6 93 fDoMatching(kFALSE),
a31b8a87 94 fNMatchJets(5),
7b822bc5 95 fNRPBins(3),
65c43de8 96 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
d7396b04 97 fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
a31b8a87 98 fFilterMask(0),
f2dd0695 99 fEventSelectionMask(0),
3d1dba29 100 fNTrigger(0),
101 fTriggerBit(0x0),
a31b8a87 102 fAnalysisType(0),
3b7ffecf 103 fTrackTypeRec(kTrackUndef),
104 fTrackTypeGen(kTrackUndef),
f4132e7d 105 fEventClass(0),
2eded560 106 fRPMethod(0),
3b7ffecf 107 fAvgTrials(1),
108 fExternalWeight(1),
a31b8a87 109 fJetRecEtaWindow(0.5),
110 fTrackRecEtaWindow(0.5),
9280dfa6 111 fMinJetPt(0),
a31b8a87 112 fMinTrackPt(0.15),
113 fDeltaPhiWindow(90./180.*TMath::Pi()),
742ee86c 114 fCentrality(100),
7b822bc5 115 fRPAngle(0),
37eb26ea 116 fMultRec(0),
117 fMultGen(0),
3d1dba29 118 fTriggerName(0x0),
3b7ffecf 119 fh1Xsec(0x0),
120 fh1Trials(0x0),
121 fh1PtHard(0x0),
122 fh1PtHardNoW(0x0),
123 fh1PtHardTrials(0x0),
57ca1193 124 fh1ZVtx(0x0),
7b822bc5 125 fh1RP(0x0),
742ee86c 126 fh1Centrality(0x0),
c8eabe24 127 fh1TmpRho(0x0),
37eb26ea 128 fh2MultRec(0x0),
129 fh2MultGen(0x0),
742ee86c 130 fh2RPCentrality(0x0),
a31b8a87 131 fh2PtFGen(0x0),
9a83d4af 132 fh2RelPtFGen(0x0),
3b7ffecf 133 fHistList(0x0)
134{
c17f867b 135
a31b8a87 136 for(int ij = 0;ij <kJetTypes;++ij){
742ee86c 137 fFlagJetType[ij] = 1; // default = on
a31b8a87 138 fh1NJets[ij] = 0;
139 fh1SumPtTrack[ij] = 0;
140 fh1PtJetsIn[ij] = 0;
65c43de8 141 fh1PtJetsInRej[ij] = 0;
b01d3ab7 142 fh1PtJetsInBest[ij] = 0;
a31b8a87 143 fh1PtTracksIn[ij] = 0;
cb883243 144 fh1PtTracksInLow[ij] = 0;
a31b8a87 145 fh2NJetsPt[ij] = 0;
146 fh2NTracksPt[ij] = 0;
2eded560 147 fp2MultRPPhiTrackPt[ij] = 0;
148 fp2CentRPPhiTrackPt[ij] = 0;
90a006c2 149 fhnJetPt[ij] = 0;
150 fhnJetPtQA[ij] = 0;
151 fhnTrackPt[ij] = 0;
152 fhnTrackPtQA[ij] = 0;
37eb26ea 153 for(int i = 0;i <= kMaxJets;++i){
d7117cbd 154 fh2LTrackPtJetPt[ij][i] = 0;
a31b8a87 155 fh1PtIn[ij][i] = 0;
a31b8a87 156 }
157
158 fh1DijetMinv[ij] = 0;
159 fh2DijetDeltaPhiPt[ij] = 0;
160 fh2DijetAsymPt[ij] = 0;
161 fh2DijetPt2vsPt1[ij] = 0;
162 fh2DijetDifvsSum[ij] = 0;
163
164 }
3b7ffecf 165}
166
167AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
168 AliAnalysisTaskSE(name),
169 fJetHeaderRec(0x0),
170 fJetHeaderGen(0x0),
a31b8a87 171 fAODIn(0x0),
172 fAODOut(0x0),
173 fAODExtension(0x0),
b99018c6 174 fhnJetContainer(0x0),
3b7ffecf 175 fhnCorrelation(0x0),
90a006c2 176 fhnEvent(0x0),
c8eabe24 177 f1PtScale(0x0),
3b7ffecf 178 fBranchRec("jets"),
179 fBranchGen(""),
37eb26ea 180 fBranchBkgRec(""),
181 fBranchBkgGen(""),
a31b8a87 182 fNonStdFile(""),
742ee86c 183 fRandomizer(0x0),
565584e8 184 fUseAODJetInput(kFALSE),
185 fUseAODTrackInput(kFALSE),
186 fUseAODMCInput(kFALSE),
b5cc0c6d 187 fUseGlobalSelection(kFALSE),
3b7ffecf 188 fUseExternalWeightOnly(kFALSE),
189 fLimitGenJetEta(kFALSE),
b99018c6 190 fDoMatching(kFALSE),
a31b8a87 191 fNMatchJets(5),
7b822bc5 192 fNRPBins(3),
65c43de8 193 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
d7396b04 194 fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
3b7ffecf 195 fFilterMask(0),
f2dd0695 196 fEventSelectionMask(0),
3d1dba29 197 fNTrigger(0),
198 fTriggerBit(0x0),
3b7ffecf 199 fAnalysisType(0),
200 fTrackTypeRec(kTrackUndef),
201 fTrackTypeGen(kTrackUndef),
f4132e7d 202 fEventClass(0),
2eded560 203 fRPMethod(0),
3b7ffecf 204 fAvgTrials(1),
205 fExternalWeight(1),
a31b8a87 206 fJetRecEtaWindow(0.5),
207 fTrackRecEtaWindow(0.5),
9280dfa6 208 fMinJetPt(0),
a31b8a87 209 fMinTrackPt(0.15),
210 fDeltaPhiWindow(90./180.*TMath::Pi()),
742ee86c 211 fCentrality(100),
7b822bc5 212 fRPAngle(0),
37eb26ea 213 fMultRec(0),
214 fMultGen(0),
3d1dba29 215 fTriggerName(0x0),
3b7ffecf 216 fh1Xsec(0x0),
217 fh1Trials(0x0),
218 fh1PtHard(0x0),
219 fh1PtHardNoW(0x0),
220 fh1PtHardTrials(0x0),
57ca1193 221 fh1ZVtx(0x0),
7b822bc5 222 fh1RP(0x0),
742ee86c 223 fh1Centrality(0x0),
c8eabe24 224 fh1TmpRho(0x0),
37eb26ea 225 fh2MultRec(0x0),
226 fh2MultGen(0x0),
742ee86c 227 fh2RPCentrality(0x0),
a31b8a87 228 fh2PtFGen(0x0),
9a83d4af 229 fh2RelPtFGen(0x0),
3b7ffecf 230 fHistList(0x0)
231{
232
a31b8a87 233 for(int ij = 0;ij <kJetTypes;++ij){
742ee86c 234 fFlagJetType[ij] = 1; // default = on
a31b8a87 235 fh1NJets[ij] = 0;
236 fh1SumPtTrack[ij] = 0;
237 fh1PtJetsIn[ij] = 0;
65c43de8 238 fh1PtJetsInRej[ij] = 0;
b01d3ab7 239 fh1PtJetsInBest[ij] = 0;
a31b8a87 240 fh1PtTracksIn[ij] = 0;
cb883243 241 fh1PtTracksInLow[ij] = 0;
a31b8a87 242 fh2NJetsPt[ij] = 0;
243 fh2NTracksPt[ij] = 0;
2eded560 244 fp2MultRPPhiTrackPt[ij] = 0;
245 fp2CentRPPhiTrackPt[ij] = 0;
90a006c2 246 fhnJetPt[ij] = 0;
247 fhnJetPtQA[ij] = 0;
248 fhnTrackPt[ij] = 0;
249 fhnTrackPtQA[ij] = 0;
37eb26ea 250 for(int i = 0;i <= kMaxJets;++i){
d7117cbd 251 fh2LTrackPtJetPt[ij][i] = 0;
a31b8a87 252 fh1PtIn[ij][i] = 0;
a31b8a87 253 }
254
255 fh1DijetMinv[ij] = 0;
256 fh2DijetDeltaPhiPt[ij] = 0;
257 fh2DijetAsymPt[ij] = 0;
258 fh2DijetPt2vsPt1[ij] = 0;
259 fh2DijetDifvsSum[ij] = 0;
260 }
261
37eb26ea 262 DefineOutput(1, TList::Class());
3b7ffecf 263}
264
265
266
267Bool_t AliAnalysisTaskJetSpectrum2::Notify()
268{
a31b8a87 269
270
271
3b7ffecf 272 //
273 // Implemented Notify() to read the cross sections
274 // and number of trials from pyxsec.root
275 //
a31b8a87 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());
e855f5c5 288 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
a31b8a87 289 if(!fAODExtension){
290 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
291 }
292 }
293
3b7ffecf 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();
a221560c 314 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
3b7ffecf 315 }
a31b8a87 316
d7117cbd 317 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
a31b8a87 318
3b7ffecf 319 return kTRUE;
320}
321
322void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
323{
324
742ee86c 325
3b7ffecf 326 // Connect the AOD
327
3b7ffecf 328 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
3b7ffecf 329 OpenFile(1);
37eb26ea 330 if(!fHistList)fHistList = new TList();
40440b28 331 PostData(1, fHistList); // post data in any case once
332
742ee86c 333 if(!fRandomizer)fRandomizer = new TRandom3(0);
334
5cb0f438 335 fHistList->SetOwner(kTRUE);
37eb26ea 336 Bool_t oldStatus = TH1::AddDirectoryStatus();
3b7ffecf 337 TH1::AddDirectory(kFALSE);
338
90a006c2 339
340
341 // event npsparse cent, mult
3d1dba29 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};
90a006c2 346
347
348 fhnEvent = new THnSparseF("fhnEvent",";cent;mult",nBinsSparse0,
349 nBins0,xmin0,xmax0);
350 fHistList->Add(fhnEvent);
351
b99018c6 352 if(fDoMatching){
353 MakeJetContainer();
354 fHistList->Add(fhnCorrelation);
355 fHistList->Add(fhnJetContainer);
356 }
3b7ffecf 357 //
358 // Histogram
359
90a006c2 360
361
ed5c9f3f 362 const Int_t nBinPt = 120;
3b7ffecf 363 Double_t binLimitsPt[nBinPt+1];
364 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
365 if(iPt == 0){
ed5c9f3f 366 binLimitsPt[iPt] = -50.0;
3b7ffecf 367 }
368 else {// 1.0
e24d63c0 369 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5;
3b7ffecf 370 }
371 }
edfbe476 372 const Int_t nBinPhi = 90;
3b7ffecf 373 Double_t binLimitsPhi[nBinPhi+1];
374 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
375 if(iPhi==0){
37eb26ea 376 binLimitsPhi[iPhi] = 0;
3b7ffecf 377 }
378 else{
379 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
380 }
381 }
382
edfbe476 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{
cf049e94 390 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
edfbe476 391 }
392 }
393
a31b8a87 394
3b7ffecf 395 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
396 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
90a006c2 397 fHistList->Add(fh1Xsec);
3b7ffecf 398 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
399 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
a31b8a87 400 fHistList->Add(fh1Trials);
3b7ffecf 401 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 402 fHistList->Add(fh1PtHard);
3b7ffecf 403 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 404 fHistList->Add(fh1PtHardNoW);
3b7ffecf 405 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 406 fHistList->Add(fh1PtHardTrials);
a31b8a87 407
57ca1193 408 fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
a31b8a87 409 fHistList->Add(fh1ZVtx);
7b822bc5 410
2eded560 411
9939142f 412 fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
7b822bc5 413 fHistList->Add(fh1RP);
414
eb1739b7 415 fh1Centrality = new TH1F("fh1Centrality","cent;cent (%)",103,-1,102);
742ee86c 416 fHistList->Add(fh1Centrality);
417
3170a3f8 418 fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",500,0,5000,500,0.,5000);
37eb26ea 419 fHistList->Add(fh2MultRec);
3170a3f8 420 fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,0,5000,500,0.,5000);
37eb26ea 421 fHistList->Add(fh2MultGen);
a31b8a87 422
a31b8a87 423
742ee86c 424 fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
425 fHistList->Add(fh2RPCentrality);
a31b8a87 426
742ee86c 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);
a31b8a87 429
742ee86c 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);
742ee86c 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 }
a31b8a87 452
742ee86c 453 if(cJetBranch.Length()==0)fFlagJetType[ij] = 0;
454 if(!fFlagJetType[ij])continue;
7b822bc5 455
742ee86c 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]);
65c43de8 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]);
d7396b04 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]);
742ee86c 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
2eded560 471 fh1PtTracksInLow[ij] = new TH1F(Form("fh1PtTracks%sInLow",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),100,0.,5.);
742ee86c 472 fHistList->Add(fh1PtTracksInLow[ij]);
473
3170a3f8 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.);
742ee86c 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
3170a3f8 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);
742ee86c 481 fHistList->Add(fh2NTracksPt[ij]);
90a006c2 482
2eded560 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
556ccab2 489 // Bins: Jet number: pTJet, cent, mult, RP, Area, trigger, leading track pT total bins = 4.5M
7bee4353 490 const Int_t nBinsSparse1 = 8;
556ccab2 491 const Int_t nBinsLeadingTrackPt = 10;
7bee4353 492 Int_t nBins1[nBinsSparse1] = { kMaxJets+1,120, 10, 25, fNRPBins, 1,fNTrigger,nBinsLeadingTrackPt};
d01931c3 493 if(cJetBranch.Contains("RandomCone")){
494 nBins1[1] = 600;
495 nBins1[5] = 1;
496 }
7bee4353 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.};
2eded560 499
556ccab2 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
7bee4353 501
3d1dba29 502 fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area;trigger",nBinsSparse1,nBins1,xmin1,xmax1);
7bee4353 503 fhnJetPt[ij]->SetBinEdges(7,binArrayLeadingTrackPt);
2eded560 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;
d01931c3 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
90a006c2 523 // change the binning ot the pT axis:
d01931c3 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
2eded560 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
d01931c3 533 }
534
90a006c2 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;
30f922ad 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()};
90a006c2 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.;
e24d63c0 553 else if(xPt4[i-1]<30)xPt4[i] = xPt4[i-1] + 2.5;
90a006c2 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;
a31b8a87 560
37eb26ea 561 for(int i = 0;i <= kMaxJets;++i){
a31b8a87 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
d7117cbd 565
90a006c2 566 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
d7117cbd 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]);
a31b8a87 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
b5d90baf 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);
a31b8a87 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]);
a31b8a87 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 }
3b7ffecf 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
3b7ffecf 608 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
609
610}
611
a31b8a87 612void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
613
b5cc0c6d 614
f2dd0695 615 Bool_t selected = kTRUE;
616
617 if(fUseGlobalSelection&&fEventSelectionMask==0){
618 selected = AliAnalysisHelperJetTasks::Selected();
619 }
620 else if(fUseGlobalSelection&&fEventSelectionMask>0){
45a11aef 621 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
f2dd0695 622 }
623
f4132e7d 624 if(fEventClass>0){
625 selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
626 }
627
f2dd0695 628 if(!selected){
b5cc0c6d 629 // no selection by the service task, we continue
f4132e7d 630 if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
b5cc0c6d 631 PostData(1, fHistList);
632 return;
633 }
634
f2dd0695 635
a31b8a87 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
3b7ffecf 645 //
646 // Execute analysis for current event
647 //
a31b8a87 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){
d95fc15a 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;
a31b8a87 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());
3b7ffecf 704 return;
705 }
3b7ffecf 706 }
a31b8a87 707
37eb26ea 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
a31b8a87 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
742ee86c 760
a31b8a87 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
742ee86c 780 // Getting some global properties
781 fCentrality = GetCentrality();
aa3ba8d2 782 if(fCentrality<=0)fCentrality = 0;
742ee86c 783 fh1Centrality->Fill(fCentrality);
784
785
a31b8a87 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__);
3b7ffecf 792 return;
793 }
a31b8a87 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
d3a3f33d 821 // CalculateReactionPlaneAngle(&recParticles);
2eded560 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);
a31b8a87 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
784b1747 837 if(aod->GetPrimaryVertex()){// No vtx for pure MC
838 fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
839 }
a31b8a87 840
a2522483 841
37eb26ea 842 Int_t recMult1 = recParticles.GetEntries();
843 Int_t genMult1 = genParticles.GetEntries();
844
742ee86c 845 Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
846 Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
37eb26ea 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;
a2522483 854
441afcb9 855 Double_t var0[3] = {0,};
90a006c2 856 var0[0] = fCentrality;
857 var0[1] = fMultRec;
3d1dba29 858 for(int it=0;it<fNTrigger;it++){
859 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
860 var0[2] = it;
861 fhnEvent->Fill(var0);
862 }
863 }
a31b8a87 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);
a2522483 872 FillJetHistos(recJetsList,recParticles,kJetRecFull);
a31b8a87 873 FillTrackHistos(recParticles,kJetRec);
874
875 FillJetHistos(genJetsListCut,genParticles,kJetGen);
a2522483 876 FillJetHistos(genJetsList,genParticles,kJetGenFull);
a31b8a87 877 FillTrackHistos(genParticles,kJetGen);
878
879 // Here follows the jet matching part
880 // delegate to a separated method?
881
b99018c6 882 if(fDoMatching){
883 FillMatchHistos(recJetsList,genJetsList);
884 }
885
a31b8a87 886 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
887 PostData(1, fHistList);
90a006c2 888}
a31b8a87 889
890void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
891
892 if(iType>=kJetTypes){
893 return;
3b7ffecf 894 }
742ee86c 895 if(!fFlagJetType[iType])return;
a31b8a87 896
37eb26ea 897 Int_t refMult = fMultRec;
898 if(iType==kJetGen||iType==kJetGenFull){
899 refMult = fMultGen;
900 }
901
a31b8a87 902 Int_t nJets = jetsList.GetEntries();
903 fh1NJets[iType]->Fill(nJets);
904
905 if(nJets<=0)return;
3b7ffecf 906
a31b8a87 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
7bee4353 915 Double_t var1[8] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leadingTrackPt
916
90a006c2 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;
a31b8a87 922
923 for(int ij = 0;ij < nJets;ij++){
924 AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
925 Float_t ptJet = jet->Pt();
4a828dd2 926 if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0);
d7396b04 927 if(jet->Trigger()&fJetTriggerBestMask){
928 fh1PtJetsInBest[iType]->Fill(ptJet);
929 }
65c43de8 930 if(jet->Trigger()&fJetTriggerExcludeMask){
931 fh1PtJetsInRej[iType]->Fill(ptJet);
932 continue;
933 }
a31b8a87 934 fh1PtJetsIn[iType]->Fill(ptJet);
a31b8a87 935 ptOld = ptJet;
936
b5d90baf 937 // find the dijets assume sorting and acceptance cut...
a31b8a87 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();
b5d90baf 947 if(phi1<0)phi1+=TMath::Pi()*2.;
a31b8a87 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
37eb26ea 957
a31b8a87 958 Float_t phiJet = jet->Phi();
f12de05f 959 Float_t etaJet = jet->Eta();
a31b8a87 960 if(phiJet<0)phiJet+=TMath::Pi()*2.;
961 fh1TmpRho->Reset();
37eb26ea 962 if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
90a006c2 963
37eb26ea 964 fh1PtIn[iType][kMaxJets]->Fill(ptJet);
a31b8a87 965 // fill leading jets...
d7117cbd 966 AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
967 // AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
742ee86c 968 Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
90a006c2 969
970 var1[1] = ptJet;
971 var1[4] = phiBin;
972 var1[5] = jet->EffectiveAreaCharged();
7bee4353 973 var1[7] = leadTrack->Pt();//pT of leading jet
90a006c2 974
975 var2[1] = ptJet;
976 var2[3] = etaJet;
977 var2[4] = phiJet;
978 var2[5] = jet->EffectiveAreaCharged();
37eb26ea 979 if(ij<kMaxJets){
d7117cbd 980 if(leadTrack)fh2LTrackPtJetPt[iType][ij]->Fill(leadTrack->Pt(),ptJet);
90a006c2 981 var1[0] = ij;
982 var2[0] = ij;
3d1dba29 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 }
90a006c2 989 fhnJetPtQA[iType]->Fill(var2);
a31b8a87 990 }
90a006c2 991 var1[0] = kMaxJets;// fill for all jets
992 var2[0] = kMaxJets;// fill for all jets
3d1dba29 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
90a006c2 1000 fhnJetPtQA[iType]->Fill(var2);
d7117cbd 1001 if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet);
a31b8a87 1002
37eb26ea 1003 if(particlesList.GetSize()&&ij<kMaxJets){
a31b8a87 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
a31b8a87 1011 }// if we have particles
a31b8a87 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);
b5d90baf 1031 fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
a31b8a87 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());
a31b8a87 1054 if(tmpJet){
e855f5c5 1055 Float_t pt = tmpJet->Pt();
a31b8a87 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
90a006c2 1075 if(fFlagJetType[iType]<=0)return;
7b822bc5 1076 Int_t refMult = fMultRec;
1077 if(iType==kJetGen||iType==kJetGenFull){
1078 refMult = fMultGen;
742ee86c 1079
7b822bc5 1080 }
1081
90a006c2 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;
a31b8a87 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);
cb883243 1115 fh1PtTracksInLow[iType]->Fill(tmpPt);
742ee86c 1116
a31b8a87 1117 sumPt += tmpPt;
1118 Float_t tmpPhi = tmpTrack->Phi();
1119 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
2eded560 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
dae7dd67 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 }
742ee86c 1135 Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
90a006c2 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
a31b8a87 1149 if(tmpTrack==leading){
90a006c2 1150 var3[0] = 0;
1151 var4[0] = 0;
1152 fhnTrackPt[iType]->Fill(var3);
1153 fhnTrackPtQA[iType]->Fill(var4);
a31b8a87 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());
3b7ffecf 1174
a31b8a87 1175 static TArrayI aRecIndex(genJetsList.GetEntries());
1176 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
3b7ffecf 1177
a31b8a87 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];
a31b8a87 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;
b99018c6 1208 fhnJetContainer->Fill(&container[3],kStep0);
a31b8a87 1209 if(JetSelected(genJet)){
b99018c6 1210 fhnJetContainer->Fill(&container[3],kStep1);
a31b8a87 1211 Int_t ir = aRecIndex[ig];
1212 if(ir>=0&&ir<recJetsList.GetEntries()){
b99018c6 1213 fhnJetContainer->Fill(&container[3],kStep2);
a31b8a87 1214 AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
b99018c6 1215 if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep4);
1216 if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep3);
a31b8a87 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;
a31b8a87 1235
b99018c6 1236 fhnJetContainer->Fill(container,kStep0+kMaxStep);
a31b8a87 1237 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1238
1239 if(JetSelected(recJet)){
b99018c6 1240 fhnJetContainer->Fill(container,kStep1+kMaxStep);
a31b8a87 1241 // Fill Correlation
1242 Int_t ig = aGenIndex[ir];
1243 if(ig>=0 && ig<genJetsList.GetEntries()){
b99018c6 1244 fhnJetContainer->Fill(container,kStep2+kMaxStep);
a31b8a87 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;
a31b8a87 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 //
b99018c6 1258 if(JetSelected(genJet))fhnJetContainer->Fill(container,kStep4+kMaxStep);
1259 fhnJetContainer->Fill(container,kStep3+kMaxStep);
1260 fhnCorrelation->Fill(container,0);
a31b8a87 1261 if(ptGen>0){
1262 Float_t delta = (ptRec-ptGen)/ptGen;
1263 fh2RelPtFGen->Fill(ptGen,delta);
1264 fh2PtFGen->Fill(ptGen,ptRec);
1265 }
a31b8a87 1266 }
a31b8a87 1267 }// loop over reconstructed jets
1268 }
1269 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1270}
1271
1272
3b7ffecf 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
d95fc15a 1279 const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
3b7ffecf 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
3b7ffecf 1288 //arrays for the number of bins in each dimension
1289 Int_t iBin[kNvar];
d95fc15a 1290 iBin[0] = 125; //bins in pt
3b7ffecf 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);
a923bd34 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;
3b7ffecf 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
b99018c6 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]);
3b7ffecf 1309 }
b99018c6 1310
3b7ffecf 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
b99018c6 1320 fhnCorrelation = new AliTHn("fahnCorrelation","AliTHn with correlations",1,2*kNvar,thnDim);
3b7ffecf 1321 for (int k=0; k<kNvar; k++) {
b99018c6 1322 fhnCorrelation->SetBinLimits(k,binEdges[k]);
1323 fhnCorrelation->SetBinLimits(k+kNvar,binEdges[k]);
3b7ffecf 1324 }
742ee86c 1325
1326 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1327 delete [] binEdges[ivar];
1328
1329
3b7ffecf 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
565584e8 1342 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
3b7ffecf 1343
1344 Int_t iCount = 0;
565584e8 1345 if(type==kTrackAOD){
3b7ffecf 1346 AliAODEvent *aod = 0;
565584e8 1347 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1348 else aod = AODEvent();
3b7ffecf 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;
a31b8a87 1355 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1356 if(tr->Pt()<fMinTrackPt)continue;
38ecb6a5 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());
fd5db301 1360 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
38ecb6a5 1361 tr->Print();
1362 // tr->Dump();
1363 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1364 if(fESD){
a2522483 1365 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
38ecb6a5 1366 esdTr->Print("");
bb3d13aa 1367 // esdTr->Dump();
38ecb6a5 1368 }
1369 }
1370 }
3b7ffecf 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);
a31b8a87 1382 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 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 }
565584e8 1394 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1395 AliAODEvent *aod = 0;
5301f754 1396 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
565584e8 1397 else aod = AODEvent();
5010a3f7 1398 if(!aod)return iCount;
3b7ffecf 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));
e855f5c5 1403 if(!part)continue;
a31b8a87 1404 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 1405 if(!part->IsPhysicalPrimary())continue;
c4631cdb 1406 if(type == kTrackAODMCAll){
3b7ffecf 1407 list->Add(part);
1408 iCount++;
1409 }
565584e8 1410 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
3b7ffecf 1411 if(part->Charge()==0)continue;
565584e8 1412 if(kTrackAODMCCharged){
1413 list->Add(part);
1414 }
1415 else {
a31b8a87 1416 if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
565584e8 1417 list->Add(part);
1418 }
3b7ffecf 1419 iCount++;
1420 }
1421 }
1422 }// AODMCparticle
cc0649e4 1423 list->Sort();
c4631cdb 1424 return iCount;
3b7ffecf 1425
1426}
a31b8a87 1427
1428
742ee86c 1429Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1430 AliAODEvent *aod = 0;
1431 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1432 else aod = AODEvent();
1433 if(!aod){
3170a3f8 1434 return 101;
742ee86c 1435 }
1436 return aod->GetHeader()->GetCentrality();
1437}
1438
1439
a31b8a87 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
37eb26ea 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}
d7117cbd 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}
742ee86c 1537
742ee86c 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
3d1dba29 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}