]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/JetTasks/AliAnalysisTaskJetSpectrum2.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
updated jet name
[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
ed5c9f3f 352 const Int_t nBinPt = 120;
3b7ffecf 353 Double_t binLimitsPt[nBinPt+1];
354 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
355 if(iPt == 0){
ed5c9f3f 356 binLimitsPt[iPt] = -50.0;
3b7ffecf 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;
ed5c9f3f 478 const Int_t nBins1[nBinsSparse1] = { kMaxJets+1,120, 10, 25, fNRPBins, 10};
479 const Double_t xmin1[nBinsSparse1] = { -0.5,-50, 0, 0, -0.5, 0.};
2eded560 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();
ed5c9f3f 897 if(ptJet<1.)jet->GetPtSubtracted(0);
65c43de8 898 if(jet->Trigger()&fJetTriggerExcludeMask){
899 fh1PtJetsInRej[iType]->Fill(ptJet);
900 continue;
901 }
a31b8a87 902 fh1PtJetsIn[iType]->Fill(ptJet);
903 if(ptJet>ptOld){
904 Printf("%s:%d Jets Type %d Not Sorted !! %d:%.3E %d:%.3E",(char*)__FILE__,__LINE__,iType,ij,ptJet,ij-1,ptOld);
905 }
906 ptOld = ptJet;
907
b5d90baf 908 // find the dijets assume sorting and acceptance cut...
a31b8a87 909 if(ij==0){
910 ij0 = ij;
911 pT0 = ptJet;
912 phi0 = jet->Phi();
913 if(phi0<0)phi0+=TMath::Pi()*2.;
914 }
915 else if(ptJet>pT1){
916 // take only the backward hemisphere??
917 phi1 = jet->Phi();
b5d90baf 918 if(phi1<0)phi1+=TMath::Pi()*2.;
a31b8a87 919 Float_t dPhi = phi1 - phi0;
920 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
921 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
922 if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
923 ij1 = ij;
924 pT1 = ptJet;
925 }
926 }
927 // fill jet histos for kmax jets
37eb26ea 928
a31b8a87 929 Float_t phiJet = jet->Phi();
f12de05f 930 Float_t etaJet = jet->Eta();
a31b8a87 931 if(phiJet<0)phiJet+=TMath::Pi()*2.;
932 fh1TmpRho->Reset();
37eb26ea 933 if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
90a006c2 934
37eb26ea 935 fh1PtIn[iType][kMaxJets]->Fill(ptJet);
a31b8a87 936 // fill leading jets...
d7117cbd 937 AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
938 // AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
742ee86c 939 Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
90a006c2 940
941 var1[1] = ptJet;
942 var1[4] = phiBin;
943 var1[5] = jet->EffectiveAreaCharged();
944
945 var2[1] = ptJet;
946 var2[3] = etaJet;
947 var2[4] = phiJet;
948 var2[5] = jet->EffectiveAreaCharged();
37eb26ea 949 if(ij<kMaxJets){
d7117cbd 950 if(leadTrack)fh2LTrackPtJetPt[iType][ij]->Fill(leadTrack->Pt(),ptJet);
90a006c2 951 var1[0] = ij;
952 var2[0] = ij;
953 fhnJetPt[iType]->Fill(var1);
954 fhnJetPtQA[iType]->Fill(var2);
a31b8a87 955 }
90a006c2 956 var1[0] = kMaxJets;// fill for all jets
957 var2[0] = kMaxJets;// fill for all jets
958 fhnJetPt[iType]->Fill(var1);
959 fhnJetPtQA[iType]->Fill(var2);
d7117cbd 960 if(leadTrack)fh2LTrackPtJetPt[iType][kMaxJets]->Fill(leadTrack->Pt(),ptJet);
a31b8a87 961
37eb26ea 962 if(particlesList.GetSize()&&ij<kMaxJets){
a31b8a87 963 // Particles... correlated with jets...
964 for(int it = 0;it<particlesList.GetEntries();++it){
965 AliVParticle *part = (AliVParticle*)particlesList.At(it);
966 Float_t deltaR = jet->DeltaR(part);
967 if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
968 }
969 // fill the jet shapes
a31b8a87 970 }// if we have particles
a31b8a87 971 }// Jet Loop
972
973
974 // do something with dijets...
975 if(ij0>=0&&ij1>0){
976 AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
977 Double_t ptJet0 = jet0->Pt();
978 Double_t phiJet0 = jet0->Phi();
979 if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
980
981 AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
982 Double_t ptJet1 = jet1->Pt();
983 Double_t phiJet1 = jet1->Phi();
984 if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
985
986 Float_t deltaPhi = phiJet0 - phiJet1;
987 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
988 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
989 deltaPhi = TMath::Abs(deltaPhi);
b5d90baf 990 fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
a31b8a87 991
992 Float_t asym = 9999;
993 if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
994 fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
995 fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
996 fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
997 Float_t minv = 2.*(jet0->P()*jet1->P()-
998 jet0->Px()*jet1->Px()-
999 jet0->Py()*jet1->Py()-
1000 jet0->Pz()*jet1->Pz()); // assume mass == 0;
1001 if(minv<0)minv=0; // prevent numerical instabilities
1002 minv = TMath::Sqrt(minv);
1003 fh1DijetMinv[iType]->Fill(minv);
1004 }
1005
1006
1007
1008 // count down the jets above thrueshold
1009 Int_t nOver = nJets;
1010 if(nOver>0){
1011 TIterator *jetIter = jetsList.MakeIterator();
1012 AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
a31b8a87 1013 if(tmpJet){
e855f5c5 1014 Float_t pt = tmpJet->Pt();
a31b8a87 1015 for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1016 Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1017 while(pt<ptCut&&tmpJet){
1018 nOver--;
1019 tmpJet = (AliAODJet*)(jetIter->Next());
1020 if(tmpJet){
1021 pt = tmpJet->Pt();
1022 }
1023 }
1024 if(nOver<=0)break;
1025 fh2NJetsPt[iType]->Fill(ptCut,nOver);
1026 }
1027 }
1028 delete jetIter;
1029 }
1030}
1031
1032void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1033
90a006c2 1034 if(fFlagJetType[iType]<=0)return;
7b822bc5 1035 Int_t refMult = fMultRec;
1036 if(iType==kJetGen||iType==kJetGenFull){
1037 refMult = fMultGen;
742ee86c 1038
7b822bc5 1039 }
1040
90a006c2 1041 //
1042 Double_t var3[5]; // track number;p_{T};cent;#tracks;RP
1043 var3[2] = fCentrality;
1044 var3[3] = refMult;
1045 Double_t var4[5]; // track number;p_{T};cent;#eta;#phi
1046 var4[2] = fCentrality;
a31b8a87 1047 Int_t nTrackOver = particlesList.GetSize();
1048 // do the same for tracks and jets
1049 if(nTrackOver>0){
1050 TIterator *trackIter = particlesList.MakeIterator();
1051 AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1052 Float_t pt = tmpTrack->Pt();
1053 for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1054 Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1055 while(pt<ptCut&&tmpTrack){
1056 nTrackOver--;
1057 tmpTrack = (AliVParticle*)(trackIter->Next());
1058 if(tmpTrack){
1059 pt = tmpTrack->Pt();
1060 }
1061 }
1062 if(nTrackOver<=0)break;
1063 fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1064 }
1065
1066
1067 trackIter->Reset();
1068 AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1069 Float_t sumPt = 0;
1070
1071 while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1072 Float_t tmpPt = tmpTrack->Pt();
1073 fh1PtTracksIn[iType]->Fill(tmpPt);
cb883243 1074 fh1PtTracksInLow[iType]->Fill(tmpPt);
742ee86c 1075
a31b8a87 1076 sumPt += tmpPt;
1077 Float_t tmpPhi = tmpTrack->Phi();
1078 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
2eded560 1079
1080
1081 Float_t phiRP = tmpPhi-fRPAngle;
1082 if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1083 if(phiRP<0)phiRP += TMath::Pi();
1084 if(phiRP<0)phiRP += TMath::Pi();
1085 const float allPhi = -1./180.*TMath::Pi();
1086
dae7dd67 1087 if(tmpPt<100){
1088 fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1089 fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1090
1091 fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1092 fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1093 }
742ee86c 1094 Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
90a006c2 1095 var3[0] = 1;
1096 var3[1] = tmpPt;
1097 var3[4] = phiBin;
1098
1099 var4[0] = 1;
1100 var4[1] = tmpPt;
1101 var4[3] = tmpTrack->Eta();
1102 var4[4] = tmpPhi;
1103
1104
1105 fhnTrackPt[iType]->Fill(var3);
1106 fhnTrackPtQA[iType]->Fill(var4);
1107
a31b8a87 1108 if(tmpTrack==leading){
90a006c2 1109 var3[0] = 0;
1110 var4[0] = 0;
1111 fhnTrackPt[iType]->Fill(var3);
1112 fhnTrackPtQA[iType]->Fill(var4);
a31b8a87 1113 continue;
1114 }
1115 }
1116 fh1SumPtTrack[iType]->Fill(sumPt);
1117 delete trackIter;
1118 }
1119
1120}
1121
1122
1123void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1124
1125
1126 // Fill al the matching histos
1127 // It is important that the acceptances for the mathing are as large as possible
1128 // to avoid false matches on the edge of acceptance
1129 // therefore we add some extra matching jets as overhead
1130
1131 static TArrayI aGenIndex(recJetsList.GetEntries());
1132 if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
3b7ffecf 1133
a31b8a87 1134 static TArrayI aRecIndex(genJetsList.GetEntries());
1135 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
3b7ffecf 1136
a31b8a87 1137 if(fDebug){
1138 Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1139 Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1140 }
1141 AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1142 &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1143 aGenIndex,aRecIndex,fDebug);
1144
1145 if(fDebug){
1146 for(int i = 0;i< aGenIndex.GetSize();++i){
1147 if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1148 }
1149 for(int i = 0;i< aRecIndex.GetSize();++i){
1150 if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1151 }
1152 }
1153
1154 Double_t container[6];
a31b8a87 1155
1156 // loop over generated jets
1157 // consider the
1158 for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1159 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1160 Double_t ptGen = genJet->Pt();
1161 Double_t phiGen = genJet->Phi();
1162 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1163 Double_t etaGen = genJet->Eta();
1164 container[3] = ptGen;
1165 container[4] = etaGen;
1166 container[5] = phiGen;
b99018c6 1167 fhnJetContainer->Fill(&container[3],kStep0);
a31b8a87 1168 if(JetSelected(genJet)){
b99018c6 1169 fhnJetContainer->Fill(&container[3],kStep1);
a31b8a87 1170 Int_t ir = aRecIndex[ig];
1171 if(ir>=0&&ir<recJetsList.GetEntries()){
b99018c6 1172 fhnJetContainer->Fill(&container[3],kStep2);
a31b8a87 1173 AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
b99018c6 1174 if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep4);
1175 if(JetSelected(recJet))fhnJetContainer->Fill(&container[3],kStep3);
a31b8a87 1176 }
1177 }
1178 }// loop over generated jets used for matching...
1179
1180
1181
1182 // loop over reconstructed jets
1183 for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1184 AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1185 Double_t etaRec = recJet->Eta();
1186 Double_t ptRec = recJet->Pt();
1187 Double_t phiRec = recJet->Phi();
1188 if(phiRec<0)phiRec+=TMath::Pi()*2.;
1189 // do something with dijets...
1190
1191 container[0] = ptRec;
1192 container[1] = etaRec;
1193 container[2] = phiRec;
a31b8a87 1194
b99018c6 1195 fhnJetContainer->Fill(container,kStep0+kMaxStep);
a31b8a87 1196 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1197
1198 if(JetSelected(recJet)){
b99018c6 1199 fhnJetContainer->Fill(container,kStep1+kMaxStep);
a31b8a87 1200 // Fill Correlation
1201 Int_t ig = aGenIndex[ir];
1202 if(ig>=0 && ig<genJetsList.GetEntries()){
b99018c6 1203 fhnJetContainer->Fill(container,kStep2+kMaxStep);
a31b8a87 1204 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1205 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1206 Double_t ptGen = genJet->Pt();
1207 Double_t phiGen = genJet->Phi();
1208 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1209 Double_t etaGen = genJet->Eta();
1210
1211 container[3] = ptGen;
1212 container[4] = etaGen;
1213 container[5] = phiGen;
a31b8a87 1214 //
1215 // we accept only jets which are detected within a smaller window, to avoid ambigious pair association at the edges of the acceptance
1216 //
b99018c6 1217 if(JetSelected(genJet))fhnJetContainer->Fill(container,kStep4+kMaxStep);
1218 fhnJetContainer->Fill(container,kStep3+kMaxStep);
1219 fhnCorrelation->Fill(container,0);
a31b8a87 1220 if(ptGen>0){
1221 Float_t delta = (ptRec-ptGen)/ptGen;
1222 fh2RelPtFGen->Fill(ptGen,delta);
1223 fh2PtFGen->Fill(ptGen,ptRec);
1224 }
a31b8a87 1225 }
a31b8a87 1226 }// loop over reconstructed jets
1227 }
1228 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1229}
1230
1231
3b7ffecf 1232void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1233 //
1234 // Create the particle container for the correction framework manager and
1235 // link it
1236 //
1237 const Int_t kNvar = 3 ; //number of variables on the grid:pt,eta, phi
d95fc15a 1238 const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
3b7ffecf 1239 const Double_t kEtamin = -3.0, kEtamax = 3.0;
1240 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
1241
1242 // can we neglect migration in eta and phi?
1243 // phi should be no problem since we cover full phi and are phi symmetric
1244 // eta migration is more difficult due to needed acceptance correction
1245 // in limited eta range
1246
3b7ffecf 1247 //arrays for the number of bins in each dimension
1248 Int_t iBin[kNvar];
d95fc15a 1249 iBin[0] = 125; //bins in pt
3b7ffecf 1250 iBin[1] = 1; //bins in eta
1251 iBin[2] = 1; // bins in phi
1252
1253 //arrays for lower bounds :
1254 Double_t* binEdges[kNvar];
1255 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1256 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1257
1258 //values for bin lower bounds
1259 // 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 1260 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 1261 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1262 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
1263
1264
b99018c6 1265 fhnJetContainer = new AliTHn(Form("fahnJetContainer"),Form("AliTHn jet info"),kMaxStep*2,kNvar,iBin);
1266 for (int k=0; k<kNvar; k++) {
1267 fhnJetContainer->SetBinLimits(k,binEdges[k]);
3b7ffecf 1268 }
b99018c6 1269
3b7ffecf 1270 //create correlation matrix for unfolding
1271 Int_t thnDim[2*kNvar];
1272 for (int k=0; k<kNvar; k++) {
1273 //first half : reconstructed
1274 //second half : MC
1275 thnDim[k] = iBin[k];
1276 thnDim[k+kNvar] = iBin[k];
1277 }
1278
b99018c6 1279 fhnCorrelation = new AliTHn("fahnCorrelation","AliTHn with correlations",1,2*kNvar,thnDim);
3b7ffecf 1280 for (int k=0; k<kNvar; k++) {
b99018c6 1281 fhnCorrelation->SetBinLimits(k,binEdges[k]);
1282 fhnCorrelation->SetBinLimits(k+kNvar,binEdges[k]);
3b7ffecf 1283 }
742ee86c 1284
1285 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1286 delete [] binEdges[ivar];
1287
1288
3b7ffecf 1289}
1290
1291void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1292{
1293// Terminate analysis
1294//
1295 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
1296}
1297
1298
1299Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1300
565584e8 1301 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
3b7ffecf 1302
1303 Int_t iCount = 0;
565584e8 1304 if(type==kTrackAOD){
3b7ffecf 1305 AliAODEvent *aod = 0;
565584e8 1306 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1307 else aod = AODEvent();
3b7ffecf 1308 if(!aod){
1309 return iCount;
1310 }
1311 for(int it = 0;it < aod->GetNumberOfTracks();++it){
1312 AliAODTrack *tr = aod->GetTrack(it);
1313 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
a31b8a87 1314 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1315 if(tr->Pt()<fMinTrackPt)continue;
38ecb6a5 1316 if(fDebug>0){
1317 if(tr->Pt()>20){
1318 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
fd5db301 1319 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
38ecb6a5 1320 tr->Print();
1321 // tr->Dump();
1322 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1323 if(fESD){
a2522483 1324 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
38ecb6a5 1325 esdTr->Print("");
bb3d13aa 1326 // esdTr->Dump();
38ecb6a5 1327 }
1328 }
1329 }
3b7ffecf 1330 list->Add(tr);
1331 iCount++;
1332 }
1333 }
1334 else if (type == kTrackKineAll||type == kTrackKineCharged){
1335 AliMCEvent* mcEvent = MCEvent();
1336 if(!mcEvent)return iCount;
1337 // we want to have alivpartilces so use get track
1338 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1339 if(!mcEvent->IsPhysicalPrimary(it))continue;
1340 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
a31b8a87 1341 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 1342 if(type == kTrackKineAll){
1343 list->Add(part);
1344 iCount++;
1345 }
1346 else if(type == kTrackKineCharged){
1347 if(part->Particle()->GetPDG()->Charge()==0)continue;
1348 list->Add(part);
1349 iCount++;
1350 }
1351 }
1352 }
565584e8 1353 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1354 AliAODEvent *aod = 0;
5301f754 1355 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
565584e8 1356 else aod = AODEvent();
5010a3f7 1357 if(!aod)return iCount;
3b7ffecf 1358 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1359 if(!tca)return iCount;
1360 for(int it = 0;it < tca->GetEntriesFast();++it){
1361 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
e855f5c5 1362 if(!part)continue;
a31b8a87 1363 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 1364 if(!part->IsPhysicalPrimary())continue;
c4631cdb 1365 if(type == kTrackAODMCAll){
3b7ffecf 1366 list->Add(part);
1367 iCount++;
1368 }
565584e8 1369 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
3b7ffecf 1370 if(part->Charge()==0)continue;
565584e8 1371 if(kTrackAODMCCharged){
1372 list->Add(part);
1373 }
1374 else {
a31b8a87 1375 if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
565584e8 1376 list->Add(part);
1377 }
3b7ffecf 1378 iCount++;
1379 }
1380 }
1381 }// AODMCparticle
cc0649e4 1382 list->Sort();
c4631cdb 1383 return iCount;
3b7ffecf 1384
1385}
a31b8a87 1386
1387
742ee86c 1388Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
1389 AliAODEvent *aod = 0;
1390 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1391 else aod = AODEvent();
1392 if(!aod){
3170a3f8 1393 return 101;
742ee86c 1394 }
1395 return aod->GetHeader()->GetCentrality();
1396}
1397
1398
a31b8a87 1399
1400Bool_t AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1401 Bool_t selected = false;
1402
1403 if(!jet)return selected;
1404
1405 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1406 selected = kTRUE;
1407 }
1408 return selected;
1409
1410}
1411
1412Int_t AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1413
1414 if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1415 Int_t iCount = 0;
1416
1417 if(!jarray){
1418 Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1419 return iCount;
1420 }
1421
1422
1423 for(int ij=0;ij<jarray->GetEntries();ij++){
1424 AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1425 if(!jet)continue;
1426 if(type==0){
1427 // No cut at all, main purpose here is sorting
1428 list->Add(jet);
1429 iCount++;
1430 }
1431 else if(type == 1){
1432 // eta cut
1433 if(JetSelected(jet)){
1434 list->Add(jet);
1435 iCount++;
1436 }
1437 }
1438 }
1439
1440 list->Sort();
1441 return iCount;
1442
1443}
1444
1445
37eb26ea 1446Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1447 if(!jets)return 0;
1448
1449 Int_t refMult = 0;
1450 for(int ij = 0;ij < jets->GetEntries();++ij){
1451 AliAODJet* jet = (AliAODJet*)jets->At(ij);
1452 if(!jet)continue;
1453 TRefArray *refs = jet->GetRefTracks();
1454 if(!refs)continue;
1455 refMult += refs->GetEntries();
1456 }
1457 return refMult;
1458
1459}
d7117cbd 1460
1461
1462AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1463 if(!jet)return 0;
1464 TRefArray *refs = jet->GetRefTracks();
1465 if(!refs) return 0;
1466 AliVParticle *leading = 0;
1467 Float_t fMaxPt = 0;
1468 for(int i = 0;i<refs->GetEntries();i++){
1469 AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1470 if(!tmp)continue;
1471 if(tmp->Pt()>fMaxPt){
1472 leading = tmp;
1473 fMaxPt = tmp->Pt();
1474 }
1475 }
1476 return leading;
1477}
1478
1479
1480AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1481 if(!jet)return 0;
1482 if(!list) return 0;
1483 AliVParticle *leading = 0;
1484 Float_t fMaxPt = 0;
1485 for(int i = 0;i<list->GetEntries();i++){
1486 AliVParticle *tmp = (AliVParticle*)(list->At(i));
1487 if(!tmp)continue;
1488 if(jet->DeltaR(tmp)>r)continue;
1489 if(tmp->Pt()>fMaxPt){
1490 leading = tmp;
1491 fMaxPt = tmp->Pt();
1492 }
1493 }
1494 return leading;
1495}
742ee86c 1496
742ee86c 1497Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1498{
1499 Int_t phibin=-1;
1500 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1501 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1502 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1503 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1504 return phibin;
1505}
1506