]>
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), |
ea603b64 | 130 | fkDeltaPhiCut(TMath::Pi()-0.6), |
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), |
ea603b64 | 207 | fkDeltaPhiCut(TMath::Pi()-0.6), |
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 | } | |
385 | ||
386 | } | |
387 | ||
388 | //-------------------------------------------------------------- | |
389 | ||
390 | void AliAnalysisTaskJetCorePP::UserCreateOutputObjects() | |
391 | { | |
392 | ||
ad869500 | 393 | // Create histograms |
394 | // Called once | |
80ac66f6 | 395 | fListJets = new TList(); //reconstructed level |
396 | ||
397 | fIsMC = (fJetBranchNameMC.Length()>0) ? kTRUE : kFALSE; | |
398 | ||
8e103a56 | 399 | fRandom = new TRandom3(0); |
400 | ||
80ac66f6 | 401 | if(fIsMC) fListJetsGen = new TList(); //generator level |
ad869500 | 402 | OpenFile(1); |
403 | if(!fOutputList) fOutputList = new TList(); | |
404 | fOutputList->SetOwner(kTRUE); | |
405 | ||
406 | Bool_t oldStatus = TH1::AddDirectoryStatus(); | |
407 | TH1::AddDirectory(kFALSE); | |
408 | ||
409 | fHistEvtSelection = new TH1I("fHistEvtSelection", "event selection", 6, -0.5, 5.5); | |
410 | fHistEvtSelection->GetXaxis()->SetBinLabel(1,"ACCEPTED"); | |
411 | fHistEvtSelection->GetXaxis()->SetBinLabel(2,"events IN"); | |
412 | fHistEvtSelection->GetXaxis()->SetBinLabel(3,"event selection (rejected)"); | |
413 | fHistEvtSelection->GetXaxis()->SetBinLabel(4,"vertex cut (rejected)"); | |
414 | fHistEvtSelection->GetXaxis()->SetBinLabel(5,"centrality (rejected)"); | |
8e103a56 | 415 | fHistEvtSelection->GetXaxis()->SetBinLabel(6,"event number (rejected)"); |
ad869500 | 416 | |
417 | fOutputList->Add(fHistEvtSelection); | |
418 | ||
419 | Int_t nBinsCentrality = (fSystem==0) ? 1 : 10; // pp=1 else 10 | |
80ac66f6 | 420 | //trigger pt spectrum (reconstructed) |
ad869500 | 421 | fh2Ntriggers = new TH2F("fh2Ntriggers","# of triggers", |
422 | nBinsCentrality,0.0,100.0,50,0.0,50.0); | |
7fc3d134 | 423 | fOutputList->Add(fh2Ntriggers); |
ad869500 | 424 | |
80ac66f6 | 425 | //Centrality, A, pTjet, pTtrigg |
426 | const Int_t dimSpec = 4; | |
427 | const Int_t nBinsSpec[dimSpec] = {nBinsCentrality, 100, 100, 50}; | |
428 | const Double_t lowBinSpec[dimSpec] = {0.0, 0.0, 0, 0.0}; | |
429 | const Double_t hiBinSpec[dimSpec] = {100.0, 1.0, 200.0,50.0}; | |
7fc3d134 | 430 | fHJetSpec = new THnSparseF("fHJetSpec", |
80ac66f6 | 431 | "Recoil jet spectrum [cent,A,pTjet,pTtrig]", |
7fc3d134 | 432 | dimSpec,nBinsSpec,lowBinSpec,hiBinSpec); |
ad869500 | 433 | fOutputList->Add(fHJetSpec); |
434 | ||
435 | //------------------- HISTOS FOR DIAGNOSTIC ---------------------- | |
436 | //Jet number density histos [Trk Mult, jet density, pT trigger] | |
437 | const Int_t dimJetDens = 3; | |
80ac66f6 | 438 | const Int_t nBinsJetDens[dimJetDens] = {100, 100, 10}; |
ad869500 | 439 | const Double_t lowBinJetDens[dimJetDens] = {0.0, 0.0, 0.0}; |
80ac66f6 | 440 | const Double_t hiBinJetDens[dimJetDens] = {500.0, 5.0, 50.0}; |
ad869500 | 441 | |
442 | fHJetDensity = new THnSparseF("fHJetDensity","Jet dens vs trk mult A>0.07", | |
443 | dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens); | |
444 | ||
445 | fHJetDensityA4 =new THnSparseF("fHJetDensityA4","Jet dens vs trk mult A>0.4", | |
446 | dimJetDens,nBinsJetDens,lowBinJetDens,hiBinJetDens); | |
447 | ||
ad869500 | 448 | fOutputList->Add(fHJetDensity); |
449 | fOutputList->Add(fHJetDensityA4); | |
450 | ||
451 | ||
452 | //inclusive azimuthal and pseudorapidity histograms | |
7fc3d134 | 453 | fhJetPhi = new TH2D("fhJetPhi","Azim dist jets vs pTjet", |
454 | 50, 0, 100, 50,-TMath::Pi(),TMath::Pi()); | |
455 | fhTriggerPhi= new TH2D("fhTriggerPhi","azim dist trig had vs pTtrigg", | |
456 | 25, 0, 50, 50,-TMath::Pi(),TMath::Pi()); | |
457 | fhJetEta = new TH2D("fhJetEta","Eta dist jets vs pTjet", | |
458 | 50,0, 100, 40,-0.9,0.9); | |
459 | fhTriggerEta = new TH2D("fhTriggerEta","Eta dist trig had vs pTtrigg", | |
460 | 25, 0, 50, 40,-0.9,0.9); | |
461 | ||
ad869500 | 462 | fhVertexZ = new TH1D("fhVertexZ","z vertex",40,-20,20); |
463 | fhVertexZAccept = new TH1D("fhVertexZAccept","z vertex after cut",40,-20,20); | |
464 | fhContribVtx = new TH1D("fhContribVtx","contrib to vtx",200,0,200); | |
465 | fhContribVtxAccept = new TH1D("fhContribVtxAccept","contrib to vtx after cut",200,0,200); | |
466 | fhDphiTriggerJet = new TH1D("fhDphiTriggerJet","Deltaphi trig-jet",50, -TMath::Pi(),TMath::Pi()); | |
467 | fhDphiTriggerJetAccept = new TH1D("fhDphiTriggerJetAccept","Deltaphi trig-jet after cut",50, -TMath::Pi(),TMath::Pi()); | |
468 | fhCentrality = new TH1D("fhCentrality","Centrality",20,0,100); | |
469 | fhCentralityAccept = new TH1D("fhCentralityAccept","Centrality after cut",20,0,100); | |
470 | ||
471 | fOutputList->Add(fhJetPhi); | |
472 | fOutputList->Add(fhTriggerPhi); | |
473 | fOutputList->Add(fhJetEta); | |
474 | fOutputList->Add(fhTriggerEta); | |
475 | fOutputList->Add(fhVertexZ); | |
476 | fOutputList->Add(fhVertexZAccept); | |
477 | fOutputList->Add(fhContribVtx); | |
478 | fOutputList->Add(fhContribVtxAccept); | |
479 | fOutputList->Add(fhDphiTriggerJet); | |
480 | fOutputList->Add(fhDphiTriggerJetAccept); | |
481 | fOutputList->Add(fhCentrality); | |
482 | fOutputList->Add(fhCentralityAccept); | |
483 | ||
7fc3d134 | 484 | // raw spectra of INCLUSIVE jets |
80ac66f6 | 485 | //Centrality, pTjet, A |
486 | const Int_t dimRaw = 3; | |
487 | const Int_t nBinsRaw[dimRaw] = {nBinsCentrality, 50, 100}; | |
488 | const Double_t lowBinRaw[dimRaw] = {0.0, 0.0, 0.0}; | |
489 | const Double_t hiBinRaw[dimRaw] = {100.0, 100, 1.0}; | |
7fc3d134 | 490 | fHJetPtRaw = new THnSparseF("fHJetPtRaw", |
80ac66f6 | 491 | "Incl. jet spectrum [cent,pTjet,A]", |
7fc3d134 | 492 | dimRaw,nBinsRaw,lowBinRaw,hiBinRaw); |
493 | fOutputList->Add(fHJetPtRaw); | |
494 | ||
495 | // raw spectra of LEADING jets | |
80ac66f6 | 496 | //Centrality, pTjet, A |
7fc3d134 | 497 | fHLeadingJetPtRaw = new THnSparseF("fHLeadingJetPtRaw", |
80ac66f6 | 498 | "Leading jet spectrum [cent,pTjet,A]", |
7fc3d134 | 499 | dimRaw,nBinsRaw,lowBinRaw,hiBinRaw); |
500 | fOutputList->Add(fHLeadingJetPtRaw); | |
501 | ||
502 | // Dphi versus pT jet | |
503 | //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg | |
504 | const Int_t dimDp = 4; | |
505 | const Int_t nBinsDp[dimDp] = {nBinsCentrality, 50, 50, 50}; | |
506 | const Double_t lowBinDp[dimDp] = {0.0, -TMath::Pi(), 0.0, 0.0}; | |
507 | const Double_t hiBinDp[dimDp] = {100.0, TMath::Pi(), 100.0, 100.0}; | |
508 | fHDphiVsJetPtAll = new THnSparseF("fHDphiVsJetPtAll", | |
509 | "Dphi vs jet pT [cent,Dphi,pTjet,pTtrigg]", | |
510 | dimDp,nBinsDp,lowBinDp,hiBinDp); | |
511 | fOutputList->Add(fHDphiVsJetPtAll); | |
512 | ||
7fc3d134 | 513 | |
80ac66f6 | 514 | //analyze MC generator level |
515 | if(fIsMC){ | |
516 | fhJetPtGenVsJetPtRec = new TH2D("fhJetPtGenVsJetPtRec","JetPtGenVsJetPtRec", 100,0,200, 100,0,200); | |
517 | fOutputList->Add(fhJetPtGenVsJetPtRec); //gen MC jet pt spectrum versus reconstructed pt spectrum | |
518 | ||
519 | fhJetPtGen = new TH1D("fhJetPtGen","Jet Pt (MC Gen)",100,0,200); //MC generator jet pt spectrum | |
520 | fOutputList->Add(fhJetPtGen); | |
7fc3d134 | 521 | |
80ac66f6 | 522 | fh2NtriggersGen = (TH2F*) fh2Ntriggers->Clone("fh2NtriggersGen"); |
523 | fh2NtriggersGen->SetTitle(Form("%s Gen MC",fh2Ntriggers->GetTitle())); | |
524 | fOutputList->Add(fh2NtriggersGen); | |
7fc3d134 | 525 | |
80ac66f6 | 526 | //Centrality, A, pT, pTtrigg |
527 | fHJetSpecGen = (THnSparseF*) fHJetSpec->Clone("fHJetSpecGen"); | |
528 | fHJetSpecGen->SetTitle(Form("%s Gen MC",fHJetSpec->GetTitle())); | |
529 | fOutputList->Add(fHJetSpecGen); | |
530 | ||
531 | //track efficiency/contamination histograms eta versus pT | |
532 | fhPtTrkTruePrimRec = new TH2D("fhPtTrkTruePrimRec","PtTrkTruePrimRec",100,0,20,18,-0.9,0.9); | |
533 | fOutputList->Add(fhPtTrkTruePrimRec); | |
534 | ||
535 | fhPtTrkTruePrimGen = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkTruePrimGen"); | |
536 | fhPtTrkTruePrimGen->SetTitle("PtTrkTruePrimGen"); | |
537 | fOutputList->Add(fhPtTrkTruePrimGen); | |
538 | ||
539 | fhPtTrkSecOrFakeRec = (TH2D*) fhPtTrkTruePrimRec->Clone("fhPtTrkSecOrFakeRec"); | |
540 | fhPtTrkSecOrFakeRec->SetTitle("PtTrkSecOrFakeRec"); | |
541 | fOutputList->Add(fhPtTrkSecOrFakeRec); | |
542 | } | |
ea603b64 | 543 | //------------------------------------- |
544 | // pythia histograms | |
545 | const Int_t nBinPt = 150; | |
546 | Double_t binLimitsPt[nBinPt+1]; | |
547 | for(Int_t iPt = 0;iPt <= nBinPt;iPt++){ | |
548 | if(iPt == 0){ | |
549 | binLimitsPt[iPt] = -50.0; | |
550 | }else{// 1.0 | |
551 | binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.; | |
552 | } | |
553 | } | |
554 | ||
555 | fh1Xsec = new TProfile("fh1Xsec","xsec from pyxsec.root",1,0,1); | |
556 | fh1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>"); | |
557 | fOutputList->Add(fh1Xsec); | |
558 | fh1Trials = new TH1F("fh1Trials","trials root file",1,0,1); | |
559 | fh1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}"); | |
560 | fOutputList->Add(fh1Trials); | |
561 | fh1AvgTrials = new TH1F("fh1AvgTrials","trials root file",1,0,1); | |
562 | fh1AvgTrials->GetXaxis()->SetBinLabel(1,"#sum{avg ntrials}"); | |
563 | fOutputList->Add(fh1AvgTrials); | |
564 | fh1PtHard = new TH1F("fh1PtHard","PYTHIA Pt hard;p_{T,hard}",nBinPt,binLimitsPt); | |
565 | fOutputList->Add(fh1PtHard); | |
566 | fh1PtHardNoW = new TH1F("fh1PtHardNoW","PYTHIA Pt hard no weight;p_{T,hard}",nBinPt,binLimitsPt); | |
567 | fOutputList->Add(fh1PtHardNoW); | |
568 | fh1PtHardTrials = new TH1F("fh1PtHardTrials","PYTHIA Pt hard weight with trials;p_{T,hard}",nBinPt,binLimitsPt); | |
569 | fOutputList->Add(fh1PtHardTrials); | |
570 | ||
ad869500 | 571 | |
572 | // =========== Switch on Sumw2 for all histos =========== | |
573 | for(Int_t i=0; i<fOutputList->GetEntries(); i++){ | |
574 | TH1 *h1 = dynamic_cast<TH1*>(fOutputList->At(i)); | |
575 | if(h1){ | |
576 | h1->Sumw2(); | |
577 | continue; | |
578 | } | |
579 | THnSparse *hn = dynamic_cast<THnSparse*>(fOutputList->At(i)); | |
580 | if(hn){ | |
581 | hn->Sumw2(); | |
582 | } | |
583 | } | |
584 | TH1::AddDirectory(oldStatus); | |
585 | ||
586 | PostData(1, fOutputList); | |
587 | } | |
588 | ||
589 | //-------------------------------------------------------------------- | |
590 | ||
591 | void AliAnalysisTaskJetCorePP::UserExec(Option_t *) | |
592 | { | |
ad869500 | 593 | //Event loop |
ea603b64 | 594 | Double_t eventW = 1.0; |
595 | Double_t ptHard = 0.0; | |
596 | Double_t nTrials = 1.0; // Trials for MC trigger | |
597 | if(fIsMC) fh1AvgTrials->Fill("#sum{avg ntrials}",fAvgTrials); | |
ad869500 | 598 | |
599 | if(TMath::Abs((Float_t) fJetParamR)<0.00001){ | |
80ac66f6 | 600 | AliError("Cone radius is set to zero."); |
ad869500 | 601 | return; |
602 | } | |
603 | if(!strlen(fJetBranchName.Data())){ | |
604 | AliError("Jet branch name not set."); | |
605 | return; | |
606 | } | |
607 | ||
608 | fESD = dynamic_cast<AliESDEvent*>(InputEvent()); | |
609 | if(!fESD){ | |
80ac66f6 | 610 | if(fDebug>1) AliError("ESD not available"); |
ad869500 | 611 | fAODIn = dynamic_cast<AliAODEvent*>(InputEvent()); |
612 | } | |
613 | ||
614 | fAODOut = dynamic_cast<AliAODEvent*>(AODEvent()); | |
615 | AliAODEvent* aod = NULL; | |
616 | // take all other information from the aod we take the tracks from | |
617 | if(!aod){ | |
618 | if(!fESD) aod = fAODIn; | |
619 | else aod = fAODOut; | |
620 | } | |
621 | ||
622 | if(fNonStdFile.Length()!=0){ | |
623 | // case that we have an AOD extension we can fetch the jets from the extended output | |
624 | AliAODHandler *aodH = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()); | |
625 | fAODExtension = aodH ? aodH->GetExtension(fNonStdFile.Data()) : 0; | |
626 | if(!fAODExtension){ | |
627 | if(fDebug>1) Printf("AODExtension found for %s",fNonStdFile.Data()); | |
628 | } | |
629 | } | |
630 | ||
80ac66f6 | 631 | // ----------------- EVENT SELECTION -------------------------- |
ad869500 | 632 | fHistEvtSelection->Fill(1); // number of events before event selection |
633 | ||
634 | // physics selection | |
635 | AliInputEventHandler* inputHandler = (AliInputEventHandler*) | |
636 | ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler()); | |
637 | ||
638 | if(!(inputHandler->IsEventSelected() & fOfflineTrgMask)){ | |
639 | if(fDebug) Printf(" Trigger Selection: event REJECTED ... "); | |
640 | fHistEvtSelection->Fill(2); | |
641 | PostData(1, fOutputList); | |
642 | return; | |
643 | } | |
644 | ||
645 | //check AOD pointer | |
646 | if(!aod){ | |
647 | if(fDebug) Printf("%s:%d No AOD",(char*)__FILE__,__LINE__); | |
648 | fHistEvtSelection->Fill(3); | |
649 | PostData(1, fOutputList); | |
650 | return; | |
651 | } | |
652 | ||
653 | // vertex selection | |
654 | AliAODVertex* primVtx = aod->GetPrimaryVertex(); | |
655 | ||
656 | if(!primVtx){ | |
657 | if(fDebug) Printf("%s:%d No primVtx",(char*)__FILE__,__LINE__); | |
658 | fHistEvtSelection->Fill(3); | |
659 | PostData(1, fOutputList); | |
660 | return; | |
661 | } | |
662 | ||
663 | Int_t nTracksPrim = primVtx->GetNContributors(); | |
664 | Float_t vtxz = primVtx->GetZ(); | |
665 | //Input events | |
666 | fhContribVtx->Fill(nTracksPrim); | |
667 | if( nTracksPrim > 0 ) fhVertexZ->Fill(vtxz); | |
668 | ||
669 | if((nTracksPrim < fMinContribVtx) || | |
670 | (vtxz < fVtxZMin) || | |
671 | (vtxz > fVtxZMax)){ | |
672 | if(fDebug) Printf("%s:%d primary vertex z = %f: event REJECTED...", | |
673 | (char*)__FILE__,__LINE__,vtxz); | |
674 | fHistEvtSelection->Fill(3); | |
675 | PostData(1, fOutputList); | |
676 | return; | |
677 | }else{ | |
678 | //Accepted events | |
679 | fhContribVtxAccept->Fill(nTracksPrim); | |
680 | fhVertexZAccept->Fill(vtxz); | |
681 | } | |
682 | //FK// No event class selection imposed | |
683 | // event class selection (from jet helper task) | |
684 | //Int_t eventClass = AliAnalysisHelperJetTasks::EventClass(); | |
685 | //if(fDebug) Printf("Event class %d", eventClass); | |
686 | //if(eventClass < fEvtClassMin || eventClass > fEvtClassMax){ | |
687 | // fHistEvtSelection->Fill(4); | |
688 | // PostData(1, fOutputList); | |
689 | // return; | |
690 | //} | |
691 | ||
80ac66f6 | 692 | //------------------ CENTRALITY SELECTION --------------- |
ad869500 | 693 | AliCentrality *cent = 0x0; |
694 | Double_t centValue = 0.0; | |
695 | if(fSystem){ //fSystem=0 for pp, fSystem=1 for pPb | |
696 | if(fESD){ | |
697 | cent = fESD->GetCentrality(); | |
698 | if(cent) centValue = cent->GetCentralityPercentile("V0M"); | |
699 | }else{ | |
700 | centValue = aod->GetHeader()->GetCentrality(); | |
701 | } | |
702 | if(fDebug) printf("centrality: %f\n", centValue); | |
703 | //Input events | |
704 | fhCentrality->Fill(centValue); | |
705 | ||
706 | if(centValue < fCentMin || centValue > fCentMax){ | |
707 | fHistEvtSelection->Fill(4); | |
708 | PostData(1, fOutputList); | |
709 | return; | |
710 | }else{ | |
711 | //Accepted events | |
712 | fhCentralityAccept->Fill( centValue ); | |
713 | } | |
714 | } | |
715 | ||
8e103a56 | 716 | //-----------------select disjunct event subsamples ---------------- |
717 | Int_t eventnum = aod->GetHeader()->GetEventNumberESDFile(); | |
718 | Int_t lastdigit = eventnum % 10; | |
719 | if(!(fEventNumberRangeLow<=lastdigit && lastdigit<=fEventNumberRangeHigh)){ | |
720 | fHistEvtSelection->Fill(5); | |
721 | PostData(1, fOutputList); | |
722 | return; | |
723 | } | |
724 | ||
ad869500 | 725 | if(fDebug) std::cout<<" ACCEPTED EVENT "<<endl; |
ad869500 | 726 | fHistEvtSelection->Fill(0); // accepted events |
ad869500 | 727 | // ------------------- end event selection -------------------- |
80ac66f6 | 728 | |
ad869500 | 729 | |
80ac66f6 | 730 | // fetch RECONSTRUCTED jets |
ad869500 | 731 | TClonesArray *aodJets = 0x0; |
732 | ||
733 | if(fAODOut && !aodJets){ | |
734 | aodJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchName.Data())); | |
735 | } | |
736 | if(fAODExtension && !aodJets){ | |
737 | aodJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchName.Data())); | |
738 | } | |
739 | if(fAODIn && !aodJets){ | |
740 | aodJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchName.Data())); | |
741 | } | |
742 | ||
80ac66f6 | 743 | //--------- Fill list of RECONSTRUCTED jets -------------- |
744 | fListJets->Clear(); | |
745 | if(aodJets){ | |
746 | if(fDebug) Printf("########## %s: %d jets",fJetBranchName.Data(), aodJets->GetEntriesFast()); | |
747 | for(Int_t iJet = 0; iJet < aodJets->GetEntriesFast(); iJet++) { | |
748 | AliAODJet *jet = dynamic_cast<AliAODJet*>((*aodJets)[iJet]); | |
749 | if (jet) fListJets->Add(jet); | |
750 | } | |
751 | } | |
752 | ||
753 | //--------- Fill list of MC GENERATOR LEVEL jets -------------- | |
754 | TList particleListGen; //list of tracks in MC | |
755 | ||
756 | if(fIsMC){ //analyze generator level MC | |
757 | // fetch MC generator level jets | |
758 | TClonesArray *aodGenJets = NULL; | |
759 | ||
760 | if(fAODOut&&!aodGenJets){ | |
761 | aodGenJets = dynamic_cast<TClonesArray*>(fAODOut->FindListObject(fJetBranchNameMC.Data())); | |
762 | } | |
763 | if(fAODExtension&&!aodGenJets){ | |
764 | aodGenJets = dynamic_cast<TClonesArray*>(fAODExtension->GetAOD()->FindListObject(fJetBranchNameMC.Data())); | |
765 | } | |
766 | if(fAODIn&&!aodGenJets){ | |
767 | aodGenJets = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(fJetBranchNameMC.Data())); | |
768 | } | |
769 | ||
770 | if(!aodGenJets){ | |
771 | Printf("%s:%d no generated Jet array with name %s in AOD", | |
772 | (char*)__FILE__,__LINE__, fJetBranchNameMC.Data()); | |
773 | PostData(1, fOutputList); | |
774 | return; | |
775 | } | |
776 | ||
777 | fListJetsGen->Clear(); | |
778 | ||
779 | //serarch for charged trigger at the MC generator level | |
8e103a56 | 780 | Int_t indexTriggGen = -1; |
781 | Double_t ptTriggGen = -1; | |
782 | Int_t iCounterGen = 0; | |
783 | Int_t triggersMC[200];//list of trigger candidates | |
784 | Int_t ntriggersMC = 0; //index in triggers array | |
785 | ||
80ac66f6 | 786 | if(fESD){//ESD input |
787 | ||
788 | AliMCEvent* mcEvent = MCEvent(); | |
789 | if(!mcEvent){ | |
790 | PostData(1, fOutputList); | |
791 | return; | |
ea603b64 | 792 | } |
793 | ||
794 | AliGenPythiaEventHeader *pythiaGenHeader = AliAnalysisHelperJetTasks::GetPythiaEventHeader(mcEvent); | |
795 | if(pythiaGenHeader){ | |
796 | nTrials = pythiaGenHeader->Trials(); | |
797 | ptHard = pythiaGenHeader->GetPtHard(); | |
798 | fh1PtHard->Fill(ptHard,eventW); | |
799 | fh1PtHardNoW->Fill(ptHard,1); | |
800 | fh1PtHardTrials->Fill(ptHard,nTrials); | |
801 | } | |
802 | ||
80ac66f6 | 803 | for(Int_t it = 0; it < mcEvent->GetNumberOfTracks(); it++){ |
804 | if(!mcEvent->IsPhysicalPrimary(it)) continue; | |
805 | AliMCParticle* part = (AliMCParticle*) mcEvent->GetTrack(it); | |
806 | if(!part) continue; | |
807 | if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){ | |
8e103a56 | 808 | |
809 | if(fHardest==0 && ntriggersMC<200){//single inclusive trigger | |
810 | if(indexTriggGen > -1){//trigger candidater was found | |
811 | triggersMC[ntriggersMC] = indexTriggGen; | |
812 | ntriggersMC++; | |
813 | } | |
814 | } | |
815 | ||
80ac66f6 | 816 | iCounterGen++; |
817 | } | |
818 | } | |
819 | }else{ //AOD input | |
820 | if(!fAODIn){ | |
821 | PostData(1, fOutputList); | |
822 | return; | |
823 | } | |
824 | TClonesArray *tca = dynamic_cast<TClonesArray*>(fAODIn->FindListObject(AliAODMCParticle::StdBranchName())); | |
825 | if(!tca){ | |
826 | PostData(1, fOutputList); | |
827 | return; | |
828 | } | |
829 | for(Int_t it = 0; it < tca->GetEntriesFast(); it++){ | |
830 | AliAODMCParticle *part = dynamic_cast<AliAODMCParticle*>(tca->At(it)); | |
831 | if(!part) continue; | |
832 | if(!part->IsPhysicalPrimary()) continue; | |
833 | if(SelectMCGenTracks((AliVParticle*) part, &particleListGen, ptTriggGen, indexTriggGen, iCounterGen)){ | |
834 | ||
8e103a56 | 835 | if(fHardest==0 && ntriggersMC<200){//single inclusive trigger |
836 | if(indexTriggGen > -1){ //trigger candidater was found | |
837 | triggersMC[ntriggersMC] = indexTriggGen; | |
838 | ntriggersMC++; | |
839 | } | |
840 | } | |
841 | ||
80ac66f6 | 842 | iCounterGen++; |
843 | } | |
844 | } | |
845 | } | |
8e103a56 | 846 | |
847 | if(fHardest==0){ //single inclusive trigger | |
848 | if(ntriggersMC>0){ //there is at least one trigger | |
849 | Int_t rnd = fRandom->Integer(ntriggersMC); //0 to ntriggers-1 | |
850 | indexTriggGen = triggersMC[rnd]; | |
80ac66f6 | 851 | }else{ |
8e103a56 | 852 | indexTriggGen = -1; //trigger not found |
80ac66f6 | 853 | } |
854 | } | |
8e103a56 | 855 | |
856 | //================== Fill jet list =================== | |
80ac66f6 | 857 | if(aodGenJets){ |
858 | if(fDebug) Printf("########## %s: %d jets",fJetBranchNameMC.Data(), aodGenJets->GetEntriesFast()); | |
859 | ||
860 | for(Int_t igJet = 0; igJet < aodGenJets->GetEntriesFast(); igJet++) { | |
861 | AliAODJet *jetGen = dynamic_cast<AliAODJet*>((*aodGenJets)[igJet]); | |
862 | if(jetGen) fListJetsGen->Add(jetGen); | |
863 | } | |
8e103a56 | 864 | } |
865 | ||
866 | //============ Generator trigger+jet ================== | |
867 | Int_t ilowGen = (fHardest==0 || fHardest==1) ? indexTriggGen : 0; | |
868 | Int_t ihighGen = (fHardest==0 || fHardest==1) ? indexTriggGen+1 : particleListGen.GetEntries(); | |
869 | ||
870 | for(Int_t igen= ilowGen; igen< ihighGen; igen++){ //loop over possible trigger | |
871 | indexTriggGen = igen; //trigger hadron | |
872 | ||
873 | if(indexTriggGen == -1) continue; | |
874 | AliVParticle* triggerGen=(AliVParticle*) particleListGen.At(indexTriggGen); | |
875 | if(!triggerGen) continue; | |
876 | ||
877 | if(fHardest >= 2){ | |
8f9e503d | 878 | if(triggerGen->Pt() < 6.0) continue; //all hadrons pt>10 |
8e103a56 | 879 | } |
880 | if(TMath::Abs((Float_t) triggerGen->Eta()) > fTriggerEtaCut) continue; | |
881 | ||
882 | ptTriggGen = triggerGen->Pt(); //count triggers | |
883 | fh2NtriggersGen->Fill(centValue, ptTriggGen); | |
884 | ||
80ac66f6 | 885 | //Count jets and trigger-jet pairs at MC generator level |
8e103a56 | 886 | if(!aodGenJets) continue; |
80ac66f6 | 887 | for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){ |
888 | AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij)); | |
889 | if(!jet) continue; | |
890 | Double_t etaJetGen = jet->Eta(); | |
80ac66f6 | 891 | |
892 | if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){ | |
80ac66f6 | 893 | |
8e103a56 | 894 | Double_t dphi = RelativePhi(triggerGen->Phi(), jet->Phi()); |
895 | if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue; | |
80ac66f6 | 896 | |
8e103a56 | 897 | //Centrality, A, pT, pTtrigg |
898 | Double_t fillspecgen[] = { centValue, | |
80ac66f6 | 899 | jet->EffectiveAreaCharged(), |
900 | jet->Pt(), | |
901 | ptTriggGen | |
902 | }; | |
8e103a56 | 903 | fHJetSpecGen->Fill(fillspecgen); |
904 | }//back to back jet-trigger pair | |
905 | }//jet loop | |
906 | }//trigger loop | |
907 | ||
908 | ||
909 | //================ RESPONSE MATRIX =============== | |
910 | if(aodGenJets){ | |
911 | //Count jets and trigger-jet pairs at MC generator level | |
912 | for(Int_t ij=0; ij<fListJetsGen->GetEntries(); ij++){ | |
913 | AliAODJet* jet = (AliAODJet*)(fListJetsGen->At(ij)); | |
914 | if(!jet) continue; | |
915 | Double_t etaJetGen = jet->Eta(); | |
916 | Double_t ptJetGen = jet->Pt(); | |
917 | ||
918 | if((fJetEtaMin<=etaJetGen) && (etaJetGen<=fJetEtaMax)){ | |
919 | fhJetPtGen->Fill(ptJetGen); // generator level pt spectum of jets response mx normalization | |
80ac66f6 | 920 | } |
921 | } | |
80ac66f6 | 922 | if(fListJets->GetEntries()>0 && fListJetsGen->GetEntries()>0){ //at least some reconstructed jets |
923 | Int_t ng = (Int_t) fListJetsGen->GetEntries(); | |
924 | Int_t nr = (Int_t) fListJets->GetEntries(); | |
925 | ||
926 | //Find closest MC generator - reconstructed jets | |
927 | if(faGenIndex.GetSize()<nr) faGenIndex.Set(nr); //idx of gen jet assoc to rec jet | |
928 | if(faRecIndex.GetSize()<ng) faRecIndex.Set(ng); //idx of rec jet assoc to gen jet | |
929 | ||
930 | if(fDebug){ | |
931 | Printf("New Rec List %d gen index Array %d",nr,faGenIndex.GetSize()); | |
932 | Printf("New Gen List %d rec index Array %d",ng,faRecIndex.GetSize()); | |
933 | } | |
934 | //matching of MC genrator level and reconstructed jets | |
935 | AliAnalysisHelperJetTasks::GetClosestJets(fListJetsGen,ng,fListJets,nr,faGenIndex,faRecIndex,fDebug); | |
936 | ||
937 | // Fill response matrix | |
938 | for(Int_t ir = 0; ir < nr; ir++){ | |
939 | AliAODJet *recJet = (AliAODJet*) fListJets->At(ir); | |
940 | Double_t etaJetRec = recJet->Eta(); | |
941 | Double_t ptJetRec = recJet->Pt(); | |
942 | //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc | |
943 | ||
944 | if((fJetEtaMin <= etaJetRec) && (etaJetRec <= fJetEtaMax)){ | |
945 | Int_t ig = faGenIndex[ir]; //associated generator level jet | |
946 | if(ig >= 0 && ig < ng){ | |
947 | if(fDebug > 10) Printf("%s:%d ig = %d ir = %d",(char*)__FILE__,__LINE__,ig,ir); | |
948 | AliAODJet *genJet = (AliAODJet*) fListJetsGen->At(ig); | |
949 | Double_t ptJetGen = genJet->Pt(); | |
950 | Double_t etaJetGen = genJet->Eta(); | |
951 | ||
952 | //fill response matrix if generator and reconstructed jets are within |eta|<0.9-fiduc | |
953 | if((fJetEtaMin <= etaJetGen) && (etaJetGen <= fJetEtaMax)){ | |
954 | fhJetPtGenVsJetPtRec->Fill(ptJetRec, ptJetGen); | |
955 | } | |
956 | }//ig>=0 | |
957 | }//rec jet in eta acceptance | |
958 | }//loop over reconstructed jets | |
959 | }// # of rec jets >0 | |
960 | }//pointer MC generator jets | |
961 | } //analyze generator level MC | |
962 | ||
8e103a56 | 963 | //============= RECONSTRUCTED INCLUSIVE JETS =============== |
ad869500 | 964 | |
ad869500 | 965 | Double_t etaJet = 0.0; |
966 | Double_t pTJet = 0.0; | |
967 | Double_t areaJet = 0.0; | |
968 | Double_t phiJet = 0.0; | |
7fc3d134 | 969 | Int_t indexLeadingJet = -1; |
970 | Double_t pTLeadingJet = -10.0; | |
7fc3d134 | 971 | Double_t areaLeadingJet = -10.0; |
8e103a56 | 972 | |
ad869500 | 973 | for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){ |
974 | AliAODJet* jet = (AliAODJet*)(fListJets->At(ij)); | |
975 | if(!jet) continue; | |
976 | etaJet = jet->Eta(); | |
977 | phiJet = jet->Phi(); | |
978 | pTJet = jet->Pt(); | |
979 | if(pTJet==0) continue; | |
980 | ||
981 | if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue; | |
982 | areaJet = jet->EffectiveAreaCharged(); | |
7fc3d134 | 983 | |
ad869500 | 984 | //Jet Diagnostics--------------------------------- |
8e103a56 | 985 | fhJetPhi->Fill(pTJet, RelativePhi(phiJet,0.0)); //phi -pi,pi |
986 | fhJetEta->Fill(pTJet, etaJet); | |
7fc3d134 | 987 | //search for leading jet |
988 | if(pTJet > pTLeadingJet){ | |
989 | indexLeadingJet = ij; | |
990 | pTLeadingJet = pTJet; | |
7fc3d134 | 991 | areaLeadingJet = areaJet; |
992 | } | |
993 | ||
994 | // raw spectra of INCLUSIVE jets | |
80ac66f6 | 995 | //Centrality, pTjet, A |
7fc3d134 | 996 | Double_t fillraw[] = { centValue, |
997 | pTJet, | |
7fc3d134 | 998 | areaJet |
999 | }; | |
1000 | fHJetPtRaw->Fill(fillraw); | |
8e103a56 | 1001 | } |
7fc3d134 | 1002 | |
8e103a56 | 1003 | if(indexLeadingJet > -1){ |
1004 | // raw spectra of LEADING jets | |
1005 | //Centrality, pTjet, A | |
1006 | Double_t fillleading[] = { centValue, | |
1007 | pTLeadingJet, | |
1008 | areaLeadingJet | |
1009 | }; | |
1010 | fHLeadingJetPtRaw->Fill(fillleading); | |
1011 | } | |
1012 | ||
1013 | ||
1014 | // =============== RECONSTRUCTED TRIGGER-JET PAIRS ================ | |
1015 | //Find Hadron trigger | |
1016 | TList particleList; //list of tracks | |
1017 | Int_t indexTrigg = GetListOfTracks(&particleList); //index of trigger hadron in Particle list | |
1018 | if(fIsMC) FillEffHistos(&particleList, &particleListGen); //Fill efficiency histos | |
7fc3d134 | 1019 | |
8e103a56 | 1020 | //set ranges of the trigger loop |
1021 | Int_t ilow = (fHardest==0 || fHardest==1) ? indexTrigg : 0; | |
1022 | Int_t ihigh = (fHardest==0 || fHardest==1) ? indexTrigg+1 : particleList.GetEntries(); | |
7fc3d134 | 1023 | |
8e103a56 | 1024 | for(Int_t itrk= ilow; itrk< ihigh; itrk++){ //loop over possible trigger |
1025 | indexTrigg = itrk; //trigger hadron with pT >10 GeV | |
7fc3d134 | 1026 | |
8e103a56 | 1027 | if(indexTrigg < 0) continue; |
1028 | ||
1029 | AliVParticle *triggerHadron = (AliVParticle*) particleList.At(indexTrigg); | |
1030 | if(!triggerHadron){ | |
1031 | PostData(1, fOutputList); | |
1032 | return; | |
1033 | } | |
1034 | if(fHardest >= 2){ | |
7d9fd386 | 1035 | if(triggerHadron->Pt() < 6.0) continue; //all hadrons pt>6 |
8e103a56 | 1036 | } |
1037 | if(TMath::Abs((Float_t) triggerHadron->Eta())> fTriggerEtaCut) continue; | |
7fc3d134 | 1038 | |
8e103a56 | 1039 | //Fill trigger histograms |
1040 | fh2Ntriggers->Fill(centValue,triggerHadron->Pt()); //trigger pT | |
1041 | fhTriggerPhi->Fill(triggerHadron->Pt(),RelativePhi(triggerHadron->Phi(),0.0)); //phi -pi,pi | |
1042 | fhTriggerEta->Fill(triggerHadron->Pt(),triggerHadron->Eta()); | |
7fc3d134 | 1043 | |
1044 | ||
8e103a56 | 1045 | //---------- make trigger-jet pairs --------- |
1046 | Int_t injet4 = 0; | |
1047 | Int_t injet = 0; | |
1048 | ||
1049 | for(Int_t ij=0; ij<fListJets->GetEntries(); ij++){ | |
1050 | AliAODJet* jet = (AliAODJet*)(fListJets->At(ij)); | |
1051 | if(!jet) continue; | |
1052 | etaJet = jet->Eta(); | |
1053 | phiJet = jet->Phi(); | |
1054 | pTJet = jet->Pt(); | |
1055 | if(pTJet==0) continue; | |
1056 | ||
1057 | if((etaJet<fJetEtaMin) || (etaJet>fJetEtaMax)) continue; | |
1058 | areaJet = jet->EffectiveAreaCharged(); | |
1059 | if(areaJet >= 0.07) injet++; | |
1060 | if(areaJet >= 0.4) injet4++; | |
1061 | ||
1062 | Double_t dphi = RelativePhi(triggerHadron->Phi(), phiJet); | |
1063 | fhDphiTriggerJet->Fill(dphi); //Input | |
1064 | ||
1065 | //Dphi versus jet pT | |
1066 | //Centrality, Dphi=phiTrig-phiJet, pTjet, pTtrigg | |
1067 | Double_t filldp[] = { centValue, | |
1068 | dphi, | |
1069 | pTJet, | |
1070 | triggerHadron->Pt() | |
1071 | }; | |
1072 | fHDphiVsJetPtAll->Fill(filldp); | |
1073 | ||
1074 | // Select back to back trigger - jet pairs | |
1075 | if(TMath::Abs((Double_t) dphi) < fkDeltaPhiCut) continue; | |
1076 | fhDphiTriggerJetAccept->Fill(dphi); //Accepted | |
1077 | ||
1078 | ||
1079 | //Centrality, A, pTjet, pTtrigg | |
1080 | Double_t fillspec[] = { centValue, | |
1081 | areaJet, | |
1082 | pTJet, | |
1083 | triggerHadron->Pt() | |
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 |