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