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