]>
Commit | Line | Data |
---|---|---|
ad869500 | 1 | |
2 | // ****************************************** | |
3 | // This task computes several jet observables like | |
4 | // the fraction of energy in inner and outer coronnas, | |
5 | // jet-track correlations,triggered jet shapes and | |
6 | // correlation strength distribution of particles inside jets. | |
7 | // Author: lcunquei@cern.ch | |
8 | // ******************************************* | |
9 | ||
10 | ||
11 | /************************************************************************** | |
12 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
13 | * * | |
14 | * Author: The ALICE Off-line Project. * | |
15 | * Contributors are mentioned in the code where appropriate. * | |
16 | * * | |
17 | * Permission to use, copy, modify and distribute this software and its * | |
18 | * documentation strictly for non-commercial purposes is hereby granted * | |
19 | * without fee, provided that the above copyright notice appears in all * | |
20 | * copies and that both the copyright notice and this permission notice * | |
21 | * appear in the supporting documentation. The authors make no claims * | |
22 | * about the suitability of this software for any purpose. It is * | |
23 | * provided "as is" without express or implied warranty. * | |
24 | **************************************************************************/ | |
25 | ||
26 | ||
27 | #include "TChain.h" | |
28 | #include "TTree.h" | |
80ac66f6 | 29 | #include "TList.h" |
ad869500 | 30 | #include "TMath.h" |
31 | #include "TH1F.h" | |
32 | #include "TH1D.h" | |
33 | #include "TH2F.h" | |
34 | #include "TH3F.h" | |
35 | #include "THnSparse.h" | |
36 | #include "TCanvas.h" | |
80ac66f6 | 37 | #include "TArrayI.h" |
ea603b64 | 38 | #include "TProfile.h" |
39 | #include "TFile.h" | |
40 | #include "TKey.h" | |
8e103a56 | 41 | #include "TRandom3.h" |
ad869500 | 42 | |
43 | #include "AliLog.h" | |
44 | ||
45 | #include "AliAnalysisTask.h" | |
46 | #include "AliAnalysisManager.h" | |
47 | ||
48 | #include "AliVEvent.h" | |
49 | #include "AliESDEvent.h" | |
80ac66f6 | 50 | #include "AliMCEvent.h" |
ad869500 | 51 | #include "AliESDInputHandler.h" |
52 | #include "AliCentrality.h" | |
53 | #include "AliAnalysisHelperJetTasks.h" | |
54 | #include "AliInputEventHandler.h" | |
55 | #include "AliAODJetEventBackground.h" | |
ea603b64 | 56 | #include "AliGenPythiaEventHeader.h" |
ad869500 | 57 | #include "AliAODMCParticle.h" |
80ac66f6 | 58 | #include "AliMCParticle.h" |
ad869500 | 59 | #include "AliAODEvent.h" |
60 | #include "AliAODHandler.h" | |
61 | #include "AliAODJet.h" | |
62 | ||
63 | #include "AliAnalysisTaskJetCorePP.h" | |
64 | ||
65 | using std::cout; | |
66 | using std::endl; | |
67 | ||
68 | ClassImp(AliAnalysisTaskJetCorePP) | |
69 | ||
70 | //Filip Krizek 1st March 2013 | |
71 | ||
72 | //--------------------------------------------------------------------- | |
73 | AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP() : | |
74 | AliAnalysisTaskSE(), | |
75 | fESD(0x0), | |
76 | fAODIn(0x0), | |
77 | fAODOut(0x0), | |
78 | fAODExtension(0x0), | |
79 | fJetBranchName(""), | |
80ac66f6 | 80 | fJetBranchNameMC(""), |
7fc3d134 | 81 | fListJets(0x0), |
80ac66f6 | 82 | fListJetsGen(0x0), |
ad869500 | 83 | fNonStdFile(""), |
84 | fSystem(0), //pp=0 pPb=1 | |
85 | fJetParamR(0.4), | |
86 | fOfflineTrgMask(AliVEvent::kAny), | |
87 | fMinContribVtx(1), | |
88 | fVtxZMin(-10.0), | |
89 | fVtxZMax(10.0), | |
90 | fFilterMask(0), | |
91 | fCentMin(0.0), | |
92 | fCentMax(100.0), | |
93 | fJetEtaMin(-0.5), | |
94 | fJetEtaMax(0.5), | |
95 | fTriggerEtaCut(0.9), | |
96 | fTrackEtaCut(0.9), | |
97 | fTrackLowPtCut(0.15), | |
98 | fOutputList(0x0), | |
99 | fHistEvtSelection(0x0), | |
100 | fh2Ntriggers(0x0), | |
101 | fHJetSpec(0x0), | |
102 | fHJetDensity(0x0), | |
103 | fHJetDensityA4(0x0), | |
104 | fhJetPhi(0x0), | |
105 | fhTriggerPhi(0x0), | |
106 | fhJetEta(0x0), | |
107 | fhTriggerEta(0x0), | |
108 | fhVertexZ(0x0), | |
109 | fhVertexZAccept(0x0), | |
110 | fhContribVtx(0x0), | |
111 | fhContribVtxAccept(0x0), | |
112 | fhDphiTriggerJet(0x0), | |
113 | fhDphiTriggerJetAccept(0x0), | |
114 | fhCentrality(0x0), | |
115 | fhCentralityAccept(0x0), | |
7fc3d134 | 116 | fHJetPtRaw(0x0), |
117 | fHLeadingJetPtRaw(0x0), | |
118 | fHDphiVsJetPtAll(0x0), | |
80ac66f6 | 119 | fhJetPtGenVsJetPtRec(0x0), |
120 | fhJetPtGen(0x0), | |
121 | fh2NtriggersGen(0x0), | |
122 | fHJetSpecGen(0x0), | |
123 | fhPtTrkTruePrimRec(0x0), | |
124 | fhPtTrkTruePrimGen(0x0), | |
125 | fhPtTrkSecOrFakeRec(0x0), | |
126 | fIsMC(0), | |
127 | faGenIndex(0), | |
128 | faRecIndex(0), | |
7fc3d134 | 129 | fkAcceptance(2.0*TMath::Pi()*1.8), |
6168afc3 | 130 | fkDeltaPhiCut(TMath::Pi()-0.8), |
ea603b64 | 131 | fh1Xsec(0x0), |
132 | fh1Trials(0x0), | |
133 | fh1AvgTrials(0x0), | |
134 | fh1PtHard(0x0), | |
135 | fh1PtHardNoW(0x0), | |
136 | fh1PtHardTrials(0x0), | |
8e103a56 | 137 | fAvgTrials(1), |
138 | fHardest(0), | |
139 | fEventNumberRangeLow(0), | |
140 | fEventNumberRangeHigh(99), | |
141 | fTriggerPtRangeLow(0.0), | |
142 | fTriggerPtRangeHigh(10000.0), | |
143 | fRandom(0x0) | |
ad869500 | 144 | { |
145 | // default Constructor | |
ad869500 | 146 | } |
147 | ||
148 | //--------------------------------------------------------------------- | |
149 | ||
150 | AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const char *name) : | |
151 | AliAnalysisTaskSE(name), | |
152 | fESD(0x0), | |
153 | fAODIn(0x0), | |
154 | fAODOut(0x0), | |
155 | fAODExtension(0x0), | |
156 | fJetBranchName(""), | |
80ac66f6 | 157 | fJetBranchNameMC(""), |
7fc3d134 | 158 | fListJets(0x0), |
80ac66f6 | 159 | fListJetsGen(0x0), |
ad869500 | 160 | fNonStdFile(""), |
161 | fSystem(0), //pp=0 pPb=1 | |
162 | fJetParamR(0.4), | |
163 | fOfflineTrgMask(AliVEvent::kAny), | |
164 | fMinContribVtx(1), | |
165 | fVtxZMin(-10.0), | |
166 | fVtxZMax(10.0), | |
167 | fFilterMask(0), | |
168 | fCentMin(0.0), | |
169 | fCentMax(100.0), | |
170 | fJetEtaMin(-0.5), | |
171 | fJetEtaMax(0.5), | |
172 | fTriggerEtaCut(0.9), | |
173 | fTrackEtaCut(0.9), | |
174 | fTrackLowPtCut(0.15), | |
175 | fOutputList(0x0), | |
176 | fHistEvtSelection(0x0), | |
177 | fh2Ntriggers(0x0), | |
178 | fHJetSpec(0x0), | |
179 | fHJetDensity(0x0), | |
180 | fHJetDensityA4(0x0), | |
181 | fhJetPhi(0x0), | |
182 | fhTriggerPhi(0x0), | |
183 | fhJetEta(0x0), | |
184 | fhTriggerEta(0x0), | |
185 | fhVertexZ(0x0), | |
186 | fhVertexZAccept(0x0), | |
187 | fhContribVtx(0x0), | |
188 | fhContribVtxAccept(0x0), | |
189 | fhDphiTriggerJet(0x0), | |
190 | fhDphiTriggerJetAccept(0x0), | |
191 | fhCentrality(0x0), | |
192 | fhCentralityAccept(0x0), | |
7fc3d134 | 193 | fHJetPtRaw(0x0), |
194 | fHLeadingJetPtRaw(0x0), | |
195 | fHDphiVsJetPtAll(0x0), | |
80ac66f6 | 196 | fhJetPtGenVsJetPtRec(0x0), |
197 | fhJetPtGen(0x0), | |
198 | fh2NtriggersGen(0x0), | |
199 | fHJetSpecGen(0x0), | |
200 | fhPtTrkTruePrimRec(0x0), | |
201 | fhPtTrkTruePrimGen(0x0), | |
202 | fhPtTrkSecOrFakeRec(0x0), | |
203 | fIsMC(0), | |
204 | faGenIndex(0), | |
205 | faRecIndex(0), | |
7fc3d134 | 206 | fkAcceptance(2.0*TMath::Pi()*1.8), |
6168afc3 | 207 | fkDeltaPhiCut(TMath::Pi()-0.8), |
ea603b64 | 208 | fh1Xsec(0x0), |
209 | fh1Trials(0x0), | |
210 | fh1AvgTrials(0x0), | |
211 | fh1PtHard(0x0), | |
212 | fh1PtHardNoW(0x0), | |
213 | fh1PtHardTrials(0x0), | |
8e103a56 | 214 | fAvgTrials(1), |
215 | fHardest(0), | |
216 | fEventNumberRangeLow(0), | |
217 | fEventNumberRangeHigh(99), | |
218 | fTriggerPtRangeLow(0.0), | |
219 | fTriggerPtRangeHigh(10000.0), | |
220 | fRandom(0x0) | |
ad869500 | 221 | { |
7fc3d134 | 222 | // Constructor |
ad869500 | 223 | |
224 | DefineOutput(1, TList::Class()); | |
225 | } | |
226 | ||
227 | //-------------------------------------------------------------- | |
228 | AliAnalysisTaskJetCorePP::AliAnalysisTaskJetCorePP(const AliAnalysisTaskJetCorePP& a): | |
229 | AliAnalysisTaskSE(a.GetName()), | |
230 | fESD(a.fESD), | |
231 | fAODIn(a.fAODIn), | |
232 | fAODOut(a.fAODOut), | |
233 | fAODExtension(a.fAODExtension), | |
234 | fJetBranchName(a.fJetBranchName), | |
80ac66f6 | 235 | fJetBranchNameMC(a.fJetBranchNameMC), |
ad869500 | 236 | fListJets(a.fListJets), |
80ac66f6 | 237 | fListJetsGen(a.fListJetsGen), |
ad869500 | 238 | fNonStdFile(a.fNonStdFile), |
239 | fSystem(a.fSystem), | |
240 | fJetParamR(a.fJetParamR), | |
241 | fOfflineTrgMask(a.fOfflineTrgMask), | |
242 | fMinContribVtx(a.fMinContribVtx), | |
243 | fVtxZMin(a.fVtxZMin), | |
244 | fVtxZMax(a.fVtxZMax), | |
245 | fFilterMask(a.fFilterMask), | |
246 | fCentMin(a.fCentMin), | |
247 | fCentMax(a.fCentMax), | |
248 | fJetEtaMin(a.fJetEtaMin), | |
249 | fJetEtaMax(a.fJetEtaMax), | |
250 | fTriggerEtaCut(a.fTriggerEtaCut), | |
251 | fTrackEtaCut(a.fTrackEtaCut), | |
252 | fTrackLowPtCut(a.fTrackLowPtCut), | |
253 | fOutputList(a.fOutputList), | |
254 | fHistEvtSelection(a.fHistEvtSelection), | |
255 | fh2Ntriggers(a.fh2Ntriggers), | |
256 | fHJetSpec(a.fHJetSpec), | |
257 | fHJetDensity(a.fHJetDensity), | |
258 | fHJetDensityA4(a.fHJetDensityA4), | |
259 | fhJetPhi(a.fhJetPhi), | |
260 | fhTriggerPhi(a.fhTriggerPhi), | |
261 | fhJetEta(a.fhJetEta), | |
262 | fhTriggerEta(a.fhTriggerEta), | |
263 | fhVertexZ(a.fhVertexZ), | |
264 | fhVertexZAccept(a.fhVertexZAccept), | |
265 | fhContribVtx(a.fhContribVtx), | |
266 | fhContribVtxAccept(a.fhContribVtxAccept), | |
267 | fhDphiTriggerJet(a.fhDphiTriggerJet), | |
268 | fhDphiTriggerJetAccept(a.fhDphiTriggerJetAccept), | |
269 | fhCentrality(a.fhCentrality), | |
270 | fhCentralityAccept(a.fhCentralityAccept), | |
7fc3d134 | 271 | fHJetPtRaw(a.fHJetPtRaw), |
272 | fHLeadingJetPtRaw(a.fHLeadingJetPtRaw), | |
273 | fHDphiVsJetPtAll(a.fHDphiVsJetPtAll), | |
80ac66f6 | 274 | fhJetPtGenVsJetPtRec(a.fhJetPtGenVsJetPtRec), |
275 | fhJetPtGen(a.fhJetPtGen), | |
276 | fh2NtriggersGen(a.fh2NtriggersGen), | |
277 | fHJetSpecGen(a.fHJetSpecGen), | |
278 | fhPtTrkTruePrimRec(a.fhPtTrkTruePrimRec), | |
279 | fhPtTrkTruePrimGen(a.fhPtTrkTruePrimGen), | |
280 | fhPtTrkSecOrFakeRec(a.fhPtTrkSecOrFakeRec), | |
281 | fIsMC(a.fIsMC), | |
282 | faGenIndex(a.faGenIndex), | |
283 | faRecIndex(a.faRecIndex), | |
7fc3d134 | 284 | fkAcceptance(a.fkAcceptance), |
ea603b64 | 285 | fkDeltaPhiCut(a.fkDeltaPhiCut), |
286 | fh1Xsec(a.fh1Xsec), | |
287 | fh1Trials(a.fh1Trials), | |
288 | fh1AvgTrials(a.fh1AvgTrials), | |
289 | fh1PtHard(a.fh1PtHard), | |
290 | fh1PtHardNoW(a.fh1PtHardNoW), | |
291 | fh1PtHardTrials(a.fh1PtHardTrials), | |
8e103a56 | 292 | fAvgTrials(a.fAvgTrials), |
293 | fHardest(a.fHardest), | |
294 | fEventNumberRangeLow(a.fEventNumberRangeLow), | |
295 | fEventNumberRangeHigh(a.fEventNumberRangeHigh), | |
296 | fTriggerPtRangeLow(a.fTriggerPtRangeLow), | |
43457761 | 297 | fTriggerPtRangeHigh(a.fTriggerPtRangeHigh), |
298 | fRandom(a.fRandom) | |
ad869500 | 299 | { |
300 | //Copy Constructor | |
43457761 | 301 | fRandom->SetSeed(0); |
ad869500 | 302 | } |
303 | //-------------------------------------------------------------- | |
304 | ||
305 | AliAnalysisTaskJetCorePP& AliAnalysisTaskJetCorePP::operator = (const AliAnalysisTaskJetCorePP& a){ | |
306 | // assignment operator | |
307 | this->~AliAnalysisTaskJetCorePP(); | |
308 | new(this) AliAnalysisTaskJetCorePP(a); | |
309 | return *this; | |
310 | } | |
311 | //-------------------------------------------------------------- | |
312 | ||
313 | AliAnalysisTaskJetCorePP::~AliAnalysisTaskJetCorePP() | |
314 | { | |
315 | //Destructor | |
316 | delete fListJets; | |
80ac66f6 | 317 | delete fListJetsGen; |
ad869500 | 318 | delete fOutputList; // ???? |
8e103a56 | 319 | delete fRandom; |
ad869500 | 320 | } |
321 | ||
322 | //-------------------------------------------------------------- | |
323 | ||
ea603b64 | 324 | |
325 | Bool_t AliAnalysisTaskJetCorePP::Notify() | |
326 | { | |
327 | //Implemented Notify() to read the cross sections | |
328 | //and number of trials from pyxsec.root | |
329 | //inspired by AliAnalysisTaskJetSpectrum2::Notify() | |
330 | if(!fIsMC) return kFALSE; | |
331 | ||
332 | fESD = dynamic_cast<AliESDEvent*>(InputEvent()); | |
333 | if(!fESD){ | |
334 | if(fDebug>1) AliError("ESD not available"); | |
335 | fAODIn = dynamic_cast<AliAODEvent*>(InputEvent()); | |
336 | } | |
337 | ||
338 | fAODOut = dynamic_cast<AliAODEvent*>(AODEvent()); | |
339 | ||
340 | ||
341 | if(fNonStdFile.Length()!=0){ | |
342 | // case that we have an AOD extension we can fetch the jets from the extended output | |
343 | AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); | |
344 | fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0; | |
345 | if(!fAODExtension){ | |
346 | if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data()); | |
347 | } | |
348 | } | |
349 | ||
350 | TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree(); | |
351 | Float_t xsection = 0; | |
352 | Float_t ftrials = 1; | |
353 | ||
354 | fAvgTrials = 1; | |
355 | if(tree){ | |
356 | TFile *curfile = tree->GetCurrentFile(); | |
357 | if(!curfile) { | |
358 | Error("Notify","No current file"); | |
359 | return kFALSE; | |
360 | } | |
361 | if(!fh1Xsec || !fh1Trials){ | |
362 | Printf("%s%d No Histogram fh1Xsec",(char*)__FILE__,__LINE__); | |
363 | return kFALSE; | |
364 | } | |
365 | AliAnalysisHelperJetTasks::PythiaInfoFromFile(curfile->GetName(),xsection,ftrials); | |
366 | fh1Xsec->Fill("<#sigma>",xsection); | |
367 | // construct a poor man average trials | |
368 | Float_t nEntries = (Float_t)tree->GetTree()->GetEntries(); | |
369 | if(ftrials>=nEntries && nEntries>0.) fAvgTrials = ftrials/nEntries; | |
370 | fh1Trials->Fill("#sum{ntrials}",ftrials); | |
371 | } | |
372 | ||
373 | if(fDebug)Printf("Reading File %s",fInputHandler->GetTree()->GetCurrentFile()->GetName()); | |
374 | ||
375 | return kTRUE; | |
376 | } | |
377 | //-------------------------------------------------------------- | |
378 | ||
ad869500 | 379 | void AliAnalysisTaskJetCorePP::Init() |
380 | { | |
381 | // check for jet branches | |
382 | if(!strlen(fJetBranchName.Data())){ | |
383 | AliError("Jet branch name not set."); | |
384 | } | |
ad869500 | 385 | } |
386 | ||
387 | //-------------------------------------------------------------- | |
388 | ||
389 | void AliAnalysisTaskJetCorePP::UserCreateOutputObjects() | |
390 | { | |
ad869500 | 391 | // Create histograms |
392 | // Called once | |
80ac66f6 | 393 | fListJets = new TList(); //reconstructed level |
394 | ||
395 | fIsMC = (fJetBranchNameMC.Length()>0) ? kTRUE : kFALSE; | |
396 | ||
8e103a56 | 397 | fRandom = new TRandom3(0); |
398 | ||
80ac66f6 | 399 | if(fIsMC) fListJetsGen = new TList(); //generator level |
ad869500 | 400 | OpenFile(1); |
401 | if(!fOutputList) fOutputList = new TList(); | |
402 | fOutputList->SetOwner(kTRUE); | |
403 | ||
404 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
405 | TH1::AddDirectory(kFALSE); | |
406 | ||
407 | fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5); | |
408 | fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED"); | |
409 | fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN"); | |
410 | fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)"); | |
411 | fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)"); | |
412 | fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)"); | |
8e103a56 | 413 | fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)"); |
ad869500 | 414 | |
415 | fOutputList->Add(fHistEvtSelection); | |
416 | ||
417 | Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10 | |
80ac66f6 | 418 | //trigger pt spectrum (reconstructed) |
ad869500 | 419 | fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers", |
420 | nBinsCentrality,0.0,100.0,50,0.0,50.0); | |
7fc3d134 | 421 | fOutputList->Add(fh2Ntriggers); |
ad869500 | 422 | |
6168afc3 | 423 | //Centrality, A, pTjet, pTtrigg, dphi |
424 | const Int_t dimSpec = 5; | |
425 | const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 100, 100, 50, TMath::Nint(10*(TMath::Pi()-fkDeltaPhiCut))}; | |
426 | const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, 0, 0.0, fkDeltaPhiCut}; | |
427 | const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0, TMath::Pi()}; | |
7fc3d134 | 428 | fHJetSpec = new THnSparseF("fHJetSpec", |
6168afc3 | 429 | "Recoil jet spectrum [cent,A,pTjet,pTtrig,dphi]", |
7fc3d134 | 430 | dimSpec,nBinsSpec,lowBinSpec,hiBinSpec); |
ad869500 | 431 | fOutputList->Add(fHJetSpec); |
432 | ||
433 | //------------------- HISTOS FOR DIAGNOSTIC ---------------------- | |
434 | //Jet number density histos [Trk Mult, jet density, pT trigger] | |
435 | const Int_t dimJetDens = 3; | |
80ac66f6 | 436 | const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10}; |
ad869500 | 437 | const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0}; |
80ac66f6 | 438 | const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0}; |
ad869500 | 439 | |
440 | fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07", | |
441 | dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens); | |
442 | ||
443 | fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4", | |
444 | dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens); | |
445 | ||
ad869500 | 446 | fOutputList->Add(fHJetDensity); |
447 | fOutputList->Add(fHJetDensityA4); | |
448 | ||
449 | ||
450 | //inclusive azimuthal and pseudorapidity histograms | |
7fc3d134 | 451 | fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet", |
452 | 50, 0, 100, 50,-TMath::Pi(),TMath::Pi()); | |
453 | fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg", | |
6168afc3 | 454 | 50, 0, 50, 50,-TMath::Pi(),TMath::Pi()); |
7fc3d134 | 455 | fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet", |
456 | 50,0, 100, 40,-0.9,0.9); | |
457 | fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg", | |
6168afc3 | 458 | 50, 0, 50, 40,-0.9,0.9); |
7fc3d134 | 459 | |
ad869500 | 460 | fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20); |
461 | fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20); | |
462 | fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200); | |
463 | fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200); | |
464 | fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi()); | |
465 | fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi()); | |
466 | fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100); | |
467 | fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100); | |
468 | ||
469 | fOutputList->Add(fhJetPhi); | |
470 | fOutputList->Add(fhTriggerPhi); | |
471 | fOutputList->Add(fhJetEta); | |
472 | fOutputList->Add(fhTriggerEta); | |
473 | fOutputList->Add(fhVertexZ); | |
474 | fOutputList->Add(fhVertexZAccept); | |
475 | fOutputList->Add(fhContribVtx); | |
476 | fOutputList->Add(fhContribVtxAccept); | |
477 | fOutputList->Add(fhDphiTriggerJet); | |
478 | fOutputList->Add(fhDphiTriggerJetAccept); | |
479 | fOutputList->Add(fhCentrality); | |
480 | fOutputList->Add(fhCentralityAccept); | |
481 | ||
7fc3d134 | 482 | // raw spectra of INCLUSIVE jets |
80ac66f6 | 483 | //Centrality, pTjet, A |
484 | const Int_t dimRaw = 3; | |
485 | const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100}; | |
486 | const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0}; | |
487 | const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0}; | |
7fc3d134 | 488 | fHJetPtRaw = new THnSparseF("fHJetPtRaw", |
80ac66f6 | 489 | "Incl. jet spectrum [cent,pTjet,A]", |
7fc3d134 | 490 | dimRaw,nBinsRaw,lowBinRaw,hiBinRaw); |
491 | fOutputList->Add(fHJetPtRaw); | |
492 | ||
493 | // raw spectra of LEADING jets | |
80ac66f6 | 494 | //Centrality, pTjet, A |
7fc3d134 | 495 | fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw", |
80ac66f6 | 496 | "Leading jet spectrum [cent,pTjet,A]", |
7fc3d134 | 497 | dimRaw,nBinsRaw,lowBinRaw,hiBinRaw); |
498 | fOutputList->Add(fHLeadingJetPtRaw); | |
499 | ||
500 | // Dphi versus pT jet | |
501 | //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg | |
502 | const Int_t dimDp = 4; | |
503 | const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50}; | |
504 | const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0}; | |
505 | const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0}; | |
506 | fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll", | |
507 | "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]", | |
508 | dimDp,nBinsDp,lowBinDp,hiBinDp); | |
509 | fOutputList->Add(fHDphiVsJetPtAll); | |
510 | ||
7fc3d134 | 511 | |
80ac66f6 | 512 | //analyze MC generator level |
513 | if(fIsMC){ | |
514 | fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 100,0,200, 100,0,200); | |
515 | fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC jet pt spectrum versus reconstructed pt spectrum | |
516 | ||
517 | fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",100,0,200); //MC generator jet pt spectrum | |
518 | fOutputList->Add(fhJetPtGen); | |
7fc3d134 | 519 | |
80ac66f6 | 520 | fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen"); |
521 | fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle())); | |
522 | fOutputList->Add(fh2NtriggersGen); | |
7fc3d134 | 523 | |
80ac66f6 | 524 | //Centrality, A, pT, pTtrigg |
525 | fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen"); | |
526 | fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle())); | |
527 | fOutputList->Add(fHJetSpecGen); | |
528 | ||
529 | //track efficiency/contamination histograms eta versus pT | |
530 | fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9); | |
531 | fOutputList->Add(fhPtTrkTruePrimRec); | |
532 | ||
533 | fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen"); | |
534 | fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen"); | |
535 | fOutputList->Add(fhPtTrkTruePrimGen); | |
536 | ||
537 | fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec"); | |
538 | fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec"); | |
539 | fOutputList->Add(fhPtTrkSecOrFakeRec); | |
540 | } | |
ea603b64 | 541 | //------------------------------------- |
542 | // pythia histograms | |
543 | const Int_t nBinPt = 150; | |
544 | Double_t binLimitsPt[nBinPt+1]; | |
545 | for(Int_t iPt = 0;iPt <= nBinPt;iPt++){ | |
546 | if(iPt == 0){ | |
547 | binLimitsPt[iPt] = -50.0; | |
548 | }else{// 1.0 | |
549 | binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.; | |
550 | } | |
551 | } | |
552 | ||
553 | fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1); | |
554 | fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>"); | |
555 | fOutputList->Add(fh1Xsec); | |
556 | fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1); | |
557 | fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}"); | |
558 | fOutputList->Add(fh1Trials); | |
559 | fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1); | |
560 | fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}"); | |
561 | fOutputList->Add(fh1AvgTrials); | |
562 | fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt); | |
563 | fOutputList->Add(fh1PtHard); | |
564 | fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt); | |
565 | fOutputList->Add(fh1PtHardNoW); | |
566 | fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt); | |
567 | fOutputList->Add(fh1PtHardTrials); | |
568 | ||
ad869500 | 569 | |
570 | // =========== Switch on Sumw2 for all histos =========== | |
571 | for(Int_t i=0; i<fOutputList->GetEntries(); i++){ | |
572 | TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i)); | |
573 | if(h1){ | |
574 | h1->Sumw2(); | |
575 | continue; | |
576 | } | |
577 | THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i)); | |
578 | if(hn){ | |
579 | hn->Sumw2(); | |
580 | } | |
581 | } | |
582 | TH1::AddDirectory(oldStatus); | |
583 | ||
584 | PostData(1, fOutputList); | |
585 | } | |
586 | ||
587 | //-------------------------------------------------------------------- | |
588 | ||
589 | void AliAnalysisTaskJetCorePP::UserExec(Option_t *) | |
590 | { | |
ad869500 | 591 | //Event loop |
ea603b64 | 592 | Double_t eventW = 1.0; |
593 | Double_t ptHard = 0.0; | |
594 | Double_t nTrials = 1.0; // Trials for MC trigger | |
595 | if(fIsMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials); | |
ad869500 | 596 | |
597 | if(TMath::Abs((Float_t) fJetParamR)<0.00001){ | |
80ac66f6 | 598 | AliError("Cone radius is set to zero."); |
ad869500 | 599 | return; |
600 | } | |
601 | if(!strlen(fJetBranchName.Data())){ | |
602 | AliError("Jet branch name not set."); | |
603 | return; | |
604 | } | |
605 | ||
606 | fESD = dynamic_cast<AliESDEvent*>(InputEvent()); | |
607 | if(!fESD){ | |
80ac66f6 | 608 | if(fDebug>1) AliError("ESD not available"); |
ad869500 | 609 | fAODIn = dynamic_cast<AliAODEvent*>(InputEvent()); |
610 | } | |
611 | ||
612 | fAODOut = dynamic_cast<AliAODEvent*>(AODEvent()); | |
613 | AliAODEvent* aod = NULL; | |
614 | // take all other information from the aod we take the tracks from | |
615 | if(!aod){ | |
616 | if(!fESD) aod = fAODIn; | |
617 | else aod = fAODOut; | |
618 | } | |
619 | ||
620 | if(fNonStdFile.Length()!=0){ | |
621 | // case that we have an AOD extension we can fetch the jets from the extended output | |
622 | AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); | |
623 | fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0; | |
624 | if(!fAODExtension){ | |
625 | if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data()); | |
626 | } | |
627 | } | |
628 | ||
80ac66f6 | 629 | // ----------------- EVENT SELECTION -------------------------- |
ad869500 | 630 | fHistEvtSelection->Fill(1); // number of events before event selection |
631 | ||
632 | // physics selection | |
633 | AliInputEventHandler* inputHandler = (AliInputEventHandler*) | |
634 | ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); | |
635 | ||
636 | if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){ | |
637 | if(fDebug) Printf(" Trigger Selection: event REJECTED ... "); | |
638 | fHistEvtSelection->Fill(2); | |
639 | PostData(1, fOutputList); | |
640 | return; | |
641 | } | |
642 | ||
643 | //check AOD pointer | |
644 | if(!aod){ | |
645 | if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__); | |
646 | fHistEvtSelection->Fill(3); | |
647 | PostData(1, fOutputList); | |
648 | return; | |
649 | } | |
650 | ||
651 | // vertex selection | |
652 | AliAODVertex* primVtx = aod->GetPrimaryVertex(); | |
653 | ||
654 | if(!primVtx){ | |
655 | if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__); | |
656 | fHistEvtSelection->Fill(3); | |
657 | PostData(1, fOutputList); | |
658 | return; | |
659 | } | |
660 | ||
661 | Int_t nTracksPrim = primVtx->GetNContributors(); | |
662 | Float_t vtxz = primVtx->GetZ(); | |
663 | //Input events | |
664 | fhContribVtx->Fill(nTracksPrim); | |
665 | if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz); | |
666 | ||
667 | if((nTracksPrim < fMinContribVtx) || | |
668 | (vtxz < fVtxZMin) || | |
669 | (vtxz > fVtxZMax)){ | |
670 | if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...", | |
671 | (char*)__FILE__,__LINE__,vtxz); | |
672 | fHistEvtSelection->Fill(3); | |
673 | PostData(1, fOutputList); | |
674 | return; | |
675 | }else{ | |
676 | //Accepted events | |
677 | fhContribVtxAccept->Fill(nTracksPrim); | |
678 | fhVertexZAccept->Fill(vtxz); | |
679 | } | |
680 | //FK// No event class selection imposed | |
681 | // event class selection (from jet helper task) | |
682 | //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass(); | |
683 | //if(fDebug) Printf("Event class %d", eventClass); | |
684 | //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){ | |
685 | // fHistEvtSelection->Fill(4); | |
686 | // PostData(1, fOutputList); | |
687 | // return; | |
688 | //} | |
689 | ||
80ac66f6 | 690 | //------------------ CENTRALITY SELECTION --------------- |
ad869500 | 691 | AliCentrality *cent = 0x0; |
692 | Double_t centValue = 0.0; | |
693 | if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb | |
694 | if(fESD){ | |
695 | cent = fESD->GetCentrality(); | |
696 | if(cent) centValue = cent->GetCentralityPercentile("V0M"); | |
697 | }else{ | |
698 | centValue = aod->GetHeader()->GetCentrality(); | |
699 | } | |
700 | if(fDebug) printf("centrality: %f\n", centValue); | |
701 | //Input events | |
702 | fhCentrality->Fill(centValue); | |
703 | ||
704 | if(centValue < fCentMin || centValue > fCentMax){ | |
705 | fHistEvtSelection->Fill(4); | |
706 | PostData(1, fOutputList); | |
707 | return; | |
708 | }else{ | |
709 | //Accepted events | |
710 | fhCentralityAccept->Fill( centValue ); | |
711 | } | |
712 | } | |
713 | ||
8e103a56 | 714 | //-----------------select disjunct event subsamples ---------------- |
715 | Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile(); | |
716 | Int_t lastdigit = eventnum % 10; | |
717 | if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){ | |
718 | fHistEvtSelection->Fill(5); | |
719 | PostData(1, fOutputList); | |
720 | return; | |
721 | } | |
722 | ||
ad869500 | 723 | if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl; |
ad869500 | 724 | fHistEvtSelection->Fill(0); // accepted events |
ad869500 | 725 | // ------------------- end event selection -------------------- |
80ac66f6 | 726 | |
ad869500 | 727 | |
80ac66f6 | 728 | // fetch RECONSTRUCTED jets |
ad869500 | 729 | TClonesArray *aodJets = 0x0; |
730 | ||
731 | if(fAODOut && !aodJets){ | |
732 | aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data())); | |
733 | } | |
734 | if(fAODExtension && !aodJets){ | |
735 | aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data())); | |
736 | } | |
737 | if(fAODIn && !aodJets){ | |
738 | aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data())); | |
739 | } | |
740 | ||
80ac66f6 | 741 | //--------- Fill list of RECONSTRUCTED jets -------------- |
742 | fListJets->Clear(); | |
743 | if(aodJets){ | |
744 | if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(), aodJets->GetEntriesFast()); | |
745 | for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) { | |
746 | AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]); | |
747 | if (jet) fListJets->Add(jet); | |
748 | } | |
749 | } | |
750 | ||
751 | //--------- Fill list of MC GENERATOR LEVEL jets -------------- | |
752 | TList particleListGen; //list of tracks in MC | |
753 | ||
754 | if(fIsMC){ //analyze generator level MC | |
755 | // fetch MC generator level jets | |
756 | TClonesArray *aodGenJets = NULL; | |
757 | ||
758 | if(fAODOut&&!aodGenJets){ | |
759 | aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMC.Data())); | |
760 | } | |
761 | if(fAODExtension&&!aodGenJets){ | |
762 | aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMC.Data())); | |
763 | } | |
764 | if(fAODIn&&!aodGenJets){ | |
765 | aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMC.Data())); | |
766 | } | |
767 | ||
768 | if(!aodGenJets){ | |
769 | Printf("%s:%d no generated Jet array with name %s in AOD", | |
770 | (char*)__FILE__,__LINE__, fJetBranchNameMC.Data()); | |
771 | PostData(1, fOutputList); | |
772 | return; | |
773 | } | |
774 | ||
775 | fListJetsGen->Clear(); | |
776 | ||
777 | //serarch for charged trigger at the MC generator level | |
8e103a56 | 778 | Int_t indexTriggGen = -1; |
779 | Double_t ptTriggGen = -1; | |
780 | Int_t iCounterGen = 0; | |
781 | Int_t triggersMC[200];//list of trigger candidates | |
782 | Int_t ntriggersMC = 0; //index in triggers array | |
783 | ||
80ac66f6 | 784 | if(fESD){//ESD input |
785 | ||
786 | AliMCEvent* mcEvent = MCEvent(); | |
787 | if(!mcEvent){ | |
788 | PostData(1, fOutputList); | |
789 | return; | |
ea603b64 | 790 | } |
791 | ||
792 | AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent); | |
793 | if(pythiaGenHeader){ | |
794 | nTrials = pythiaGenHeader->Trials(); | |
795 | ptHard = pythiaGenHeader->GetPtHard(); | |
796 | fh1PtHard->Fill(ptHard,eventW); | |
797 | fh1PtHardNoW->Fill(ptHard,1); | |
798 | fh1PtHardTrials->Fill(ptHard,nTrials); | |
799 | } | |
800 | ||
80ac66f6 | 801 | for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){ |
802 | if(!mcEvent->IsPhysicalPrimary(it)) continue; | |
803 | AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it); | |
804 | if(!part) continue; | |
805 | if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){ | |
8e103a56 | 806 | |
807 | if(fHardest==0 && ntriggersMC<200){//single inclusive trigger | |
808 | if(indexTriggGen > -1){//trigger candidater was found | |
809 | triggersMC[ntriggersMC] = indexTriggGen; | |
810 | ntriggersMC++; | |
811 | } | |
812 | } | |
813 | ||
80ac66f6 | 814 | iCounterGen++; |
815 | } | |
816 | } | |
817 | }else{ //AOD input | |
818 | if(!fAODIn){ | |
819 | PostData(1, fOutputList); | |
820 | return; | |
821 | } | |
822 | TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName())); | |
823 | if(!tca){ | |
824 | PostData(1, fOutputList); | |
825 | return; | |
826 | } | |
827 | for(Int_t it = 0; it < tca->GetEntriesFast(); it++){ | |
828 | AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it)); | |
829 | if(!part) continue; | |
830 | if(!part->IsPhysicalPrimary()) continue; | |
831 | if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){ | |
832 | ||
8e103a56 | 833 | if(fHardest==0 && ntriggersMC<200){//single inclusive trigger |
834 | if(indexTriggGen > -1){ //trigger candidater was found | |
835 | triggersMC[ntriggersMC] = indexTriggGen; | |
836 | ntriggersMC++; | |
837 | } | |
838 | } | |
839 | ||
80ac66f6 | 840 | iCounterGen++; |
841 | } | |
842 | } | |
843 | } | |
8e103a56 | 844 | |
845 | if(fHardest==0){ //single inclusive trigger | |
846 | if(ntriggersMC>0){ //there is at least one trigger | |
847 | Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1 | |
848 | indexTriggGen = triggersMC[rnd]; | |
80ac66f6 | 849 | }else{ |
8e103a56 | 850 | indexTriggGen = -1; //trigger not found |
80ac66f6 | 851 | } |
852 | } | |
8e103a56 | 853 | |
854 | //================== Fill jet list =================== | |
80ac66f6 | 855 | if(aodGenJets){ |
856 | if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMC.Data(), aodGenJets->GetEntriesFast()); | |
857 | ||
858 | for(Int_t igJet = 0; igJet < aodGenJets->GetEntriesFast(); igJet++) { | |
859 | AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJets)[igJet]); | |
860 | if(jetGen) fListJetsGen->Add(jetGen); | |
861 | } | |
8e103a56 | 862 | } |
863 | ||
864 | //============ Generator trigger+jet ================== | |
865 | Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0; | |
866 | Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries(); | |
867 | ||
868 | for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger | |
869 | indexTriggGen = igen; //trigger hadron | |
870 | ||
871 | if(indexTriggGen == -1) continue; | |
872 | AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen); | |
873 | if(!triggerGen) continue; | |
874 | ||
875 | if(fHardest >= 2){ | |
6168afc3 | 876 | if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>6 |
8e103a56 | 877 | } |
878 | if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue; | |
879 | ||
880 | ptTriggGen = triggerGen->Pt(); //count triggers | |
881 | fh2NtriggersGen->Fill(centValue, ptTriggGen); | |
882 | ||
80ac66f6 | 883 | //Count jets and trigger-jet pairs at MC generator level |
8e103a56 | 884 | if(!aodGenJets) continue; |
80ac66f6 | 885 | for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){ |
886 | AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij)); | |
887 | if(!jet) continue; | |
888 | Double_t etaJetGen = jet->Eta(); | |
80ac66f6 | 889 | |
890 | if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){ | |
80ac66f6 | 891 | |
8e103a56 | 892 | Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi()); |
893 | if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue; | |
80ac66f6 | 894 | |
8e103a56 | 895 | //Centrality, A, pT, pTtrigg |
896 | Double_t fillspecgen[] = { centValue, | |
80ac66f6 | 897 | jet->EffectiveAreaCharged(), |
898 | jet->Pt(), | |
6168afc3 | 899 | ptTriggGen, |
900 | TMath::Abs((Double_t) dphi) | |
80ac66f6 | 901 | }; |
8e103a56 | 902 | fHJetSpecGen->Fill(fillspecgen); |
903 | }//back to back jet-trigger pair | |
904 | }//jet loop | |
905 | }//trigger loop | |
906 | ||
907 | ||
908 | //================ RESPONSE MATRIX =============== | |
909 | if(aodGenJets){ | |
910 | //Count jets and trigger-jet pairs at MC generator level | |
911 | for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){ | |
912 | AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij)); | |
913 | if(!jet) continue; | |
914 | Double_t etaJetGen = jet->Eta(); | |
915 | Double_t ptJetGen = jet->Pt(); | |
916 | ||
917 | if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){ | |
918 | fhJetPtGen->Fill(ptJetGen); // generator level pt spectum of jets response mx normalization | |
80ac66f6 | 919 | } |
920 | } | |
80ac66f6 | 921 | if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets |
922 | Int_t ng = (Int_t) fListJetsGen->GetEntries(); | |
923 | Int_t nr = (Int_t) fListJets->GetEntries(); | |
924 | ||
925 | //Find closest MC generator - reconstructed jets | |
926 | if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet | |
927 | if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet | |
928 | ||
929 | if(fDebug){ | |
930 | Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize()); | |
931 | Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize()); | |
932 | } | |
933 | //matching of MC genrator level and reconstructed jets | |
934 | AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug); | |
935 | ||
936 | // Fill response matrix | |
937 | for(Int_t ir = 0; ir < nr; ir++){ | |
938 | AliAODJet *recJet = (AliAODJet*) fListJets->At(ir); | |
939 | Double_t etaJetRec = recJet->Eta(); | |
940 | Double_t ptJetRec = recJet->Pt(); | |
941 | //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc | |
942 | ||
943 | if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){ | |
944 | Int_t ig = faGenIndex[ir]; //associated generator level jet | |
945 | if(ig >= 0 && ig < ng){ | |
946 | if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir); | |
947 | AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig); | |
948 | Double_t ptJetGen = genJet->Pt(); | |
949 | Double_t etaJetGen = genJet->Eta(); | |
950 | ||
951 | //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc | |
952 | if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){ | |
953 | fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen); | |
954 | } | |
955 | }//ig>=0 | |
956 | }//rec jet in eta acceptance | |
957 | }//loop over reconstructed jets | |
958 | }// # of rec jets >0 | |
959 | }//pointer MC generator jets | |
960 | } //analyze generator level MC | |
961 | ||
8e103a56 | 962 | //============= RECONSTRUCTED INCLUSIVE JETS =============== |
ad869500 | 963 | |
ad869500 | 964 | Double_t etaJet = 0.0; |
965 | Double_t pTJet = 0.0; | |
966 | Double_t areaJet = 0.0; | |
967 | Double_t phiJet = 0.0; | |
7fc3d134 | 968 | Int_t indexLeadingJet = -1; |
969 | Double_t pTLeadingJet = -10.0; | |
7fc3d134 | 970 | Double_t areaLeadingJet = -10.0; |
8e103a56 | 971 | |
ad869500 | 972 | for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){ |
973 | AliAODJet* jet = (AliAODJet*)(fListJets->At(ij)); | |
974 | if(!jet) continue; | |
975 | etaJet = jet->Eta(); | |
976 | phiJet = jet->Phi(); | |
977 | pTJet = jet->Pt(); | |
978 | if(pTJet==0) continue; | |
979 | ||
980 | if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue; | |
981 | areaJet = jet->EffectiveAreaCharged(); | |
7fc3d134 | 982 | |
ad869500 | 983 | //Jet Diagnostics--------------------------------- |
8e103a56 | 984 | fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi |
985 | fhJetEta->Fill(pTJet, etaJet); | |
7fc3d134 | 986 | //search for leading jet |
987 | if(pTJet > pTLeadingJet){ | |
988 | indexLeadingJet = ij; | |
989 | pTLeadingJet = pTJet; | |
7fc3d134 | 990 | areaLeadingJet = areaJet; |
991 | } | |
992 | ||
993 | // raw spectra of INCLUSIVE jets | |
80ac66f6 | 994 | //Centrality, pTjet, A |
7fc3d134 | 995 | Double_t fillraw[] = { centValue, |
996 | pTJet, | |
7fc3d134 | 997 | areaJet |
998 | }; | |
999 | fHJetPtRaw->Fill(fillraw); | |
8e103a56 | 1000 | } |
7fc3d134 | 1001 | |
8e103a56 | 1002 | if(indexLeadingJet > -1){ |
1003 | // raw spectra of LEADING jets | |
1004 | //Centrality, pTjet, A | |
1005 | Double_t fillleading[] = { centValue, | |
1006 | pTLeadingJet, | |
1007 | areaLeadingJet | |
1008 | }; | |
1009 | fHLeadingJetPtRaw->Fill(fillleading); | |
1010 | } | |
1011 | ||
1012 | ||
1013 | // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================ | |
1014 | //Find Hadron trigger | |
1015 | TList particleList; //list of tracks | |
1016 | Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list | |
1017 | if(fIsMC) FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos | |
7fc3d134 | 1018 | |
8e103a56 | 1019 | //set ranges of the trigger loop |
1020 | Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0; | |
1021 | Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries(); | |
7fc3d134 | 1022 | |
8e103a56 | 1023 | for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger |
1024 | indexTrigg = itrk; //trigger hadron with pT >10 GeV | |
7fc3d134 | 1025 | |
8e103a56 | 1026 | if(indexTrigg < 0) continue; |
1027 | ||
1028 | AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg); | |
1029 | if(!triggerHadron){ | |
1030 | PostData(1, fOutputList); | |
1031 | return; | |
1032 | } | |
1033 | if(fHardest >= 2){ | |
7d9fd386 | 1034 | if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6 |
8e103a56 | 1035 | } |
1036 | if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue; | |
7fc3d134 | 1037 | |
8e103a56 | 1038 | //Fill trigger histograms |
1039 | fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT | |
1040 | fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi | |
1041 | fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta()); | |
7fc3d134 | 1042 | |
1043 | ||
8e103a56 | 1044 | //---------- make trigger-jet pairs --------- |
1045 | Int_t injet4 = 0; | |
1046 | Int_t injet = 0; | |
1047 | ||
1048 | for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){ | |
1049 | AliAODJet* jet = (AliAODJet*)(fListJets->At(ij)); | |
1050 | if(!jet) continue; | |
1051 | etaJet = jet->Eta(); | |
1052 | phiJet = jet->Phi(); | |
1053 | pTJet = jet->Pt(); | |
1054 | if(pTJet==0) continue; | |
1055 | ||
1056 | if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue; | |
1057 | areaJet = jet->EffectiveAreaCharged(); | |
1058 | if(areaJet >= 0.07) injet++; | |
1059 | if(areaJet >= 0.4) injet4++; | |
1060 | ||
1061 | Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet); | |
1062 | fhDphiTriggerJet->Fill(dphi); //Input | |
1063 | ||
1064 | //Dphi versus jet pT | |
1065 | //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg | |
1066 | Double_t filldp[] = { centValue, | |
1067 | dphi, | |
1068 | pTJet, | |
1069 | triggerHadron->Pt() | |
1070 | }; | |
1071 | fHDphiVsJetPtAll->Fill(filldp); | |
1072 | ||
1073 | // Select back to back trigger - jet pairs | |
1074 | if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue; | |
1075 | fhDphiTriggerJetAccept->Fill(dphi); //Accepted | |
1076 | ||
1077 | ||
1078 | //Centrality, A, pTjet, pTtrigg | |
1079 | Double_t fillspec[] = { centValue, | |
1080 | areaJet, | |
1081 | pTJet, | |
6168afc3 | 1082 | triggerHadron->Pt(), |
1083 | TMath::Abs((Double_t) dphi) | |
8e103a56 | 1084 | }; |
1085 | fHJetSpec->Fill(fillspec); | |
1086 | }//jet loop | |
1087 | ||
1088 | //Fill Jet Density In the Event A>0.07 | |
1089 | if(injet>0){ | |
1090 | Double_t filldens[]={ (Double_t) particleList.GetEntries(), | |
ad869500 | 1091 | injet/fkAcceptance, |
1092 | triggerHadron->Pt() | |
1093 | }; | |
8e103a56 | 1094 | fHJetDensity->Fill(filldens); |
1095 | } | |
ad869500 | 1096 | |
8e103a56 | 1097 | //Fill Jet Density In the Event A>0.4 |
1098 | if(injet4>0){ | |
1099 | Double_t filldens4[]={ (Double_t) particleList.GetEntries(), | |
1100 | injet4/fkAcceptance, | |
1101 | triggerHadron->Pt() | |
1102 | }; | |
1103 | fHJetDensityA4->Fill(filldens4); | |
1104 | } | |
ad869500 | 1105 | } |
1106 | ||
8e103a56 | 1107 | |
ad869500 | 1108 | PostData(1, fOutputList); |
1109 | } | |
1110 | ||
1111 | //---------------------------------------------------------------------------- | |
1112 | void AliAnalysisTaskJetCorePP::Terminate(const Option_t *) | |
1113 | { | |
1114 | // Draw result to the screen | |
1115 | // Called once at the end of the query | |
1116 | ||
1117 | if(fDebug) printf("AliAnalysisTaskJetCorePP DONE\n"); | |
1118 | if(!GetOutputData(1)) return; | |
1119 | } | |
1120 | ||
1121 | ||
1122 | //---------------------------------------------------------------------------- | |
1123 | Int_t AliAnalysisTaskJetCorePP::GetListOfTracks(TList *list){ | |
1124 | //Fill the list of accepted tracks (passed track cut) | |
1125 | //return consecutive index of the hardest ch hadron in the list | |
1126 | Int_t iCount = 0; | |
1127 | AliAODEvent *aodevt = NULL; | |
1128 | ||
1129 | if(!fESD) aodevt = fAODIn; | |
1130 | else aodevt = fAODOut; | |
1131 | ||
1132 | if(!aodevt) return -1; | |
1133 | ||
1134 | Int_t index = -1; //index of the highest particle in the list | |
1135 | Double_t ptmax = -10; | |
8e103a56 | 1136 | Int_t triggers[200]; |
1137 | Int_t ntriggers = 0; //index in triggers array | |
ad869500 | 1138 | |
80ac66f6 | 1139 | for(Int_t it = 0; it < aodevt->GetNumberOfTracks(); it++){ |
ad869500 | 1140 | AliAODTrack *tr = aodevt->GetTrack(it); |
1141 | ||
1142 | //if((fFilterMask > 0) && !(tr->TestFilterBit(fFilterMask))) continue; | |
1143 | if((fFilterMask > 0) && !(tr->IsHybridGlobalConstrainedGlobal())) continue; | |
1144 | if(TMath::Abs((Float_t) tr->Eta()) > fTrackEtaCut) continue; | |
1145 | if(tr->Pt() < fTrackLowPtCut) continue; | |
1146 | list->Add(tr); | |
8e103a56 | 1147 | |
1148 | //Search trigger: | |
1149 | if(fHardest>0){ //leading particle | |
1150 | if(tr->Pt()>ptmax){ | |
1151 | ptmax = tr->Pt(); | |
1152 | index = iCount; | |
1153 | } | |
1154 | } | |
1155 | ||
1156 | if(fHardest==0 && ntriggers<200){ //single inclusive | |
1157 | if(fTriggerPtRangeLow <= tr->Pt() && | |
1158 | tr->Pt() < fTriggerPtRangeHigh){ | |
1159 | triggers[ntriggers] = iCount; | |
1160 | ntriggers++; | |
1161 | } | |
ad869500 | 1162 | } |
8e103a56 | 1163 | |
ad869500 | 1164 | iCount++; |
1165 | } | |
1166 | ||
8e103a56 | 1167 | if(fHardest==0 && ntriggers>0){ //select random inclusive trigger |
1168 | Int_t rnd = fRandom->Integer(ntriggers); //0 to ntriggers-1 | |
1169 | index = triggers[rnd]; | |
ad869500 | 1170 | } |
8e103a56 | 1171 | |
1172 | return index; | |
ad869500 | 1173 | } |
1174 | ||
1175 | //---------------------------------------------------------------------------- | |
80ac66f6 | 1176 | /* |
ad869500 | 1177 | Double_t AliAnalysisTaskJetCorePP::GetBackgroundInPerpCone(Float_t jetR, Double_t jetPhi, Double_t jetEta, TList* trkList){ |
1178 | //calculate sum of track pT in the cone perpendicular in phi to the jet | |
1179 | //jetR = cone radius | |
1180 | // jetPhi, jetEta = direction of the jet | |
1181 | Int_t numberOfTrks = trkList->GetEntries(); | |
1182 | Double_t pTsum = 0.0; | |
1183 | Double_t perpConePhi = jetPhi + TMath::Pi()/2;//perp cone w.r.t. jet in phi | |
1184 | for(Int_t it=0; it<numberOfTrks; it++){ | |
1185 | AliVParticle *trk = (AliVParticle*) trkList->At(it); | |
1186 | Double_t dphi = RelativePhi(perpConePhi,trk->Phi()); | |
1187 | Double_t deta = trk->Eta()-jetEta; | |
1188 | ||
1189 | if( (dphi*dphi + deta*deta)< (jetR*jetR)){ | |
1190 | pTsum += trk->Pt(); | |
1191 | } | |
1192 | } | |
1193 | ||
1194 | return pTsum; | |
1195 | } | |
80ac66f6 | 1196 | */ |
ad869500 | 1197 | //---------------------------------------------------------------------------- |
1198 | ||
1199 | Double_t AliAnalysisTaskJetCorePP::RelativePhi(Double_t mphi,Double_t vphi){ | |
1200 | //Get relative azimuthal angle of two particles -pi to pi | |
1201 | if (vphi < -TMath::Pi()) vphi += TMath::TwoPi(); | |
1202 | else if (vphi > TMath::Pi()) vphi -= TMath::TwoPi(); | |
1203 | ||
1204 | if (mphi < -TMath::Pi()) mphi += TMath::TwoPi(); | |
1205 | else if (mphi > TMath::Pi()) mphi -= TMath::TwoPi(); | |
1206 | ||
1207 | Double_t dphi = mphi - vphi; | |
1208 | if (dphi < -TMath::Pi()) dphi += TMath::TwoPi(); | |
1209 | else if (dphi > TMath::Pi()) dphi -= TMath::TwoPi(); | |
1210 | ||
1211 | return dphi;//dphi in [-Pi, Pi] | |
1212 | } | |
1213 | ||
1214 | ||
80ac66f6 | 1215 | //---------------------------------------------------------------------------- |
1216 | Bool_t AliAnalysisTaskJetCorePP::SelectMCGenTracks(AliVParticle *trk, TList *trkList, Double_t &ptLeading, Int_t &index, Int_t counter){ | |
1217 | //fill track efficiency denominator | |
1218 | if(trk->Pt() < fTrackLowPtCut) return kFALSE; | |
1219 | if(trk->Charge()==0) return kFALSE; | |
1220 | if(TMath::Abs((Float_t) trk->Eta()) > fTrackEtaCut) return kFALSE; | |
1221 | trkList->Add(trk); | |
1222 | fhPtTrkTruePrimGen->Fill(trk->Pt(),trk->Eta()); //Efficiency denominator | |
8e103a56 | 1223 | |
1224 | //Search trigger: | |
1225 | if(fHardest>0){ //leading particle | |
1226 | if(ptLeading < trk->Pt()){ | |
1227 | index = counter; | |
1228 | ptLeading = trk->Pt(); | |
1229 | } | |
80ac66f6 | 1230 | } |
8e103a56 | 1231 | |
1232 | if(fHardest==0){ //single inclusive | |
1233 | index = -1; | |
1234 | if(fTriggerPtRangeLow <= trk->Pt() && | |
1235 | trk->Pt() < fTriggerPtRangeHigh){ | |
1236 | index = counter; | |
1237 | } | |
1238 | } | |
1239 | ||
80ac66f6 | 1240 | return kTRUE; |
1241 | } | |
1242 | ||
1243 | //---------------------------------------------------------------------------- | |
1244 | void AliAnalysisTaskJetCorePP::FillEffHistos(TList *recList, TList *genList){ | |
1245 | //fill track efficiency numerator | |
1246 | if(!recList) return; | |
1247 | if(!genList) return; | |
1248 | Int_t nRec = recList->GetEntries(); | |
1249 | Int_t nGen = genList->GetEntries(); | |
1250 | for(Int_t ir=0; ir<nRec; ir++){ | |
1251 | AliAODTrack *trkRec = (AliAODTrack*) recList->At(ir); | |
1252 | if(!trkRec) continue; | |
1253 | Int_t recLabel = TMath::Abs(trkRec->GetLabel()); | |
1254 | Bool_t matched = kFALSE; | |
1255 | ||
1256 | for(Int_t ig=0; ig<nGen; ig++){ | |
1257 | AliVParticle *trkGen = (AliVParticle*) genList->At(ig); | |
1258 | if(!trkGen) continue; | |
1259 | Int_t genLabel = TMath::Abs(trkGen->GetLabel()); | |
1260 | if(recLabel==genLabel){ | |
1261 | fhPtTrkTruePrimRec->Fill(trkGen->Pt(),trkGen->Eta()); | |
1262 | matched = kTRUE; | |
1263 | break; | |
1264 | } | |
1265 | } | |
1266 | if(!matched) fhPtTrkSecOrFakeRec->Fill(trkRec->Pt(),trkRec->Eta()); | |
1267 | } | |
1268 | ||
1269 | return; | |
1270 | } | |
1271 | ||
1272 |