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