Adapt add macro and particle/cluster containers for the analysis on jets
[u/mrichter/AliRoot.git] / PWGJE / AliAnalysisTaskJetSpectrum2.cxx
CommitLineData
3b7ffecf 1// **************************************
3170a3f8 2// used for the correction of determiantion of reconstructed jet spectra
3b7ffecf 3// Compares input (gen) and output (rec) jets
4// *******************************************
5
6
7/**************************************************************************
8 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
9 * *
10 * Author: The ALICE Off-line Project. *
11 * Contributors are mentioned in the code where appropriate. *
12 * *
13 * Permission to use, copy, modify and distribute this software and its *
14 * documentation strictly for non-commercial purposes is hereby granted *
15 * without fee, provided that the above copyright notice appears in all *
16 * copies and that both the copyright notice and this permission notice *
17 * appear in the supporting documentation. The authors make no claims *
18 * about the suitability of this software for any purpose. It is *
19 * provided "as is" without express or implied warranty. *
20 **************************************************************************/
3b7ffecf 21
22#include <TROOT.h>
23#include <TRandom.h>
742ee86c 24#include <TRandom3.h>
3b7ffecf 25#include <TSystem.h>
26#include <TInterpreter.h>
27#include <TChain.h>
28#include <TFile.h>
29#include <TKey.h>
30#include <TH1F.h>
31#include <TH2F.h>
32#include <TH3F.h>
33#include <TProfile.h>
2eded560 34#include <TProfile2D.h>
3b7ffecf 35#include <TList.h>
36#include <TLorentzVector.h>
37#include <TClonesArray.h>
37eb26ea 38#include <TRefArray.h>
3b7ffecf 39#include "TDatabasePDG.h"
40
41#include "AliAnalysisTaskJetSpectrum2.h"
42#include "AliAnalysisManager.h"
43#include "AliJetFinder.h"
b99018c6 44#include "AliTHn.h"
3b7ffecf 45#include "AliJetHeader.h"
46#include "AliJetReader.h"
47#include "AliJetReaderHeader.h"
48#include "AliUA1JetHeaderV1.h"
49#include "AliESDEvent.h"
50#include "AliAODEvent.h"
51#include "AliAODHandler.h"
52#include "AliAODTrack.h"
53#include "AliAODJet.h"
c2785065 54#include "AliAODJetEventBackground.h"
3b7ffecf 55#include "AliAODMCParticle.h"
56#include "AliMCEventHandler.h"
57#include "AliMCEvent.h"
58#include "AliStack.h"
59#include "AliGenPythiaEventHeader.h"
60#include "AliJetKineReaderHeader.h"
61#include "AliGenCocktailEventHeader.h"
62#include "AliInputEventHandler.h"
152a0744 63#include "AliCFContainer.h"
3b7ffecf 64
65#include "AliAnalysisHelperJetTasks.h"
66
67ClassImp(AliAnalysisTaskJetSpectrum2)
68
37eb26ea 69AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2():
adf3920d 70AliAnalysisTaskSE(),
a31b8a87 71 fJetHeaderRec(0x0),
72 fJetHeaderGen(0x0),
73 fAODIn(0x0),
74 fAODOut(0x0),
75 fAODExtension(0x0),
b99018c6 76 fhnJetContainer(0x0),
a31b8a87 77 fhnCorrelation(0x0),
90a006c2 78 fhnEvent(0x0),
c8eabe24 79 f1PtScale(0x0),
a31b8a87 80 fBranchRec("jets"),
81 fBranchGen(""),
37eb26ea 82 fBranchBkgRec(""),
83 fBranchBkgGen(""),
a31b8a87 84 fNonStdFile(""),
742ee86c 85 fRandomizer(0x0),
a31b8a87 86 fUseAODJetInput(kFALSE),
87 fUseAODTrackInput(kFALSE),
88 fUseAODMCInput(kFALSE),
89 fUseGlobalSelection(kFALSE),
90 fUseExternalWeightOnly(kFALSE),
91 fLimitGenJetEta(kFALSE),
b99018c6 92 fDoMatching(kFALSE),
a31b8a87 93 fNMatchJets(5),
7b822bc5 94 fNRPBins(3),
4f6d1cc1 95 fTRP(0),
96 fDebug(0),
65c43de8 97 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
d7396b04 98 fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
a31b8a87 99 fFilterMask(0),
f2dd0695 100 fEventSelectionMask(0),
3d1dba29 101 fNTrigger(0),
102 fTriggerBit(0x0),
db745d29 103 fNAcceptance(0),
85b5b73e 104 fNBinsLeadingTrackPt(10),
105 fNBinsMult(20),
a31b8a87 106 fAnalysisType(0),
3b7ffecf 107 fTrackTypeRec(kTrackUndef),
108 fTrackTypeGen(kTrackUndef),
f4132e7d 109 fEventClass(0),
2eded560 110 fRPMethod(0),
3b7ffecf 111 fAvgTrials(1),
112 fExternalWeight(1),
a31b8a87 113 fJetRecEtaWindow(0.5),
114 fTrackRecEtaWindow(0.5),
9280dfa6 115 fMinJetPt(0),
a31b8a87 116 fMinTrackPt(0.15),
117 fDeltaPhiWindow(90./180.*TMath::Pi()),
69a0d320 118 fAcceptancePhiMin(0x0),
119 fAcceptancePhiMax(0x0),
120 fAcceptanceEtaMin(0x0),
121 fAcceptanceEtaMax(0x0),
742ee86c 122 fCentrality(100),
7b822bc5 123 fRPAngle(0),
37eb26ea 124 fMultRec(0),
125 fMultGen(0),
3d1dba29 126 fTriggerName(0x0),
3b7ffecf 127 fh1Xsec(0x0),
128 fh1Trials(0x0),
432abb2b 129 fh1AvgTrials(0x0),
3b7ffecf 130 fh1PtHard(0x0),
131 fh1PtHardNoW(0x0),
132 fh1PtHardTrials(0x0),
57ca1193 133 fh1ZVtx(0x0),
7b822bc5 134 fh1RP(0x0),
742ee86c 135 fh1Centrality(0x0),
c8eabe24 136 fh1TmpRho(0x0),
37eb26ea 137 fh2MultRec(0x0),
138 fh2MultGen(0x0),
742ee86c 139 fh2RPCentrality(0x0),
a31b8a87 140 fh2PtFGen(0x0),
e39a4f22 141 fh2deltaPt1Pt2(0x0),
9a83d4af 142 fh2RelPtFGen(0x0),
346e5a93 143 fh3RelPtFGenLeadTrkPt(0x0),
0c964c83 144 fh1EvtSelection(0),
145 fMaxVertexZ(100.),
146 fMinNcontributors(0),
a2d2ff89 147 fRejectPileup(0),
148 fHistList(0x0)
3b7ffecf 149{
c17f867b 150
a31b8a87 151 for(int ij = 0;ij <kJetTypes;++ij){
742ee86c 152 fFlagJetType[ij] = 1; // default = on
a31b8a87 153 fh1NJets[ij] = 0;
154 fh1SumPtTrack[ij] = 0;
155 fh1PtJetsIn[ij] = 0;
65c43de8 156 fh1PtJetsInRej[ij] = 0;
b01d3ab7 157 fh1PtJetsInBest[ij] = 0;
a31b8a87 158 fh1PtTracksIn[ij] = 0;
cb883243 159 fh1PtTracksInLow[ij] = 0;
a31b8a87 160 fh2NJetsPt[ij] = 0;
161 fh2NTracksPt[ij] = 0;
2eded560 162 fp2MultRPPhiTrackPt[ij] = 0;
163 fp2CentRPPhiTrackPt[ij] = 0;
90a006c2 164 fhnJetPt[ij] = 0;
17f4943e 165 fhnJetPtBest[ij] = 0;
166 fhnJetPtRej[ij] = 0;
90a006c2 167 fhnJetPtQA[ij] = 0;
168 fhnTrackPt[ij] = 0;
169 fhnTrackPtQA[ij] = 0;
37eb26ea 170 for(int i = 0;i <= kMaxJets;++i){
d7117cbd 171 fh2LTrackPtJetPt[ij][i] = 0;
a31b8a87 172 fh1PtIn[ij][i] = 0;
a31b8a87 173 }
174
175 fh1DijetMinv[ij] = 0;
176 fh2DijetDeltaPhiPt[ij] = 0;
177 fh2DijetAsymPt[ij] = 0;
178 fh2DijetPt2vsPt1[ij] = 0;
179 fh2DijetDifvsSum[ij] = 0;
180
181 }
3b7ffecf 182}
183
184AliAnalysisTaskJetSpectrum2::AliAnalysisTaskJetSpectrum2(const char* name):
185 AliAnalysisTaskSE(name),
186 fJetHeaderRec(0x0),
187 fJetHeaderGen(0x0),
a31b8a87 188 fAODIn(0x0),
189 fAODOut(0x0),
190 fAODExtension(0x0),
b99018c6 191 fhnJetContainer(0x0),
3b7ffecf 192 fhnCorrelation(0x0),
90a006c2 193 fhnEvent(0x0),
c8eabe24 194 f1PtScale(0x0),
3b7ffecf 195 fBranchRec("jets"),
196 fBranchGen(""),
37eb26ea 197 fBranchBkgRec(""),
198 fBranchBkgGen(""),
a31b8a87 199 fNonStdFile(""),
742ee86c 200 fRandomizer(0x0),
565584e8 201 fUseAODJetInput(kFALSE),
202 fUseAODTrackInput(kFALSE),
203 fUseAODMCInput(kFALSE),
b5cc0c6d 204 fUseGlobalSelection(kFALSE),
3b7ffecf 205 fUseExternalWeightOnly(kFALSE),
206 fLimitGenJetEta(kFALSE),
b99018c6 207 fDoMatching(kFALSE),
a31b8a87 208 fNMatchJets(5),
7b822bc5 209 fNRPBins(3),
4f6d1cc1 210 fTRP(0),
211 fDebug(0),
65c43de8 212 fJetTriggerExcludeMask(AliAODJet::kHighTrackPtTriggered),
d7396b04 213 fJetTriggerBestMask(AliAODJet::kHighTrackPtBest),
3b7ffecf 214 fFilterMask(0),
f2dd0695 215 fEventSelectionMask(0),
3d1dba29 216 fNTrigger(0),
217 fTriggerBit(0x0),
db745d29 218 fNAcceptance(0),
85b5b73e 219 fNBinsLeadingTrackPt(10),
220 fNBinsMult(20),
3b7ffecf 221 fAnalysisType(0),
222 fTrackTypeRec(kTrackUndef),
223 fTrackTypeGen(kTrackUndef),
f4132e7d 224 fEventClass(0),
2eded560 225 fRPMethod(0),
3b7ffecf 226 fAvgTrials(1),
227 fExternalWeight(1),
a31b8a87 228 fJetRecEtaWindow(0.5),
229 fTrackRecEtaWindow(0.5),
9280dfa6 230 fMinJetPt(0),
a31b8a87 231 fMinTrackPt(0.15),
232 fDeltaPhiWindow(90./180.*TMath::Pi()),
69a0d320 233 fAcceptancePhiMin(0x0),
234 fAcceptancePhiMax(0x0),
235 fAcceptanceEtaMin(0x0),
236 fAcceptanceEtaMax(0x0),
742ee86c 237 fCentrality(100),
7b822bc5 238 fRPAngle(0),
37eb26ea 239 fMultRec(0),
240 fMultGen(0),
3d1dba29 241 fTriggerName(0x0),
3b7ffecf 242 fh1Xsec(0x0),
243 fh1Trials(0x0),
432abb2b 244 fh1AvgTrials(0x0),
3b7ffecf 245 fh1PtHard(0x0),
246 fh1PtHardNoW(0x0),
247 fh1PtHardTrials(0x0),
57ca1193 248 fh1ZVtx(0x0),
7b822bc5 249 fh1RP(0x0),
742ee86c 250 fh1Centrality(0x0),
c8eabe24 251 fh1TmpRho(0x0),
37eb26ea 252 fh2MultRec(0x0),
253 fh2MultGen(0x0),
742ee86c 254 fh2RPCentrality(0x0),
a31b8a87 255 fh2PtFGen(0x0),
e39a4f22 256 fh2deltaPt1Pt2(0x0),
9a83d4af 257 fh2RelPtFGen(0x0),
346e5a93 258 fh3RelPtFGenLeadTrkPt(0x0),
0c964c83 259 fh1EvtSelection(0),
260 fMaxVertexZ(100.),
261 fMinNcontributors(0),
a2d2ff89 262 fRejectPileup(0),
263 fHistList(0x0)
3b7ffecf 264{
265
a31b8a87 266 for(int ij = 0;ij <kJetTypes;++ij){
742ee86c 267 fFlagJetType[ij] = 1; // default = on
a31b8a87 268 fh1NJets[ij] = 0;
269 fh1SumPtTrack[ij] = 0;
270 fh1PtJetsIn[ij] = 0;
65c43de8 271 fh1PtJetsInRej[ij] = 0;
b01d3ab7 272 fh1PtJetsInBest[ij] = 0;
a31b8a87 273 fh1PtTracksIn[ij] = 0;
cb883243 274 fh1PtTracksInLow[ij] = 0;
a31b8a87 275 fh2NJetsPt[ij] = 0;
276 fh2NTracksPt[ij] = 0;
2eded560 277 fp2MultRPPhiTrackPt[ij] = 0;
278 fp2CentRPPhiTrackPt[ij] = 0;
90a006c2 279 fhnJetPt[ij] = 0;
17f4943e 280 fhnJetPtBest[ij] = 0;
281 fhnJetPtRej[ij] = 0;
90a006c2 282 fhnJetPtQA[ij] = 0;
283 fhnTrackPt[ij] = 0;
284 fhnTrackPtQA[ij] = 0;
37eb26ea 285 for(int i = 0;i <= kMaxJets;++i){
d7117cbd 286 fh2LTrackPtJetPt[ij][i] = 0;
a31b8a87 287 fh1PtIn[ij][i] = 0;
a31b8a87 288 }
289
290 fh1DijetMinv[ij] = 0;
291 fh2DijetDeltaPhiPt[ij] = 0;
292 fh2DijetAsymPt[ij] = 0;
293 fh2DijetPt2vsPt1[ij] = 0;
294 fh2DijetDifvsSum[ij] = 0;
295 }
296
37eb26ea 297 DefineOutput(1, TList::Class());
3b7ffecf 298}
299
300
301
302Bool_t AliAnalysisTaskJetSpectrum2::Notify()
303{
a31b8a87 304
305
306
3b7ffecf 307 //
308 // Implemented Notify() to read the cross sections
309 // and number of trials from pyxsec.root
310 //
a31b8a87 311
312 // Fetch the aod also from the input in,
313 // have todo it in notify
314
315
316 fAODIn = dynamic_cast<AliAODEvent*>(InputEvent());
317 // assume that the AOD is in the general output...
318 fAODOut = AODEvent();
319
320 if(fNonStdFile.Length()!=0){
321 // case that we have an AOD extension we need can fetch the jets from the extended output
322 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
e855f5c5 323 fAODExtension = (aodH?aodH->GetExtension(fNonStdFile.Data()):0);
a31b8a87 324 if(!fAODExtension){
325 if(fDebug>1)Printf("AODExtension found for %s",fNonStdFile.Data());
326 }
327 }
328
3b7ffecf 329
330 TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
331 Float_t xsection = 0;
332 Float_t ftrials = 1;
333
334 fAvgTrials = 1;
335 if(tree){
336 TFile *curfile = tree->GetCurrentFile();
337 if (!curfile) {
338 Error("Notify","No current file");
339 return kFALSE;
340 }
341 if(!fh1Xsec||!fh1Trials){
346cfe4a 342 Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__);
3b7ffecf 343 return kFALSE;
344 }
345 AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials);
346 fh1Xsec->Fill("<#sigma>",xsection);
347 // construct a poor man average trials
348 Float_t nEntries = (Float_t)tree->GetTree()->GetEntries();
a221560c 349 if(ftrials>=nEntries && nEntries>0.)fAvgTrials = ftrials/nEntries;
432abb2b 350 fh1Trials->Fill("#sum{ntrials}",ftrials);
3b7ffecf 351 }
a31b8a87 352
d7117cbd 353 if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName());
a31b8a87 354
3b7ffecf 355 return kTRUE;
356}
357
358void AliAnalysisTaskJetSpectrum2::UserCreateOutputObjects()
359{
360
742ee86c 361
3b7ffecf 362 // Connect the AOD
363
3b7ffecf 364 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::UserCreateOutputObjects() \n");
3b7ffecf 365 OpenFile(1);
37eb26ea 366 if(!fHistList)fHistList = new TList();
40440b28 367 PostData(1, fHistList); // post data in any case once
368
742ee86c 369 if(!fRandomizer)fRandomizer = new TRandom3(0);
370
5cb0f438 371 fHistList->SetOwner(kTRUE);
37eb26ea 372 Bool_t oldStatus = TH1::AddDirectoryStatus();
3b7ffecf 373 TH1::AddDirectory(kFALSE);
374
90a006c2 375
376
377 // event npsparse cent, mult
e5ee1c19 378 const Int_t nBinsSparse0 = 4;
379 const Int_t nBins0[nBinsSparse0] = { 100, 500,fNTrigger,125};
380 const Double_t xmin0[nBinsSparse0] = { 0, 0, -0.5,-2};
381 const Double_t xmax0[nBinsSparse0] = { 100,5000,fNTrigger-0.5,248};
90a006c2 382
383
e5ee1c19 384 fhnEvent = new THnSparseF("fhnEvent",";cent;mult:trigger;#rho",nBinsSparse0,
90a006c2 385 nBins0,xmin0,xmax0);
386 fHistList->Add(fhnEvent);
387
b99018c6 388 if(fDoMatching){
389 MakeJetContainer();
390 fHistList->Add(fhnCorrelation);
391 fHistList->Add(fhnJetContainer);
392 }
3b7ffecf 393 //
394 // Histogram
395
90a006c2 396
0c964c83 397 fh1EvtSelection = new TH1F("fh1EvtSelection", "Event Selection", 7, -0.5, 6.5);
398 fh1EvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED");
399 fh1EvtSelection->GetXaxis()->SetBinLabel(2,"event selection: rejected");
400 fh1EvtSelection->GetXaxis()->SetBinLabel(3,"event class: rejected");
401 fh1EvtSelection->GetXaxis()->SetBinLabel(4,"vertex Ncontr: rejected");
402 fh1EvtSelection->GetXaxis()->SetBinLabel(5,"vertex z: rejected");
403 fh1EvtSelection->GetXaxis()->SetBinLabel(6,"vertex type: rejected");
404 fh1EvtSelection->GetXaxis()->SetBinLabel(7,"pileup: rejected");
90a006c2 405
0c964c83 406 fHistList->Add(fh1EvtSelection);
407
66fd4b2b 408 const Int_t nBinPt = 150;
3b7ffecf 409 Double_t binLimitsPt[nBinPt+1];
410 for(Int_t iPt = 0;iPt <= nBinPt;iPt++){
411 if(iPt == 0){
ed5c9f3f 412 binLimitsPt[iPt] = -50.0;
3b7ffecf 413 }
414 else {// 1.0
66fd4b2b 415 binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.;
3b7ffecf 416 }
417 }
edfbe476 418 const Int_t nBinPhi = 90;
3b7ffecf 419 Double_t binLimitsPhi[nBinPhi+1];
420 for(Int_t iPhi = 0;iPhi<=nBinPhi;iPhi++){
421 if(iPhi==0){
37eb26ea 422 binLimitsPhi[iPhi] = 0;
3b7ffecf 423 }
424 else{
425 binLimitsPhi[iPhi] = binLimitsPhi[iPhi-1] + 1/(Float_t)nBinPhi * TMath::Pi()*2;
426 }
427 }
428
edfbe476 429 const Int_t nBinEta = 40;
430 Double_t binLimitsEta[nBinEta+1];
431 for(Int_t iEta = 0;iEta<=nBinEta;iEta++){
432 if(iEta==0){
433 binLimitsEta[iEta] = -2.0;
434 }
435 else{
cf049e94 436 binLimitsEta[iEta] = binLimitsEta[iEta-1] + 0.1;
edfbe476 437 }
438 }
439
a31b8a87 440
3b7ffecf 441 fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1);
442 fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
90a006c2 443 fHistList->Add(fh1Xsec);
3b7ffecf 444 fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1);
445 fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
a31b8a87 446 fHistList->Add(fh1Trials);
432abb2b 447 fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1);
448 fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}");
449 fHistList->Add(fh1AvgTrials);
3b7ffecf 450 fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 451 fHistList->Add(fh1PtHard);
3b7ffecf 452 fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 453 fHistList->Add(fh1PtHardNoW);
3b7ffecf 454 fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt);
a31b8a87 455 fHistList->Add(fh1PtHardTrials);
a31b8a87 456
57ca1193 457 fh1ZVtx = new TH1F("fh1ZVtx","z vtx;z_{vtx} (cm)",400,-20,20);
a31b8a87 458 fHistList->Add(fh1ZVtx);
7b822bc5 459
2eded560 460
9939142f 461 fh1RP = new TH1F("fh1RP","RP;#Psi",440, -1.*TMath::Pi(), 2.*TMath::Pi());
7b822bc5 462 fHistList->Add(fh1RP);
463
eb1739b7 464 fh1Centrality = new TH1F("fh1Centrality","cent;cent (%)",103,-1,102);
742ee86c 465 fHistList->Add(fh1Centrality);
466
3170a3f8 467 fh2MultRec = new TH2F("fh2MultRec","multiplicity rec;# tracks;# jetrefs",500,0,5000,500,0.,5000);
37eb26ea 468 fHistList->Add(fh2MultRec);
3170a3f8 469 fh2MultGen = new TH2F("fh2MultGen","multiplicity gen;# tracks;# jetrefs",400,0,5000,500,0.,5000);
37eb26ea 470 fHistList->Add(fh2MultGen);
a31b8a87 471
a31b8a87 472
742ee86c 473 fh2RPCentrality = new TH2F("fh2RPCentrality" ,"Reaction Plane Angle" , 20, 0.,100., 180, 0, TMath::Pi());
474 fHistList->Add(fh2RPCentrality);
a31b8a87 475
742ee86c 476 fh2PtFGen = new TH2F("fh2PtFGen",Form("%s vs. %s;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt);
477 fHistList->Add(fh2PtFGen);
a31b8a87 478
e39a4f22 479 // fh2deltaPt1Pt2(0x0),
480 fh2deltaPt1Pt2 = new TH3F("fh2deltaPt1Pt2",Form("%s vs. %s;delta;p_{T,gen};p_{T,rec}",fBranchRec.Data(),fBranchGen.Data()),nBinPt,binLimitsPt,nBinPt,binLimitsPt,nBinPt,binLimitsPt);
481 fHistList->Add(fh2deltaPt1Pt2);
482
742ee86c 483 fh2RelPtFGen = new TH2F("fh2RelPtFGen",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen}",nBinPt,binLimitsPt,241,-2.41,2.41);
484 fHistList->Add(fh2RelPtFGen);
346e5a93 485
486 const Int_t nBinsLeadingTrackPt = 10;
487 const Double_t binArrayLeadingTrackPt[nBinsLeadingTrackPt+1] = {0.,1.,2.,3.,4.,5.,6.,8.,10.,12.,200.}; //store pT of leading track in jet
488
489 const Int_t nBinDeltaPt = 241;
490 Double_t binLimitsDeltaPt[nBinDeltaPt+1];
491 for(Int_t iDeltaPt = 0;iDeltaPt <= nBinDeltaPt;iDeltaPt++){
492 if(iDeltaPt == 0){
493 binLimitsDeltaPt[iDeltaPt] = -2.41;
494 }
495 else {
496 binLimitsDeltaPt[iDeltaPt] = binLimitsDeltaPt[iDeltaPt-1] + 0.02;
497 }
498 }
499
500 fh3RelPtFGenLeadTrkPt = new TH3F("fh3RelPtFGenLeadTrkPt",";p_{T,gen};p_{T,rec}-p_{T,gen}/p_{T,Gen};p_{T}^{leading track}",nBinPt,binLimitsPt,nBinDeltaPt,binLimitsDeltaPt,nBinsLeadingTrackPt,binArrayLeadingTrackPt);
501 fHistList->Add(fh3RelPtFGenLeadTrkPt);
502
742ee86c 503 for(int ij = 0;ij <kJetTypes;++ij){
504 TString cAdd = "";
505 TString cJetBranch = "";
506 if(ij==kJetRec){
507 cAdd = "Rec";
508 cJetBranch = fBranchRec.Data();
509 }
510 else if (ij==kJetGen){
511 cAdd = "Gen";
512 cJetBranch = fBranchGen.Data();
513 }
514 else if (ij==kJetRecFull){
515 cAdd = "RecFull";
516 cJetBranch = fBranchRec.Data();
517 }
518 else if (ij==kJetGenFull){
519 cAdd = "GenFull";
520 cJetBranch = fBranchGen.Data();
521 }
a31b8a87 522
742ee86c 523 if(cJetBranch.Length()==0)fFlagJetType[ij] = 0;
524 if(!fFlagJetType[ij])continue;
7b822bc5 525
742ee86c 526 fh1NJets[ij] =new TH1F(Form("fh1N%sJets",cAdd.Data()),Form("N %s jets",cAdd.Data()),50,-0.5,49.5);
527 fHistList->Add(fh1NJets[ij]);
528
529 fh1PtJetsIn[ij] = new TH1F(Form("fh1PtJets%sIn",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
530 fHistList->Add(fh1PtJetsIn[ij]);
65c43de8 531
532 fh1PtJetsInRej[ij] = new TH1F(Form("fh1PtJets%sInRej",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
533 fHistList->Add(fh1PtJetsInRej[ij]);
d7396b04 534
535 fh1PtJetsInBest[ij] = new TH1F(Form("fh1PtJets%sInBest",cAdd.Data()),Form("%s jets p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
536 fHistList->Add(fh1PtJetsInBest[ij]);
742ee86c 537
538 fh1PtTracksIn[ij] = new TH1F(Form("fh1PtTracks%sIn",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),nBinPt,binLimitsPt);
539 fHistList->Add(fh1PtTracksIn[ij]);
540
2eded560 541 fh1PtTracksInLow[ij] = new TH1F(Form("fh1PtTracks%sInLow",cAdd.Data()),Form("%s track p_T;p_{T} (GeV/c)",cAdd.Data()),100,0.,5.);
742ee86c 542 fHistList->Add(fh1PtTracksInLow[ij]);
543
3170a3f8 544 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 545 fHistList->Add(fh1SumPtTrack[ij]);
546
547 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);
548 fHistList->Add(fh2NJetsPt[ij]);
549
3170a3f8 550 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 551 fHistList->Add(fh2NTracksPt[ij]);
90a006c2 552
2eded560 553
554 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");
555 fHistList->Add(fp2MultRPPhiTrackPt[ij]);
556 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");
557 fHistList->Add(fp2CentRPPhiTrackPt[ij]);
558
556ccab2 559 // Bins: Jet number: pTJet, cent, mult, RP, Area, trigger, leading track pT total bins = 4.5M
db745d29 560 const Int_t nBinsSparse1 = 9;
152a0744 561 const Int_t nBinsArea = 10;
4f6d1cc1 562 Int_t nbinr5=120;
563 Int_t nbinr5min=-50;
564
565 if(cJetBranch.Contains("KT05")){
566 nbinr5=132;
567 nbinr5min=-80;}
568
569 Int_t nBins1[nBinsSparse1] = { kMaxJets+1,nbinr5, 10, fNBinsMult, fNRPBins, nBinsArea,fNTrigger,fNBinsLeadingTrackPt,fNAcceptance+1};
d01931c3 570 if(cJetBranch.Contains("RandomCone")){
571 nBins1[1] = 600;
572 nBins1[5] = 1;
573 }
2942f542 574 const Double_t xmin1[nBinsSparse1] = { -0.5, static_cast<Double_t>(nbinr5min), 0, 0, -0.5, 0., -0.5, 0., -0.5,};
db745d29 575 const Double_t xmax1[nBinsSparse1] = {kMaxJets+0.5,250,100,4000,fNRPBins-0.5,1.0,fNTrigger-0.5,200.,fNAcceptance+0.5};
4f6d1cc1 576
adf3920d 577 const Double_t binArrayArea[nBinsArea+1] = {xmin1[5],0.07,0.1,0.15,0.2,0.25,0.3,0.4,0.5,0.6,xmax1[5]};
346e5a93 578
7bee4353 579
db745d29 580 fhnJetPt[ij] = new THnSparseF(Form("fhnJetPt%s",cAdd.Data()),";jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leading track p_{T};acceptance bin",nBinsSparse1,nBins1,xmin1,xmax1);
152a0744 581 fhnJetPt[ij]->SetBinEdges(5,binArrayArea);
85b5b73e 582 if(fNBinsLeadingTrackPt>1) fhnJetPt[ij]->SetBinEdges(7,binArrayLeadingTrackPt);
2eded560 583 fHistList->Add(fhnJetPt[ij]);
17f4943e 584
585
586 // Bins: pTJet, cent, trigger,
587 const Int_t nBinsSparse1b = 3;
588 Int_t nBins1b[nBinsSparse1b] = {120, 10,fNTrigger};
589 const Double_t xmin1b[nBinsSparse1b] = {-50, 0,-0.5};
590 const Double_t xmax1b[nBinsSparse1b] = {250,100,fNTrigger-0.5};
591
592 fhnJetPtBest[ij] = new THnSparseF(Form("fhnJetPtBest%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b);
593 fHistList->Add(fhnJetPtBest[ij]);
594
595 fhnJetPtRej[ij] = new THnSparseF(Form("fhnJetPtRej%s",cAdd.Data()),";p_{T,jet};cent;trigger",nBinsSparse1b,nBins1b,xmin1b,xmax1b);
596 fHistList->Add(fhnJetPtRej[ij]);
2eded560 597
adf3920d 598 // Bins: Jet number: pTJet, cent, eta, phi, Area, trigger, acceptance, signed pT leading
5662640e 599 const Int_t nBinsSparse2 = 9;
fedb5a47 600 Int_t nBins2[nBinsSparse2] = { kMaxJets+1, 60, 8, 6, 90, 100,fNTrigger,fNAcceptance+1,20};
d01931c3 601 if(cJetBranch.Contains("RandomCone")){
602 nBins2[5] = 1;
603 }
adf3920d 604 const Double_t xmin2[nBinsSparse2] = { -0.5, -50, 0,-0.9, 0, 0., -0.5, -0.5, -100};
605 const Double_t xmax2[nBinsSparse2] = {kMaxJets+0.5, 250, 80, 0.9, 2.*TMath::Pi(),1.0,fNTrigger-0.5,fNAcceptance+0.5, 100};
5662640e 606 fhnJetPtQA[ij] = new THnSparseF(Form("fhnJetPtQA%s",cAdd.Data()),";jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin;signed pt leading",nBinsSparse2,nBins2,xmin2,xmax2);
e39a4f22 607 // fhnJetPt[ij]->SetBinEdges(5,binArrayArea);
d01931c3 608 fHistList->Add(fhnJetPtQA[ij]);
609
5a7c81dc 610 // Bins:track number pTtrack, cent, mult, RP. charge total bins = 224 k
611 const Int_t nBinsSparse3 = 7;
85b5b73e 612 const Int_t nRPBinsSparse3 = 3;
5a7c81dc 613 const Int_t nBins3[nBinsSparse3] = { 2, 100, 10, 1, nRPBinsSparse3,fNTrigger,2};
614 const Double_t xmin3[nBinsSparse3] = { -0.5, 0, 0, 0, -0.5,-0.5,-1.5};
615 const Double_t xmax3[nBinsSparse3] = { 1.5, 200, 100, 4000,nRPBinsSparse3-0.5,fNTrigger-0.5,1.5};
d01931c3 616
adf3920d 617 // change the binning of the pT axis:
d01931c3 618 Double_t *xPt3 = new Double_t[nBins3[1]+1];
619 xPt3[0] = 0.;
620 for(int i = 1; i<=nBins3[1];i++){
621 if(xPt3[i-1]<2)xPt3[i] = xPt3[i-1] + 0.05; // 1 - 40
622 else if(xPt3[i-1]<4)xPt3[i] = xPt3[i-1] + 0.2; // 41 - 50
623 else if(xPt3[i-1]<10)xPt3[i] = xPt3[i-1] + 0.5; // 50 - 62
17f4943e 624 else if(xPt3[i-1]<15)xPt3[i] = xPt3[i-1] + 1.; // 62 - 67
625 else if(xPt3[i-1]<20)xPt3[i] = xPt3[i-1] + 2.; // 67 - 72
626 else if(xPt3[i-1]<60)xPt3[i] = xPt3[i-1] + 5; // 72 - 76
627 else xPt3[i] = xPt3[i-1] + 10; // 76 - 100 = 140
d01931c3 628 }
629
adf3920d 630 fhnTrackPt[ij] = new THnSparseF(Form("fhnTrackPt%s",cAdd.Data()),";track number;p_{T};cent;#tracks;RP;trigger",nBinsSparse3,nBins3,xmin3,xmax3);
631 fhnTrackPt[ij]->SetBinEdges(1,xPt3);
632 fHistList->Add(fhnTrackPt[ij]);
633 delete [] xPt3;
634
635 // Track QA bins track nr, pTrack, cent, eta, phi bins 5.4 M
5a7c81dc 636 const Int_t nBinsSparse4 = 6;
637 const Int_t nBins4[nBinsSparse4] = { 2, 50, 10, 20, 180,2};
638 const Double_t xmin4[nBinsSparse4] = { -0.5, 0, 0, -1.0, 0.,-1.5};
639 const Double_t xmax4[nBinsSparse4] = { 1.5,150, 100, 1.0,2.*TMath::Pi(),1.5};
adf3920d 640
641 // change the binning ot the pT axis:
642 Double_t *xPt4 = new Double_t[nBins4[1]+1];
643 xPt4[0] = 0.;
644 for(int i = 1; i<=nBins4[1];i++){
645 if(xPt4[i-1]<2)xPt4[i] = xPt4[i-1] + 0.1;
646 else if(xPt4[i-1]<10)xPt4[i] = xPt4[i-1] + 0.5;
647 else if(xPt4[i-1]<20)xPt4[i] = xPt4[i-1] + 1.;
648 else if(xPt4[i-1]<30)xPt4[i] = xPt4[i-1] + 2.5;
649 else xPt4[i] = xPt4[i-1] + 5.;
650 }
5a7c81dc 651 fhnTrackPtQA[ij] = new THnSparseF(Form("fhnTrackPtQA%s",cAdd.Data()),";track number;p_{T};cent;#eta;#phi;sign",nBinsSparse4,nBins4,xmin4,xmax4);
adf3920d 652 fhnTrackPtQA[ij]->SetBinEdges(1,xPt4);
653 fHistList->Add(fhnTrackPtQA[ij]);
654 delete [] xPt4;
655
656 for(int i = 0;i <= kMaxJets;++i){
657 fh1PtIn[ij][i] = new TH1F(Form("fh1Pt%sIn_j%d",cAdd.Data(),i),Form("%s p_T input ;p_{T}",cAdd.Data()),nBinPt,binLimitsPt);
658 fHistList->Add(fh1PtIn[ij][i]);
a31b8a87 659
660
adf3920d 661 if(!fh1TmpRho)fh1TmpRho = new TH1F("fh1TmpRho","tmp histo for jet shape",40,0.,2);
662 fh2LTrackPtJetPt[ij][i] = new TH2F(Form("fh2LTrackPtJetPt%s_j%d",cAdd.Data(),i),
663 Form("pt of leadin track within a jet vs jet %s;p_{T,lead in jet};p_{T.jet};",
664 cAdd.Data()),
665 200,0.,200.,nBinPt,binLimitsPt);
666 fHistList->Add(fh2LTrackPtJetPt[ij][i]);
667 }
668
a31b8a87 669
adf3920d 670 fh1DijetMinv[ij] = new TH1F(Form("fh1Dijet%sMinv",cAdd.Data()),"Dijet invariant mass;m_{JJ}",nBinPt,binLimitsPt);
671 fHistList->Add(fh1DijetMinv[ij]);
a31b8a87 672
adf3920d 673 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);
674 fHistList->Add(fh2DijetDeltaPhiPt[ij]);
a31b8a87 675
adf3920d 676 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);
677 fHistList->Add(fh2DijetAsymPt[ij]);
678
679 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.);
680 fHistList->Add(fh2DijetPt2vsPt1[ij]);
681 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.);
682 fHistList->Add( fh2DijetDifvsSum[ij]);
683 }
3b7ffecf 684 // =========== Switch on Sumw2 for all histos ===========
685 for (Int_t i=0; i<fHistList->GetEntries(); ++i) {
686 TH1 *h1 = dynamic_cast<TH1*>(fHistList->At(i));
687 if (h1){
688 h1->Sumw2();
689 continue;
690 }
691 THnSparse *hn = dynamic_cast<THnSparse*>(fHistList->At(i));
692 if(hn)hn->Sumw2();
693 }
694 TH1::AddDirectory(oldStatus);
695}
696
697void AliAnalysisTaskJetSpectrum2::Init()
698{
699 //
700 // Initialization
701 //
702
3b7ffecf 703 if (fDebug > 1) printf("AnalysisTaskJetSpectrum2::Init() \n");
704
705}
706
a31b8a87 707void AliAnalysisTaskJetSpectrum2::UserExec(Option_t */*option*/){
708
b5cc0c6d 709
f2dd0695 710 Bool_t selected = kTRUE;
711
712 if(fUseGlobalSelection&&fEventSelectionMask==0){
713 selected = AliAnalysisHelperJetTasks::Selected();
714 }
715 else if(fUseGlobalSelection&&fEventSelectionMask>0){
45a11aef 716 selected = AliAnalysisHelperJetTasks::TestSelectInfo(fEventSelectionMask);
f2dd0695 717 }
718
0c964c83 719 if(!selected){
720 // no selection by the service task, we continue
721 if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
722 fh1EvtSelection->Fill(1.);
723 PostData(1, fHistList);
724 return;
725 }
726
f4132e7d 727 if(fEventClass>0){
728 selected = selected&&(AliAnalysisHelperJetTasks::EventClass()==fEventClass);
729 }
730
f2dd0695 731 if(!selected){
b5cc0c6d 732 // no selection by the service task, we continue
f4132e7d 733 if (fDebug > 1)Printf("Not selected %s:%d SelectInfo %d Class %d",(char*)__FILE__,__LINE__, AliAnalysisHelperJetTasks::Selected(),AliAnalysisHelperJetTasks::EventClass());
0c964c83 734 fh1EvtSelection->Fill(2.);
b5cc0c6d 735 PostData(1, fHistList);
736 return;
737 }
738
f2dd0695 739
a31b8a87 740 static AliAODEvent* aod = 0;
741
742 // take all other information from the aod we take the tracks from
743 if(!aod){
adf3920d 744 if(fUseAODTrackInput)aod = fAODIn;
745 else aod = fAODOut;
a31b8a87 746 }
747
748
3b7ffecf 749 //
750 // Execute analysis for current event
751 //
a31b8a87 752 if (fDebug > 1)printf("Analysing event # %5d\n", (Int_t) fEntry);
753 AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
754 if(!aodH){
755 Printf("%s:%d no output aodHandler found Jet",(char*)__FILE__,__LINE__);
756 return;
757 }
758
759 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
760 TClonesArray *aodRecJets = 0;
761
762 if(fAODOut&&!aodRecJets){
763 aodRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchRec.Data()));
764 }
765 if(fAODExtension&&!aodRecJets){
766 aodRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchRec.Data()));
767 }
768 if(fAODIn&&!aodRecJets){
769 aodRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchRec.Data()));
770 }
771
772
773
774 if(!aodRecJets){
d95fc15a 775 if(fDebug){
776
adf3920d 777 Printf("%s:%d no reconstructed Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchRec.Data());
778 if(fAODIn){
779 Printf("Input AOD >>>>");
780 fAODIn->Print();
781 }
782 if(fAODExtension){
783 Printf("AOD Extension >>>>");
784 fAODExtension->Print();
785 }
786 if(fAODOut){
787 Printf("Output AOD >>>>");
788 fAODOut->Print();
789 }
d95fc15a 790 }
adf3920d 791 return;
a31b8a87 792 }
793
794 TClonesArray *aodGenJets = 0;
795 if(fBranchGen.Length()>0){
796 if(fAODOut&&!aodGenJets){
797 aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchGen.Data()));
798 }
799 if(fAODExtension&&!aodGenJets){
800 aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchGen.Data()));
801 }
802 if(fAODIn&&!aodGenJets){
803 aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchGen.Data()));
804 }
805
806 if(!aodGenJets){
807 Printf("%s:%d no generated Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchGen.Data());
3b7ffecf 808 return;
809 }
3b7ffecf 810 }
a31b8a87 811
37eb26ea 812 TClonesArray *aodBackRecJets = 0;
813 if(fBranchBkgRec.Length()>0){
814 if(fAODOut&&!aodBackRecJets){
815 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgRec.Data()));
816 }
817 if(fAODExtension&&!aodBackRecJets){
818 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgRec.Data()));
819 }
820 if(fAODIn&&!aodBackRecJets){
821 aodBackRecJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgRec.Data()));
822 }
823
824 if(!aodBackRecJets){
825 Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgRec.Data());
826 return;
827 }
828 }
829
830
831 TClonesArray *aodBackGenJets = 0;
832
833 if(fBranchBkgGen.Length()>0){
834 if(fAODOut&&!aodBackGenJets){
835 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fBranchBkgGen.Data()));
836 }
837 if(fAODExtension&&!aodBackGenJets){
838 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fBranchBkgGen.Data()));
839 }
840 if(fAODIn&&!aodBackGenJets){
841 aodBackGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fBranchBkgGen.Data()));
842 }
843
844 if(!aodBackGenJets){
845 Printf("%s:%d no background rec Jet array with name %s in AOD",(char*)__FILE__,__LINE__,fBranchBkgGen.Data());
846 return;
847 }
848 }
849
0c964c83 850 AliAODVertex* primVtx = aod->GetPrimaryVertex();
851 Int_t nTracksPrim = primVtx->GetNContributors();
852
853
854 if (fDebug > 1) Printf("%s:%d primary vertex selection: %d", (char*)__FILE__,__LINE__,nTracksPrim);
855 if(nTracksPrim < fMinNcontributors){
856 if (fDebug > 1) Printf("%s:%d primary vertex selection: event REJECTED...",(char*)__FILE__,__LINE__);
857 fh1EvtSelection->Fill(3.);
858 PostData(1, fHistList);
859 return;
860 }
861
862 if(TMath::Abs(primVtx->GetZ())>fMaxVertexZ){
863 if (fDebug > 1) Printf("%s:%d primary vertex z = %f: event REJECTED...",(char*)__FILE__,__LINE__,primVtx->GetZ());
864 fh1EvtSelection->Fill(4.);
865 PostData(1, fHistList);
866 return;
867 }
868
869 TString primVtxName(primVtx->GetName());
870
871 if(primVtxName.CompareTo("TPCVertex",TString::kIgnoreCase) == 1){
872 if (fDebug > 1) Printf("%s:%d primary vertex selection: TPC vertex, event REJECTED...",(char*)__FILE__,__LINE__);
873 fh1EvtSelection->Fill(5.);
874 PostData(1, fHistList);
875 return;
876 }
877
878 if(fRejectPileup && AliAnalysisHelperJetTasks::IsPileUp()){
879 if (fDebug > 1) Printf("%s:%d SPD pileup: event REJECTED...",(char*)__FILE__,__LINE__);
880 fh1EvtSelection->Fill(6.);
881 PostData(1, fHistList);
882 return;
883 }
884
885 if (fDebug > 1) Printf("%s:%d event ACCEPTED ...",(char*)__FILE__,__LINE__);
886 fh1EvtSelection->Fill(0.);
887
a31b8a87 888
889 // new Scheme
890 // first fill all the pure histograms, i.e. full jets
891 // in case of the correaltion limit to the n-leading jets
892
893 // reconstructed
894
895
896 // generated
897
898
899 // second fill the correlation histograms, here we limit to the n-th leading jet in case of the reconstructed
900
901
742ee86c 902
a31b8a87 903 TList genJetsList; // full acceptance
904 TList genJetsListCut; // acceptance cut
905 TList recJetsList; // full acceptance
906 TList recJetsListCut; // acceptance cut
907
908
909 GetListOfJets(&recJetsList,aodRecJets,0);
910 GetListOfJets(&recJetsListCut,aodRecJets,1);
911
912 if(fBranchGen.Length()>0){
913 GetListOfJets(&genJetsList,aodGenJets,0);
914 GetListOfJets(&genJetsListCut,aodGenJets,1);
915 }
916
917 Double_t eventW = 1;
918 Double_t ptHard = 0;
919 Double_t nTrials = 1; // Trials for MC trigger
432abb2b 920 fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials);
a31b8a87 921
742ee86c 922 // Getting some global properties
923 fCentrality = GetCentrality();
aa3ba8d2 924 if(fCentrality<=0)fCentrality = 0;
742ee86c 925 fh1Centrality->Fill(fCentrality);
926
927
a31b8a87 928 if((fAnalysisType&kAnaMCESD)==kAnaMCESD){
929 // this is the part we only use when we have MC information
930 AliMCEvent* mcEvent = MCEvent();
931 // AliStack *pStack = 0;
932 if(!mcEvent){
933 Printf("%s:%d no mcEvent",(char*)__FILE__,__LINE__);
3b7ffecf 934 return;
935 }
a31b8a87 936 AliGenPythiaEventHeader* pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent);
937 if(pythiaGenHeader){
938 nTrials = pythiaGenHeader->Trials();
939 ptHard = pythiaGenHeader->GetPtHard();
940 int iProcessType = pythiaGenHeader->ProcessType();
941 // 11 f+f -> f+f
942 // 12 f+barf -> f+barf
943 // 13 f+barf -> g+g
944 // 28 f+g -> f+g
945 // 53 g+g -> f+barf
946 // 68 g+g -> g+g
947 if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType);
948 if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent);
949 }
950 }// (fAnalysisType&kMCESD)==kMCESD)
951
952
953 // we simply fetch the tracks/mc particles as a list of AliVParticles
954
955 TList recParticles;
956 TList genParticles;
957
958 Int_t nT = GetListOfTracks(&recParticles,fTrackTypeRec);
4f6d1cc1 959
a31b8a87 960 if(fDebug>2)Printf("%s:%d Selected Rec tracks: %d %d",(char*)__FILE__,__LINE__,nT,recParticles.GetEntries());
961 nT = GetListOfTracks(&genParticles,fTrackTypeGen);
962 if(fDebug>2)Printf("%s:%d Selected Gen tracks: %d %d",(char*)__FILE__,__LINE__,nT,genParticles.GetEntries());
963
d3a3f33d 964 // CalculateReactionPlaneAngle(&recParticles);
2eded560 965 fRPAngle = 0;
966
0a918d8d 967 if(fRPMethod==0)fRPAngle = ((AliVAODHeader*)aod->GetHeader())->GetEventplane();
2eded560 968 else if(fRPMethod==1||fRPMethod==2){
0a918d8d 969 AliAODHeader * aodheader = dynamic_cast<AliAODHeader*>(aod->GetHeader());
970 if(!aodheader) AliFatal("Not a standard AOD");
971 fRPAngle = aodheader->GetQTheta(fRPMethod);
2eded560 972 }
973 fh1RP->Fill(fRPAngle);
974 fh2RPCentrality->Fill(fCentrality,fRPAngle);
a31b8a87 975 // Event control and counting ...
976 // MC
977 fh1PtHard->Fill(ptHard,eventW);
978 fh1PtHardNoW->Fill(ptHard,1);
979 fh1PtHardTrials->Fill(ptHard,nTrials);
980
981 // Real
784b1747 982 if(aod->GetPrimaryVertex()){// No vtx for pure MC
983 fh1ZVtx->Fill(aod->GetPrimaryVertex()->GetZ());
984 }
a31b8a87 985
a2522483 986
37eb26ea 987 Int_t recMult1 = recParticles.GetEntries();
988 Int_t genMult1 = genParticles.GetEntries();
989
742ee86c 990 Int_t recMult2 = MultFromJetRefs(aodBackRecJets);
991 Int_t genMult2 = MultFromJetRefs(aodBackGenJets);
37eb26ea 992
993 fh2MultRec->Fill(recMult1,recMult2);
994 fh2MultGen->Fill(genMult1,genMult2);
995 fMultRec = recMult1;
996 if(fMultRec<=0)fMultRec = recMult2;
997 fMultGen = genMult1;
998 if(fMultGen<=0)fMultGen = genMult2;
a2522483 999
e5ee1c19 1000 Double_t var0[4] = {0,};
90a006c2 1001 var0[0] = fCentrality;
1002 var0[1] = fMultRec;
3d1dba29 1003 for(int it=0;it<fNTrigger;it++){
1004 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1005 var0[2] = it;
e5ee1c19 1006 var0[3] = GetRho(recJetsList);
3d1dba29 1007 fhnEvent->Fill(var0);
1008 }
1009 }
a31b8a87 1010 // the loops for rec and gen should be indentical... pass it to a separate
1011 // function ...
1012 // Jet Loop
1013 // Track Loop
1014 // Jet Jet Loop
1015 // Jet Track loop
1016
1017 FillJetHistos(recJetsListCut,recParticles,kJetRec);
a2522483 1018 FillJetHistos(recJetsList,recParticles,kJetRecFull);
a31b8a87 1019 FillTrackHistos(recParticles,kJetRec);
1020
1021 FillJetHistos(genJetsListCut,genParticles,kJetGen);
a2522483 1022 FillJetHistos(genJetsList,genParticles,kJetGenFull);
a31b8a87 1023 FillTrackHistos(genParticles,kJetGen);
1024
1025 // Here follows the jet matching part
1026 // delegate to a separated method?
1027
b99018c6 1028 if(fDoMatching){
1029 FillMatchHistos(recJetsList,genJetsList);
1030 }
1031
a31b8a87 1032 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1033 PostData(1, fHistList);
90a006c2 1034}
a31b8a87 1035
1036void AliAnalysisTaskJetSpectrum2::FillJetHistos(TList &jetsList,TList &particlesList,Int_t iType){
1037
1038 if(iType>=kJetTypes){
1039 return;
3b7ffecf 1040 }
742ee86c 1041 if(!fFlagJetType[iType])return;
a31b8a87 1042
37eb26ea 1043 Int_t refMult = fMultRec;
1044 if(iType==kJetGen||iType==kJetGenFull){
1045 refMult = fMultGen;
1046 }
1047
a31b8a87 1048 Int_t nJets = jetsList.GetEntries();
1049 fh1NJets[iType]->Fill(nJets);
1050
1051 if(nJets<=0)return;
3b7ffecf 1052
a31b8a87 1053 Float_t pT1 = 0;
1054 Float_t phi0 = 0;
1055 Float_t phi1 = 0;
1056 Int_t ij0 = -1;
1057 Int_t ij1 = -1;
1058
db745d29 1059 Double_t var1[9] = {0,}; // jet number;p_{T,jet};cent;# tracks;RP;area;trigger;leadingTrackPt
17f4943e 1060 Double_t var1b[3] = {0,}; // p_{T,jet};cent;trigger;
7bee4353 1061
90a006c2 1062 var1[2] = fCentrality;
17f4943e 1063 var1b[1] = fCentrality;
90a006c2 1064 var1[3] = refMult;
1065
17f4943e 1066
1067
5662640e 1068 Double_t var2[9] = {0,}; // jet number;p_{T,jet};cent;#eta;#phi;area;trigger
90a006c2 1069 var2[2] = fCentrality;
a31b8a87 1070
1071 for(int ij = 0;ij < nJets;ij++){
1072 AliAODJet *jet = (AliAODJet*)jetsList.At(ij);
1073 Float_t ptJet = jet->Pt();
17f4943e 1074
1075
4a828dd2 1076 if(ptJet<0.150)ptJet = jet->GetPtSubtracted(0);
17f4943e 1077
1078 var1b[0] = ptJet;
d7396b04 1079 if(jet->Trigger()&fJetTriggerBestMask){
1080 fh1PtJetsInBest[iType]->Fill(ptJet);
17f4943e 1081 for(int it = 0;it <fNTrigger;it++){
1082 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1083 var1b[2] = it;
1084 fhnJetPtBest[iType]->Fill(var1b);
1085 }
1086 }
d7396b04 1087 }
65c43de8 1088 if(jet->Trigger()&fJetTriggerExcludeMask){
1089 fh1PtJetsInRej[iType]->Fill(ptJet);
17f4943e 1090 for(int it = 0;it <fNTrigger;it++){
1091 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1092 var1b[2] = it;
1093 fhnJetPtRej[iType]->Fill(var1b);
1094 }
1095 }
65c43de8 1096 continue;
1097 }
a31b8a87 1098 fh1PtJetsIn[iType]->Fill(ptJet);
a31b8a87 1099
b5d90baf 1100 // find the dijets assume sorting and acceptance cut...
a31b8a87 1101 if(ij==0){
1102 ij0 = ij;
a31b8a87 1103 phi0 = jet->Phi();
1104 if(phi0<0)phi0+=TMath::Pi()*2.;
1105 }
1106 else if(ptJet>pT1){
1107 // take only the backward hemisphere??
1108 phi1 = jet->Phi();
b5d90baf 1109 if(phi1<0)phi1+=TMath::Pi()*2.;
a31b8a87 1110 Float_t dPhi = phi1 - phi0;
1111 if(dPhi>TMath::Pi())dPhi = dPhi - 2.*TMath::Pi();
1112 if(dPhi<(-1.*TMath::Pi()))dPhi = dPhi + 2.*TMath::Pi();
1113 if(TMath::Abs(TMath::Pi()-dPhi)<fDeltaPhiWindow){
1114 ij1 = ij;
1115 pT1 = ptJet;
1116 }
1117 }
1118 // fill jet histos for kmax jets
adf3920d 1119 Float_t phiJet = jet->Phi();
1120 Float_t etaJet = jet->Eta();
1121 if(phiJet<0)phiJet+=TMath::Pi()*2.;
1122 fh1TmpRho->Reset();
1123 if(ij<kMaxJets)fh1PtIn[iType][ij]->Fill(ptJet);
1124
1125 fh1PtIn[iType][kMaxJets]->Fill(ptJet);
1126 // fill leading jets...
1127 AliVParticle *leadTrack = LeadingTrackFromJetRefs(jet);
4f6d1cc1 1128 //AliVParticle *leadTrack = LeadingTrackInCone(jet,&particlesList);
adf3920d 1129 Int_t phiBin = GetPhiBin(phiJet-fRPAngle);
1130 Double_t ptLead = jet->GetPtLeading(); //pT of leading jet
4f6d1cc1 1131 AliVParticle *tt = (AliVParticle*)particlesList.At(0);
1132 Float_t ttphi=tt->Phi();
1133 if(ttphi<0)ttphi+=TMath::Pi()*2.;
1134 Float_t ttpt=tt->Pt();
1135 Int_t phiBintt = GetPhiBin(ttphi-fRPAngle);
fedb5a47 1136 Double_t dphitrigjet=RelativePhi(ttphi,phiJet);
4f6d1cc1 1137 if(fTRP==1){
fedb5a47 1138 if(TMath::Abs(dphitrigjet)<TMath::Pi()-0.6) continue;
1139 var1[1] = ptJet;
1140 var1[4] = phiBintt;
1141 var1[5] = jet->EffectiveAreaCharged();
1142 var1[7] = ttpt;
1143 var1[8] = CheckAcceptance(phiJet,etaJet);
1144 }
4f6d1cc1 1145
fedb5a47 1146 if(fTRP==0){
1147 var1[1] = ptJet;
1148 var1[4] = phiBin;
1149 var1[5] = jet->EffectiveAreaCharged();
1150 var1[7] = ptLead;
1151 var1[8] = CheckAcceptance(phiJet,etaJet);
1152 }
1153
1154 //jet number;p_{T,jet};cent;#eta;#phi;area;trigger;acceptance bin;signed pt leading
adf3920d 1155 var2[1] = ptJet;
1156 var2[3] = etaJet;
1157 var2[4] = phiJet;
1158 var2[5] = jet->EffectiveAreaCharged();
1159 var2[7] = var1[8];
1160 var2[8] = (leadTrack?leadTrack->Charge()*ptLead:0);//pT of leading jet x charge
1161
1162 if(ij<kMaxJets){
1163 fh2LTrackPtJetPt[iType][ij]->Fill(jet->GetPtLeading(),ptJet);
1164 var1[0] = ij;
1165 var2[0] = ij;
3d1dba29 1166 for(int it = 0;it <fNTrigger;it++){
1167 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1168 var1[6] = it;
adf3920d 1169 var2[6] = it;
3d1dba29 1170 fhnJetPt[iType]->Fill(var1);
1d478650 1171 fhnJetPtQA[iType]->Fill(var2);
3d1dba29 1172 }
1173 }
adf3920d 1174 }
1175 var1[0] = kMaxJets;// fill for all jets
1176 var2[0] = kMaxJets;// fill for all jets
1177 for(int it = 0;it <fNTrigger;it++){
1178 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1179 var1[6] = it;
94c3d5e4 1180 var2[6] = it;
adf3920d 1181 fhnJetPt[iType]->Fill(var1);
1182 fhnJetPtQA[iType]->Fill(var2);
1183 }
1184 }
3d1dba29 1185
adf3920d 1186 fh2LTrackPtJetPt[iType][kMaxJets]->Fill(jet->GetPtLeading(),ptJet);
a31b8a87 1187
adf3920d 1188 if(particlesList.GetSize()&&ij<kMaxJets){
1189 // Particles... correlated with jets...
1190 for(int it = 0;it<particlesList.GetEntries();++it){
1191 AliVParticle *part = (AliVParticle*)particlesList.At(it);
1192 Float_t deltaR = jet->DeltaR(part);
1193 if(ptJet>0)fh1TmpRho->Fill(deltaR,part->Pt()/ptJet);
1194 }
1195 // fill the jet shapes
1196 }// if we have particles
a31b8a87 1197 }// Jet Loop
1198
1199
1200 // do something with dijets...
1201 if(ij0>=0&&ij1>0){
1202 AliAODJet *jet0 = (AliAODJet*)jetsList.At(ij0);
1203 Double_t ptJet0 = jet0->Pt();
1204 Double_t phiJet0 = jet0->Phi();
1205 if(phiJet0<0)phiJet0+=TMath::Pi()*2.;
1206
1207 AliAODJet *jet1 = (AliAODJet*)jetsList.At(ij1);
1208 Double_t ptJet1 = jet1->Pt();
1209 Double_t phiJet1 = jet1->Phi();
1210 if(phiJet1<0)phiJet1+=TMath::Pi()*2.;
1211
1212 Float_t deltaPhi = phiJet0 - phiJet1;
1213 if(deltaPhi>TMath::Pi())deltaPhi = deltaPhi - 2.*TMath::Pi();
1214 if(deltaPhi<(-1.*TMath::Pi()))deltaPhi = deltaPhi + 2.*TMath::Pi();
1215 deltaPhi = TMath::Abs(deltaPhi);
b5d90baf 1216 fh2DijetDeltaPhiPt[iType]->Fill(deltaPhi,ptJet1);
a31b8a87 1217
1218 Float_t asym = 9999;
1219 if((ptJet0+ptJet1)>0)asym = (ptJet0-ptJet1)/(ptJet0+ptJet1);
adf3920d 1220 fh2DijetAsymPt[iType]->Fill(asym,ptJet0);
1221 fh2DijetPt2vsPt1[iType]->Fill(ptJet0,ptJet1);
1222 fh2DijetDifvsSum[iType]->Fill(ptJet0+ptJet1,ptJet0-ptJet1);
1223 Float_t minv = 2.*(jet0->P()*jet1->P()-
1224 jet0->Px()*jet1->Px()-
1225 jet0->Py()*jet1->Py()-
1226 jet0->Pz()*jet1->Pz()); // assume mass == 0;
1227 if(minv<0)minv=0; // prevent numerical instabilities
1228 minv = TMath::Sqrt(minv);
1229 fh1DijetMinv[iType]->Fill(minv);
a31b8a87 1230 }
1231
1232
1233
1234 // count down the jets above thrueshold
1235 Int_t nOver = nJets;
1236 if(nOver>0){
1237 TIterator *jetIter = jetsList.MakeIterator();
1238 AliAODJet *tmpJet = (AliAODJet*)(jetIter->Next());
a31b8a87 1239 if(tmpJet){
e855f5c5 1240 Float_t pt = tmpJet->Pt();
a31b8a87 1241 for(int i = 1;i <= fh2NJetsPt[iType]->GetNbinsX();i++){
1242 Float_t ptCut = fh2NJetsPt[iType]->GetXaxis()->GetBinCenter(i);
1243 while(pt<ptCut&&tmpJet){
1244 nOver--;
1245 tmpJet = (AliAODJet*)(jetIter->Next());
1246 if(tmpJet){
1247 pt = tmpJet->Pt();
1248 }
1249 }
1250 if(nOver<=0)break;
1251 fh2NJetsPt[iType]->Fill(ptCut,nOver);
1252 }
1253 }
1254 delete jetIter;
1255 }
1256}
1257
1258void AliAnalysisTaskJetSpectrum2::FillTrackHistos(TList &particlesList,int iType){
1259
90a006c2 1260 if(fFlagJetType[iType]<=0)return;
7b822bc5 1261 Int_t refMult = fMultRec;
1262 if(iType==kJetGen||iType==kJetGenFull){
1263 refMult = fMultGen;
742ee86c 1264
7b822bc5 1265 }
1266
90a006c2 1267 //
5a7c81dc 1268 Double_t var3[7]; // track number;p_{T};cent;#tracks;RP;cahrge
90a006c2 1269 var3[2] = fCentrality;
1270 var3[3] = refMult;
5a7c81dc 1271 Double_t var4[6]; // track number;p_{T};cent;#eta;#phi;sign
90a006c2 1272 var4[2] = fCentrality;
a31b8a87 1273 Int_t nTrackOver = particlesList.GetSize();
1274 // do the same for tracks and jets
1275 if(nTrackOver>0){
1276 TIterator *trackIter = particlesList.MakeIterator();
1277 AliVParticle *tmpTrack = (AliVParticle*)(trackIter->Next());
1278 Float_t pt = tmpTrack->Pt();
1279 for(int i = 1;i <= fh2NTracksPt[iType]->GetNbinsX();i++){
1280 Float_t ptCut = fh2NTracksPt[iType]->GetXaxis()->GetBinCenter(i);
1281 while(pt<ptCut&&tmpTrack){
1282 nTrackOver--;
1283 tmpTrack = (AliVParticle*)(trackIter->Next());
1284 if(tmpTrack){
1285 pt = tmpTrack->Pt();
1286 }
1287 }
1288 if(nTrackOver<=0)break;
1289 fh2NTracksPt[iType]->Fill(ptCut,nTrackOver);
1290 }
1291
1292
1293 trackIter->Reset();
1294 AliVParticle *leading = (AliVParticle*)particlesList.At(0);
1295 Float_t sumPt = 0;
1296
1297 while((tmpTrack = (AliVParticle*)(trackIter->Next()))){
1298 Float_t tmpPt = tmpTrack->Pt();
1299 fh1PtTracksIn[iType]->Fill(tmpPt);
cb883243 1300 fh1PtTracksInLow[iType]->Fill(tmpPt);
742ee86c 1301
a31b8a87 1302 sumPt += tmpPt;
1303 Float_t tmpPhi = tmpTrack->Phi();
1304 if(tmpPhi<0)tmpPhi+=TMath::Pi()*2.;
2eded560 1305
1306
1307 Float_t phiRP = tmpPhi-fRPAngle;
1308 if(phiRP>TMath::Pi())phiRP -= TMath::Pi();
1309 if(phiRP<0)phiRP += TMath::Pi();
1310 if(phiRP<0)phiRP += TMath::Pi();
1311 const float allPhi = -1./180.*TMath::Pi();
1312
dae7dd67 1313 if(tmpPt<100){
1314 fp2MultRPPhiTrackPt[iType]->Fill(refMult,phiRP,tmpPt);
1315 fp2MultRPPhiTrackPt[iType]->Fill(refMult,allPhi,tmpPt);
1316
1317 fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,phiRP,tmpPt);
1318 fp2CentRPPhiTrackPt[iType]->Fill(fCentrality,allPhi,tmpPt);
1319 }
742ee86c 1320 Int_t phiBin = GetPhiBin(tmpPhi-fRPAngle);
90a006c2 1321 var3[0] = 1;
1322 var3[1] = tmpPt;
1323 var3[4] = phiBin;
5a7c81dc 1324 var3[6] = tmpTrack->Charge();
90a006c2 1325
1326 var4[0] = 1;
1327 var4[1] = tmpPt;
1328 var4[3] = tmpTrack->Eta();
1329 var4[4] = tmpPhi;
5662640e 1330 var4[5] = tmpTrack->Charge();
90a006c2 1331
17f4943e 1332
1333 for(int it = 0;it <fNTrigger;it++){
1334 if(fInputHandler->IsEventSelected()&fTriggerBit[it]){
1335 var3[0] = 1;
1336 var3[5] = it;
1337 fhnTrackPt[iType]->Fill(var3);
1338 if(tmpTrack==leading){
1339 var3[0] = 0;
1340 fhnTrackPt[iType]->Fill(var3);
1341 }
1342 }
1343 }
1344
1345
1346
5a7c81dc 1347
1348 fhnTrackPtQA[iType]->Fill(var4);
1349 if(tmpTrack==leading){
1350 var4[0] = 0;
1617846c 1351 fhnTrackPtQA[iType]->Fill(var4);
5a7c81dc 1352 continue;
1617846c 1353 }
17f4943e 1354
1355
90a006c2 1356
1617846c 1357
1358
a31b8a87 1359 }
1360 fh1SumPtTrack[iType]->Fill(sumPt);
1361 delete trackIter;
1362 }
1363
1364}
1365
1366
1367void AliAnalysisTaskJetSpectrum2::FillMatchHistos(TList &recJetsList,TList &genJetsList){
1368
1369
1370 // Fill al the matching histos
1371 // It is important that the acceptances for the mathing are as large as possible
1372 // to avoid false matches on the edge of acceptance
1373 // therefore we add some extra matching jets as overhead
1374
1375 static TArrayI aGenIndex(recJetsList.GetEntries());
1376 if(aGenIndex.GetSize()<recJetsList.GetEntries())aGenIndex.Set(recJetsList.GetEntries());
3b7ffecf 1377
a31b8a87 1378 static TArrayI aRecIndex(genJetsList.GetEntries());
1379 if(aRecIndex.GetSize()<genJetsList.GetEntries())aRecIndex.Set(genJetsList.GetEntries());
3b7ffecf 1380
a31b8a87 1381 if(fDebug){
1382 Printf("New Gens List %d rec index Array %d",genJetsList.GetEntries(),aRecIndex.GetSize());
1383 Printf("New Rec List %d gen indey Array %d",recJetsList.GetEntries(),aGenIndex.GetSize());
1384 }
1385 AliAnalysisHelperJetTasks::GetClosestJets(&genJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)genJetsList.GetEntries()),
1386 &recJetsList,TMath::Min((Int_t)fNMatchJets,(Int_t)recJetsList.GetEntries()),
1387 aGenIndex,aRecIndex,fDebug);
1388
1389 if(fDebug){
1390 for(int i = 0;i< aGenIndex.GetSize();++i){
1391 if(aGenIndex[i]>=0)Printf("iGenFound: %d -> %d",i,aGenIndex[i]);
1392 }
1393 for(int i = 0;i< aRecIndex.GetSize();++i){
1394 if(aRecIndex[i]>=0)Printf("iRecFound: %d -> %d",i,aRecIndex[i]);
1395 }
1396 }
1397
152a0744 1398 Double_t container[8];
1399 Double_t containerRec[4];
1400 Double_t containerGen[4];
a31b8a87 1401
1402 // loop over generated jets
1403 // consider the
1404 for(int ig = 0;ig < genJetsList.GetEntries();++ig){
1405 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1406 Double_t ptGen = genJet->Pt();
1407 Double_t phiGen = genJet->Phi();
1408 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1409 Double_t etaGen = genJet->Eta();
152a0744 1410 containerGen[0] = ptGen;
1411 containerGen[1] = etaGen;
1412 containerGen[2] = phiGen;
1413 containerGen[3] = genJet->GetPtLeading();
1414
1415 fhnJetContainer->Fill(containerGen,kStep0); //all generated jets
1416
346e5a93 1417 if(JetSelected(genJet))
1418 fhnJetContainer->Fill(containerGen,kStep1); // all generated jets in eta window
1419
152a0744 1420 Int_t ir = aRecIndex[ig];
1421 if(ir>=0&&ir<recJetsList.GetEntries()){
346e5a93 1422 AliAODJet* recJet = (AliAODJet*)recJetsList.At(ir);
152a0744 1423
346e5a93 1424 fhnJetContainer->Fill(containerGen,kStep2); // all generated jets with reconstructed partner
152a0744 1425
346e5a93 1426 if(JetSelected(genJet)){
152a0744 1427 fhnJetContainer->Fill(containerGen,kStep3); // all generated jets in eta window with reconstructed partner
346e5a93 1428
152a0744 1429 if(JetSelected(recJet)) {
1430 fhnJetContainer->Fill(containerGen,kStep4); // all generated jets in eta window with reconstructed partner in eta window
346e5a93 1431
152a0744 1432 }
346e5a93 1433 containerGen[3] = recJet->GetPtLeading();
1434 fhnJetContainer->Fill(containerGen,kStep5); // all generated jets in eta window with reconstructed partner with leading track on reconstructed level
a31b8a87 1435 }
1436 }
1437 }// loop over generated jets used for matching...
1438
1439
1440
1441 // loop over reconstructed jets
1442 for(int ir = 0;ir < recJetsList.GetEntries();++ir){
1443 AliAODJet *recJet = (AliAODJet*)recJetsList.At(ir);
1444 Double_t etaRec = recJet->Eta();
1445 Double_t ptRec = recJet->Pt();
1446 Double_t phiRec = recJet->Phi();
1447 if(phiRec<0)phiRec+=TMath::Pi()*2.;
a31b8a87 1448
152a0744 1449 containerRec[0] = ptRec;
1450 containerRec[1] = etaRec;
1451 containerRec[2] = phiRec;
1452 containerRec[3] = recJet->GetPtLeading();
1453
a31b8a87 1454 container[0] = ptRec;
1455 container[1] = etaRec;
1456 container[2] = phiRec;
152a0744 1457 container[3] = recJet->GetPtLeading();
a31b8a87 1458
a31b8a87 1459 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1460
1461 if(JetSelected(recJet)){
152a0744 1462 fhnJetContainer->Fill(containerRec,kStep7); //all rec jets in eta window
a31b8a87 1463 // Fill Correlation
1464 Int_t ig = aGenIndex[ir];
1465 if(ig>=0 && ig<genJetsList.GetEntries()){
a31b8a87 1466 if (fDebug > 10)Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir);
1467 AliAODJet *genJet = (AliAODJet*)genJetsList.At(ig);
1468 Double_t ptGen = genJet->Pt();
1469 Double_t phiGen = genJet->Phi();
1470 if(phiGen<0)phiGen+=TMath::Pi()*2.;
1471 Double_t etaGen = genJet->Eta();
152a0744 1472
1473 containerGen[0] = ptGen;
1474 containerGen[1] = etaGen;
1475 containerGen[2] = phiGen;
1476 containerGen[3] = genJet->GetPtLeading();
1477
1478 container[4] = ptGen;
1479 container[5] = etaGen;
1480 container[6] = phiGen;
1481 container[7] = genJet->GetPtLeading();
1482
1483 fhnJetContainer->Fill(containerGen,kStep6); // all rec jets in eta window with generated partner
1484
1485 fhnCorrelation->Fill(container);
e39a4f22 1486
1487 Float_t delta = (ptRec-ptGen)/ptGen;
1488 fh2RelPtFGen->Fill(ptGen,delta);
1489 fh3RelPtFGenLeadTrkPt->Fill(ptGen,delta,recJet->GetPtLeading());
1490 fh2PtFGen->Fill(ptGen,ptRec);
1491
1492 fh2deltaPt1Pt2->Fill(ptRec-ptGen,ptRec,ptGen);
1493
a31b8a87 1494 }
a31b8a87 1495 }// loop over reconstructed jets
1496 }
1497 if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__);
1498}
1499
1500
3b7ffecf 1501void AliAnalysisTaskJetSpectrum2::MakeJetContainer(){
1502 //
1503 // Create the particle container for the correction framework manager and
1504 // link it
1505 //
152a0744 1506 const Int_t kNvar = 4 ; //number of variables on the grid:pt,eta, phi
d95fc15a 1507 const Double_t kPtmin = 0.0, kPtmax = 250.; // we do not want to have empty bins at the beginning...
152a0744 1508 const Double_t kEtamin = -1.0, kEtamax = 1.0;
3b7ffecf 1509 const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
152a0744 1510 const Double_t kPtLeadingTrackPtMin = 0., kPtLeadingTrackPtMax = 200.;
3b7ffecf 1511
1512 // can we neglect migration in eta and phi?
1513 // phi should be no problem since we cover full phi and are phi symmetric
1514 // eta migration is more difficult due to needed acceptance correction
1515 // in limited eta range
1516
3b7ffecf 1517 //arrays for the number of bins in each dimension
1518 Int_t iBin[kNvar];
d95fc15a 1519 iBin[0] = 125; //bins in pt
152a0744 1520 iBin[1] = 4; //bins in eta
1521 iBin[2] = 1; // bins in phi
1522 iBin[3] = 10; // bins in leading track Pt
3b7ffecf 1523
1524 //arrays for lower bounds :
1525 Double_t* binEdges[kNvar];
1526 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1527 binEdges[ivar] = new Double_t[iBin[ivar] + 1];
1528
1529 //values for bin lower bounds
1530 // 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 1531 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 1532 for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
1533 for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
152a0744 1534 binEdges[3][0]= kPtLeadingTrackPtMin;
1535 binEdges[3][1]= 1.;
1536 binEdges[3][2]= 2.;
1537 binEdges[3][3]= 3.;
1538 binEdges[3][4]= 4.;
1539 binEdges[3][5]= 5.;
1540 binEdges[3][6]= 6.;
1541 binEdges[3][7]= 8.;
1542 binEdges[3][8]= 10.;
1543 binEdges[3][9]= 12.;
1544 binEdges[3][10]= kPtLeadingTrackPtMax;
1545
1546 fhnJetContainer = new AliCFContainer(Form("fhnJetContainer"),Form("AliCFContainer jet info"),kMaxStep,kNvar,iBin);
b99018c6 1547 for (int k=0; k<kNvar; k++) {
1548 fhnJetContainer->SetBinLimits(k,binEdges[k]);
3b7ffecf 1549 }
b99018c6 1550
3b7ffecf 1551 //create correlation matrix for unfolding
1552 Int_t thnDim[2*kNvar];
1553 for (int k=0; k<kNvar; k++) {
1554 //first half : reconstructed
1555 //second half : MC
1556 thnDim[k] = iBin[k];
1557 thnDim[k+kNvar] = iBin[k];
1558 }
1559
152a0744 1560 fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparseF with correlations",2*kNvar,thnDim);
3b7ffecf 1561 for (int k=0; k<kNvar; k++) {
152a0744 1562 fhnCorrelation->SetBinEdges(k,binEdges[k]);
1563 fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
3b7ffecf 1564 }
742ee86c 1565
1566 for(Int_t ivar = 0; ivar < kNvar; ivar++)
1567 delete [] binEdges[ivar];
1568
1569
3b7ffecf 1570}
1571
1572void AliAnalysisTaskJetSpectrum2::Terminate(Option_t */*option*/)
1573{
adf3920d 1574 // Terminate analysis
1575 //
1576 if (fDebug > 1) printf("AnalysisJetSpectrum2: Terminate() \n");
3b7ffecf 1577}
1578
1579
1580Int_t AliAnalysisTaskJetSpectrum2::GetListOfTracks(TList *list,Int_t type){
1581
565584e8 1582 if(fDebug>2)Printf("%s:%d Selecting tracks with %d",(char*)__FILE__,__LINE__,type);
3b7ffecf 1583
1584 Int_t iCount = 0;
565584e8 1585 if(type==kTrackAOD){
3b7ffecf 1586 AliAODEvent *aod = 0;
565584e8 1587 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1588 else aod = AODEvent();
3b7ffecf 1589 if(!aod){
1590 return iCount;
1591 }
1592 for(int it = 0;it < aod->GetNumberOfTracks();++it){
f15c1f69 1593 AliAODTrack *tr = dynamic_cast<AliAODTrack*>(aod->GetTrack(it));
1594 if(!tr) AliFatal("Not a standard AOD");
3b7ffecf 1595 if((fFilterMask>0)&&!(tr->TestFilterBit(fFilterMask)))continue;
a31b8a87 1596 if(TMath::Abs(tr->Eta())>fTrackRecEtaWindow)continue;
1597 if(tr->Pt()<fMinTrackPt)continue;
38ecb6a5 1598 if(fDebug>0){
1599 if(tr->Pt()>20){
1600 Printf("High pT track found in Event %d with p_T, %E",(int)Entry(),tr->Pt());
fd5db301 1601 Printf("%s read event, %d",fInputHandler->GetTree()->GetCurrentFile()->GetName(),(Int_t)fInputHandler->GetTree()->GetReadEntry());
38ecb6a5 1602 tr->Print();
1603 // tr->Dump();
1604 AliESDEvent *fESD = dynamic_cast<AliESDEvent*> (InputEvent());
1605 if(fESD){
a2522483 1606 AliESDtrack *esdTr = (AliESDtrack*)fESD->GetTrack(TMath::Abs(tr->GetID()+1));
38ecb6a5 1607 esdTr->Print("");
bb3d13aa 1608 // esdTr->Dump();
38ecb6a5 1609 }
1610 }
1611 }
3b7ffecf 1612 list->Add(tr);
1613 iCount++;
1614 }
1615 }
1616 else if (type == kTrackKineAll||type == kTrackKineCharged){
1617 AliMCEvent* mcEvent = MCEvent();
1618 if(!mcEvent)return iCount;
1619 // we want to have alivpartilces so use get track
1620 for(int it = 0;it < mcEvent->GetNumberOfTracks();++it){
1621 if(!mcEvent->IsPhysicalPrimary(it))continue;
1622 AliMCParticle* part = (AliMCParticle*)mcEvent->GetTrack(it);
a31b8a87 1623 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 1624 if(type == kTrackKineAll){
1625 list->Add(part);
1626 iCount++;
1627 }
1628 else if(type == kTrackKineCharged){
1629 if(part->Particle()->GetPDG()->Charge()==0)continue;
1630 list->Add(part);
1631 iCount++;
1632 }
1633 }
1634 }
565584e8 1635 else if (type == kTrackAODMCCharged || type == kTrackAODMCAll || type == kTrackAODMCChargedAcceptance) {
1636 AliAODEvent *aod = 0;
5301f754 1637 if(fUseAODMCInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
565584e8 1638 else aod = AODEvent();
5010a3f7 1639 if(!aod)return iCount;
3b7ffecf 1640 TClonesArray *tca = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
1641 if(!tca)return iCount;
1642 for(int it = 0;it < tca->GetEntriesFast();++it){
1643 AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it));
e855f5c5 1644 if(!part)continue;
a31b8a87 1645 if(part->Pt()<fMinTrackPt)continue;
3b7ffecf 1646 if(!part->IsPhysicalPrimary())continue;
c4631cdb 1647 if(type == kTrackAODMCAll){
3b7ffecf 1648 list->Add(part);
1649 iCount++;
1650 }
565584e8 1651 else if (type == kTrackAODMCCharged || type == kTrackAODMCChargedAcceptance ){
3b7ffecf 1652 if(part->Charge()==0)continue;
565584e8 1653 if(kTrackAODMCCharged){
1654 list->Add(part);
1655 }
1656 else {
a31b8a87 1657 if(TMath::Abs(part->Eta())>fTrackRecEtaWindow)continue;
565584e8 1658 list->Add(part);
1659 }
3b7ffecf 1660 iCount++;
1661 }
1662 }
1663 }// AODMCparticle
cc0649e4 1664 list->Sort();
c4631cdb 1665 return iCount;
3b7ffecf 1666
1667}
a31b8a87 1668
1669
742ee86c 1670Float_t AliAnalysisTaskJetSpectrum2::GetCentrality(){
adf3920d 1671 AliAODEvent *aod = 0;
1672 if(fUseAODTrackInput)aod = dynamic_cast<AliAODEvent*>(InputEvent());
1673 else aod = AODEvent();
1674 if(!aod){
1675 return 101;
1676 }
0a918d8d 1677 return ((AliVAODHeader*)aod->GetHeader())->GetCentrality();
742ee86c 1678}
1679
1680
a31b8a87 1681
1682Bool_t AliAnalysisTaskJetSpectrum2::JetSelected(AliAODJet *jet){
1683 Bool_t selected = false;
1684
1685 if(!jet)return selected;
1686
1687 if(fabs(jet->Eta())<fJetRecEtaWindow&&jet->Pt()>fMinJetPt){
1688 selected = kTRUE;
1689 }
1690 return selected;
1691
1692}
1693
1694Int_t AliAnalysisTaskJetSpectrum2::GetListOfJets(TList *list,TClonesArray* jarray,Int_t type){
1695
1696 if(fDebug>2)Printf("%s:%d Selecting jets with cuts %d",(char*)__FILE__,__LINE__,type);
1697 Int_t iCount = 0;
1698
1699 if(!jarray){
1700 Printf("%s:%d no Jet array",(char*)__FILE__,__LINE__);
1701 return iCount;
1702 }
1703
1704
1705 for(int ij=0;ij<jarray->GetEntries();ij++){
1706 AliAODJet* jet = (AliAODJet*)jarray->At(ij);
1707 if(!jet)continue;
1708 if(type==0){
1709 // No cut at all, main purpose here is sorting
1710 list->Add(jet);
1711 iCount++;
1712 }
1713 else if(type == 1){
1714 // eta cut
1715 if(JetSelected(jet)){
1716 list->Add(jet);
1717 iCount++;
1718 }
1719 }
1720 }
1721
1722 list->Sort();
1723 return iCount;
1724
1725}
1726
1727
37eb26ea 1728Int_t AliAnalysisTaskJetSpectrum2::MultFromJetRefs(TClonesArray* jets){
1729 if(!jets)return 0;
1730
1731 Int_t refMult = 0;
1732 for(int ij = 0;ij < jets->GetEntries();++ij){
1733 AliAODJet* jet = (AliAODJet*)jets->At(ij);
1734 if(!jet)continue;
1735 TRefArray *refs = jet->GetRefTracks();
1736 if(!refs)continue;
1737 refMult += refs->GetEntries();
1738 }
1739 return refMult;
1740
1741}
d7117cbd 1742
1743
1744AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackFromJetRefs(AliAODJet* jet){
1745 if(!jet)return 0;
1746 TRefArray *refs = jet->GetRefTracks();
1747 if(!refs) return 0;
1748 AliVParticle *leading = 0;
1749 Float_t fMaxPt = 0;
1750 for(int i = 0;i<refs->GetEntries();i++){
1751 AliVParticle *tmp = dynamic_cast<AliVParticle*>(refs->At(i));
1752 if(!tmp)continue;
1753 if(tmp->Pt()>fMaxPt){
1754 leading = tmp;
1755 fMaxPt = tmp->Pt();
1756 }
1757 }
1758 return leading;
1759}
1760
1761
1762AliVParticle *AliAnalysisTaskJetSpectrum2::LeadingTrackInCone(AliAODJet* jet,TList *list,Float_t r){
1763 if(!jet)return 0;
1764 if(!list) return 0;
1765 AliVParticle *leading = 0;
1766 Float_t fMaxPt = 0;
1767 for(int i = 0;i<list->GetEntries();i++){
1768 AliVParticle *tmp = (AliVParticle*)(list->At(i));
1769 if(!tmp)continue;
1770 if(jet->DeltaR(tmp)>r)continue;
1771 if(tmp->Pt()>fMaxPt){
1772 leading = tmp;
1773 fMaxPt = tmp->Pt();
1774 }
1775 }
1776 return leading;
1777}
742ee86c 1778
4f6d1cc1 1779
1780Double_t AliAnalysisTaskJetSpectrum2::RelativePhi(Double_t mphi,Double_t vphi){
1781
1782 if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
1783 else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
1784 if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
1785 else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
1786 double dphi = mphi-vphi;
1787 if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
1788 else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
1789 return dphi;//dphi in [-Pi, Pi]
1790}
1791
1792
1793
1794
1795
1796
742ee86c 1797Int_t AliAnalysisTaskJetSpectrum2::GetPhiBin(Double_t phi)
1798{
adf3920d 1799 Int_t phibin=-1;
1800 if(!(TMath::Abs(phi)<=2*TMath::Pi())){AliError("phi w.r.t. RP out of defined range");return -1;}
1801 Double_t phiwrtrp=TMath::ACos(TMath::Abs(TMath::Cos(phi)));
1802 phibin=Int_t(fNRPBins*phiwrtrp/(0.5*TMath::Pi()));
1803 if(phibin<0||phibin>=fNRPBins){AliError("Phi Bin not defined");}
1804 return phibin;
742ee86c 1805}
1806
3d1dba29 1807void AliAnalysisTaskJetSpectrum2::SetNTrigger(Int_t n){
db745d29 1808 //
1809 // set number of trigger bins
1810 //
3d1dba29 1811 if(n>0){
adf3920d 1812 fNTrigger = n;
1813 delete [] fTriggerName;
1814 fTriggerName = new TString [fNTrigger];
1815 delete [] fTriggerBit;fTriggerBit = 0;
1816 fTriggerBit = new UInt_t [fNTrigger];
3d1dba29 1817 }
1818 else{
1819 fNTrigger = 0;
1820 }
1821}
1822
db745d29 1823
3d1dba29 1824void AliAnalysisTaskJetSpectrum2::SetTrigger(Int_t i,UInt_t it,const char* c){
db745d29 1825 //
1826 // set trigger bin
1827 //
3d1dba29 1828 if(i<fNTrigger){
1829 fTriggerBit[i] = it;
1830 fTriggerName[i] = c;
1831 }
1832}
1833
db745d29 1834
1835void AliAnalysisTaskJetSpectrum2::SetNAcceptance(Int_t n){
1836 //
1837 // Set number of acceptance bins
1838 //
1839
1840 if(n>0){
adf3920d 1841 fNAcceptance = n;
1842 delete [] fAcceptancePhiMin;
1843 delete [] fAcceptancePhiMax;
1844 delete [] fAcceptanceEtaMin;
1845 delete [] fAcceptanceEtaMax;
1846
1847 fAcceptancePhiMin = new Float_t[fNAcceptance];
1848 fAcceptancePhiMax = new Float_t[fNAcceptance];
1849 fAcceptanceEtaMin = new Float_t[fNAcceptance];
1850 fAcceptanceEtaMax = new Float_t[fNAcceptance];
db745d29 1851 }
1852 else{
1853 fNTrigger = 0;
1854 }
1855}
1856
1857void AliAnalysisTaskJetSpectrum2::SetAcceptance(Int_t i,Float_t phiMin,Float_t phiMax,Float_t etaMin,Float_t etaMax){
1858 //
1859 // Set acceptance bins
1860 //
1861
1862 if(i<fNAcceptance){
24cc3482 1863 Printf("%s:%d Setting acceptance %d phi %3.3f-%3.3f eta %3.3f-%3.3f",(char*)__FILE__,__LINE__,i,phiMin,phiMax,etaMin,etaMax);
1864
1865 fAcceptancePhiMin[i] = phiMin;
1866 fAcceptancePhiMax[i] = phiMax;
1867 fAcceptanceEtaMin[i] = etaMin;
1868 fAcceptanceEtaMax[i] = etaMax;
db745d29 1869 }
1870 else{
1871 fNTrigger = 0;
1872 }
1873}
1874
1875
1876Int_t AliAnalysisTaskJetSpectrum2::CheckAcceptance(Float_t phi,Float_t eta){
1877
1878 //
1879 // return aceptnace bin do no use ovelapping bins
1880 //
1881
24cc3482 1882 for(int i = 0;i<fNAcceptance;i++){
db745d29 1883 if(phi<fAcceptancePhiMin[i])continue;
1884 if(phi>fAcceptancePhiMax[i])continue;
1885 if(eta<fAcceptanceEtaMin[i])continue;
1886 if(eta>fAcceptanceEtaMax[i])continue;
1887 return i;
1888 }
1889 // catch the rest
1890 return fNAcceptance;
1891}
1892
e5ee1c19 1893Float_t AliAnalysisTaskJetSpectrum2::GetRho(TList &list){
db745d29 1894
e5ee1c19 1895 // invert the correction
1896 AliAODJet *jet = (AliAODJet*)list.At(0); // highest pt jet
1897 if(!jet)return -1;
1898 if(jet->EffectiveAreaCharged()<=0)return -1;
1899 Float_t rho = jet->ChargedBgEnergy()/jet->EffectiveAreaCharged();
1900 return rho;
1901
1902}
db745d29 1903
1904
1905
3d1dba29 1906AliAnalysisTaskJetSpectrum2::~AliAnalysisTaskJetSpectrum2(){
1907 //
1908 delete [] fTriggerBit;
1909 delete [] fTriggerName;
db745d29 1910
1911 delete [] fAcceptancePhiMin;
1912 delete [] fAcceptancePhiMax;
1913 delete [] fAcceptanceEtaMin;
1914 delete [] fAcceptanceEtaMax;
1915
1916
1917
3d1dba29 1918}