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