]>
Commit | Line | Data |
---|---|---|
b6e50056 | 1 | // ************************************** |
2 | // Full pp jet task - ESD input only | |
7c915297 | 3 | // Extract the jet spectrum and all the |
4 | // systematic uncertainties | |
b6e50056 | 5 | // -R. Ma, Mar 2011 |
6 | // ************************************** | |
7 | ||
8 | #include <TCanvas.h> | |
9 | #include <TChain.h> | |
10 | #include <TFormula.h> | |
11 | #include <TH1.h> | |
12 | #include <TH2.h> | |
13 | #include <TH3.h> | |
14 | #include <TProfile2D.h> | |
15 | #include <THnSparse.h> | |
16 | #include <TROOT.h> | |
17 | #include <TTree.h> | |
18 | #include <TArrayI.h> | |
19 | #include <TClonesArray.h> | |
20 | #include <TRandom3.h> | |
21 | #include <TFile.h> | |
22 | #include <TF1.h> | |
23 | ||
24 | #include "AliAODEvent.h" | |
25 | #include "AliAODInputHandler.h" | |
26 | #include "AliAnalysisManager.h" | |
27 | #include "AliAnalysisTask.h" | |
28 | #include "AliCentrality.h" | |
29 | #include "AliAnalysisTaskFullppJet.h" | |
30 | #include "AliESDEvent.h" | |
31 | #include "AliESDInputHandler.h" | |
32 | #include "AliESDtrackCuts.h" | |
33 | #include "AliVParticle.h" | |
34 | #include "AliInputEventHandler.h" | |
35 | #include "AliEMCALGeometry.h" | |
36 | #include "AliEMCALRecoUtils.h" | |
37 | #include "TGeoManager.h" | |
38 | #include "AliMCEvent.h" | |
7c915297 | 39 | #include "AliMCParticle.h" |
b6e50056 | 40 | #include "AliStack.h" |
41 | #include "AliGenEventHeader.h" | |
42 | #include "AliGenPythiaEventHeader.h" | |
43 | #include "TGeoGlobalMagField.h" | |
44 | ||
45 | #include "AliAODJet.h" | |
46 | #include "AliFJWrapper.h" | |
47 | ||
48 | #include <iostream> | |
49 | using std::cout; | |
50 | using std::cerr; | |
51 | using std::endl; | |
a7eebe5b | 52 | using std::vector; |
b6e50056 | 53 | |
54 | ClassImp(AliAnalysisTaskFullppJet) | |
55 | ||
7c915297 | 56 | const Float_t kRadius[3] = {0.4,0.2,0.3}; |
b6e50056 | 57 | const Double_t kPI = TMath::Pi(); |
7c915297 | 58 | const Double_t kdRCut[3] = {0.25,0.1,0.15}; |
b6e50056 | 59 | |
60 | //________________________________________________________________________ | |
61 | AliAnalysisTaskFullppJet::AliAnalysisTaskFullppJet() : | |
62 | AliAnalysisTaskSE("default"), | |
7c915297 | 63 | fVerbosity(0), fEDSFileCounter(0), fNTracksPerChunk(0), fRejectPileup(kFALSE), fRejectExoticTrigger(kTRUE), |
64 | fAnaType(0), fPeriod("lhc11a"), fESD(0), fAOD(0), fMC(0), fStack(0), fTrackArray(0x0), fClusterArray(0x0), fMcPartArray(0x0), | |
65 | fIsMC(kFALSE), fPhySelForMC(kFALSE), fChargedMC(kFALSE), fXsecScale(0.), | |
b6e50056 | 66 | fCentrality(99), fZVtxMax(10), |
7c915297 | 67 | fTriggerType(-1), fCheckTriggerMask(kTRUE), fIsTPCOnlyVtx(kFALSE), |
b6e50056 | 68 | fIsExoticEvent3GeV(kFALSE), fIsExoticEvent5GeV(kFALSE), fIsEventTriggerBit(kFALSE), fOfflineTrigger(kFALSE), fTriggerMask(0x0), |
69 | fGeom(0x0), fRecoUtil(0x0), | |
70 | fEsdTrackCuts(0x0), fHybridTrackCuts1(0x0), fHybridTrackCuts2(0x0),fTrackCutsType(0), fKinCutType(0), fTrkEtaMax(0.9), | |
71 | fdEdxMin(73), fdEdxMax(90), fEoverPMin(0.8), fEoverPMax(1.2), | |
72 | fMatchType(0), fRejectExoticCluster(kTRUE), fRemoveBadChannel(kFALSE), fUseGoodSM(kFALSE), | |
73 | fStudySubEInHC(kFALSE), fStudyMcOverSubE(kFALSE), | |
7c915297 | 74 | fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0), |
b6e50056 | 75 | fJetNEFMin(0.01), fJetNEFMax(0.98), |
76 | fSpotGoodJet(kFALSE), fFindChargedOnlyJet(kFALSE), fFindNeutralOnlyJet(kFALSE), | |
7c915297 | 77 | fCheckTrkEffCorr(kFALSE), fTrkEffCorrCutZ(0.3), fRandomGen(0x0), fRunUE(0), fCheckTPCOnlyVtx(0), fRunSecondaries(0), |
78 | fSysJetTrigEff(kFALSE), fVaryJetTrigEff(0), | |
79 | fSysTrkPtRes(kFALSE), fVaryTrkPtRes(0), fSysTrkEff(kFALSE), fVaryTrkEff(0.), | |
b6e50056 | 80 | fSysTrkClsMth(kFALSE), fCutdEta(0.05), fCutdPhi(0.05), |
7c915297 | 81 | fSysNonLinearity(kFALSE), fNonLinear(0x0), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0), |
82 | fSysClusterizer(0), | |
b6e50056 | 83 | fNonStdBranch(""), fNonStdFile(""), fAlgorithm("aKt"), fRadius("0.4"), fRecombinationScheme(5), |
84 | fConstrainChInEMCal(kFALSE), fRejectNK(kFALSE), fRejectWD(kFALSE), fSmearMC(kFALSE), fTrkPtResData(0x0), | |
7c915297 | 85 | fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventCount(0x0), fhJetEventStat(0x0), fhEventStatTPCVtx(0x0), fhChunkQA(0x0) |
b6e50056 | 86 | { |
87 | // Constructor | |
88 | ||
89 | // Define input and output slots here | |
90 | // Input slot #0 works with a TChain | |
91 | DefineInput(0, TChain::Class()); | |
92 | DefineOutput(1, TList::Class()); | |
93 | // Output slot #0 id reserved by the base class for AOD | |
94 | for(Int_t i=0; i<2; i++) | |
95 | { | |
96 | fTrkPtMin[i] = 0.15; | |
97 | fTrkPtMax[i] = 200; | |
98 | fClsEtMin[i] = 0.15; | |
99 | fClsEtMax[i] = 200; | |
7c915297 | 100 | |
101 | fVertexGenZ[i] = 0x0; | |
102 | fEventZ[i] = 0x0; | |
103 | fhNTrials[i] = 0x0; | |
104 | fhNMatchedTrack[i] = 0x0; | |
105 | fhClsE[i] = 0x0; | |
106 | for(Int_t k=0; k<2; k++) | |
107 | { | |
108 | fhSysClusterE[i][k] = 0x0; | |
109 | fhSysNCellVsClsE[i][k] = 0x0; | |
110 | } | |
b6e50056 | 111 | } |
112 | ||
7c915297 | 113 | for(Int_t r=0; r<kNRadii; r++) |
b6e50056 | 114 | { |
b6e50056 | 115 | for(Int_t j=0; j<5; j++) |
116 | { | |
117 | fHCOverSubE[r][j] = 0x0; | |
118 | fHCOverSubEFrac[r][j] = 0x0; | |
119 | for(Int_t k=0; k<2; k++) | |
120 | { | |
121 | fhSubClsEVsJetPt[k][r][j] = 0x0; | |
122 | fhHCTrkPtClean[k][r][j] = 0x0; | |
123 | fhHCTrkPtAmbig[k][r][j] = 0x0; | |
124 | } | |
7c915297 | 125 | |
126 | if(j<4) fhSubEVsTrkPt[r][j] = 0x0; | |
127 | ||
128 | if(j<3) | |
129 | { | |
130 | fJetCount[j][r] = 0x0; | |
131 | fhNeutralPtInJet[j][r] = 0x0; | |
132 | fhTrigNeuPtInJet[j][r] = 0x0; | |
133 | fhChargedPtInJet[j][r] = 0x0; | |
134 | fhLeadNePtInJet[j][r] = 0x0; | |
135 | fhLeadChPtInJet[j][r] = 0x0; | |
136 | fJetEnergyFraction[j][r] = 0x0; | |
137 | fJetNPartFraction[j][r] = 0x0; | |
138 | } | |
139 | if(j<2) | |
b6e50056 | 140 | { |
7c915297 | 141 | fRelTrkCon[j][r] = 0x0; |
142 | fhFcrossVsZleading[j][r] = 0x0; | |
143 | fhChLeadZVsJetPt[j][r] = 0x0; | |
144 | fhJetPtWithTrkThres[j][r]= 0x0; | |
145 | fhJetPtWithClsThres[j][r]= 0x0; | |
146 | for(Int_t k=0; k<2; k++) | |
147 | { | |
148 | fhJetPtVsLowPtCons[j][r][k] = 0x0; | |
149 | } | |
150 | fhJetPtInExoticEvent[j][r] = 0x0; | |
151 | fhJetInTPCOnlyVtx[j][r] = 0x0; | |
152 | fhCorrTrkEffPtBin[j][r] = 0x0; | |
153 | for(Int_t i=0; i<kNBins; i++) | |
154 | { | |
155 | fhCorrTrkEffSample[j][r][i] = 0x0; | |
156 | } | |
157 | } | |
158 | ||
159 | if(r==0 && j<3) | |
160 | { | |
161 | for(Int_t k=0; k<2; k++) | |
162 | { | |
163 | for(Int_t l=0; l<2; l++) | |
164 | { | |
165 | fhUEJetPtNorm[j][k][l] = 0x0; | |
166 | fhUEJetPtVsSumPt[j][k][l] = 0x0; | |
167 | fhUEJetPtVsConsPt[j][k][l] = 0x0; | |
168 | } | |
169 | } | |
b6e50056 | 170 | } |
171 | } | |
7c915297 | 172 | for(Int_t i=0; i<2; i++) |
173 | { | |
174 | fhNKFracVsJetPt[i][r] = 0x0; | |
175 | fhWeakFracVsJetPt[i][r] = 0x0; | |
176 | fhJetResponseNK[i][r] = 0x0; | |
177 | fhJetResponseWP[i][r] = 0x0; | |
178 | fhJetResolutionNK[i][r] = 0x0; | |
179 | fhJetResolutionWP[i][r] = 0x0; | |
180 | fhJetResponseNKSM[i][r] = 0x0; | |
181 | fhJetResponseWPSM[i][r] = 0x0; | |
182 | fhJetResolutionNKSM[i][r] = 0x0; | |
183 | fhJetResolutionWPSM[i][r] = 0x0; | |
184 | } | |
b6e50056 | 185 | } |
186 | ||
187 | for(Int_t s=0; s<3; s++) | |
188 | { | |
189 | for(Int_t a=0; a<2; a++) | |
190 | { | |
7c915297 | 191 | for(Int_t r=0; r<kNRadii; r++) |
b6e50056 | 192 | { |
193 | fDetJetFinder[s][a][r] = 0x0; | |
194 | fJetTCA[s][a][r] = 0x0; | |
195 | if(s==0 && a==0) | |
196 | { | |
197 | fTrueJetFinder[r] = 0x0; | |
198 | fMcTruthAntikt[r] = 0x0; | |
199 | } | |
200 | } | |
201 | } | |
202 | } | |
203 | ||
204 | for(Int_t j=0; j<3; j++) | |
205 | { | |
206 | fTrkEffFunc[j] = 0x0; | |
7c915297 | 207 | fhSecondaryResponse[j] = 0x0; |
b6e50056 | 208 | } |
209 | for(Int_t i=0; i<10; i++) | |
210 | { | |
211 | fTriggerCurve[i] = 0x0; | |
212 | fTriggerEfficiency[i] = 0x0; | |
213 | } | |
214 | ||
215 | } | |
216 | ||
217 | //________________________________________________________________________ | |
218 | AliAnalysisTaskFullppJet::AliAnalysisTaskFullppJet(const char *name) : | |
219 | AliAnalysisTaskSE(name), | |
7c915297 | 220 | fVerbosity(0), fEDSFileCounter(0), fNTracksPerChunk(0), fRejectPileup(kFALSE), fRejectExoticTrigger(kTRUE), |
221 | fAnaType(0), fPeriod("lhc11a"), fESD(0), fAOD(0), fMC(0), fStack(0), fTrackArray(0x0), fClusterArray(0x0), fMcPartArray(0x0), | |
222 | fIsMC(kFALSE), fPhySelForMC(kFALSE), fChargedMC(kFALSE), fXsecScale(0.), | |
b6e50056 | 223 | fCentrality(99), fZVtxMax(10), |
7c915297 | 224 | fTriggerType(-1), fCheckTriggerMask(kTRUE), fIsTPCOnlyVtx(kFALSE), |
b6e50056 | 225 | fIsExoticEvent3GeV(kFALSE), fIsExoticEvent5GeV(kFALSE), fIsEventTriggerBit(kFALSE), fOfflineTrigger(kFALSE), fTriggerMask(0x0), |
226 | fGeom(0x0), fRecoUtil(0x0), | |
227 | fEsdTrackCuts(0x0), fHybridTrackCuts1(0x0), fHybridTrackCuts2(0x0), fTrackCutsType(0), fKinCutType(0), fTrkEtaMax(0.9), | |
228 | fdEdxMin(73), fdEdxMax(90), fEoverPMin(0.8), fEoverPMax(1.2), | |
229 | fMatchType(0), fRejectExoticCluster(kTRUE), fRemoveBadChannel(kFALSE), fUseGoodSM(kFALSE), | |
230 | fStudySubEInHC(kFALSE), fStudyMcOverSubE(kFALSE), | |
7c915297 | 231 | fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0), |
b6e50056 | 232 | fJetNEFMin(0.01), fJetNEFMax(0.98), |
233 | fSpotGoodJet(kFALSE), fFindChargedOnlyJet(kFALSE), fFindNeutralOnlyJet(kFALSE), | |
7c915297 | 234 | fCheckTrkEffCorr(kFALSE), fTrkEffCorrCutZ(0.3), fRandomGen(0x0), fRunUE(0), fCheckTPCOnlyVtx(0), fRunSecondaries(0), |
235 | fSysJetTrigEff(kFALSE), fVaryJetTrigEff(0), | |
236 | fSysTrkPtRes(kFALSE), fVaryTrkPtRes(0), fSysTrkEff(kFALSE), fVaryTrkEff(0.), | |
b6e50056 | 237 | fSysTrkClsMth(kFALSE), fCutdEta(0.05), fCutdPhi(0.05), |
7c915297 | 238 | fSysNonLinearity(kFALSE), fNonLinear(0x0), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0), |
239 | fSysClusterizer(0), | |
b6e50056 | 240 | fNonStdBranch(""), fNonStdFile(""), fAlgorithm("aKt"), fRadius("0.4"), fRecombinationScheme(5), |
241 | fConstrainChInEMCal(kFALSE), fRejectNK(kFALSE), fRejectWD(kFALSE), fSmearMC(kFALSE), fTrkPtResData(0x0), | |
7c915297 | 242 | fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventCount(0x0), fhJetEventStat(0x0), fhEventStatTPCVtx(0x0), fhChunkQA(0x0) |
b6e50056 | 243 | { |
244 | // Constructor | |
245 | ||
246 | // Define input and output slots here | |
247 | // Input slot #0 works with a TChain | |
248 | DefineInput(0, TChain::Class()); | |
249 | DefineOutput(1, TList::Class()); | |
250 | // Output slot #0 id reserved by the base class for AOD | |
251 | for(Int_t i=0; i<2; i++) | |
252 | { | |
253 | fTrkPtMin[i] = 0.15; | |
254 | fTrkPtMax[i] = 200; | |
255 | fClsEtMin[i] = 0.15; | |
256 | fClsEtMax[i] = 200; | |
7c915297 | 257 | |
258 | fVertexGenZ[i] = 0x0; | |
259 | fEventZ[i] = 0x0; | |
260 | fhNTrials[i] = 0x0; | |
261 | fhNMatchedTrack[i] = 0x0; | |
262 | fhClsE[i] = 0x0; | |
263 | for(Int_t k=0; k<2; k++) | |
264 | { | |
265 | fhSysClusterE[i][k] = 0x0; | |
266 | fhSysNCellVsClsE[i][k] = 0x0; | |
267 | } | |
b6e50056 | 268 | } |
269 | ||
7c915297 | 270 | for(Int_t r=0; r<kNRadii; r++) |
b6e50056 | 271 | { |
b6e50056 | 272 | for(Int_t j=0; j<5; j++) |
273 | { | |
274 | fHCOverSubE[r][j] = 0x0; | |
275 | fHCOverSubEFrac[r][j] = 0x0; | |
276 | for(Int_t k=0; k<2; k++) | |
277 | { | |
278 | fhSubClsEVsJetPt[k][r][j] = 0x0; | |
279 | fhHCTrkPtClean[k][r][j] = 0x0; | |
280 | fhHCTrkPtAmbig[k][r][j] = 0x0; | |
281 | } | |
7c915297 | 282 | |
283 | if(j<4) fhSubEVsTrkPt[r][j] = 0x0; | |
284 | ||
285 | if(j<3) | |
286 | { | |
287 | fJetCount[j][r] = 0x0; | |
288 | fhNeutralPtInJet[j][r] = 0x0; | |
289 | fhTrigNeuPtInJet[j][r] = 0x0; | |
290 | fhChargedPtInJet[j][r] = 0x0; | |
291 | fhLeadNePtInJet[j][r] = 0x0; | |
292 | fhLeadChPtInJet[j][r] = 0x0; | |
293 | fJetEnergyFraction[j][r] = 0x0; | |
294 | fJetNPartFraction[j][r] = 0x0; | |
295 | } | |
296 | if(j<2) | |
b6e50056 | 297 | { |
7c915297 | 298 | fRelTrkCon[j][r] = 0x0; |
299 | fhFcrossVsZleading[j][r] = 0x0; | |
300 | fhChLeadZVsJetPt[j][r] = 0x0; | |
301 | fhJetPtWithTrkThres[j][r]= 0x0; | |
302 | fhJetPtWithClsThres[j][r]= 0x0; | |
303 | for(Int_t k=0; k<2; k++) | |
304 | { | |
305 | fhJetPtVsLowPtCons[j][r][k] = 0x0; | |
306 | } | |
307 | fhJetPtInExoticEvent[j][r] = 0x0; | |
308 | fhJetInTPCOnlyVtx[j][r] = 0x0; | |
309 | fhCorrTrkEffPtBin[j][r] = 0x0; | |
310 | for(Int_t i=0; i<kNBins; i++) | |
311 | { | |
312 | fhCorrTrkEffSample[j][r][i] = 0x0; | |
313 | } | |
314 | } | |
315 | ||
316 | if(r==0 && j<3) | |
317 | { | |
318 | for(Int_t k=0; k<2; k++) | |
319 | { | |
320 | for(Int_t l=0; l<2; l++) | |
321 | { | |
322 | fhUEJetPtNorm[j][k][l] = 0x0; | |
323 | fhUEJetPtVsSumPt[j][k][l] = 0x0; | |
324 | fhUEJetPtVsConsPt[j][k][l] = 0x0; | |
325 | } | |
326 | } | |
b6e50056 | 327 | } |
328 | } | |
7c915297 | 329 | for(Int_t i=0; i<2; i++) |
330 | { | |
331 | fhNKFracVsJetPt[i][r] = 0x0; | |
332 | fhWeakFracVsJetPt[i][r] = 0x0; | |
333 | fhJetResponseNK[i][r] = 0x0; | |
334 | fhJetResponseWP[i][r] = 0x0; | |
335 | fhJetResolutionNK[i][r] = 0x0; | |
336 | fhJetResolutionWP[i][r] = 0x0; | |
337 | fhJetResponseNKSM[i][r] = 0x0; | |
338 | fhJetResponseWPSM[i][r] = 0x0; | |
339 | fhJetResolutionNKSM[i][r] = 0x0; | |
340 | fhJetResolutionWPSM[i][r] = 0x0; | |
341 | } | |
b6e50056 | 342 | } |
343 | ||
344 | for(Int_t s=0; s<3; s++) | |
345 | { | |
346 | for(Int_t a=0; a<2; a++) | |
347 | { | |
7c915297 | 348 | for(Int_t r=0; r<kNRadii; r++) |
b6e50056 | 349 | { |
350 | fDetJetFinder[s][a][r] = 0x0; | |
351 | fJetTCA[s][a][r] = 0x0; | |
352 | if(s==0 && a==0) | |
353 | { | |
354 | fTrueJetFinder[r] = 0x0; | |
355 | fMcTruthAntikt[r] = 0x0; | |
356 | } | |
357 | } | |
358 | } | |
359 | } | |
7c915297 | 360 | |
b6e50056 | 361 | for(Int_t j=0; j<3; j++) |
362 | { | |
363 | fTrkEffFunc[j] = 0x0; | |
7c915297 | 364 | fhSecondaryResponse[j] = 0x0; |
b6e50056 | 365 | } |
366 | for(Int_t i=0; i<10; i++) | |
367 | { | |
368 | fTriggerCurve[i] = 0x0; | |
369 | fTriggerEfficiency[i] = 0x0; | |
370 | } | |
371 | } | |
372 | ||
373 | //________________________________________________________________________ | |
374 | AliAnalysisTaskFullppJet::~AliAnalysisTaskFullppJet() | |
375 | { | |
376 | //Destructor | |
377 | ||
378 | if(fEsdTrackCuts) delete fEsdTrackCuts; | |
379 | if(fHybridTrackCuts1) delete fHybridTrackCuts1; | |
380 | if(fHybridTrackCuts2) delete fHybridTrackCuts2; | |
381 | if(fOutputList) | |
382 | { fOutputList->Delete(); delete fOutputList;} | |
383 | for(Int_t s=0; s<3; s++) | |
384 | { | |
385 | for(Int_t a=0; a<2; a++) | |
386 | { | |
7c915297 | 387 | for(Int_t r=0; r<kNRadii; r++) |
b6e50056 | 388 | { |
389 | if(fDetJetFinder[s][a][r]) delete fDetJetFinder[s][a][r]; | |
390 | if(fJetTCA[s][a][r]) { fJetTCA[s][a][r]->Delete(); delete fJetTCA[s][a][r]; } | |
391 | if(s==0 && a==0) | |
392 | { | |
393 | if(fTrueJetFinder[r]) delete fTrueJetFinder[r]; | |
394 | if(fMcTruthAntikt[r]) { fMcTruthAntikt[r]->Delete(); delete fMcTruthAntikt[r]; } | |
395 | } | |
396 | } | |
397 | } | |
398 | } | |
399 | if(fRandomGen) delete fRandomGen; | |
400 | for(Int_t i=0; i<3; i++) | |
401 | { | |
402 | if(fTrkEffFunc[i]) delete fTrkEffFunc[i]; | |
7c915297 | 403 | if(fhSecondaryResponse[i]) delete fhSecondaryResponse[i]; |
b6e50056 | 404 | } |
405 | for(Int_t i=0; i<10; i++) | |
406 | { | |
407 | if(fTriggerEfficiency[i]) delete fTriggerEfficiency[i]; | |
408 | if(fTriggerCurve[i]) delete fTriggerCurve[i]; | |
409 | } | |
7c915297 | 410 | for(Int_t r=0; r<kNRadii; r++) |
b6e50056 | 411 | { |
412 | for(Int_t j=0; j<2; j++) | |
413 | { | |
414 | if(fhCorrTrkEffPtBin[j][r]) delete fhCorrTrkEffPtBin[j][r]; | |
415 | for(Int_t i=0; i<kNBins; i++) | |
416 | { | |
417 | if(fhCorrTrkEffSample[j][r][i]) delete fhCorrTrkEffSample[j][r][i]; | |
418 | } | |
419 | } | |
420 | } | |
421 | if(fTrackArray) { fTrackArray->Delete(); delete fTrackArray; } | |
422 | if(fClusterArray) { fClusterArray->Delete(); delete fClusterArray; } | |
423 | if(fMcPartArray) { fMcPartArray->Reset(); delete fMcPartArray; } | |
424 | ||
425 | if(fRecoUtil) delete fRecoUtil; | |
426 | if(fClusterEResolution) delete fClusterEResolution; | |
7c915297 | 427 | if(fNonLinear) delete fNonLinear; |
b6e50056 | 428 | if(fTrkPtResData) delete fTrkPtResData; |
429 | } | |
430 | ||
431 | //________________________________________________________________________ | |
432 | Bool_t AliAnalysisTaskFullppJet::Notify() | |
433 | { | |
434 | // | |
435 | // Fill the number of tracks per chunk | |
436 | // | |
437 | ||
438 | fhChunkQA->SetBinContent(fEDSFileCounter,fNTracksPerChunk); | |
439 | fNTracksPerChunk = 0; | |
440 | fEDSFileCounter++; | |
441 | return kTRUE; | |
442 | } | |
443 | ||
444 | //________________________________________________________________________ | |
445 | void AliAnalysisTaskFullppJet::UserCreateOutputObjects() | |
446 | { | |
447 | // Create histograms | |
448 | // Called once | |
449 | // | |
7c915297 | 450 | |
451 | if(fRunUE) fFindChargedOnlyJet = kTRUE; | |
b6e50056 | 452 | |
453 | const Int_t nTrkPtBins = 100; | |
454 | const Float_t lowTrkPtBin=0, upTrkPtBin=100.; | |
455 | ||
456 | const Int_t nbins = 220; | |
457 | Double_t xbins[221]; | |
458 | for(Int_t i=0; i<nbins+1; i++) | |
459 | xbins[i] = i; | |
460 | ||
461 | OpenFile(1); | |
462 | fOutputList = new TList(); | |
463 | fOutputList->SetOwner(1); | |
464 | ||
7c915297 | 465 | fhJetEventStat = new TH1F("fhJetEventStat","Event statistics for jet analysis",12,0,12); |
b6e50056 | 466 | fhJetEventStat->GetXaxis()->SetBinLabel(1,"ALL"); |
467 | fhJetEventStat->GetXaxis()->SetBinLabel(2,"MB"); | |
468 | fhJetEventStat->GetXaxis()->SetBinLabel(3,"MB+vtx+10cm"); | |
469 | fhJetEventStat->GetXaxis()->SetBinLabel(4,"EMC"); | |
470 | fhJetEventStat->GetXaxis()->SetBinLabel(5,"EMC+vtx+10cm"); | |
471 | fhJetEventStat->GetXaxis()->SetBinLabel(6,"MB+vtx"); | |
472 | fhJetEventStat->GetXaxis()->SetBinLabel(7,"EMC+vtx"); | |
473 | fhJetEventStat->GetXaxis()->SetBinLabel(8,"TriggerBit"); | |
474 | fhJetEventStat->GetXaxis()->SetBinLabel(9,"LED"); | |
7c915297 | 475 | fhJetEventStat->GetXaxis()->SetBinLabel(10,"MB+TVtx+10cm"); |
476 | fhJetEventStat->GetXaxis()->SetBinLabel(11,"EMC+TVtx+10cm"); | |
477 | fhJetEventStat->GetXaxis()->SetBinLabel(12,"ALL-Pileup"); | |
b6e50056 | 478 | fOutputList->Add(fhJetEventStat); |
479 | ||
7c915297 | 480 | fhJetEventCount = new TH1F("fhJetEventCount","Event statistics for jet analysis",12,0,12); |
481 | fhJetEventCount->GetXaxis()->SetBinLabel(1,"ALL"); | |
482 | fhJetEventCount->GetXaxis()->SetBinLabel(2,"MB"); | |
483 | fhJetEventCount->GetXaxis()->SetBinLabel(3,"MB-pileup"); | |
484 | fhJetEventCount->GetXaxis()->SetBinLabel(4,"MB+Vtx"); | |
485 | fhJetEventCount->GetXaxis()->SetBinLabel(5,"MB+Vtx+10cm"); | |
486 | fhJetEventCount->GetXaxis()->SetBinLabel(6,"EMC"); | |
487 | fhJetEventCount->GetXaxis()->SetBinLabel(7,"EMC-pileup"); | |
488 | fhJetEventCount->GetXaxis()->SetBinLabel(8,"EMC+Vtx"); | |
489 | fhJetEventCount->GetXaxis()->SetBinLabel(9,"EMC+Vtx+10cm"); | |
490 | fhJetEventCount->GetXaxis()->SetBinLabel(10,"Good EMC"); | |
491 | fOutputList->Add(fhJetEventCount); | |
492 | ||
493 | if(fCheckTPCOnlyVtx) | |
494 | { | |
495 | fhEventStatTPCVtx = new TH1F("fhEventStatTPCVtx","Event statistics for TPC only vertex",9,0,9); | |
496 | fhEventStatTPCVtx->GetXaxis()->SetBinLabel(1,"FastOnly"); | |
497 | fhEventStatTPCVtx->GetXaxis()->SetBinLabel(2,"FastOnly+PVtx"); | |
498 | fhEventStatTPCVtx->GetXaxis()->SetBinLabel(3,"FastOnly+TVtx"); | |
499 | fhEventStatTPCVtx->GetXaxis()->SetBinLabel(4,"MB"); | |
500 | fhEventStatTPCVtx->GetXaxis()->SetBinLabel(5,"MB+PVtx"); | |
501 | fhEventStatTPCVtx->GetXaxis()->SetBinLabel(6,"MB+TVtx"); | |
502 | fhEventStatTPCVtx->GetXaxis()->SetBinLabel(7,"EMC"); | |
503 | fhEventStatTPCVtx->GetXaxis()->SetBinLabel(8,"EMC+PVtx"); | |
504 | fhEventStatTPCVtx->GetXaxis()->SetBinLabel(9,"EMC+TVtx"); | |
505 | fOutputList->Add(fhEventStatTPCVtx); | |
506 | } | |
b6e50056 | 507 | |
508 | fhChunkQA = new TH1F("fhChunkQA","# of hybrid tracks per chunk",200,0,200); | |
509 | fOutputList->Add(fhChunkQA); | |
510 | ||
511 | const Int_t dim1 = 3; | |
512 | Int_t nBins1[dim1] = {200,50,110}; | |
513 | Double_t lowBin1[dim1] = {0,0,0,}; | |
514 | Double_t upBin1[dim1] = {200,0.5,1.1}; | |
515 | ||
516 | const Int_t dim2 = 3; | |
517 | Int_t nBins2[dim2] = {200,50,50}; | |
518 | Double_t lowBin2[dim2] = {0,0,0,}; | |
519 | Double_t upBin2[dim2] = {200,0.5,50}; | |
520 | ||
b6e50056 | 521 | const char* triggerName[3] = {"MB","EMC","MC"}; |
522 | const char* triggerTitle[3] = {"MB","EMCal-trigger","MC true"}; | |
523 | const char* fraction[5] = {"MIP","30","50","70","100"}; | |
524 | const char* exotic[2] = {"3GeV","5GeV"}; | |
525 | const char* vertexType[2] = {"All MB","MB with vertex"}; | |
526 | const char *vertexName[2] = {"All","Vertex"}; | |
7c915297 | 527 | const char *clusterizerName[2] = {"before","after"}; |
528 | const char *UEName[2] = {"charged","charged_neutral"}; | |
529 | const char *UETitle[2] = {"charged","charged+neutral"}; | |
530 | const char *UEEventName[2] = {"LeadingJet","Back-To-Back"}; | |
b6e50056 | 531 | |
532 | if(fIsMC) | |
533 | { | |
534 | for(Int_t i=0; i<2; i++) | |
535 | { | |
536 | fhNTrials[i] = new TH1F(Form("MC_%s_fhNTrials",triggerName[i]),Form("MC-%s: # of trials",triggerName[i]),1,0,1); | |
537 | fOutputList->Add(fhNTrials[i]); | |
538 | ||
539 | fVertexGenZ[i] = new TH1F(Form("%s_fVertexGenZ",vertexName[i]),Form("Distribution of vertex z (%s); z (cm)",vertexType[i]),60,-30,30); | |
540 | fOutputList->Add(fVertexGenZ[i]); | |
541 | } | |
542 | } | |
543 | ||
544 | for(Int_t i=0; i<3; i++) | |
545 | { | |
546 | if(!fIsMC && i==2) continue; | |
547 | ||
548 | if(fSaveQAHistos) | |
549 | { | |
550 | if(i<2) | |
551 | { | |
552 | fEventZ[i] = new TH1F(Form("%s_fEventZ",triggerName[i]),Form("%s: Distribution of vertex z; z (cm)",triggerTitle[i]),60,-30,30); | |
553 | fOutputList->Add(fEventZ[i]); | |
554 | } | |
555 | ||
7c915297 | 556 | for(Int_t r=0; r<kNRadii; r++) |
b6e50056 | 557 | { |
7c915297 | 558 | if(!fRadius.Contains("0.4") && r==0) continue; |
559 | if(!fRadius.Contains("0.2") && r==1) continue; | |
560 | if(!fRadius.Contains("0.3") && r==2) continue; | |
561 | ||
b6e50056 | 562 | fJetCount[i][r] = new TH1F(Form("%s_fJetCount_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins); |
563 | fOutputList->Add(fJetCount[i][r]); | |
564 | ||
7c915297 | 565 | fhNeutralPtInJet[i][r] = new TH2F(Form("%s_fhNeutralPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of neutral constituents vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,nTrkPtBins*10,lowTrkPtBin,upTrkPtBin); |
b6e50056 | 566 | fOutputList->Add(fhNeutralPtInJet[i][r]); |
7c915297 | 567 | |
568 | fhTrigNeuPtInJet[i][r] = new TH2F(Form("%s_fhTrigNeuPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of triggered neutral constituents vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,nTrkPtBins*10,lowTrkPtBin,upTrkPtBin); | |
569 | fOutputList->Add(fhTrigNeuPtInJet[i][r]); | |
b6e50056 | 570 | |
7c915297 | 571 | fhChargedPtInJet[i][r] = new TH2F(Form("%s_fhChargedPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of charged constituents vs jet p_{T} in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,nTrkPtBins*10,lowTrkPtBin,upTrkPtBin); |
b6e50056 | 572 | fOutputList->Add(fhChargedPtInJet[i][r]); |
573 | ||
574 | fhLeadNePtInJet[i][r] = new TH2F(Form("%s_fhLeadNePtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of leading neutral constituent vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T}^{ne} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,100); | |
575 | fOutputList->Add(fhLeadNePtInJet[i][r]); | |
576 | ||
577 | fhLeadChPtInJet[i][r] = new TH2F(Form("%s_fhLeadChPtInJet_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of leading charged constituent vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); p_{T}^{ch} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,100); | |
578 | fOutputList->Add(fhLeadChPtInJet[i][r]); | |
579 | ||
580 | fhJetPtVsZ[i][r] = new TH3F(Form("%s_fhJetPtVsZ_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} vs Z_{h} vs constituent type in EMCal (R=%1.1f);p_{T}^{jet} (GeV/c); Z_{h};constituent type",triggerTitle[i],kRadius[r]),200,0,200,110,-0.05,1.05,4,0,4); | |
581 | fOutputList->Add(fhJetPtVsZ[i][r]); | |
582 | ||
583 | fJetEnergyFraction[i][r] = new THnSparseF(Form("%s_fJetEnergyFraction_%1.1f",triggerName[i],kRadius[r]),Form("%s: Jet p_{T} vs radius vs energy fraction in EMCal (R=%1.1f,Z<0.98); p_{T};R;Fraction",triggerName[i],kRadius[r]),dim1,nBins1,lowBin1,upBin1); | |
584 | fOutputList->Add(fJetEnergyFraction[i][r]); | |
585 | ||
586 | fJetNPartFraction[i][r] = new THnSparseF(Form("%s_fJetNPartFraction_%1.1f",triggerName[i],kRadius[r]),Form("%s: Jet p_{T} vs radius vs NPart in EMCal (R=%1.1f,Z<0.98); p_{T};R;NPart",triggerName[i],kRadius[r]),dim2,nBins2,lowBin2,upBin2); | |
587 | fOutputList->Add(fJetNPartFraction[i][r]); | |
588 | ||
589 | if(i<2) | |
590 | { | |
591 | fhJetPtInExoticEvent[i][r] = new TH1F(Form("EMC_fhJetPtInExoticEvent_%1.1f_%s",kRadius[r],exotic[i]),Form("EMC: jet p_{T} in events with exotic cluster > %s (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c)",exotic[i],kRadius[r]),nbins,xbins); | |
592 | fOutputList->Add(fhJetPtInExoticEvent[i][r]); | |
593 | ||
594 | fRelTrkCon[i][r] = new TH3F(Form("%s_fRelTrkCon_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} vs (sum p_{T}^{ch})/p_{T}^{jet} vs track class in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); (sum p_{T}^{ch})/p_{T}^{jet}; track class",triggerTitle[i],kRadius[r]),200,0,200,110,-0.05,1.05,3,0,3); | |
595 | fOutputList->Add(fRelTrkCon[i][r]); | |
596 | ||
597 | fhFcrossVsZleading[i][r] = new TH3F(Form("%s_fhFcrossVsZleading_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet p_{T} vs F_{cross} vs Z_{leading}^{ne} in EMCal (R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c);F_{cross};Z_{leading}^{ne}",triggerTitle[i],kRadius[r]),200,0,200,55,-0.05,1.05,55,-0.05,1.05); | |
598 | fOutputList->Add(fhFcrossVsZleading[i][r]); | |
599 | ||
600 | fhChLeadZVsJetPt[i][r] = new TH2F(Form("%s_fhChLeadZVsJetPt_%1.1f",triggerName[i],kRadius[r]),Form("%s: Z of leading charged constituent vs jet p_{T} in EMCal(R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c); Z_{leading}^{ch}",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,1); | |
601 | fOutputList->Add(fhChLeadZVsJetPt[i][r]); | |
602 | ||
603 | fhJetPtWithTrkThres[i][r] = new TH2F(Form("%s_fhJetPtWithTrkThres_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of jets containing tracks above certain threshold (15,25,40 GeV/c) (EMCal, R=%1.1f,Z<0.98);Threshold type;p_{T}^{jet} (GeV/c)",triggerTitle[i],kRadius[r]),3,0,3,nbins,xbins); | |
604 | fOutputList->Add(fhJetPtWithTrkThres[i][r]); | |
605 | ||
606 | fhJetPtWithClsThres[i][r] = new TH2F(Form("%s_fhJetPtWithClsThres_%1.1f",triggerName[i],kRadius[r]),Form("%s: p_{T} of jets containing clusters above certain threshold (15,25,40 GeV) (EMCal, R=%1.1f,Z<0.98);Threshold type;p_{T}^{jet} (GeV/c)",triggerTitle[i],kRadius[r]),3,0,3,nbins,xbins); | |
607 | fOutputList->Add(fhJetPtWithClsThres[i][r]); | |
608 | ||
7c915297 | 609 | fhJetPtVsLowPtCons[i][r][0] = new TH2F(Form("%s_fhJetPtVsLowPtCons_150-300MeV_%1.1f",triggerName[i],kRadius[r]),Form("%s: energy carried by constituents in 150-300MeV (EMCal, R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c);p_{T,em}^{low} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,1); |
610 | fOutputList->Add(fhJetPtVsLowPtCons[i][r][0]); | |
611 | ||
612 | fhJetPtVsLowPtCons[i][r][1] = new TH2F(Form("%s_fhJetPtVsLowPtCons_300-500MeV_%1.1f",triggerName[i],kRadius[r]),Form("%s: energy carried by constituents in 300-500MeV (EMCal, R=%1.1f,Z<0.98);p_{T}^{jet} (GeV/c);p_{T,em}^{low} (GeV/c)",triggerTitle[i],kRadius[r]),nbins,xbins,100,0,1); | |
613 | fOutputList->Add(fhJetPtVsLowPtCons[i][r][1]); | |
614 | ||
615 | if(fCheckTPCOnlyVtx) | |
616 | { | |
617 | fhJetInTPCOnlyVtx[i][r] = new TH3F(Form("%s_fhJetInTPCOnlyVtx_%1.1f",triggerName[i],kRadius[r]),Form("%s: jet pt in events with only TPC vertex (Full, R=%1.1f);p_{T}^{jet} (GeV/c);#phi;#eta",triggerTitle[i],kRadius[r]),20,0,100,36,0,360,20,-1,1); | |
618 | fOutputList->Add(fhJetInTPCOnlyVtx[i][r]); | |
619 | } | |
b6e50056 | 620 | |
621 | if(fStudySubEInHC) | |
622 | { | |
623 | for(Int_t k=0; k<5; k++) | |
624 | { | |
625 | fhSubClsEVsJetPt[i][r][k] = new TH2F(Form("%s_fhSubClsEVsJetPt_%s_%1.1f",triggerName[i],fraction[k],kRadius[r]),Form("%s: relative %s%% subtracted cluster Et vs jet pt (R=%1.1f);jet p_{T} (GeV/c);E_{t}^{sub}/p_{T,jet}",triggerTitle[i],fraction[k], kRadius[r]),nbins,xbins,50,0,0.5); | |
626 | fOutputList->Add(fhSubClsEVsJetPt[i][r][k]); | |
627 | ||
628 | fhHCTrkPtClean[i][r][k] = new TH2F(Form("%s_fhHCTrkPtClean_%s_%1.1f",triggerName[i],fraction[k],kRadius[r]),Form("%s: sum of track p_{T} that are cleanly subtracted vs jet pt (%s%%, R=%1.1f);jet p_{T} (GeV/c);#sum(p_{T,trk}^{clean})/p_{T,jet}",triggerTitle[i],fraction[k], kRadius[r]),nbins,xbins,100,-0.005,0.995); | |
629 | fOutputList->Add(fhHCTrkPtClean[i][r][k]); | |
630 | ||
631 | fhHCTrkPtAmbig[i][r][k] = new TH2F(Form("%s_fhHCTrkPtAmbig_%s_%1.1f",triggerName[i],fraction[k],kRadius[r]),Form("%s: sum of track p_{T} that are ambiguously subtracted vs jet pt (%s%%, R=%1.1f);jet p_{T} (GeV/c);#sum(p_{T,trk}^{ambig})/p_{T,jet}",triggerTitle[i],fraction[k], kRadius[r]),nbins,xbins,100,-0.005,0.995); | |
632 | fOutputList->Add(fhHCTrkPtAmbig[i][r][k]); | |
633 | } | |
634 | } | |
635 | } | |
636 | } | |
637 | if(i<2) | |
638 | { | |
639 | fhNMatchedTrack[i] = new TH1F(Form("%s_fhNMatchedTrack",triggerName[i]),Form("%s: # of matched tracks per cluster; N_{mth}",triggerTitle[i]),5,0,5); | |
640 | fOutputList->Add(fhNMatchedTrack[i]); | |
641 | ||
642 | for(Int_t j=0; j<4; j++) | |
643 | { | |
644 | fhSubEVsTrkPt[i][j] = new TH2F(Form("%s_fhSubEVsTrkPt_%s",triggerName[i],fraction[j+1]),Form("%s: fractional subtracted energy (%s%% HC);#sum(p_{ch,T}^{mth}) (GeV/c);E_{sub}/#sum(P_{ch}^{mth})",triggerTitle[i],fraction[j+1]),50,0,50,110,0,1.1); | |
645 | fOutputList->Add(fhSubEVsTrkPt[i][j]); | |
646 | } | |
647 | } | |
648 | } | |
7c915297 | 649 | |
650 | if(fRunUE) | |
b6e50056 | 651 | { |
7c915297 | 652 | for(Int_t k=0; k<2; k++) |
653 | { | |
654 | for(Int_t l=0; l<2; l++) | |
655 | { | |
656 | fhUEJetPtNorm[i][k][l] = new TH1F(Form("%s_fhUEJetPtNorm_%s_%s",triggerName[i],UEName[k],UEEventName[l]),Form("%s: leading jet p_{T} in TPC (%s in %s event);p_{T,jet}^{ch} (GeV/c)",triggerTitle[i],UETitle[k], UEEventName[l]),nbins,xbins); | |
657 | fOutputList->Add(fhUEJetPtNorm[i][k][l]); | |
658 | ||
659 | fhUEJetPtVsSumPt[i][k][l] = new TH2F(Form("%s_fhUEJetPtVsSumPt_%s_%s",triggerName[i],UEName[k],UEEventName[l]),Form("%s: leading jet p_{T} vs underlying event contribution (R=0.4,%s in %s event);p_{T,jet}^{ch} (GeV/c);p_{T,UE}^{ch} (GeV/c)",triggerTitle[i],UETitle[k], UEEventName[l]),nbins,xbins,500,0,50); | |
660 | fOutputList->Add(fhUEJetPtVsSumPt[i][k][l]); | |
661 | ||
662 | fhUEJetPtVsConsPt[i][k][l] = new TH2F(Form("%s_fhUEJetPtVsConsPt_%s_%s",triggerName[i],UEName[k],UEEventName[l]),Form("%s: leading jet p_{T} vs constituent pt in UE (R=0.4,%s in %s event);p_{T,jet}^{ch} (GeV/c);p_{T,cons}^{ch} (GeV/c)",triggerTitle[i],UETitle[k], UEEventName[l]),nbins,xbins,500,0,50); | |
663 | fOutputList->Add(fhUEJetPtVsConsPt[i][k][l]); | |
664 | } | |
665 | } | |
b6e50056 | 666 | } |
7c915297 | 667 | |
b6e50056 | 668 | if(fSysJetTrigEff) |
669 | { | |
670 | fhClsE[i] = new TH1F(Form("%s_fhClsE",triggerName[i]),Form("%s: cluster E;E (GeV)",triggerTitle[i]),1000,0,100); | |
671 | fOutputList->Add(fhClsE[i]); | |
672 | } | |
7c915297 | 673 | |
674 | if(fSysClusterizer) | |
675 | { | |
676 | for(Int_t k=0; k<2; k++) | |
677 | { | |
678 | fhSysClusterE[i][k] = new TH1F(Form("%s_fhSysClusterE_%sHC",triggerName[i],clusterizerName[k]),Form("%s: cluster E %s hadronic correction;E (GeV)",triggerTitle[i],clusterizerName[k]),100,0,100); | |
679 | fOutputList->Add(fhSysClusterE[i][k]); | |
680 | ||
681 | fhSysNCellVsClsE[i][k] = new TH2F(Form("%s_fhSysNCellVsClsE_%sHC",triggerName[i],clusterizerName[k]),Form("%s: NCell vs cluster E %s hadronic correction;E (GeV);NCell",triggerTitle[i],clusterizerName[k]),100,0,100,50,0,50); | |
682 | fOutputList->Add(fhSysNCellVsClsE[i][k]); | |
683 | } | |
684 | } | |
b6e50056 | 685 | } |
686 | ||
687 | ||
688 | if(fIsMC) | |
689 | { | |
690 | if(fStudyMcOverSubE) | |
691 | { | |
7c915297 | 692 | for(Int_t r=0; r<kNRadii; r++) |
b6e50056 | 693 | { |
7c915297 | 694 | if(!fRadius.Contains("0.4") && r==0) continue; |
695 | if(!fRadius.Contains("0.2") && r==1) continue; | |
696 | if(!fRadius.Contains("0.3") && r==2) continue; | |
b6e50056 | 697 | for(Int_t i=0; i<5; i++) |
698 | { | |
699 | fHCOverSubE[r][i] = new TH2F(Form("%s_HC_over_sub_e_%s_%1.1f",triggerName[2],fraction[i],kRadius[r]),Form("%s: oversubtracted neutral Et by %s%% HC (R=%1.1f);particle jet p_{T} (GeV/c);#DeltaE_{t} (GeV)",triggerName[2],fraction[i],kRadius[r]),nbins,xbins,200,-49.75,50.25); | |
700 | fOutputList->Add(fHCOverSubE[r][i]); | |
701 | fHCOverSubEFrac[r][i] = new TH2F(Form("%s_HC_over_sub_e_frac_%s_%1.1f",triggerName[2],fraction[i],kRadius[r]),Form("%s: relative oversubtracted neutral Et fraction by %s%% HC (R=%1.1f);jet p_{T} (GeV/c);#DeltaE_{t}/p_{T}^{jet}",triggerName[2],fraction[i],kRadius[r]),nbins,xbins,200,-0.995,1.005); | |
702 | fOutputList->Add(fHCOverSubEFrac[r][i]); | |
703 | } | |
704 | } | |
705 | } | |
7c915297 | 706 | if(fRunSecondaries && fAnaType==0) |
707 | { | |
708 | for(Int_t i=0; i<2; i++) | |
709 | { | |
710 | for(Int_t r=0; r<3; r++) | |
711 | { | |
712 | fhNKFracVsJetPt[i][r] = new TH2F(Form("%s_fhNKFracVsJetPt_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: energy fraction carried by n,k^{0}_{L} vs jet p_{T} (R=%1.1f,|#eta|<%1.1f);p_{T,jet} (GeV/c);fraction",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,0,1); | |
713 | fOutputList->Add(fhNKFracVsJetPt[i][r]); | |
714 | ||
715 | fhWeakFracVsJetPt[i][r] = new TH2F(Form("%s_fhWeakFracVsJetPt_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: energy fraction carried by k^{0}_{S},hyperon vs jet p_{T} (R=%1.1f,|#eta|<%1.1f);p_{T,jet} (GeV/c);fraction",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,0,1); | |
716 | fOutputList->Add(fhWeakFracVsJetPt[i][r]); | |
717 | ||
718 | fhJetResponseNK[i][r] = new TH2F(Form("%s_fhJetResponseNK_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to missing n and k^{0}_{L} (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins); | |
719 | fOutputList->Add(fhJetResponseNK[i][r]); | |
720 | ||
721 | fhJetResponseWP[i][r] = new TH2F(Form("%s_fhJetResponseWP_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to k^{0}_{S}, #Lambda and hyperon (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins); | |
722 | fOutputList->Add(fhJetResponseWP[i][r]); | |
723 | ||
724 | fhJetResolutionNK[i][r] = new TH2F(Form("%s_fhJetResolutionNK_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to missing n and k^{0}_{L}: (p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{rec} vs p_{T,jet}^{rec} (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T}/p_{T}",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,-0.995,1.005); | |
725 | fOutputList->Add(fhJetResolutionNK[i][r]); | |
726 | ||
727 | fhJetResolutionWP[i][r] = new TH2F(Form("%s_fhJetResolutionWP_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to k^{0}_{S}, #Lambda and hyperon: (p_{T,jet}^{rec}-p_{T,jet}^{gen})/p_{T,jet}^{rec} vs p_{T,jet}^{rec} (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T}/p_{T}",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,200,-0.995,1.005); | |
728 | fOutputList->Add(fhJetResolutionWP[i][r]); | |
729 | ||
730 | fhJetResponseNKSM[i][r] = new TH2F(Form("%s_fhJetResponseNKSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to missing n and k^{0}_{L} via matching (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins); | |
731 | fOutputList->Add(fhJetResponseNKSM[i][r]); | |
732 | ||
733 | fhJetResponseWPSM[i][r] = new TH2F(Form("%s_fhJetResponseWPSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet response due to k^{0}_{S}, #Lambda and hyperon via matching(R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{gen} (GeV/c);p_{T,jet}^{rec} (GeV/c)",triggerName[2],kRadius[r],i*0.5+0.5),nbins,xbins,nbins,xbins); | |
734 | fOutputList->Add(fhJetResponseWPSM[i][r]); | |
735 | ||
736 | fhJetResolutionNKSM[i][r] = new TH3F(Form("%s_fhJetResolutionNKSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet resolution due to missing n and k^{0}_{L} via matching (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T,jet}/p_{T,jet}^{rec};dR",triggerName[2],kRadius[r],i*0.5+0.5),220,0,220,200,-0.995,1.005,100,0,1); | |
737 | fOutputList->Add(fhJetResolutionNKSM[i][r]); | |
738 | ||
739 | fhJetResolutionWPSM[i][r] = new TH3F(Form("%s_fhJetResolutionWPSM_%1.1f_EtaCut%1.1f",triggerName[2],kRadius[r],i*0.5+0.5),Form("%s: jet resolution due tok^{0}_{S}, #Lambda and hyperon via matching (R=%1.1f,|#eta|<%1.1f);p_{T,jet}^{rec} (GeV/c);#Deltap_{T,jet}/p_{T,jet}^{rec};dR",triggerName[2],kRadius[r],i*0.5+0.5),220,0,220,200,-0.995,1.005,100,0,1); | |
740 | fOutputList->Add(fhJetResolutionWPSM[i][r]); | |
741 | } | |
742 | } | |
743 | } | |
b6e50056 | 744 | } |
745 | ||
746 | printf("\n=======================================\n"); | |
747 | printf("===== Jet task configuration ==========\n"); | |
7c915297 | 748 | |
b6e50056 | 749 | if(fNonStdBranch.Length()!=0) |
750 | { | |
751 | const char* species[3] = {"in","ch","ne"}; | |
752 | const char* algorithm[2] = {"akt","kt"}; | |
7c915297 | 753 | const char* radii[kNRadii] = {"04","02","03"}; |
b6e50056 | 754 | for(Int_t s=0; s<3; s++) |
755 | { | |
756 | if(!fFindChargedOnlyJet && s==1) continue; | |
757 | if(!fFindNeutralOnlyJet && s==2) continue; | |
758 | for(Int_t a=0; a<2; a++) | |
759 | { | |
760 | if(!fAlgorithm.Contains("aKt") && a==0) continue; | |
761 | if(!fAlgorithm.Contains("kt") && a==1) continue; | |
7c915297 | 762 | for(Int_t r=0; r<kNRadii; r++) |
b6e50056 | 763 | { |
764 | if(!fRadius.Contains("0.4") && r==0) continue; | |
765 | if(!fRadius.Contains("0.2") && r==1) continue; | |
7c915297 | 766 | if(!fRadius.Contains("0.3") && r==2) continue; |
767 | if(fAnaType==0) | |
768 | { | |
769 | fJetTCA[s][a][r] = new TClonesArray("AliAODJet",0); | |
770 | fJetTCA[s][a][r]->SetName(Form("Jet_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data())); | |
771 | AddAODBranch("TClonesArray",&fJetTCA[s][a][r],fNonStdFile.Data()); | |
772 | printf("Add branch: Jet_%s_%s_%s_%s\n",species[s],algorithm[a],radii[r],fNonStdBranch.Data()); | |
773 | } | |
b6e50056 | 774 | |
775 | fDetJetFinder[s][a][r] = new AliFJWrapper(Form("DetJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()),Form("DetJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data())); | |
776 | if(a==0) fDetJetFinder[s][a][r]->SetAlgorithm(fastjet::antikt_algorithm); | |
777 | if(a==1) fDetJetFinder[s][a][r]->SetAlgorithm(fastjet::kt_algorithm); | |
778 | if(fRecombinationScheme==0) fDetJetFinder[s][a][r]->SetRecombScheme(fastjet::E_scheme); | |
779 | fDetJetFinder[s][a][r]->SetR(kRadius[r]); | |
780 | fDetJetFinder[s][a][r]->SetMaxRap(0.9); | |
781 | fDetJetFinder[s][a][r]->Clear(); | |
782 | ||
783 | if(s==0 && a==0) | |
784 | { | |
785 | if(fIsMC && fNonStdBranch.Contains("Baseline",TString::kIgnoreCase)) | |
786 | { | |
7c915297 | 787 | if(fAnaType==0) |
788 | { | |
789 | fMcTruthAntikt[r] = new TClonesArray("AliAODJet",0); | |
790 | fMcTruthAntikt[r]->SetName(Form("Jet_mc_truth_in_akt_%s_%s",radii[r],fNonStdBranch.Data())); | |
791 | AddAODBranch("TClonesArray",&fMcTruthAntikt[r],fNonStdFile.Data()); | |
792 | printf("Add branch: Jet_mc_truth_in_akt_%s_%s\n",radii[r],fNonStdBranch.Data()); | |
793 | } | |
b6e50056 | 794 | |
795 | fTrueJetFinder[r] = new AliFJWrapper(Form("TrueJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data()),Form("TrueJetFinder_%s_%s_%s_%s",species[s],algorithm[a],radii[r],fNonStdBranch.Data())); | |
796 | fTrueJetFinder[r]->SetAlgorithm(fastjet::antikt_algorithm); | |
797 | fTrueJetFinder[r]->SetR(kRadius[r]); | |
798 | fTrueJetFinder[r]->SetMaxRap(0.9); | |
799 | if(fRecombinationScheme==0) fTrueJetFinder[r]->SetRecombScheme(fastjet::E_scheme); | |
800 | fTrueJetFinder[r]->Clear(); | |
801 | } | |
802 | } | |
803 | } | |
804 | } | |
805 | } | |
806 | } | |
807 | ||
808 | fRandomGen = new TRandom3(0); | |
7c915297 | 809 | if(fCheckTrkEffCorr && fAnaType==0) |
b6e50056 | 810 | { |
811 | TFile f ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TrkEffFit.root","read"); | |
812 | for(Int_t i=0; i<3; i++) | |
813 | { | |
814 | fTrkEffFunc[i] = new TF1(*((TF1*)f.Get(Form("Trk_eff_fit_%d",i+1)))); | |
815 | } | |
816 | f.Close(); | |
817 | ||
818 | if(fPeriod.CompareTo("lhc11a",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc11aJet",TString::kIgnoreCase)==0) | |
819 | { | |
820 | TFile f1 ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TrkEffSampling.root","read"); | |
821 | for(Int_t j=0; j<2; j++) | |
822 | { | |
823 | Int_t tmp = j; | |
824 | if(fPeriod.CompareTo("lhc11aJet",TString::kIgnoreCase)==0) tmp = 2; | |
5543d2a2 | 825 | for(Int_t r=0; r<kNRadii; r++) |
b6e50056 | 826 | { |
827 | fhCorrTrkEffPtBin[j][r] = new TH1F(*((TH1F*)f1.Get(Form("%s_%s_NTrackPerPtBin_%1.1f",fPeriod.Data(),triggerName[tmp],kRadius[r])))); | |
828 | for(Int_t i=0; i<kNBins; i++) | |
829 | { | |
830 | fhCorrTrkEffSample[j][r][i] = new TH1F(*((TH1F*)f1.Get(Form("%s_%s_ChTrkPt_Bin%d_%1.1f",fPeriod.Data(),triggerName[tmp],i+1,kRadius[r])))); | |
831 | } | |
832 | } | |
833 | } | |
834 | f1.Close(); | |
835 | } | |
836 | } | |
837 | ||
7c915297 | 838 | if(fRunSecondaries && fAnaType==0) |
839 | { | |
840 | const char *secondaryName[3] = {"k0S","lamda","hyperon"}; | |
841 | TFile f2 ("/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/SecondaryResponse.root","read"); | |
842 | for(Int_t j=0; j<3; j++) | |
843 | { | |
844 | fhSecondaryResponse[j] = new TH2F(*((TH2F*)f2.Get(Form("DetectorReponse_%s",secondaryName[j])))); | |
845 | } | |
846 | } | |
b6e50056 | 847 | |
7c915297 | 848 | if(fCheckTriggerMask) |
b6e50056 | 849 | { |
5543d2a2 | 850 | TString name = "TriggerCurve.root"; |
7c915297 | 851 | if(fAnaType==0) name = "/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TriggerCurve.root"; |
5543d2a2 | 852 | TFile f (name.Data(),"read"); |
7c915297 | 853 | fTriggerMask = new TH2F(*((TH2F*)f.Get("lhc11a_TriggerMask"))); |
854 | if(fOfflineTrigger) | |
b6e50056 | 855 | { |
7c915297 | 856 | for(Int_t i=0; i<10; i++) |
857 | { | |
858 | fTriggerEfficiency[i] = new TF1(*((TF1*)f.Get(Form("lhc11a_TriggerEfficiency_SM%d_fit",i)))); | |
859 | fTriggerCurve[i] = new TH1D(*((TH1D*)f.Get(Form("lhc11a_TriggerCurve_SM%d",i)))); | |
860 | } | |
861 | } | |
862 | f.Close(); | |
b6e50056 | 863 | } |
b6e50056 | 864 | |
865 | fClusterEResolution = new TF1("fClusterEResolution","sqrt([0]^2+[1]^2*x+([2]*x)^2)*0.01"); | |
866 | fClusterEResolution->SetParameters(4.35,9.07,1.63); | |
867 | ||
b6e50056 | 868 | fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1"); |
7c915297 | 869 | if(!fGeom) |
870 | { | |
871 | AliError("No EMCal geometry is available!"); | |
872 | return; | |
873 | } | |
b6e50056 | 874 | fRecoUtil = new AliEMCALRecoUtils(); |
7c915297 | 875 | if(fRejectExoticCluster || fRejectExoticTrigger) |
b6e50056 | 876 | fRecoUtil->SwitchOnRejectExoticCluster(); |
877 | else | |
878 | { | |
879 | fRecoUtil->SwitchOffRejectExoticCluster(); | |
880 | fRecoUtil->SwitchOffRejectExoticCell(); | |
881 | } | |
882 | if(fRemoveBadChannel) | |
883 | { | |
884 | fRecoUtil->SwitchOnBadChannelsRemoval(); | |
885 | // Remove the problematic SM4 region due to wrong event sequence | |
886 | for(Int_t ieta=0; ieta<36; ieta++) | |
887 | for(Int_t iphi=0; iphi<8; iphi++) | |
888 | fRecoUtil->SetEMCALChannelStatus(4,ieta,iphi); | |
889 | } | |
890 | else | |
891 | fRecoUtil->SwitchOffBadChannelsRemoval(); | |
892 | ||
893 | fRecoUtil->SetNonLinearityFunction(AliEMCALRecoUtils::kBeamTestCorrected); | |
894 | if(fSysNonLinearity) | |
895 | { | |
7c915297 | 896 | fNonLinear = new TF1("TB_oldBest","([0])*(1./(1.+[1]*exp(-x/[2]))*1./(1.+[3]*exp((x-[4])/[5])))*1/([6])",0.1,110.); |
897 | fNonLinear->SetParameters(0.99078, 0.161499, 0.655166, 0.134101, 163.282, 23.6904, 0.978); | |
b6e50056 | 898 | } |
899 | if(fSysTrkClsMth) | |
900 | { | |
901 | fRecoUtil->SwitchOnCutEtaPhiSeparate(); | |
902 | fRecoUtil->SetCutEta(fCutdEta); | |
903 | fRecoUtil->SetCutPhi(fCutdPhi); | |
904 | } | |
905 | ||
906 | fTrackArray = new TObjArray(); | |
907 | fTrackArray->SetOwner(1); | |
908 | fClusterArray = new TObjArray(); | |
909 | fClusterArray->SetOwner(1); | |
910 | fMcPartArray = new TArrayI(); | |
911 | ||
912 | ||
913 | //error calculation in THnSparse | |
914 | Int_t nObj = fOutputList->GetEntries(); | |
915 | for(Int_t i=0; i<nObj; i++) | |
916 | { | |
917 | TObject *obj = (TObject*) fOutputList->At(i); | |
918 | if (obj->IsA()->InheritsFrom( "THnSparse" )) | |
919 | { | |
920 | THnSparseF *hn = (THnSparseF*)obj; | |
921 | hn->Sumw2(); | |
922 | } | |
923 | } | |
924 | ||
925 | ||
926 | PrintConfig(); | |
927 | BookHistos(); | |
b6e50056 | 928 | PostData(1, fOutputList); |
929 | } | |
930 | ||
931 | ||
932 | ||
933 | //________________________________________________________________________ | |
934 | void AliAnalysisTaskFullppJet::BookHistos() | |
935 | { | |
936 | // Book histograms. | |
937 | ||
938 | if(fVerbosity>10) printf("[i] Booking histograms \n"); | |
939 | return; | |
940 | } | |
941 | ||
942 | //________________________________________________________________________ | |
943 | void AliAnalysisTaskFullppJet::UserExec(Option_t *) | |
944 | { | |
945 | // | |
946 | // Main loop, called for each event. | |
947 | // | |
948 | ||
b6e50056 | 949 | // Get event pointers, check for signs of life |
950 | Double_t vtxTrueZ = -100; | |
951 | if(fIsMC) | |
952 | { | |
953 | fMC = MCEvent(); | |
954 | if (!fMC) | |
955 | { | |
7c915297 | 956 | printf("ERROR: Could not retrieve MC event"); |
b6e50056 | 957 | return; |
958 | } | |
959 | fStack = fMC->Stack(); | |
960 | TParticle *particle = fStack->Particle(0); | |
5543d2a2 | 961 | if(particle) |
962 | { | |
963 | vtxTrueZ = particle->Vz(); // True vertex z in MC events | |
964 | if(fVerbosity>10) printf("Generated vertex coordinate: (x,y,z) = (%4.2f, %4.2f, %4.2f)\n", particle->Vx(), particle->Vy(), particle->Vz()); | |
965 | } | |
b6e50056 | 966 | } |
967 | ||
968 | fESD = dynamic_cast<AliESDEvent*>(InputEvent()); | |
969 | if (!fESD) | |
970 | { | |
5543d2a2 | 971 | AliError("fESD is not available"); |
b6e50056 | 972 | return; |
973 | } | |
974 | ||
975 | fhJetEventStat->Fill(0.5); | |
7c915297 | 976 | fhJetEventCount->Fill(0.5); |
977 | if(fIsMC) | |
978 | { | |
979 | GetMCInfo(); | |
980 | if(fVertexGenZ[0]) fVertexGenZ[0]->Fill(vtxTrueZ); | |
981 | } | |
982 | ||
b6e50056 | 983 | // Centrality, vertex, other event variables... |
5543d2a2 | 984 | UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected(); |
985 | if (trigger==0) return; | |
986 | if(fCheckTPCOnlyVtx) CheckTPCOnlyVtx(trigger); | |
987 | if(!fIsMC) | |
b6e50056 | 988 | { |
5543d2a2 | 989 | if (trigger & AliVEvent::kFastOnly) return; // Reject fast trigger cluster |
990 | if(fPeriod.CompareTo("lhc11a",TString::kIgnoreCase)==0) | |
b6e50056 | 991 | { |
5543d2a2 | 992 | if (trigger & AliVEvent::kMB) fTriggerType = 0; |
993 | else if(trigger & AliVEvent::kEMC1) fTriggerType = 1; | |
994 | else fTriggerType = -1; | |
b6e50056 | 995 | } |
5543d2a2 | 996 | else if (fPeriod.CompareTo("lhc11c",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc11d",TString::kIgnoreCase)==0) |
b6e50056 | 997 | { |
5543d2a2 | 998 | if (trigger & AliVEvent::kINT7) fTriggerType = 0; |
999 | else if(trigger & AliVEvent::kEMC7) fTriggerType = 1; | |
7c915297 | 1000 | else fTriggerType = -1; |
b6e50056 | 1001 | } |
5543d2a2 | 1002 | else |
b6e50056 | 1003 | { |
b6e50056 | 1004 | return; |
1005 | } | |
1006 | } | |
5543d2a2 | 1007 | else |
1008 | { | |
1009 | if(!fPhySelForMC) fTriggerType = 0; | |
1010 | else if (trigger & AliVEvent::kAnyINT) fTriggerType = 0; | |
1011 | else fTriggerType = -1; | |
1012 | ||
1013 | if(fOfflineTrigger) | |
1014 | { | |
1015 | RunOfflineTrigger(); | |
1016 | if(fIsEventTriggerBit) fTriggerType = 1; | |
1017 | } | |
1018 | } | |
1019 | ||
1020 | if(fTriggerType==-1) | |
1021 | { | |
1022 | if(fVerbosity>10) printf("Error: worng trigger type %s\n",(fESD->GetFiredTriggerClasses()).Data()); | |
1023 | return; | |
1024 | } | |
7c915297 | 1025 | |
1026 | fhJetEventCount->Fill(1.5+fTriggerType*4); | |
1027 | if(fRejectPileup && fESD->IsPileupFromSPD() ) return; // reject pileup | |
1028 | fhJetEventStat->Fill(11.5); | |
1029 | fhJetEventCount->Fill(2.5+fTriggerType*4); | |
1030 | ||
1031 | ||
1032 | fIsTPCOnlyVtx = 0; | |
1033 | // Reject LED events | |
b6e50056 | 1034 | if (IsLEDEvent()) |
1035 | { | |
1036 | fhJetEventStat->Fill(8.5); | |
1037 | return; | |
1038 | } | |
1039 | ||
1040 | fhJetEventStat->Fill(1.5+fTriggerType*2); | |
1041 | ||
7c915297 | 1042 | // Check if primary vertex exists |
1043 | if (!HasPrimaryVertex()) return; | |
1044 | fhJetEventStat->Fill(5.5+fTriggerType); | |
1045 | fhJetEventCount->Fill(3.5+fTriggerType*4); | |
1046 | ||
1047 | const AliESDVertex* vtx = fESD->GetPrimaryVertex(); | |
1048 | Double_t zVertex = vtx->GetZ(); | |
1049 | if(fEventZ[fTriggerType]) fEventZ[fTriggerType]->Fill(zVertex); | |
1050 | if(fVertexGenZ[1]) fVertexGenZ[1]->Fill(vtxTrueZ); | |
1051 | ||
1052 | // Check if |Z_vtx|<10cm | |
1053 | if( TMath::Abs(zVertex) > fZVtxMax ) return; | |
1054 | fhJetEventCount->Fill(4.5+fTriggerType*4); | |
1055 | ||
1056 | // Check if only TPC vertex exists primitive | |
b6e50056 | 1057 | fIsTPCOnlyVtx = IsTPCOnlyVtx(); |
1058 | ||
1059 | if(!fIsMC && fTriggerType==1) | |
1060 | { | |
7c915297 | 1061 | // Check if event has valid trigger bit |
b6e50056 | 1062 | CheckEventTriggerBit(); |
7c915297 | 1063 | if(fIsEventTriggerBit) |
1064 | { | |
1065 | fhJetEventStat->Fill(7.5); | |
1066 | fhJetEventCount->Fill(9.5); | |
1067 | } | |
b6e50056 | 1068 | else return; |
1069 | } | |
1070 | fhJetEventStat->Fill(2.5+fTriggerType*2); | |
7c915297 | 1071 | if(fIsTPCOnlyVtx) fhJetEventStat->Fill(9.5+fTriggerType); |
b6e50056 | 1072 | |
7c915297 | 1073 | // Check if event contains exotic clusters |
b6e50056 | 1074 | CheckExoticEvent(); |
1075 | ||
7c915297 | 1076 | // Write jet tree into AOD output in local analysis mode |
1077 | if(fAnaType==0) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE); | |
b6e50056 | 1078 | |
7c915297 | 1079 | // Clean up arrays |
b6e50056 | 1080 | for(Int_t s=0; s<3; s++) |
1081 | { | |
1082 | for(Int_t a=0; a<2; a++) | |
1083 | { | |
7c915297 | 1084 | for(Int_t r=0; r<kNRadii; r++) |
b6e50056 | 1085 | { |
7c915297 | 1086 | if(fJetTCA[s][a][r]) fJetTCA[s][a][r]->Delete(); |
1087 | if(fDetJetFinder[s][a][r]) fDetJetFinder[s][a][r]->Clear(); | |
b6e50056 | 1088 | |
1089 | if(s==0 && a==0) | |
1090 | { | |
7c915297 | 1091 | if(fMcTruthAntikt[r]) fMcTruthAntikt[r]->Delete(); |
1092 | if(fTrueJetFinder[r]) fTrueJetFinder[r]->Clear(); | |
b6e50056 | 1093 | } |
1094 | } | |
1095 | } | |
1096 | } | |
1097 | ||
7c915297 | 1098 | if(fVerbosity>5) printf("# of jets after clear: %d\n",fJetTCA[0][0][0]->GetEntries()); |
b6e50056 | 1099 | |
1100 | if(fTrackArray) fTrackArray->Delete(); | |
1101 | if(fClusterArray) fClusterArray->Delete(); | |
1102 | fMcPartArray->Reset(); | |
7c915297 | 1103 | |
5543d2a2 | 1104 | |
1105 | // get the tracks and fill the input vector for the jet finders | |
1106 | GetESDTrax(); | |
1107 | ||
1108 | // get EMCal clusters and fill the input vector for the jet finders | |
1109 | GetESDEMCalClusters(); | |
1110 | ||
1111 | for(Int_t s=0; s<3; s++) | |
1112 | { | |
1113 | for(Int_t a=0; a<2; a++) | |
1114 | { | |
1115 | for(Int_t r=0; r<kNRadii; r++) | |
1116 | { | |
1117 | //Detector jets | |
1118 | if(fDetJetFinder[s][a][r]) | |
1119 | { | |
1120 | FindDetJets(s,a,r); | |
1121 | if(fJetTCA[s][a][r]) FillAODJets(fJetTCA[s][a][r], fDetJetFinder[s][a][r], 0); | |
1122 | if(s==0 && a==0 && fNonStdBranch.Contains("Baseline",TString::kIgnoreCase)) | |
1123 | { | |
1124 | if(fSaveQAHistos) AnalyzeJets(fDetJetFinder[s][a][r],fTriggerType, r); // run analysis on jets | |
1125 | } | |
1126 | } | |
1127 | } | |
1128 | } | |
1129 | } | |
1130 | if(fRunUE && fDetJetFinder[1][0][0]) RunAnalyzeUE(fDetJetFinder[1][0][0],fTriggerType,0); // Run analysis on underlying event | |
b6e50056 | 1131 | |
1132 | if(fIsMC) | |
1133 | { | |
7c915297 | 1134 | for(Int_t r=0; r<kNRadii; r++) |
1135 | if(fTrueJetFinder[r]) ProcessMC(r); // find particle level jets | |
b6e50056 | 1136 | } |
5543d2a2 | 1137 | |
b6e50056 | 1138 | // Fast Jet calls END -------------------------------------------------------- |
1139 | ||
1140 | PostData(1, fOutputList); | |
1141 | return; | |
1142 | } | |
1143 | ||
1144 | ||
1145 | //________________________________________________________________________ | |
1146 | void AliAnalysisTaskFullppJet::FindDetJets(const Int_t s, const Int_t a, const Int_t r) | |
1147 | { | |
1148 | // | |
1149 | // Jet finding is happening here | |
1150 | // | |
1151 | ||
1152 | Bool_t isCh = kTRUE; | |
1153 | Bool_t isNe = kTRUE; | |
1154 | if(s==1) isNe = kFALSE; | |
1155 | if(s==2) isCh = kFALSE; | |
1156 | ||
1157 | if(isCh) | |
1158 | { | |
7c915297 | 1159 | // Feed in charged tracks |
b6e50056 | 1160 | Int_t countTracks = 0; |
1161 | Int_t ntracks = fTrackArray->GetEntriesFast(); | |
1162 | for(Int_t it=0; it<ntracks; it++) | |
1163 | { | |
1164 | AliESDtrack *t = (AliESDtrack*) fTrackArray->At(it); | |
7c915297 | 1165 | if(!t) continue; |
b6e50056 | 1166 | countTracks++; |
1167 | Double_t e = t->P(); | |
1168 | Double_t pt = t->Pt(); | |
7c915297 | 1169 | if(s==1 && fRunUE && pt<1) continue; |
b6e50056 | 1170 | if(fRecombinationScheme==0) e = TMath::Sqrt(t->P()*t->P()+0.139*0.139); |
1171 | if(fSysTrkPtRes) pt += GetSmearedTrackPt(t); | |
1172 | Double_t phi = t->Phi(); | |
1173 | Double_t eta = t->Eta(); | |
1174 | Double_t px = pt*TMath::Cos(phi); | |
1175 | Double_t py = pt*TMath::Sin(phi); | |
1176 | if(fConstrainChInEMCal && ( TMath::Abs(eta)>0.7 || phi>kPI || phi<TMath::DegToRad()*80) ) continue; | |
1177 | fDetJetFinder[s][a][r]->AddInputVector(px,py,t->Pz(), e, it+1); | |
7c915297 | 1178 | if(fVerbosity>10) printf("%d: m_{ch}=%f\n",it+1,t->P()*t->P()-t->Px()*t->Px()-t->Py()*t->Py()-t->Pz()*t->Pz()); |
b6e50056 | 1179 | } |
1180 | if(fVerbosity>5) printf("[i] # of tracks filled: %d\n",countTracks); | |
1181 | } | |
1182 | if(isNe) | |
1183 | { | |
7c915297 | 1184 | // Feed in EMCal clusters |
b6e50056 | 1185 | Double_t vertex[3] = {0, 0, 0}; |
1186 | fESD->GetVertex()->GetXYZ(vertex) ; | |
1187 | TLorentzVector gamma; | |
1188 | ||
1189 | Int_t countClusters = 0; | |
1190 | Int_t nclusters = fClusterArray->GetEntriesFast(); | |
1191 | for (Int_t ic = 0; ic < nclusters; ic++) | |
1192 | { | |
1193 | AliESDCaloCluster * cl = (AliESDCaloCluster *)fClusterArray->At(ic); | |
1194 | if(!cl) continue; | |
1195 | cl->GetMomentum(gamma, vertex); | |
1196 | countClusters++; | |
1197 | Double_t e = gamma.P(); | |
1198 | if(fRecombinationScheme==0) e = TMath::Sqrt(gamma.P()*gamma.P()+0.139*0.139); | |
1199 | fDetJetFinder[s][a][r]->AddInputVector(gamma.Px(), gamma.Py(), gamma.Pz(), e, (ic+1)*(-1)); | |
1200 | } | |
1201 | if(fVerbosity>5) printf("[i] # of clusters filled: %d\n",countClusters); | |
1202 | } | |
7c915297 | 1203 | // Run jet finding |
b6e50056 | 1204 | fDetJetFinder[s][a][r]->Run(); |
1205 | } | |
1206 | ||
1207 | ||
1208 | //________________________________________________________________________ | |
1209 | void AliAnalysisTaskFullppJet::FillAODJets(TClonesArray *fJetArray, AliFJWrapper *jetFinder, const Bool_t isTruth) | |
1210 | { | |
1211 | // | |
1212 | // Fill the found jets into AOD output | |
1213 | // Only consider jets pointing to EMCal and with pt above 1GeV/c | |
7c915297 | 1214 | // |
b6e50056 | 1215 | |
1216 | Int_t radiusIndex = 0; | |
1217 | TString arrayName = fJetArray->GetName(); | |
1218 | if(arrayName.Contains("_02_")) radiusIndex = 1; | |
7c915297 | 1219 | if(arrayName.Contains("_03_")) radiusIndex = 2; |
b6e50056 | 1220 | std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets(); |
7c915297 | 1221 | if(fVerbosity>5 && radiusIndex==0) printf("[i] # of jets in %s : %d\n",fJetArray->GetName(),(Int_t)jetsIncl.size()); |
b6e50056 | 1222 | AliAODJet *jet = 0; |
1223 | Int_t jetCount = 0; | |
1224 | for(UInt_t ij=0; ij<jetsIncl.size(); ij++) | |
1225 | { | |
7c915297 | 1226 | if(fVerbosity>10) printf("fastjet: eta=%f, phi=%f\n",jetsIncl[ij].eta(),jetsIncl[ij].phi()*TMath::RadToDeg()); |
b6e50056 | 1227 | if(!IsGoodJet(jetsIncl[ij],0)) continue; |
1228 | ||
1229 | AliAODJet tmpJet (jetsIncl[ij].px(), jetsIncl[ij].py(), jetsIncl[ij].pz(), jetsIncl[ij].E()); | |
1230 | jet = new ((*fJetArray)[jetCount]) AliAODJet(tmpJet); | |
1231 | jetCount++; | |
7c915297 | 1232 | if(fVerbosity>10 && radiusIndex==0) printf("AOD jet: ij=%d, pt=%f, eta=%f, phi=%f\n",ij, jet->Pt(), jet->Eta(),jet->Phi()*TMath::RadToDeg()); |
b6e50056 | 1233 | |
1234 | jet->GetRefTracks()->Clear(); | |
7c915297 | 1235 | std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij); |
b6e50056 | 1236 | Double_t totE=0, totPt=0, totChPt=0, leadChPt=0, neE=0, totNePt=0, leadNePt=0, leadPt=0; |
1237 | Double_t leadTrkType=0, nChPart = 0, nNePart = 0; | |
1238 | Bool_t isHighPtTrigger = kFALSE, isTriggering = kFALSE, isHyperTrack = kFALSE; | |
7c915297 | 1239 | Int_t leadIndex = -1; |
b6e50056 | 1240 | for(UInt_t ic=0; ic<constituents.size(); ic++) |
1241 | { | |
7c915297 | 1242 | if(fVerbosity>10 && radiusIndex==0) printf("ic=%d: user_index=%d, E=%f\n",ic,constituents[ic].user_index(),constituents[ic].E()); |
b6e50056 | 1243 | if(constituents[ic].user_index()<0) |
1244 | { | |
1245 | nNePart ++; | |
1246 | totNePt += constituents[ic].perp(); | |
1247 | if(constituents[ic].perp()>leadNePt) | |
1248 | leadNePt = constituents[ic].perp(); | |
1249 | ||
1250 | neE += constituents[ic].E(); | |
1251 | if(constituents[ic].perp()>fClsEtMax[fTriggerType]) | |
1252 | isHighPtTrigger = kTRUE; | |
1253 | if(!isTruth) | |
1254 | { | |
1255 | AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1); | |
7c915297 | 1256 | if(cluster->Chi2()>0.5) isTriggering = kTRUE; |
b6e50056 | 1257 | } |
1258 | } | |
1259 | else | |
1260 | { | |
1261 | totChPt += constituents[ic].perp(); | |
1262 | nChPart ++; | |
1263 | if(constituents[ic].perp()>leadChPt) | |
1264 | { | |
1265 | leadChPt = constituents[ic].perp(); | |
1266 | if(!isTruth) | |
1267 | { | |
1268 | AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1); | |
1269 | leadTrkType = track->GetTRDQuality(); | |
1270 | if(track->GetIntegratedLength()>500) isHyperTrack = kTRUE; | |
1271 | } | |
1272 | } | |
1273 | if(constituents[ic].perp()>fTrkPtMax[fTriggerType]) | |
1274 | isHighPtTrigger = kTRUE; | |
1275 | } | |
1276 | TParticle part; | |
1277 | part.SetMomentum(constituents[ic].px(),constituents[ic].py(),constituents[ic].pz(),constituents[ic].E()); | |
1278 | jet->AddTrack(&part); //The references are not usable, this line is aimed to count the number of contituents | |
1279 | totE += constituents[ic].E(); | |
1280 | totPt += constituents[ic].perp(); | |
1281 | if(constituents[ic].perp()>leadPt) | |
7c915297 | 1282 | { |
1283 | leadPt = constituents[ic].perp(); | |
1284 | leadIndex = ic; | |
1285 | } | |
b6e50056 | 1286 | } |
1287 | if(fIsEventTriggerBit) jet->SetTrigger(AliAODJet::kTRDTriggered); | |
1288 | if(isHighPtTrigger) jet->SetTrigger(AliAODJet::kHighTrackPtTriggered); | |
1289 | if(fTriggerType==1) jet->SetTrigger(AliAODJet::kEMCALTriggered); | |
1290 | if(GetLeadingZ(ij,jetFinder) > 0.98 ) jet->SetTrigger(AliAnalysisTaskFullppJet::kHighZ); | |
1291 | if(isTriggering) jet->SetTrigger(AliAnalysisTaskFullppJet::kTrigger); | |
1292 | if(isHyperTrack) jet->SetTrigger(AliAnalysisTaskFullppJet::kSuspicious); | |
1293 | if(fIsTPCOnlyVtx) jet->SetTrigger(AliAnalysisTaskFullppJet::kTPCOnlyVtx); | |
7c915297 | 1294 | if(constituents[leadIndex].user_index()>0) jet->SetTrigger(AliAnalysisTaskFullppJet::kLeadCh); |
1295 | ||
b6e50056 | 1296 | if(jetsIncl[ij].E()>0) jet->SetNEF(neE/jetsIncl[ij].E()); |
7c915297 | 1297 | jet->SetPtLeading(leadPt); |
1298 | jet->SetPtSubtracted(leadTrkType,0); | |
b6e50056 | 1299 | |
1300 | Double_t effAErrCh = 0, effAErrNe = leadTrkType; | |
7c915297 | 1301 | Double_t chBgEnergy = 10, neBgEnergy = nNePart; |
b6e50056 | 1302 | if(!isTruth) |
1303 | { | |
1304 | effAErrCh = GetMeasuredJetPtResolution(ij,jetFinder); | |
7c915297 | 1305 | if(fIsMC && constituents[leadIndex].user_index()>0) |
1306 | { | |
1307 | AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[leadIndex].user_index()-1); | |
1308 | Double_t pt = track->Pt(); | |
1309 | Int_t ipart = track->GetLabel(); | |
1310 | if(ipart>-1 && ipart<fMC->GetNumberOfTracks()) | |
1311 | { | |
1312 | AliVParticle* vParticle = fMC->GetTrack(ipart); | |
1313 | if(vParticle) | |
1314 | { | |
1315 | Double_t truePt = vParticle->Pt(); | |
1316 | chBgEnergy = (truePt-pt)/truePt; | |
1317 | } | |
1318 | } | |
1319 | } | |
b6e50056 | 1320 | } |
7c915297 | 1321 | |
b6e50056 | 1322 | jet->SetEffArea(leadPt, jetFinder->GetJetArea(ij),effAErrCh,effAErrNe); |
1323 | jet->SetBgEnergy(chBgEnergy,neBgEnergy); | |
7c915297 | 1324 | if(fVerbosity>10 && isTruth) cout<<jet->ErrorEffectiveAreaCharged()<<" "<<jet->ErrorEffectiveAreaNeutral()<<endl; |
1325 | if(fVerbosity>10) cout<<"jet pt="<<jetsIncl[ij].perp()<<", nef="<<jet->GetNEF()<<", trk eff corr="<<chBgEnergy<<endl; | |
b6e50056 | 1326 | |
7c915297 | 1327 | if(fVerbosity>5) |
b6e50056 | 1328 | printf("# of ref tracks: %d = %d, and nef=%f\n",jet->GetRefTracks()->GetEntries(), (Int_t)constituents.size(), jet->GetNEF()); |
1329 | ||
1330 | // For catch good high-E jets | |
1331 | if(fSpotGoodJet && !isTruth) | |
1332 | { | |
1333 | if(jetsIncl[ij].perp()>100) | |
1334 | { | |
1335 | printf("\n\n--- HAHAHA: High pt jets ---\n"); | |
1336 | printf("File: %s, event = %d, pt=%f\n",CurrentFileName(),(Int_t)Entry(),jetsIncl[ij].perp()); | |
1337 | printf("%s , pt < %f\n", fJetArray->GetName(), fTrkPtMax[1]); | |
1338 | printf("Jet: E=%f, eta=%f, phi=%f, # of constituents=%d, nef=%f\n",jetsIncl[ij].E(),jetsIncl[ij].eta(), jetsIncl[ij].phi()*TMath::RadToDeg(), (Int_t)constituents.size(),jet->GetNEF()); | |
1339 | for(UInt_t ic=0; ic<constituents.size(); ic++) | |
1340 | { | |
1341 | if(constituents[ic].user_index()<0) | |
1342 | { | |
1343 | AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1); | |
1344 | printf("id = %d, cluster with pt = %f, ncell = %d\n",ic,constituents[ic].perp(), cluster->GetNCells()); | |
1345 | } | |
1346 | else | |
1347 | { | |
1348 | AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1); | |
1349 | printf("id = %d, track with pt = %f, track class = %d, originalPt = %f\n",ic,constituents[ic].perp(),(Int_t)track->GetTRDQuality(), track->GetIntegratedLength()); | |
1350 | ||
1351 | } | |
1352 | } | |
1353 | printf("==============================\n\n"); | |
1354 | } | |
1355 | } | |
1356 | // End of catching good high-E jets | |
7c915297 | 1357 | } |
1358 | if(fVerbosity>5) printf("%s has %d jets\n",fJetArray->GetName(), fJetArray->GetEntries()); | |
b6e50056 | 1359 | } |
1360 | ||
5543d2a2 | 1361 | //________________________________________________________________________ |
1362 | Int_t AliAnalysisTaskFullppJet::GetParentParton(const Int_t ipart) | |
1363 | { | |
1364 | // | |
1365 | // Find the parent parton for a given MC particle by tracing back the parent | |
1366 | // DUMMY; NEED improvement | |
1367 | ||
1368 | return ipart; | |
1369 | } | |
1370 | ||
b6e50056 | 1371 | |
1372 | //________________________________________________________________________ | |
1373 | Bool_t AliAnalysisTaskFullppJet::IsGoodJet(fastjet::PseudoJet jet, Double_t rad) | |
1374 | { | |
1375 | // | |
1376 | // Check if the jet pt and direction fulfill the requirement | |
1377 | // | |
1378 | ||
1379 | if(jet.perp()<1) return kFALSE; | |
1380 | if(TMath::Abs(jet.eta())>(0.7-rad)) return kFALSE; | |
1381 | if(jet.phi() < (80*TMath::DegToRad()+rad) || jet.phi() > (180*TMath::DegToRad()-rad) ) return kFALSE; | |
1382 | return kTRUE; | |
1383 | } | |
1384 | ||
1385 | ||
1386 | //________________________________________________________________________ | |
1387 | void AliAnalysisTaskFullppJet::ProcessMC(const Int_t r) | |
1388 | { | |
1389 | // | |
1390 | // Main function for jet finding in MC | |
1391 | // | |
1392 | ||
1393 | Int_t npart = fMC->GetNumberOfTracks(); | |
1394 | fMcPartArray->Set(npart); | |
1395 | Int_t countPart = 0; | |
1396 | for(Int_t ipart=0; ipart<npart; ipart++) | |
1397 | { | |
1398 | AliVParticle* vParticle = fMC->GetTrack(ipart); | |
1399 | if(!IsGoodMcPartilce(vParticle,ipart)) continue; | |
1400 | ||
1401 | Int_t pdgCode = vParticle->PdgCode(); | |
1402 | if( fRejectNK && (pdgCode==130 || pdgCode==2112) ) continue; | |
1403 | ||
1404 | if( fRejectWD && (pdgCode==310 || pdgCode==3112 || pdgCode==3122 || pdgCode==3222 || pdgCode==3312 || pdgCode==3322 || pdgCode==3334)) continue; | |
1405 | ||
1406 | if( fChargedMC && vParticle->Charge()==0 ) continue; | |
1407 | ||
1408 | Int_t index = 1; | |
7c915297 | 1409 | if(vParticle->Charge()==0) { index=-1; } |
b6e50056 | 1410 | fMcPartArray->AddAt(ipart, countPart); |
1411 | countPart++; | |
7c915297 | 1412 | |
1413 | fTrueJetFinder[r]->AddInputVector(vParticle->Px(), vParticle->Py(), vParticle->Pz(), vParticle->P(), (ipart+1)*index); | |
1414 | if(fVerbosity>10) | |
1415 | printf("Input particle: ipart=%d, pdg=%d, species=%s, charge=%d, E=%4.3f\n",(ipart+1)*index,pdgCode,vParticle->GetName(), vParticle->Charge(),vParticle->P()); | |
b6e50056 | 1416 | } |
1417 | fMcPartArray->Set(countPart); | |
1418 | fTrueJetFinder[r]->Run(); | |
1419 | ||
7c915297 | 1420 | if(fMcTruthAntikt[r]) FillAODJets(fMcTruthAntikt[r], fTrueJetFinder[r], 1); |
1421 | if(fSaveQAHistos) AnalyzeJets(fTrueJetFinder[r], 2, r); | |
1422 | if(fRunUE && r==0) RunAnalyzeUE(fTrueJetFinder[r], 2, 1); | |
1423 | ||
1424 | // Run analysis on secondary particles | |
1425 | if(fRunSecondaries && fAnaType==0) | |
1426 | { | |
1427 | for(Int_t i=0; i<2; i++) | |
1428 | { | |
1429 | AnalyzeSecondaryContribution(fTrueJetFinder[r],r,i); | |
1430 | AnalyzeSecondaryContributionViaMatching(fTrueJetFinder[r],r,0,i); | |
1431 | AnalyzeSecondaryContributionViaMatching(fTrueJetFinder[r],r,1,i); | |
1432 | } | |
1433 | } | |
b6e50056 | 1434 | } |
1435 | ||
1436 | ||
1437 | //________________________________________________________________________ | |
1438 | void AliAnalysisTaskFullppJet::AnalyzeJets(AliFJWrapper *jetFinder, const Int_t type, const Int_t r) | |
1439 | { | |
1440 | // | |
1441 | // Fill all the QA plots for jets | |
1442 | // Especailly all the constituents plot must be filled here since the information | |
1443 | // is not available in the output AOD | |
1444 | ||
1445 | Double_t vertex[3] = {0, 0, 0}; | |
1446 | fESD->GetVertex()->GetXYZ(vertex) ; | |
1447 | TLorentzVector gamma; | |
1448 | ||
1449 | const Int_t nBins = fJetEnergyFraction[type][r]->GetAxis(1)->GetNbins(); | |
1450 | Float_t radiusCut[nBins]; | |
1451 | Float_t eFraction[nBins]; | |
1452 | Int_t nPart[nBins]; | |
1453 | for(Int_t i=0; i<nBins; i++) | |
1454 | { | |
1455 | radiusCut[i] = (fJetEnergyFraction[type][r]->GetAxis(1)->GetBinWidth(1)) * (i+1); | |
1456 | eFraction[i] = 0.; | |
1457 | nPart[nBins] = 0; | |
1458 | } | |
1459 | std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets(); | |
1460 | ||
1461 | for(UInt_t ij=0; ij<jetsIncl.size(); ij++) | |
1462 | { | |
7c915297 | 1463 | if(type<2 && fCheckTPCOnlyVtx && fIsTPCOnlyVtx) |
1464 | { | |
1465 | fhJetInTPCOnlyVtx[type][r]->Fill(jetsIncl[ij].perp(),jetsIncl[ij].phi()*TMath::RadToDeg(),jetsIncl[ij].eta()); | |
1466 | } | |
b6e50056 | 1467 | if(!IsGoodJet(jetsIncl[ij],kRadius[r])) continue; // Fidiual cut |
1468 | Float_t jetEta = jetsIncl[ij].eta(); | |
1469 | Float_t jetPhi = jetsIncl[ij].phi(); | |
1470 | Float_t jetE = jetsIncl[ij].E(); | |
1471 | Float_t jetPt = jetsIncl[ij].perp(); | |
1472 | ||
7c915297 | 1473 | std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij); |
b6e50056 | 1474 | Double_t neE=0, leadingZ = 0, maxPt = 0; |
1475 | Int_t constituentType = -1, leadingIndex = 0; | |
1476 | for(UInt_t ic=0; ic<constituents.size(); ic++) | |
1477 | { | |
1478 | if(constituents[ic].perp()>maxPt) | |
1479 | { | |
1480 | maxPt = constituents[ic].perp(); | |
1481 | leadingIndex = constituents[ic].user_index(); | |
1482 | } | |
1483 | ||
1484 | if(constituents[ic].user_index()<0) | |
1485 | { | |
1486 | neE += constituents[ic].E(); | |
1487 | constituentType = 3; | |
1488 | } | |
1489 | else | |
1490 | { | |
1491 | if(type==2) constituentType = 0; | |
1492 | else | |
1493 | { | |
1494 | AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1); | |
1495 | constituentType = (Int_t)track->GetTRDQuality(); | |
1496 | } | |
1497 | } | |
1498 | Double_t cz = GetZ(constituents[ic].px(),constituents[ic].py(),constituents[ic].pz(),jetsIncl[ij].px(),jetsIncl[ij].py(),jetsIncl[ij].pz()); | |
1499 | fhJetPtVsZ[type][r]->Fill(jetPt,cz, constituentType); | |
1500 | if(cz>leadingZ) leadingZ = cz; | |
1501 | } | |
1502 | ||
1503 | if(type<2 && leadingIndex<0) | |
1504 | { | |
1505 | AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(leadingIndex*(-1)-1); | |
1506 | fhFcrossVsZleading[type][r]->Fill(jetPt,GetExoticEnergyFraction(cluster),leadingZ); | |
1507 | } | |
1508 | ||
1509 | if(leadingZ > 0.98) continue; // Z cut | |
1510 | ||
1511 | fJetCount[type][r]->Fill(jetPt); | |
7c915297 | 1512 | if(type==1 && fIsExoticEvent3GeV) fhJetPtInExoticEvent[0][r]->Fill(jetPt); |
1513 | if(type==1 && fIsExoticEvent5GeV) fhJetPtInExoticEvent[1][r]->Fill(jetPt); | |
b6e50056 | 1514 | |
1515 | Double_t totTrkPt[3] = {0.,0.,0.}; | |
1516 | Double_t chPt = 0; | |
1517 | for(Int_t i=0; i<nBins; i++) { eFraction[i] = 0.; nPart[i] = 0;} | |
1518 | Double_t mcSubE=0.; | |
1519 | Double_t subClsE[5] = {0,0,0,0,0}; | |
1520 | Double_t subTrkPtClean[5] = {0,0,0,0,0}; | |
1521 | Double_t subTrkPtAmbig[5] = {0,0,0,0,0}; | |
1522 | Double_t fraction[5] = {270,0.3,0.5,0.7,1}; | |
1523 | Double_t leadChPt=0., leadNePt=0.; | |
1524 | Int_t leadChIndex = -1; | |
1525 | for(UInt_t ic=0; ic<constituents.size(); ic++) | |
1526 | { | |
1527 | Float_t partEta = constituents[ic].eta(); | |
1528 | Float_t partPhi = constituents[ic].phi(); | |
1529 | Float_t partE = constituents[ic].E(); | |
1530 | Float_t partPt = constituents[ic].perp(); | |
1531 | ||
1532 | if(constituents[ic].user_index()<0) | |
1533 | { | |
1534 | fhNeutralPtInJet[type][r]->Fill(jetPt, partPt); | |
1535 | if(partPt>leadNePt) | |
1536 | leadNePt = partPt; | |
7c915297 | 1537 | if(type<2) |
b6e50056 | 1538 | { |
1539 | AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1); | |
1540 | Double_t clsE = cluster->E(); | |
7c915297 | 1541 | if(cluster->Chi2()>0.5) fhTrigNeuPtInJet[type][r]->Fill(jetPt, partPt); |
1542 | if(fStudySubEInHC || fStudyMcOverSubE) | |
b6e50056 | 1543 | { |
7c915297 | 1544 | cluster->GetMomentum(gamma, vertex); |
1545 | Double_t sinTheta = TMath::Sqrt(1-TMath::Power(gamma.CosTheta(),2)); | |
1546 | ||
1547 | Double_t subEtmp = cluster->GetDispersion(); | |
1548 | Double_t mipETmp = cluster->GetEmcCpvDistance(); | |
1549 | mcSubE += cluster->GetDistanceToBadChannel()*sinTheta; | |
1550 | subClsE[0] += ((mipETmp>clsE)?clsE:mipETmp)*sinTheta; | |
1551 | for(Int_t j=1; j<5; j++) | |
1552 | { | |
1553 | subClsE[j] += ((fraction[j]*subEtmp>clsE)?clsE:fraction[j]*subEtmp)*sinTheta; | |
1554 | } | |
b6e50056 | 1555 | } |
1556 | } | |
1557 | } | |
1558 | else | |
1559 | { | |
1560 | if(partPt>leadChPt) {leadChPt = partPt;leadChIndex=ic;} | |
1561 | fhChargedPtInJet[type][r]->Fill(jetPt, partPt); | |
1562 | chPt += constituents[ic].perp(); | |
1563 | if(type<2) | |
1564 | { | |
1565 | AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1); | |
1566 | Int_t trkType = (Int_t)track->GetTRDQuality(); | |
1567 | totTrkPt[trkType] += partPt; | |
1568 | Int_t clusterIndex = track->GetEMCALcluster(); | |
1569 | Int_t clusterPos = -1; | |
1570 | Bool_t isExist = kFALSE; | |
1571 | for(Int_t j=0; j<fClusterArray->GetEntriesFast(); j++) | |
1572 | { | |
1573 | AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(j); | |
1574 | if( clusterIndex == cluster->GetID() ) | |
1575 | { | |
1576 | isExist = kTRUE; | |
1577 | clusterPos = j; | |
1578 | break; | |
1579 | } | |
1580 | } | |
1581 | if(isExist) | |
1582 | { | |
1583 | AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(clusterPos); | |
1584 | Double_t subEtmp = cluster->GetDispersion(); | |
1585 | Double_t mipETmp = cluster->GetEmcCpvDistance(); | |
1586 | for(Int_t k=0; k<5; k++) | |
1587 | { | |
1588 | if(k==0) | |
1589 | { | |
1590 | if(mipETmp>cluster->E()) subTrkPtClean[k] += partPt; | |
1591 | else subTrkPtAmbig[k] += partPt; | |
1592 | } | |
1593 | else | |
1594 | { | |
1595 | if(fraction[k]*subEtmp>cluster->E()) subTrkPtClean[k] += partPt; | |
1596 | else subTrkPtAmbig[k] += partPt; | |
1597 | } | |
1598 | } | |
1599 | } | |
1600 | } | |
1601 | } | |
1602 | ||
1603 | Float_t dR = TMath::Sqrt( (partEta-jetEta)*(partEta-jetEta) + (partPhi-jetPhi)*(partPhi-jetPhi) ); | |
1604 | for(Int_t i=0; i<nBins; i++) | |
1605 | { | |
1606 | if( dR < radiusCut[i] ) | |
1607 | { | |
1608 | eFraction[i] += partE; | |
1609 | nPart[i]++; | |
1610 | } | |
1611 | } | |
1612 | } | |
1613 | ||
1614 | fhLeadNePtInJet[type][r]->Fill(jetPt, leadNePt); | |
1615 | fhLeadChPtInJet[type][r]->Fill(jetPt, leadChPt); | |
1616 | if(type<2 && leadChIndex>-1) | |
1617 | { | |
1618 | Double_t cz = GetZ(constituents[leadChIndex].px(),constituents[leadChIndex].py(),constituents[leadChIndex].pz(),jetsIncl[ij].px(),jetsIncl[ij].py(),jetsIncl[ij].pz()); | |
1619 | fhChLeadZVsJetPt[type][r]->Fill(jetPt, cz); | |
1620 | } | |
1621 | ||
1622 | if((fStudySubEInHC||fStudyMcOverSubE) && type<2) | |
1623 | { | |
1624 | for(Int_t i=0; i<5; i++) | |
1625 | { | |
1626 | if(fStudySubEInHC) | |
1627 | { | |
1628 | fhSubClsEVsJetPt[type][r][i]->Fill(jetPt-subClsE[4],subClsE[i]/(jetPt-subClsE[4])); | |
1629 | fhHCTrkPtClean[type][r][i]->Fill(jetPt-subClsE[4],subTrkPtClean[i]/(jetPt-subClsE[4])); | |
1630 | fhHCTrkPtAmbig[type][r][i]->Fill(jetPt-subClsE[4],subTrkPtAmbig[i]/(jetPt-subClsE[4])); | |
1631 | } | |
1632 | if(type==0 && fIsMC && fStudyMcOverSubE) | |
1633 | { | |
1634 | fHCOverSubE[r][i]->Fill(jetPt-subClsE[4],subClsE[i]-mcSubE); | |
1635 | fHCOverSubEFrac[r][i]->Fill(jetPt-subClsE[4],(subClsE[i]-mcSubE)/(jetPt-subClsE[4])); | |
1636 | } | |
1637 | } | |
1638 | } | |
1639 | if(type<2 && chPt>0) | |
1640 | { | |
1641 | for(Int_t i=0; i<3; i++) | |
1642 | { | |
1643 | fRelTrkCon[type][r]->Fill(jetPt,totTrkPt[i]/chPt,i); | |
1644 | } | |
1645 | } | |
1646 | ||
1647 | ||
1648 | for(Int_t ibin=0; ibin<nBins; ibin++) | |
1649 | { | |
1650 | Double_t fill1[3] = {jetPt,radiusCut[ibin]-0.005,eFraction[ibin]/jetE}; | |
1651 | fJetEnergyFraction[type][r]->Fill(fill1); | |
1652 | Double_t fill2[3] = {jetPt,radiusCut[ibin]-0.005,nPart[ibin]}; | |
1653 | fJetNPartFraction[type][r]->Fill(fill2); | |
1654 | } | |
1655 | ||
1656 | // Get the jet pt containing tracks or clusters above some threshold | |
1657 | if(type<2) | |
1658 | { | |
1659 | Double_t thres[3] = {15,25,40}; | |
1660 | Int_t okCh[3] = {0,0,0}; | |
1661 | Int_t okNe[3] = {0,0,0}; | |
7c915297 | 1662 | Double_t lowPt[2] = {0.,0.}; |
b6e50056 | 1663 | for(UInt_t ic=0; ic<constituents.size(); ic++) |
1664 | { | |
1665 | Float_t partPt = constituents[ic].perp(); | |
7c915297 | 1666 | |
1667 | if(partPt<0.3) lowPt[0] += partPt; | |
1668 | else if(partPt<0.5) lowPt[1] += partPt; | |
1669 | ||
b6e50056 | 1670 | if(constituents[ic].user_index()>0) |
1671 | { | |
1672 | for(Int_t it=0; it<3; it++) | |
1673 | { | |
1674 | if(partPt>thres[it]) okCh[it] = 1; | |
1675 | } | |
1676 | } | |
1677 | else | |
1678 | { | |
1679 | for(Int_t icl=0; icl<3; icl++) | |
1680 | { | |
1681 | if(partPt>thres[icl]) okNe[icl] = 1; | |
1682 | } | |
1683 | } | |
1684 | } | |
1685 | for(Int_t i=0; i<3; i++) | |
1686 | { | |
1687 | if(okCh[i]==1) | |
1688 | fhJetPtWithTrkThres[type][r]->Fill(i,jetPt); | |
1689 | if(okNe[i]==1) | |
1690 | fhJetPtWithClsThres[type][r]->Fill(i,jetPt); | |
1691 | } | |
7c915297 | 1692 | for(Int_t k=0; k<2; k++) fhJetPtVsLowPtCons[type][r][k]->Fill(jetPt,lowPt[k]/jetPt); |
b6e50056 | 1693 | } |
1694 | } | |
1695 | } | |
1696 | ||
1697 | //________________________________________________________________________ | |
7c915297 | 1698 | void AliAnalysisTaskFullppJet::AnalyzeSecondaryContribution(AliFJWrapper *jetFinder, const Int_t r, const Int_t etaCut) |
1699 | { | |
1700 | // | |
1701 | // Analyze secondaries | |
1702 | // | |
1703 | ||
1704 | std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets(); | |
1705 | for(UInt_t ij=0; ij<jetsIncl.size(); ij++) | |
1706 | { | |
1707 | Float_t jetEta = jetsIncl[ij].eta(); | |
1708 | Float_t jetPt = jetsIncl[ij].perp(); | |
1709 | if(TMath::Abs(jetEta)>0.5*etaCut+0.5 || jetPt<1) continue; | |
1710 | std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij); | |
1711 | Double_t dNKPt = 0, dWPPt = 0, weakPt=0, nkPt=0; | |
1712 | Int_t type = -1; | |
1713 | for(UInt_t ic=0; ic<constituents.size(); ic++) | |
1714 | { | |
1715 | Int_t ipart = TMath::Abs(constituents[ic].user_index())-1; | |
1716 | AliMCParticle* mcParticle = (AliMCParticle*)fMC->GetTrack(ipart); | |
1717 | if(!mcParticle) continue; | |
1718 | Int_t pdg = mcParticle->PdgCode(); | |
1719 | if(pdg==310 || (pdg<=3400 && pdg>=3100)) weakPt += mcParticle->Pt(); | |
1720 | if(pdg==130 || pdg==2112) nkPt += mcParticle->Pt(); | |
1721 | ||
1722 | // Weak decaying particles | |
1723 | if(pdg==310) type=0; | |
1724 | else if (pdg==3122) type=1; | |
1725 | else if (pdg<=3400 && pdg>=3100) type=2; | |
1726 | else type=-1; | |
1727 | if(type>-1) | |
1728 | { | |
1729 | Int_t binx = fhSecondaryResponse[type]->GetXaxis()->FindFixBin(mcParticle->Pt()); | |
1730 | TH1F *htmp = (TH1F*)fhSecondaryResponse[type]->ProjectionY(Form("hpro_%d_%d_%d",(Int_t)Entry(),ij,ic),binx,binx); | |
1731 | Double_t pro = htmp->GetRandom(); | |
1732 | dWPPt += pro*mcParticle->Pt() - mcParticle->Pt(); | |
1733 | delete htmp; | |
1734 | } | |
1735 | ||
1736 | // Missing neutron and K0L | |
1737 | if(pdg==130 || pdg==2112) dNKPt -= mcParticle->Pt(); | |
1738 | } | |
1739 | ||
1740 | fhNKFracVsJetPt[etaCut][r]->Fill(jetPt,nkPt/jetPt); | |
1741 | fhWeakFracVsJetPt[etaCut][r]->Fill(jetPt,weakPt/jetPt); | |
1742 | ||
1743 | Double_t jetNewWPPt = jetPt + dWPPt; | |
1744 | fhJetResponseWP[etaCut][r]->Fill(jetPt,jetNewWPPt); | |
1745 | fhJetResolutionWP[etaCut][r]->Fill(jetPt,(jetPt-jetNewWPPt)/jetPt); | |
1746 | ||
1747 | Double_t jetNewNKPt = jetPt + dNKPt; | |
1748 | fhJetResponseNK[etaCut][r]->Fill(jetPt,jetNewNKPt); | |
1749 | fhJetResolutionNK[etaCut][r]->Fill(jetPt,(jetPt-jetNewNKPt)/jetPt); | |
1750 | } | |
1751 | } | |
1752 | ||
1753 | ||
1754 | //________________________________________________________________________ | |
1755 | void AliAnalysisTaskFullppJet::AnalyzeSecondaryContributionViaMatching(AliFJWrapper *jetFinder, const Int_t r, const Int_t type, const Int_t etaCut) | |
1756 | { | |
1757 | // | |
1758 | // Estimate the contribution of missing energies via matching | |
1759 | // type = 0 -> NK | |
1760 | // type = 1 -> WP | |
1761 | // | |
1762 | ||
1763 | if(type!=0 && type!=1) return; | |
1764 | ||
1765 | AliFJWrapper fj("fj","fj"); | |
1766 | fj.CopySettingsFrom(*jetFinder); | |
1767 | ||
1768 | Int_t npart = fMC->GetNumberOfTracks(); | |
1769 | Int_t pType = -1; | |
1770 | for(Int_t ipart=0; ipart<npart; ipart++) | |
1771 | { | |
1772 | AliVParticle* vParticle = fMC->GetTrack(ipart); | |
1773 | if(!IsGoodMcPartilce(vParticle,ipart)) continue; | |
1774 | Int_t pdgCode = vParticle->PdgCode(); | |
1775 | Double_t pt = vParticle->Pt(); | |
1776 | if( type==0 && (pdgCode==130 || pdgCode==2112) ) continue; // Reject NK | |
1777 | if(type==1) | |
1778 | { | |
1779 | if(pdgCode==310) pType=0; | |
1780 | else if (pdgCode==3122) pType=1; | |
1781 | else if (pdgCode<=3400 && pdgCode>=3100) pType=2; | |
1782 | else pType=-1; | |
1783 | if(pType>-1) | |
1784 | { | |
1785 | Int_t binx = fhSecondaryResponse[pType]->GetXaxis()->FindFixBin(vParticle->Pt()); | |
1786 | TH1F *htmp = (TH1F*)fhSecondaryResponse[pType]->ProjectionY(Form("hpro_%d_%d",(Int_t)Entry(),ipart),binx,binx); | |
1787 | Double_t pro = htmp->GetRandom(); | |
1788 | pt = pro * vParticle->Pt(); | |
1789 | delete htmp; | |
1790 | } | |
1791 | } | |
1792 | Int_t index = 1; | |
1793 | if(vParticle->Charge()==0) { index=-1; } | |
1794 | fj.AddInputVector(vParticle->Px()*pt/vParticle->Pt(), vParticle->Py()*pt/vParticle->Pt(), vParticle->Pz(), vParticle->P(), (ipart+1)*index); | |
1795 | } | |
1796 | ||
1797 | fj.Run(); | |
1798 | std::vector<fastjet::PseudoJet> jets_incl_mc = fj.GetInclusiveJets(); | |
1799 | std::vector<fastjet::PseudoJet> jets_incl_mc_std = jetFinder->GetInclusiveJets(); | |
1800 | ||
1801 | for(UInt_t ij=0; ij<jets_incl_mc.size(); ij++) | |
1802 | { | |
1803 | Float_t jetEta = jets_incl_mc[ij].eta(); | |
1804 | Float_t jetPt = jets_incl_mc[ij].perp(); | |
1805 | if(TMath::Abs(jetEta)>0.5*etaCut+0.5 || jetPt<5) continue; | |
1806 | Double_t dEta, dPhi; | |
1807 | Int_t index = FindSpatialMatchedJet(jets_incl_mc[ij], jetFinder, dEta, dPhi, 1); | |
1808 | if(index>-1) | |
1809 | { | |
1810 | Double_t jetPtStd = jets_incl_mc_std[index].perp(); | |
1811 | Double_t dR = TMath::Sqrt(dEta*dEta+dPhi*dPhi); | |
1812 | if(type==0) fhJetResolutionNKSM[etaCut][r]->Fill(jetPtStd,(jetPtStd-jetPt)/jetPtStd,dR); | |
1813 | if(type==1) fhJetResolutionWPSM[etaCut][r]->Fill(jetPtStd,(jetPtStd-jetPt)/jetPtStd,dR); | |
1814 | ||
1815 | if(dR<kdRCut[r]) | |
1816 | { | |
1817 | if(type==0) fhJetResponseNKSM[etaCut][r]->Fill(jetPtStd,jetPt); | |
1818 | if(type==1) fhJetResponseWPSM[etaCut][r]->Fill(jetPtStd,jetPt); | |
1819 | } | |
1820 | } | |
1821 | } | |
1822 | } | |
1823 | ||
1824 | //________________________________________________________________________ | |
1825 | void AliAnalysisTaskFullppJet::RunAnalyzeUE(AliFJWrapper *jetFinder, const Int_t type, const Bool_t isMCTruth) | |
b6e50056 | 1826 | { |
1827 | // | |
1828 | // Run analysis to estimate the underlying event | |
7c915297 | 1829 | // Adapted from Oliver |
b6e50056 | 1830 | |
1831 | std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets(); | |
1832 | Double_t leadpt=0, leadPhi = -999, leadEta = -999; | |
7c915297 | 1833 | Int_t leadIndex = -1; |
b6e50056 | 1834 | for(UInt_t ij=0; ij<jetsIncl.size(); ij++) |
1835 | { | |
b6e50056 | 1836 | Double_t jetEta = jetsIncl[ij].eta(); |
7c915297 | 1837 | Double_t jetPt = jetsIncl[ij].perp(); |
b6e50056 | 1838 | Double_t jetPhi = jetsIncl[ij].phi(); |
b6e50056 | 1839 | if(leadpt<jetPt) |
1840 | { | |
1841 | leadpt = jetPt; | |
1842 | leadPhi = jetPhi; | |
1843 | leadEta = jetEta; | |
7c915297 | 1844 | leadIndex = ij; |
b6e50056 | 1845 | } |
1846 | } | |
7c915297 | 1847 | if(leadpt<1 || TMath::Abs(leadEta)>0.5 ) return; |
b6e50056 | 1848 | |
7c915297 | 1849 | Int_t backIndex = -1; |
1850 | Bool_t isBackToBack = kFALSE; | |
1851 | for(UInt_t ij=0; ij<jetsIncl.size(); ij++) | |
b6e50056 | 1852 | { |
7c915297 | 1853 | Double_t dPhi = TMath::Abs(jetsIncl[ij].phi()-leadPhi); |
1854 | if(dPhi > kPI) dPhi = 2*kPI - dPhi; | |
1855 | if(dPhi>150*TMath::DegToRad() && jetsIncl[ij].perp()/leadpt > 0.8) | |
b6e50056 | 1856 | { |
7c915297 | 1857 | backIndex = ij; |
1858 | isBackToBack = kTRUE; | |
b6e50056 | 1859 | } |
1860 | } | |
b6e50056 | 1861 | |
7c915297 | 1862 | for(Int_t ij=0; ij<(Int_t)jetsIncl.size(); ij++) |
1863 | { | |
1864 | if(ij==leadIndex || ij==backIndex) continue; | |
1865 | if(TMath::Abs(jetsIncl[ij].eta())>0.5) continue; | |
1866 | if(jetsIncl[ij].perp()>15) isBackToBack = kFALSE; | |
1867 | } | |
1868 | ||
1869 | ||
1870 | Double_t axis[2] = {leadPhi+0.5 * kPI, leadPhi-0.5 * kPI}; | |
1871 | if(axis[0]>2*kPI) axis[0] -= 2*kPI; | |
1872 | if(axis[1]<0) axis[1] += 2*kPI; | |
1873 | ||
1874 | fhUEJetPtNorm[type][0][0]->Fill(leadpt); | |
1875 | if(isBackToBack) fhUEJetPtNorm[type][0][1]->Fill(leadpt); | |
1876 | ||
1877 | if(!isMCTruth) | |
1878 | { | |
1879 | Int_t ntracks = fTrackArray->GetEntriesFast(); | |
1880 | Double_t vertex[3] = {0, 0, 0}; | |
1881 | fESD->GetVertex()->GetXYZ(vertex) ; | |
1882 | TLorentzVector gamma; | |
1883 | Int_t nclusters = fClusterArray->GetEntriesFast(); | |
1884 | ||
1885 | for(Int_t j=0; j<2; j++) | |
1886 | { | |
1887 | Double_t ueCh = 0, ueChNe = 0.; | |
1888 | for(Int_t it=0; it<ntracks; it++) | |
1889 | { | |
1890 | AliESDtrack *t = (AliESDtrack*) fTrackArray->At(it); | |
1891 | if(!t) continue; | |
1892 | Double_t dPhi = TMath::Abs(axis[j]-t->Phi()); | |
1893 | Double_t dEta = TMath::Abs(leadEta-t->Eta()); | |
1894 | if(dPhi > kPI) dPhi = 2*kPI - dPhi; | |
1895 | if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4) | |
1896 | { | |
1897 | fhUEJetPtVsConsPt[type][0][0]->Fill(leadpt,t->Pt()); | |
1898 | if(isBackToBack) fhUEJetPtVsConsPt[type][0][1]->Fill(leadpt,t->Pt()); | |
1899 | ueCh += t->Pt(); | |
1900 | if(TMath::Abs(leadEta-0.4)<0.7 && axis[j]>80*TMath::DegToRad()+0.4 && axis[j]<180*TMath::DegToRad()-0.4) | |
1901 | { | |
1902 | fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,t->Pt()); | |
1903 | if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,t->Pt()); | |
1904 | ueChNe += t->Pt(); | |
1905 | } | |
1906 | } | |
1907 | } | |
1908 | fhUEJetPtVsSumPt[type][0][0]->Fill(leadpt,ueCh); | |
1909 | if(isBackToBack) fhUEJetPtVsSumPt[type][0][1]->Fill(leadpt,ueCh); | |
1910 | ||
1911 | if(!(TMath::Abs(leadEta-0.4)<0.7 && axis[j]>80*TMath::DegToRad()+0.4 && axis[j]<180*TMath::DegToRad()-0.4)) continue; | |
1912 | fhUEJetPtNorm[type][1][0]->Fill(leadpt); | |
1913 | if(isBackToBack) fhUEJetPtNorm[type][1][1]->Fill(leadpt); | |
1914 | for(Int_t ic=0; ic<nclusters; ic++) | |
1915 | { | |
1916 | AliESDCaloCluster * cl = (AliESDCaloCluster *)fClusterArray->At(ic); | |
1917 | if(!cl) continue; | |
1918 | cl->GetMomentum(gamma, vertex); | |
1919 | Double_t clsPhi = gamma.Phi(); | |
1920 | if(clsPhi<0) clsPhi += 2*kPI; | |
1921 | Double_t dPhi = TMath::Abs(axis[j]-clsPhi); | |
1922 | Double_t dEta = TMath::Abs(leadEta-gamma.Eta()); | |
1923 | if(dPhi > kPI) dPhi = 2*kPI - dPhi; | |
1924 | if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4) | |
1925 | { | |
1926 | ueChNe += gamma.Pt(); | |
1927 | fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,gamma.Pt()); | |
1928 | if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,gamma.Pt()); | |
1929 | } | |
1930 | } | |
1931 | fhUEJetPtVsSumPt[type][1][0]->Fill(leadpt,ueChNe); | |
1932 | if(isBackToBack) fhUEJetPtVsSumPt[type][1][1]->Fill(leadpt,ueChNe); | |
1933 | } | |
1934 | } | |
1935 | else | |
1936 | { | |
1937 | fhUEJetPtNorm[type][1][0]->Fill(leadpt); | |
1938 | if(isBackToBack) fhUEJetPtNorm[type][1][1]->Fill(leadpt); | |
1939 | Int_t npart = fMcPartArray->GetSize(); | |
1940 | for(Int_t j=0; j<2; j++) | |
1941 | { | |
1942 | Double_t ueCh = 0, ueChNe = 0., ueChNe2 = 0.; | |
1943 | for(Int_t ipos=0; ipos<npart; ipos++) | |
1944 | { | |
1945 | AliVParticle* vParticle = fMC->GetTrack(fMcPartArray->At(ipos)); | |
1946 | if(!vParticle) continue; | |
1947 | Double_t dPhi = TMath::Abs(axis[j]-vParticle->Phi()); | |
1948 | Double_t dEta = TMath::Abs(leadEta-vParticle->Eta()); | |
1949 | if(dPhi > kPI) dPhi = 2*kPI - dPhi; | |
1950 | if(TMath::Sqrt(dPhi*dPhi+dEta*dEta)<0.4) | |
1951 | { | |
1952 | Double_t pt = vParticle->Pt(); | |
1953 | ueChNe += pt; | |
1954 | fhUEJetPtVsConsPt[type][1][0]->Fill(leadpt,pt); | |
1955 | if(isBackToBack) fhUEJetPtVsConsPt[type][1][1]->Fill(leadpt,pt); | |
1956 | if( vParticle->Charge()!=0 ) | |
1957 | { | |
1958 | fhUEJetPtVsConsPt[type][0][0]->Fill(leadpt,pt); | |
1959 | if(isBackToBack) fhUEJetPtVsConsPt[type][0][1]->Fill(leadpt,pt); | |
1960 | ueCh += pt; | |
1961 | ueChNe2 += pt; | |
1962 | } | |
1963 | else | |
1964 | { | |
1965 | if(pt>0.5) ueChNe2 += pt; | |
1966 | } | |
1967 | } | |
1968 | } | |
1969 | fhUEJetPtVsSumPt[type][0][0]->Fill(leadpt,ueCh); | |
1970 | fhUEJetPtVsSumPt[type][1][0]->Fill(leadpt,ueChNe); | |
1971 | if(isBackToBack) fhUEJetPtVsSumPt[type][0][1]->Fill(leadpt,ueCh); | |
1972 | if(isBackToBack) fhUEJetPtVsSumPt[type][1][1]->Fill(leadpt,ueChNe); | |
1973 | } | |
1974 | } | |
1975 | } | |
b6e50056 | 1976 | |
1977 | //________________________________________________________________________ | |
1978 | Bool_t AliAnalysisTaskFullppJet::IsGoodJet(AliAODJet *jet, Double_t rad) | |
1979 | { | |
1980 | // | |
1981 | // Check if it is a good jet | |
1982 | // | |
1983 | ||
1984 | if(jet->Pt()<1) return kFALSE; | |
1985 | if(TMath::Abs(jet->Eta())>(0.7-rad)) return kFALSE; | |
1986 | if(jet->Phi() < (80*TMath::DegToRad()+rad) || jet->Phi() > (180*TMath::DegToRad()-rad) ) return kFALSE; | |
1987 | return kTRUE; | |
1988 | } | |
1989 | ||
1990 | ||
1991 | //________________________________________________________________________ | |
1992 | Double_t AliAnalysisTaskFullppJet::GetLeadingZ(const Int_t jetIndex, AliFJWrapper *jetFinder) | |
1993 | { | |
1994 | // | |
1995 | // Get the leading z of the input jet | |
1996 | // | |
1997 | ||
1998 | Double_t z = 0; | |
7c915297 | 1999 | std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex); |
b6e50056 | 2000 | std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets(); |
2001 | Int_t index = -1; | |
2002 | Double_t maxPt = 0; | |
2003 | for(UInt_t ic=0; ic<constituents.size(); ic++) | |
2004 | { | |
2005 | if(constituents[ic].perp()>maxPt) | |
2006 | { | |
2007 | maxPt = constituents[ic].perp(); | |
2008 | index = ic; | |
2009 | } | |
2010 | } | |
2011 | if(index>-1) | |
2012 | z = GetZ(constituents[index].px(),constituents[index].py(),constituents[index].pz(),jetsIncl[jetIndex].px(),jetsIncl[jetIndex].py(),jetsIncl[jetIndex].pz()); | |
2013 | return z; | |
2014 | } | |
2015 | ||
2016 | //________________________________________________________________________ | |
2017 | Double_t AliAnalysisTaskFullppJet::GetZ(const Double_t trkPx, const Double_t trkPy, const Double_t trkPz, const Double_t jetPx, const Double_t jetPy, const Double_t jetPz) const | |
2018 | { | |
2019 | // | |
2020 | // Get the z of a constituent inside of a jet | |
2021 | // | |
2022 | ||
2023 | return (trkPx*jetPx+trkPy*jetPy+trkPz*jetPz)/(jetPx*jetPx+jetPy*jetPy+jetPz*jetPz); | |
2024 | } | |
2025 | ||
2026 | //________________________________________________________________________ | |
2027 | Double_t AliAnalysisTaskFullppJet::GetMeasuredJetPtResolution(const Int_t jetIndex, AliFJWrapper *jetFinder) | |
2028 | { | |
2029 | // | |
2030 | // Get jet energy resoultion due to intrinsic detector effects | |
2031 | // | |
2032 | ||
2033 | Double_t jetSigma2 = 0; | |
7c915297 | 2034 | std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex); |
b6e50056 | 2035 | for(UInt_t ic=0; ic<constituents.size(); ic++) |
2036 | { | |
2037 | if(constituents[ic].user_index()>0) | |
2038 | { | |
2039 | AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1); | |
2040 | jetSigma2 += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2); | |
2041 | } | |
2042 | else | |
2043 | { | |
2044 | AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1); | |
2045 | jetSigma2 += TMath::Power(cluster->GetTOF(),2); | |
2046 | } | |
2047 | } | |
2048 | return TMath::Sqrt(jetSigma2); | |
2049 | } | |
2050 | ||
2051 | ||
2052 | ||
2053 | //________________________________________________________________________ | |
2054 | Double_t AliAnalysisTaskFullppJet::GetJetMissPtDueToTrackingEfficiency(const Int_t jetIndex, AliFJWrapper *jetFinder, const Int_t radiusIndex) | |
2055 | { | |
2056 | // | |
2057 | // Correct the tracking inefficiency explicitly | |
2058 | // | |
2059 | ||
2060 | Double_t misspt = 0; | |
2061 | std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets(); | |
2062 | Double_t jetPt = jetsIncl[jetIndex].perp(); | |
2063 | if(!fhCorrTrkEffPtBin[fTriggerType][radiusIndex]) | |
2064 | { | |
2065 | printf("Warning: can't get the mean # of tracks per jet with pt=%f in: trigger=%d, radiusIndex=%d\n",jetPt,fTriggerType,radiusIndex); | |
2066 | return 0; | |
2067 | } | |
2068 | Int_t ibin = fhCorrTrkEffPtBin[fTriggerType][radiusIndex]->FindFixBin(jetPt); | |
2069 | if(!fhCorrTrkEffSample[fTriggerType][radiusIndex][ibin-1] || fhCorrTrkEffSample[fTriggerType][radiusIndex][ibin-1]->Integral()<0.001) | |
2070 | { | |
2071 | printf("Warning: no sampling distrubtion for jet with pt=%f\n",jetPt); | |
2072 | return 0; | |
2073 | } | |
2074 | Double_t ntrack = fhCorrTrkEffPtBin[fTriggerType][radiusIndex]->GetBinContent(ibin); | |
2075 | Int_t nTrk = (Int_t) ntrack; | |
2076 | Double_t res = ntrack-nTrk; | |
2077 | Double_t pro1 = fRandomGen->Uniform(); | |
2078 | if(pro1<res) nTrk++; | |
2079 | for(Int_t itry=0; itry<nTrk; itry++) | |
2080 | { | |
2081 | Double_t trkPt = fhCorrTrkEffSample[fTriggerType][radiusIndex][ibin-1]->GetRandom(); | |
2082 | if(trkPt/jetPt>fTrkEffCorrCutZ) continue; | |
2083 | Double_t eff = GetTrkEff(trkPt); | |
2084 | Double_t pro = fRandomGen->Uniform(); | |
2085 | if(pro>eff) misspt += trkPt; | |
2086 | } | |
2087 | return misspt; | |
2088 | } | |
2089 | ||
2090 | ||
2091 | //________________________________________________________________________ | |
2092 | Bool_t AliAnalysisTaskFullppJet::IsGoodMcPartilce(const AliVParticle* vParticle, const Int_t ipart) | |
2093 | { | |
2094 | // | |
2095 | // Select good primary particles to feed into jet finder | |
2096 | // | |
2097 | ||
2098 | if(!vParticle) return kFALSE; | |
2099 | if(!fMC->IsPhysicalPrimary(ipart)) return kFALSE; | |
2100 | if (TMath::Abs(vParticle->Eta())>2) return kFALSE; | |
2101 | return kTRUE; | |
2102 | } | |
2103 | ||
2104 | //________________________________________________________________________ | |
7c915297 | 2105 | Int_t AliAnalysisTaskFullppJet::FindSpatialMatchedJet(fastjet::PseudoJet jet, AliFJWrapper *jetFinder, Double_t &dEta, Double_t &dPhi, Double_t maxR) |
b6e50056 | 2106 | { |
2107 | // | |
2108 | // Find spatially matched detector and particle jets | |
2109 | // | |
2110 | ||
b6e50056 | 2111 | Int_t index=-1; |
7c915297 | 2112 | dEta=-999, dPhi=-999; |
2113 | std::vector<fastjet::PseudoJet> jets_incl = jetFinder->GetInclusiveJets(); | |
2114 | for(UInt_t ij=0; ij<jets_incl.size(); ij++) | |
b6e50056 | 2115 | { |
7c915297 | 2116 | if(jets_incl[ij].perp()<5) continue; |
2117 | if(TMath::Abs(jets_incl[ij].eta())>1) continue; | |
2118 | Double_t tmpR = TMath::Sqrt( TMath::Power(jet.eta()-jets_incl[ij].eta(), 2) + TMath::Power(jet.phi()-jets_incl[ij].phi(), 2) ); | |
b6e50056 | 2119 | if(tmpR<maxR) |
2120 | { | |
2121 | maxR=tmpR; | |
2122 | index=ij; | |
7c915297 | 2123 | dEta = jet.eta()-jets_incl[ij].eta(); |
2124 | dPhi = jet.phi()-jets_incl[ij].phi(); | |
b6e50056 | 2125 | } |
2126 | } | |
2127 | return index; | |
2128 | } | |
2129 | ||
2130 | //________________________________________________________________________ | |
2131 | Int_t AliAnalysisTaskFullppJet::FindEnergyMatchedJet(AliFJWrapper *jetFinder1, const Int_t index1, AliFJWrapper *jetFinder2, const Double_t fraction) | |
2132 | { | |
2133 | // | |
2134 | // Find matched detector and particle jets based on shared constituents | |
2135 | // | |
2136 | ||
2137 | Int_t matchedIndex = -1; | |
2138 | std::vector<fastjet::PseudoJet> jetsIncl1 = jetFinder1->GetInclusiveJets(); | |
2139 | std::vector<fastjet::PseudoJet> jetsIncl2 = jetFinder2->GetInclusiveJets(); | |
7c915297 | 2140 | std::vector<fastjet::PseudoJet> constituents1 = jetFinder1->GetJetConstituents(index1); |
b6e50056 | 2141 | Double_t jetPt1 = jetsIncl1[index1].perp(); |
2142 | if(jetPt1<0) return matchedIndex; | |
2143 | ||
2144 | for(UInt_t ij2=0; ij2<jetsIncl2.size(); ij2++) | |
2145 | { | |
2146 | Double_t jetPt2 = jetsIncl2[ij2].perp(); | |
2147 | if(jetPt2<0) return matchedIndex; | |
7c915297 | 2148 | std::vector<fastjet::PseudoJet> constituents2 = jetFinder2->GetJetConstituents(ij2); |
b6e50056 | 2149 | Double_t sharedPt1 = 0., sharedPt2 = 0.; |
2150 | for(UInt_t ic2=0; ic2<constituents2.size(); ic2++) | |
2151 | { | |
2152 | Int_t mcLabel = constituents2[ic2].user_index()-1; | |
2153 | Double_t consPt2 = constituents2[ic2].perp(); | |
2154 | for(UInt_t ic1=0; ic1<constituents1.size(); ic1++) | |
2155 | { | |
2156 | Double_t consPt1 = constituents1[ic1].perp(); | |
2157 | if(constituents1[ic1].user_index()>0) | |
2158 | { | |
2159 | AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents1[ic1].user_index()-1); | |
2160 | if(track->GetLabel()==mcLabel) {sharedPt2 += consPt2; sharedPt1 += consPt1; cout<<"Found a matched track"<<endl;break;} | |
2161 | } | |
2162 | else | |
2163 | { | |
2164 | AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents1[ic1].user_index()*(-1)-1); | |
2165 | if(cluster->GetLabel()==mcLabel) {sharedPt2 += consPt2; sharedPt1 += consPt1; cout<<"Found a matched cluster"<<endl;break;} | |
2166 | } | |
2167 | } | |
2168 | } | |
2169 | cout<<sharedPt1/jetPt1<<" "<<sharedPt2/jetPt2<<endl; | |
2170 | if(sharedPt1/jetPt1 > fraction && sharedPt2/jetPt2 > fraction) | |
2171 | { | |
2172 | matchedIndex = ij2; | |
2173 | break; | |
2174 | } | |
2175 | } | |
2176 | return matchedIndex; | |
2177 | } | |
2178 | ||
2179 | //________________________________________________________________________ | |
2180 | void AliAnalysisTaskFullppJet::GetESDTrax() | |
2181 | { | |
2182 | // | |
2183 | // Get esd tracks. | |
2184 | // | |
2185 | ||
2186 | Int_t nTrax = fESD->GetNumberOfTracks(); | |
2187 | if(fVerbosity>5) printf("[i] # of tracks in event: %d\n",nTrax); | |
2188 | ||
2189 | for (Int_t i = 0; i < nTrax; ++i) | |
2190 | { | |
2191 | AliESDtrack* esdtrack = fESD->GetTrack(i); | |
2192 | if (!esdtrack) | |
2193 | { | |
2194 | AliError(Form("Couldn't get ESD track %d\n", i)); | |
2195 | continue; | |
2196 | } | |
2197 | ||
2198 | AliESDtrack *newtrack = GetAcceptTrack(esdtrack); | |
2199 | if(!newtrack) continue; | |
2200 | ||
2201 | Double_t pt = newtrack->Pt(); | |
2202 | Double_t eta = newtrack->Eta(); | |
7c915297 | 2203 | if (pt < fTrkPtMin[fTriggerType] || pt > fTrkPtMax[fTriggerType] || TMath::Abs(eta) > fTrkEtaMax) |
b6e50056 | 2204 | { |
2205 | delete newtrack; | |
2206 | continue; | |
2207 | } | |
2208 | ||
2209 | if(fSysTrkEff) | |
2210 | { | |
2211 | Double_t rand = fRandomGen->Uniform(); | |
2212 | if(rand<fVaryTrkEff) | |
2213 | { | |
2214 | delete newtrack; | |
2215 | continue; | |
2216 | } | |
2217 | } | |
2218 | newtrack->SetIntegratedLength(esdtrack->Pt()); | |
2219 | fTrackArray->Add(newtrack); | |
2220 | } | |
2221 | if(fVerbosity>5) printf("[i] # of tracks in event: %d\n", fTrackArray->GetEntries()); | |
2222 | fNTracksPerChunk += fTrackArray->GetEntries(); | |
2223 | } | |
2224 | ||
2225 | // | |
2226 | //________________________________________________________________________ | |
2227 | // | |
2228 | AliESDtrack *AliAnalysisTaskFullppJet::GetAcceptTrack(AliESDtrack *esdtrack) | |
2229 | { | |
2230 | // | |
2231 | // Get the hybrid tracks | |
2232 | // | |
2233 | ||
2234 | AliESDtrack *newTrack = 0x0; | |
2235 | ||
2236 | if(fTrackCutsType==0 || fTrackCutsType==3) | |
2237 | { | |
2238 | if(fEsdTrackCuts->AcceptTrack(esdtrack)) | |
2239 | { | |
2240 | newTrack = new AliESDtrack(*esdtrack); | |
2241 | newTrack->SetTRDQuality(0); | |
2242 | } | |
2243 | else if(fHybridTrackCuts1->AcceptTrack(esdtrack)) | |
2244 | { | |
2245 | if(esdtrack->GetConstrainedParam()) | |
2246 | { | |
2247 | newTrack = new AliESDtrack(*esdtrack); | |
2248 | const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam(); | |
2249 | newTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance()); | |
2250 | newTrack->SetTRDQuality(1); | |
2251 | } | |
2252 | else | |
2253 | return 0x0; | |
2254 | } | |
2255 | else if(fHybridTrackCuts2->AcceptTrack(esdtrack)) | |
2256 | { | |
2257 | if(esdtrack->GetConstrainedParam()) | |
2258 | { | |
2259 | newTrack = new AliESDtrack(*esdtrack); | |
2260 | const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam(); | |
2261 | newTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance()); | |
2262 | newTrack->SetTRDQuality(2); | |
2263 | } | |
2264 | else | |
2265 | return 0x0; | |
2266 | } | |
2267 | else | |
2268 | { | |
2269 | return 0x0; | |
2270 | } | |
2271 | } | |
2272 | else if(fTrackCutsType==1) | |
2273 | { | |
2274 | if(fEsdTrackCuts->AcceptTrack(esdtrack)) | |
2275 | { | |
2276 | newTrack = AliESDtrackCuts::GetTPCOnlyTrack(const_cast<AliESDEvent*>(fESD),esdtrack->GetID());// use TPC only tracks with non default SPD vertex | |
2277 | if(!newTrack) return 0x0; | |
2278 | AliExternalTrackParam exParam; | |
2279 | Bool_t relate = newTrack->RelateToVertexTPC(fESD->GetPrimaryVertexSPD(),fESD->GetMagneticField(),kVeryBig,&exParam); //constrain to SPD vertex | |
2280 | if( !relate ) | |
2281 | { | |
2282 | delete newTrack; | |
2283 | return 0x0; | |
2284 | } | |
2285 | newTrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance()); | |
2286 | newTrack->SetTRDQuality(1); | |
2287 | } | |
2288 | else | |
2289 | { | |
2290 | return 0x0; | |
2291 | } | |
2292 | } | |
2293 | else if (fTrackCutsType==2) | |
2294 | { | |
2295 | if(fEsdTrackCuts->AcceptTrack(esdtrack)) | |
2296 | { | |
2297 | newTrack = new AliESDtrack(*esdtrack); | |
2298 | newTrack->SetTRDQuality(0); | |
2299 | } | |
2300 | else | |
2301 | return 0x0; | |
2302 | } | |
2303 | else | |
2304 | { | |
2305 | printf("Unknown track cuts type: %d\n",fTrackCutsType); | |
2306 | return 0x0; | |
2307 | } | |
2308 | ||
2309 | return newTrack; | |
2310 | } | |
2311 | ||
2312 | //________________________________________________________________________ | |
2313 | void AliAnalysisTaskFullppJet::GetESDEMCalClusters() | |
2314 | { | |
2315 | // | |
2316 | // Get emcal clusters - selected | |
2317 | // | |
2318 | ||
2319 | if(fSysTrkClsMth) | |
2320 | { | |
2321 | if (!TGeoGlobalMagField::Instance()->GetField()) fESD->InitMagneticField(); | |
2322 | fRecoUtil->FindMatches(fESD,0x0,fGeom); | |
2323 | fRecoUtil->SetClusterMatchedToTrack(fESD); | |
2324 | fRecoUtil->SetTracksMatchedToCluster(fESD); | |
2325 | } | |
2326 | ||
b6e50056 | 2327 | const Int_t nCaloClusters = fESD->GetNumberOfCaloClusters(); |
2328 | for(Int_t i = 0 ; i < nCaloClusters; i++) | |
2329 | { | |
2330 | AliESDCaloCluster * cl = (AliESDCaloCluster *) fESD->GetCaloCluster(i); | |
2331 | if(!IsGoodCluster(cl)) continue; | |
2332 | AliESDCaloCluster *newCluster = new AliESDCaloCluster(*cl); | |
7c915297 | 2333 | if (newCluster->E() < fClsEtMin[fTriggerType] || newCluster->E() > fClsEtMax[fTriggerType]) |
2334 | {delete newCluster; continue;} | |
b6e50056 | 2335 | |
7c915297 | 2336 | // Absolute scale |
2337 | if(fPeriod.Contains("lhc11a",TString::kIgnoreCase)) | |
b6e50056 | 2338 | { |
7c915297 | 2339 | Double_t newE = newCluster->E() * 1.02; |
b6e50056 | 2340 | newCluster->SetE(newE); |
b6e50056 | 2341 | } |
2342 | ||
7c915297 | 2343 | // Trigger efficiency systematic uncertainty |
2344 | if(fSysJetTrigEff) | |
b6e50056 | 2345 | { |
7c915297 | 2346 | Double_t newE = newCluster->E() * (1+fVaryJetTrigEff); |
2347 | newCluster->SetE(newE); | |
2348 | if(fTriggerType==0) fhClsE[fTriggerType]->Fill(newCluster->E()); | |
2349 | if(fTriggerType==1) | |
2350 | { | |
2351 | if(fPeriod.Contains("lhc12a15a",TString::kIgnoreCase)) fhClsE[0]->Fill(newCluster->E()); | |
2352 | if(newCluster->Chi2()>0.5) fhClsE[fTriggerType]->Fill(newCluster->E()); | |
2353 | } | |
b6e50056 | 2354 | } |
2355 | ||
7c915297 | 2356 | // Cluster energy scale systematic uncertainty |
b6e50056 | 2357 | if(fSysClusterEScale) |
2358 | { | |
2359 | Double_t newE = newCluster->E() * (1+fVaryClusterEScale); | |
2360 | newCluster->SetE(newE); | |
2361 | } | |
2362 | ||
7c915297 | 2363 | // Cluster energy resolution systematic uncertainty |
b6e50056 | 2364 | if(fSysClusterERes) |
7c915297 | 2365 | { |
2366 | Double_t oldE = newCluster->E(); | |
2367 | Double_t resolution = fClusterEResolution->Eval(oldE); | |
2368 | Double_t smear = resolution * TMath::Sqrt((1+fVaryClusterERes)*(1+fVaryClusterERes)-1); | |
2369 | Double_t newE = oldE + fRandomGen->Gaus(0, smear); | |
2370 | newCluster->SetE(newE); | |
2371 | } | |
2372 | ||
2373 | // non-linearity systematic uncertainty | |
2374 | if(fSysNonLinearity) | |
b6e50056 | 2375 | { |
2376 | Double_t oldE = newCluster->E(); | |
7c915297 | 2377 | Double_t newE = oldE/fNonLinear->Eval(oldE) * 1.012; |
b6e50056 | 2378 | newCluster->SetE(newE); |
2379 | } | |
2380 | ||
7c915297 | 2381 | // clusterizer systematic uncertainty |
2382 | if(fSysClusterizer) | |
2383 | { | |
2384 | fhSysClusterE[fTriggerType][0]->Fill(newCluster->E()); | |
2385 | fhSysNCellVsClsE[fTriggerType][0]->Fill(newCluster->E(),newCluster->GetNCells()); | |
2386 | } | |
b6e50056 | 2387 | |
7c915297 | 2388 | // hadronic correction |
b6e50056 | 2389 | Double_t clsE = newCluster->E(); |
2390 | Double_t subE = 0., eRes = 0., mcSubE = 0, MIPE = 0.; | |
2391 | if(fElectronRejection || fHadronicCorrection) | |
2392 | subE= SubtractClusterEnergy(newCluster, eRes, MIPE, mcSubE); | |
2393 | ||
2394 | if(!fStudySubEInHC && !fStudyMcOverSubE) clsE -= subE; | |
2395 | if(clsE<0) {delete newCluster; continue;} | |
2396 | newCluster->SetE(clsE); | |
2397 | newCluster->SetTOF(eRes); | |
2398 | newCluster->SetDispersion(subE); | |
2399 | newCluster->SetEmcCpvDistance(MIPE); | |
2400 | newCluster->SetDistanceToBadChannel(mcSubE); | |
2401 | fClusterArray->Add(newCluster); | |
7c915297 | 2402 | |
2403 | // clusterizer systematic uncertainty | |
2404 | if(fSysClusterizer) | |
2405 | { | |
2406 | fhSysClusterE[fTriggerType][1]->Fill(newCluster->E()); | |
2407 | fhSysNCellVsClsE[fTriggerType][1]->Fill(newCluster->E(),newCluster->GetNCells()); | |
2408 | } | |
b6e50056 | 2409 | } |
2410 | ||
2411 | if(fVerbosity>5) printf("[i] # of EMCal clusters in event: %d\n", fClusterArray->GetEntries()); | |
2412 | ||
2413 | } | |
2414 | ||
2415 | ||
2416 | //________________________________________________________________________ | |
2417 | Bool_t AliAnalysisTaskFullppJet::IsGoodCluster(AliESDCaloCluster *cluster) | |
2418 | { | |
2419 | // | |
2420 | // Select good clusters | |
2421 | // | |
2422 | ||
2423 | if(!cluster) return kFALSE; | |
2424 | if (!cluster->IsEMCAL()) return kFALSE; | |
2425 | if(fRejectExoticCluster && fRecoUtil->IsExoticCluster(cluster, (AliVCaloCells*)fESD->GetEMCALCells())) return kFALSE; | |
2426 | if(fRemoveBadChannel && fRecoUtil->ClusterContainsBadChannel(fGeom, cluster->GetCellsAbsId(),cluster->GetNCells())) return kFALSE; | |
2427 | return kTRUE; | |
2428 | } | |
2429 | ||
2430 | ||
2431 | //________________________________________________________________________ | |
2432 | Double_t AliAnalysisTaskFullppJet::SubtractClusterEnergy(AliESDCaloCluster *cluster, Double_t &eRes, Double_t &MIPE, Double_t &mcSubE) | |
2433 | { | |
2434 | // | |
2435 | // Hadronic correction | |
2436 | // | |
2437 | ||
2438 | mcSubE = 0; | |
2439 | eRes = 0; | |
2440 | MIPE = 0; | |
2441 | ||
2442 | eRes += TMath::Power(fClusterEResolution->Eval(cluster->E()),2); | |
2443 | Double_t subE = 0., sumTrkPt = 0., sumTrkP = 0.; | |
2444 | Int_t nTrack = 0; | |
2445 | TArrayI *matched = cluster->GetTracksMatched(); | |
2446 | if(!matched) return 0; | |
2447 | for(Int_t im=0; im<matched->GetSize(); im++) | |
2448 | { | |
2449 | Int_t trkIndex = matched->At(im); | |
2450 | if(trkIndex<0 || trkIndex>=fESD->GetNumberOfTracks()) continue; | |
2451 | Bool_t isSelected = kFALSE; | |
2452 | Int_t index = -1; | |
2453 | for(Int_t j=0; j<fTrackArray->GetEntriesFast(); j++) | |
2454 | { | |
2455 | AliESDtrack *tr = (AliESDtrack*)fTrackArray->At(j); | |
2456 | if( trkIndex == tr->GetID() ) | |
2457 | { | |
2458 | isSelected = kTRUE; | |
2459 | index = j; | |
2460 | break; | |
2461 | } | |
2462 | } | |
2463 | if(!isSelected) continue; | |
2464 | nTrack++; | |
2465 | AliESDtrack *track = (AliESDtrack*)fTrackArray->At(index); | |
2466 | Double_t trkP = track->P(); | |
2467 | sumTrkPt += track->Pt(); sumTrkP += track->P(); | |
2468 | if(fSysTrkPtRes) trkP = TMath::Sqrt(TMath::Power(track->Pt()+GetSmearedTrackPt(track),2)+TMath::Power(track->Pz(),2)); | |
7c915297 | 2469 | if(IsElectron(track,cluster->E())) //electrons |
b6e50056 | 2470 | { |
2471 | if(fElectronRejection) | |
2472 | { | |
2473 | subE+= trkP; | |
2474 | MIPE+= trkP; | |
2475 | eRes += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2); | |
2476 | } | |
2477 | } | |
7c915297 | 2478 | else //hadrons |
b6e50056 | 2479 | { |
2480 | if(fHadronicCorrection) | |
2481 | { | |
7c915297 | 2482 | MIPE += (trkP>0.27)?0.27:trkP; //MIP correction |
b6e50056 | 2483 | if(fFractionHC>2) |
2484 | { | |
2485 | if(trkP>0.27) | |
2486 | { | |
2487 | subE+=0.27; | |
2488 | } | |
2489 | else | |
2490 | { | |
2491 | subE+=trkP; | |
2492 | eRes += TMath::Power(track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2); | |
2493 | } | |
2494 | } | |
2495 | else | |
2496 | { | |
2497 | if(trkP>fHCLowerPtCutMIP) subE += 0.27; | |
2498 | else subE+=fFractionHC*trkP; | |
2499 | eRes += TMath::Power(fFractionHC*track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()),2); | |
2500 | } | |
2501 | } | |
2502 | } | |
2503 | } | |
2504 | ||
2505 | if(fSaveQAHistos) fhNMatchedTrack[fTriggerType]->Fill(nTrack); | |
2506 | if(nTrack>0) | |
2507 | { | |
2508 | Double_t fraction[4] = {0.3,0.5,0.7,1.0}; | |
2509 | for(Int_t j=0; j<4; j++) | |
2510 | { | |
2511 | Double_t subETmp = sumTrkP*fraction[j]; | |
2512 | if(subETmp>cluster->E()) subETmp = cluster->E(); | |
2513 | if(fSaveQAHistos) fhSubEVsTrkPt[fTriggerType][j]->Fill(sumTrkPt,subETmp/sumTrkP); | |
2514 | } | |
2515 | } | |
2516 | ||
2517 | eRes = TMath::Sqrt(eRes); | |
2518 | ||
2519 | if(fIsMC && nTrack>0) | |
2520 | { | |
2521 | Double_t neutralE = 0; | |
2522 | TArrayI* labels = cluster->GetLabelsArray(); | |
b6e50056 | 2523 | if(labels) |
2524 | { | |
b6e50056 | 2525 | for(Int_t il=0; il<labels->GetSize(); il++) |
2526 | { | |
2527 | Int_t ipart = labels->At(il); | |
2528 | if(ipart>-1 && ipart<fMC->GetNumberOfTracks()) | |
2529 | { | |
2530 | AliVParticle* vParticle = fMC->GetTrack(ipart); | |
2531 | if(vParticle->Charge()==0) | |
2532 | { | |
2533 | neutralE += vParticle->E(); | |
2534 | } | |
2535 | } | |
2536 | } | |
2537 | } | |
2538 | mcSubE = cluster->E() - neutralE; | |
2539 | if(mcSubE<0) | |
2540 | mcSubE=0; | |
b6e50056 | 2541 | } |
2542 | ||
2543 | ||
2544 | return subE; | |
2545 | } | |
2546 | ||
2547 | ||
2548 | // | |
2549 | //________________________________________________________________________ | |
2550 | // | |
2551 | Double_t AliAnalysisTaskFullppJet::GetExoticEnergyFraction(AliESDCaloCluster *cluster) | |
2552 | { | |
2553 | // | |
2554 | // Exotic fraction: f_cross | |
7c915297 | 2555 | // Adpated from AliEMCalRecoUtils |
b6e50056 | 2556 | |
2557 | if(!cluster) return -1; | |
2558 | AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells(); | |
2559 | if(!cells) return -1; | |
2560 | ||
2561 | // Get highest energy tower | |
2562 | Int_t iSupMod = -1, absID = -1, ieta = -1, iphi = -1,iTower = -1, iIphi = -1, iIeta = -1; | |
2563 | Bool_t share = kFALSE; | |
2564 | fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, iSupMod, ieta, iphi, share); | |
2565 | fGeom->GetCellIndex(absID,iSupMod,iTower,iIphi,iIeta); | |
2566 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta); | |
2567 | ||
2568 | Int_t absID1 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi+1, ieta); | |
2569 | Int_t absID2 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi-1, ieta); | |
2570 | Int_t absID3 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi, ieta+1); | |
2571 | Int_t absID4 = fGeom-> GetAbsCellIdFromCellIndexes(iSupMod, iphi, ieta-1); | |
2572 | ||
2573 | Float_t ecell = 0, ecell1 = 0, ecell2 = 0, ecell3 = 0, ecell4 = 0; | |
2574 | Double_t tcell = 0, tcell1 = 0, tcell2 = 0, tcell3 = 0, tcell4 = 0; | |
2575 | Bool_t accept = 0, accept1 = 0, accept2 = 0, accept3 = 0, accept4 = 0; | |
2576 | const Int_t bc = 0; | |
2577 | ||
2578 | accept = fRecoUtil->AcceptCalibrateCell(absID, bc, ecell ,tcell ,cells); | |
2579 | ||
2580 | if(!accept) return -1; | |
2581 | ||
2582 | if(ecell < 0.5) return -1; | |
2583 | ||
2584 | accept1 = fRecoUtil->AcceptCalibrateCell(absID1,bc, ecell1,tcell1,cells); | |
2585 | accept2 = fRecoUtil->AcceptCalibrateCell(absID2,bc, ecell2,tcell2,cells); | |
2586 | accept3 = fRecoUtil->AcceptCalibrateCell(absID3,bc, ecell3,tcell3,cells); | |
2587 | accept4 = fRecoUtil->AcceptCalibrateCell(absID4,bc, ecell4,tcell4,cells); | |
2588 | ||
2589 | const Double_t exoticCellDiffTime = 1e6; | |
2590 | if(TMath::Abs(tcell-tcell1)*1.e9 > exoticCellDiffTime) ecell1 = 0 ; | |
2591 | if(TMath::Abs(tcell-tcell2)*1.e9 > exoticCellDiffTime) ecell2 = 0 ; | |
2592 | if(TMath::Abs(tcell-tcell3)*1.e9 > exoticCellDiffTime) ecell3 = 0 ; | |
2593 | if(TMath::Abs(tcell-tcell4)*1.e9 > exoticCellDiffTime) ecell4 = 0 ; | |
2594 | ||
2595 | Float_t eCross = ecell1+ecell2+ecell3+ecell4; | |
2596 | ||
2597 | return 1-eCross/ecell; | |
2598 | } | |
2599 | ||
2600 | ||
2601 | //________________________________________________________________________ | |
2602 | void AliAnalysisTaskFullppJet::GetMCInfo() | |
2603 | { | |
2604 | // | |
2605 | // Get # of trials per ESD event | |
2606 | // | |
2607 | ||
2608 | AliStack *stack = fMC->Stack(); | |
2609 | if(stack) | |
2610 | { | |
2611 | AliGenEventHeader *head = dynamic_cast<AliGenEventHeader*>(fMC->GenEventHeader()); | |
2612 | if (head == 0x0) | |
2613 | { | |
2614 | AliError("Could not get the event header"); | |
2615 | return; | |
2616 | } | |
2617 | ||
2618 | AliGenPythiaEventHeader *headPy = dynamic_cast<AliGenPythiaEventHeader*>(head); | |
2619 | if (headPy != 0x0) | |
2620 | { | |
2621 | if (headPy->Trials() > 0) | |
2622 | { | |
2623 | fhNTrials[0]->Fill(0.5,headPy->Trials()); | |
2624 | ||
2625 | } | |
2626 | } | |
2627 | } | |
2628 | } | |
2629 | ||
2630 | ||
2631 | ||
2632 | //________________________________________________________________________ | |
2633 | void AliAnalysisTaskFullppJet::Terminate(Option_t *) | |
2634 | { | |
2635 | // | |
2636 | // Called once at the end of the query | |
2637 | // | |
7c915297 | 2638 | Info("Terminate","Terminate"); |
2639 | AliAnalysisTaskSE::Terminate(); | |
b6e50056 | 2640 | |
2641 | } | |
2642 | ||
2643 | // | |
2644 | //________________________________________________________________________ | |
2645 | // | |
2646 | Int_t AliAnalysisTaskFullppJet::RunOfflineTrigger() | |
2647 | { | |
2648 | // | |
2649 | // Run trigger offline | |
2650 | // | |
2651 | ||
2652 | fIsEventTriggerBit = 0; | |
2653 | Int_t isTrigger = 0; | |
2654 | Int_t ncl = fESD->GetNumberOfCaloClusters(); | |
2655 | for(Int_t icl=0; icl<ncl; icl++) | |
2656 | { | |
7c915297 | 2657 | // Check every cluster |
b6e50056 | 2658 | AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl); |
2659 | if(!IsGoodCluster(cluster)) continue; | |
2660 | Double_t pro = GetOfflineTriggerProbability(cluster); | |
2661 | Double_t rand = fRandomGen->Uniform(); | |
2662 | if(rand<pro) | |
2663 | { | |
2664 | isTrigger = 1; | |
7c915297 | 2665 | fIsEventTriggerBit = 1; |
b6e50056 | 2666 | cluster->SetChi2(1); |
2667 | } | |
7c915297 | 2668 | else |
2669 | cluster->SetChi2(0); | |
b6e50056 | 2670 | } |
2671 | return isTrigger; | |
2672 | } | |
2673 | ||
2674 | ||
2675 | // | |
2676 | //________________________________________________________________________ | |
2677 | // | |
2678 | Double_t AliAnalysisTaskFullppJet::GetOfflineTriggerProbability(AliESDCaloCluster *cluster) | |
2679 | { | |
2680 | // | |
2681 | // Get the probablity of the given cluster to trigger the event | |
2682 | // | |
2683 | ||
2684 | Double_t pro = 0; | |
2685 | // Check the trigger mask | |
2686 | AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells(); | |
2687 | Int_t iSupMod = -1, absID = -1, ieta = -1, iphi = -1,iTower = -1, iIphi = -1, iIeta = -1; | |
2688 | Bool_t share = kFALSE; | |
7c915297 | 2689 | fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, iSupMod, ieta, iphi, share); // Get the position of the most energetic cell in the cluster |
b6e50056 | 2690 | fGeom->GetCellIndex(absID,iSupMod,iTower,iIphi,iIeta); |
2691 | fGeom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,iIphi, iIeta,iphi,ieta); | |
2692 | // convert co global phi eta | |
2693 | Int_t gphi = iphi + 24*(iSupMod/2); | |
2694 | Int_t geta = ieta + 48*(iSupMod%2); | |
2695 | // get corresponding FALTRO | |
2696 | Int_t fphi = gphi / 2; | |
2697 | Int_t feta = geta / 2; | |
7c915297 | 2698 | if(fCheckTriggerMask && fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5 && iSupMod>-1) // check the trigger mask |
b6e50056 | 2699 | { |
2700 | Double_t clsE = cluster->E(); | |
7c915297 | 2701 | if(fSysClusterEScale) clsE = clsE * (1+fVaryClusterEScale); // Used for systematic uncertainty. Not needed for regular analysis |
2702 | if(clsE>10) pro = 1; // Probability is 1 at high E | |
b6e50056 | 2703 | else |
2704 | { | |
2705 | Int_t bin = fTriggerCurve[iSupMod]->FindFixBin(clsE); | |
7c915297 | 2706 | pro = fTriggerCurve[iSupMod]->GetBinContent(bin)/fTriggerEfficiency[iSupMod]->Eval(10); // Read the probability from trigger turn-on curves |
b6e50056 | 2707 | } |
2708 | } | |
2709 | return pro; | |
2710 | } | |
2711 | ||
2712 | ||
2713 | ||
2714 | // | |
2715 | //________________________________________________________________________ | |
2716 | // | |
2717 | Int_t AliAnalysisTaskFullppJet::GetClusterSuperModule(AliESDCaloCluster *cluster) | |
2718 | { | |
2719 | // | |
2720 | // Return the given cluster supermodule | |
2721 | // | |
2722 | ||
2723 | Float_t pos[3]; | |
2724 | cluster->GetPosition(pos); | |
2725 | TVector3 clsVec(pos[0],pos[1],pos[2]); | |
2726 | ||
2727 | Int_t sMod=-1; | |
2728 | fGeom->SuperModuleNumberFromEtaPhi(clsVec.Eta(), clsVec.Phi(), sMod); | |
2729 | return sMod; | |
2730 | } | |
2731 | ||
2732 | ||
2733 | //________________________________________________________________________ | |
2734 | Bool_t AliAnalysisTaskFullppJet::IsElectron(AliESDtrack *track, Double_t clsE) const | |
2735 | { | |
2736 | // | |
2737 | // Check if the given track is an electron candidate based on de/dx | |
2738 | // | |
2739 | ||
2740 | if(track->GetTPCsignal()<=fdEdxMax && track->GetTPCsignal()>=fdEdxMin && (clsE/track->P())<=fEoverPMax && (clsE/track->P())>=fEoverPMin ) | |
2741 | return kTRUE; | |
2742 | else | |
2743 | return kFALSE; | |
2744 | } | |
2745 | ||
2746 | ||
2747 | //________________________________________________________________________ | |
2748 | Double_t AliAnalysisTaskFullppJet::GetTrkEff(Double_t inPt) | |
2749 | { | |
2750 | // | |
2751 | // Get tracking efficiency estimated from simulation | |
2752 | // | |
2753 | ||
2754 | Double_t eff = 1; | |
2755 | Double_t ptBound[4] = {0, 0.5, 3.8, 300}; | |
2756 | ||
2757 | for(Int_t i=0; i<3; i++) | |
2758 | { | |
2759 | if( inPt < ptBound[i+1] && inPt >= ptBound[i]) | |
2760 | { | |
2761 | eff = fTrkEffFunc[i]->Eval(inPt); | |
2762 | break; | |
2763 | } | |
2764 | } | |
2765 | return eff; | |
2766 | } | |
2767 | ||
2768 | //________________________________________________________________________ | |
2769 | Double_t AliAnalysisTaskFullppJet::GetSmearedTrackPt(AliESDtrack *track) | |
2770 | { | |
2771 | // | |
2772 | // Smear track momentum | |
2773 | // | |
2774 | ||
7c915297 | 2775 | Double_t resolution = track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2()); |
2776 | Double_t smear = resolution*TMath::Sqrt((1+fVaryTrkPtRes)*(1+fVaryTrkPtRes)-1); | |
b6e50056 | 2777 | return fRandomGen->Gaus(0, smear); |
b6e50056 | 2778 | |
b6e50056 | 2779 | } |
2780 | ||
2781 | ||
2782 | //________________________________________________________________________ | |
2783 | void AliAnalysisTaskFullppJet::CheckEventTriggerBit() | |
2784 | { | |
2785 | // | |
2786 | // Check if the triggered events have correct trigger bit | |
2787 | // | |
2788 | ||
2789 | fIsEventTriggerBit = 0; | |
2790 | // constants | |
2791 | const Int_t nColsModule = 2; | |
2792 | const Int_t nRowsModule = 5; | |
2793 | const Int_t nRowsFeeModule = 24; | |
2794 | const Int_t nColsFeeModule = 48; | |
2795 | const Int_t nColsFaltroModule = 24; | |
2796 | const Int_t nRowsFaltroModule = 12; | |
2797 | ||
2798 | // part 1, trigger extraction ------------------------------------- | |
2799 | Int_t trigtimes[30], globCol, globRow, ntimes; | |
2800 | Int_t trigger[nColsFaltroModule*nColsModule][nRowsFaltroModule*nRowsModule]; | |
2801 | ||
2802 | // erase trigger maps | |
2803 | for( Int_t i = 0; i < nColsFaltroModule*nColsModule; i++ ) | |
2804 | { | |
2805 | for( Int_t j = 0; j < nRowsFaltroModule*nRowsModule; j++ ) | |
2806 | { | |
2807 | trigger[i][j] = 0; | |
2808 | } | |
2809 | } | |
2810 | ||
2811 | Int_t fTrigCutLow = 7, fTrigCutHigh = 10, trigInCut = 0; | |
2812 | AliESDCaloTrigger * fCaloTrigger = fESD->GetCaloTrigger( "EMCAL" ); | |
2813 | // go through triggers | |
2814 | if( fCaloTrigger->GetEntries() > 0 ) | |
2815 | { | |
2816 | // needs reset | |
2817 | fCaloTrigger->Reset(); | |
2818 | while( fCaloTrigger->Next() ) | |
2819 | { | |
2820 | fCaloTrigger->GetPosition( globCol, globRow ); // get position in global 2x2 tower coordinates | |
2821 | fCaloTrigger->GetNL0Times( ntimes ); // get dimension of time arrays | |
2822 | if( ntimes < 1 ) continue; // no L0s in this channel. Presence of the channel in the iterator still does not guarantee that L0 was produced!! | |
2823 | fCaloTrigger->GetL0Times( trigtimes ); // get timing array | |
7c915297 | 2824 | if(fCheckTriggerMask && fTriggerMask->GetBinContent(globCol+1,globRow+1)<0.5) continue; |
b6e50056 | 2825 | trigInCut = 0; |
b6e50056 | 2826 | for( Int_t i = 0; i < ntimes; i++ ) |
2827 | { | |
b6e50056 | 2828 | if( trigtimes[i] > fTrigCutLow && trigtimes[i] < fTrigCutHigh ) // check if in cut |
2829 | { | |
2830 | trigInCut = 1; | |
2831 | } | |
2832 | } | |
2833 | if(trigInCut==1) trigger[globCol][globRow] = 1; | |
2834 | } // calo trigger entries | |
2835 | } // has calo trigger entries | |
2836 | ||
2837 | ||
2838 | ||
2839 | // part 2 go through the clusters here ----------------------------------- | |
2840 | Int_t nCell, iCell; | |
2841 | UShort_t *cellAddrs; | |
2842 | Int_t absID = -1, nSupMod=-1, nModule=-1, nIphi=-1, nIeta=-1, iphi=-1, ieta=-1, gphi=-1, geta=-1, feta=-1, fphi=-1; | |
2843 | Bool_t share = kFALSE; | |
2844 | AliVCaloCells *cells = (AliVCaloCells*)fESD->GetEMCALCells(); | |
2845 | for(Int_t icl=0; icl<fESD->GetNumberOfCaloClusters(); icl++) | |
2846 | { | |
2847 | AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl); | |
7c915297 | 2848 | if(!cluster || !cluster->IsEMCAL()) continue; |
b6e50056 | 2849 | cluster->SetChi2(0); |
7c915297 | 2850 | if(!IsGoodCluster(cluster)) continue; |
2851 | ||
b6e50056 | 2852 | //Clusters with most energetic cell in dead region can't be triggered |
2853 | fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, nSupMod, ieta, iphi, share); | |
2854 | fGeom->GetCellIndex(absID,nSupMod,nModule, nIphi, nIeta); | |
2855 | fGeom->GetCellPhiEtaIndexInSModule(nSupMod,nModule, nIphi, nIeta, iphi, ieta); | |
2856 | gphi = iphi + nRowsFeeModule*(nSupMod/2); | |
2857 | geta = ieta + nColsFeeModule*(nSupMod%2); | |
2858 | fphi = gphi / 2; | |
2859 | feta = geta / 2; | |
2860 | ||
7c915297 | 2861 | if(fCheckTriggerMask && fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5) |
b6e50056 | 2862 | { |
2863 | nCell = cluster->GetNCells(); // get cluster cells | |
2864 | cellAddrs = cluster->GetCellsAbsId(); // get the cell addresses | |
2865 | for( iCell = 0; iCell < nCell; iCell++ ) | |
2866 | { | |
2867 | // get cell position | |
2868 | fGeom->GetCellIndex( cellAddrs[iCell], nSupMod, nModule, nIphi, nIeta ); | |
2869 | fGeom->GetCellPhiEtaIndexInSModule( nSupMod,nModule, nIphi, nIeta, iphi, ieta); | |
2870 | ||
2871 | // convert co global phi eta | |
2872 | gphi = iphi + nRowsFeeModule*(nSupMod/2); | |
2873 | geta = ieta + nColsFeeModule*(nSupMod%2); | |
2874 | ||
2875 | // get corresponding FALTRO | |
2876 | fphi = gphi / 2; | |
2877 | feta = geta / 2; | |
2878 | // try to match with a triggered | |
2879 | if( trigger[feta][fphi] == 1) | |
2880 | { | |
2881 | cluster->SetChi2(1); | |
2882 | fIsEventTriggerBit = 1; | |
2883 | //break; | |
2884 | } | |
2885 | } // cells | |
2886 | } | |
2887 | } // clusters | |
2888 | } | |
2889 | ||
2890 | ||
2891 | //________________________________________________________________________ | |
2892 | void AliAnalysisTaskFullppJet::PrintConfig() | |
2893 | { | |
2894 | // | |
2895 | // Print configuration | |
2896 | // | |
2897 | ||
2898 | const char *trackCutName[3] = {"Hybrid tracks","TPCOnly tracks","Golden tracks"}; | |
2899 | const char *triggerType[2] = {"MB","EMC"}; | |
2900 | const char *decision[2] = {"no","yes"}; | |
2901 | const char *recombination[] = {"E_scheme","pt_scheme","pt2_scheme","Et_scheme","Et2_scheme","BIpt_scheme","BIpt2_scheme"}; | |
7c915297 | 2902 | const char *type[2] = {"local","grid"}; |
b6e50056 | 2903 | if(fStudySubEInHC || fStudyMcOverSubE) |
2904 | { | |
2905 | printf("\n\n=================================\n"); | |
2906 | printf("======WARNING: HC is ingored!======\n"); | |
2907 | printf("====== NOT for PHYSICS! ======\n\n"); | |
2908 | } | |
2909 | printf("Run period: %s\n",fPeriod.Data()); | |
7c915297 | 2910 | printf("Reject SPD pileup: %s\n",decision[fRejectPileup]); |
2911 | printf("Reject exotic triggered events: %s\n",decision[fRejectExoticTrigger]); | |
b6e50056 | 2912 | printf("Is this MC data: %s\n",decision[fIsMC]); |
2913 | printf("Only find charged jets in MC: %s\n",decision[fChargedMC]); | |
7c915297 | 2914 | printf("Analyze on local or grid? %s\n",type[fAnaType]); |
b6e50056 | 2915 | printf("Run offline trigger on MC: %s\n",decision[fOfflineTrigger]); |
2916 | if(fIsMC) | |
2917 | printf("Is K0 and n included: %s\n",decision[1-fRejectNK]); | |
2918 | printf("Constrain tracks in EMCal acceptance: %s\n",decision[fConstrainChInEMCal]); | |
2919 | printf("Track selection: %s, |eta| < %2.1f\n",trackCutName[fTrackCutsType], fTrkEtaMax); | |
2920 | for(Int_t i=0; i<2; i++) | |
2921 | { | |
2922 | printf("Track pt cut: %s -> %2.2f < pT < %2.1f\n",triggerType[i], fTrkPtMin[i], fTrkPtMax[i]); | |
2923 | } | |
2924 | for(Int_t i=0; i<2; i++) | |
2925 | { | |
2926 | printf("Cluster selection: %s -> %2.2f < Et < %2.1f\n",triggerType[i],fClsEtMin[i], fClsEtMax[i]); | |
2927 | } | |
2928 | printf("Electron selectoin: %2.0f < dE/dx < %2.0f, %1.1f < E/P < %1.1f\n",fdEdxMin, fdEdxMax, fEoverPMin, fEoverPMax); | |
2929 | printf("Reject exotic cluster: %s\n",decision[fRejectExoticCluster]); | |
2930 | printf("Remove problematic region in SM4: %s\n",decision[fRemoveBadChannel]); | |
2931 | printf("Use only good SM (1,2,6,7,8,9) for trigger: %s\n",decision[fUseGoodSM]); | |
2932 | printf("Reject electron: %s\n", decision[fElectronRejection]); | |
2933 | printf("Correct hadron: %s\n",decision[fHadronicCorrection]); | |
2934 | printf("HC fraction: %2.1f up to %2.0f GeV/c\n",fFractionHC,fHCLowerPtCutMIP); | |
2935 | printf("Find charged jets: %s\n", decision[fFindChargedOnlyJet]); | |
2936 | printf("Find netural jets: %s\n", decision[fFindNeutralOnlyJet]); | |
2937 | printf("Find good jets: %s\n",decision[fSpotGoodJet]); | |
2938 | printf("Jet radius: %s\n",fRadius.Data()); | |
2939 | printf("Jet recombination scheme: %s\n",recombination[fRecombinationScheme]); | |
2940 | printf("Correct tracking efficiency: %s\n",decision[fCheckTrkEffCorr]); | |
2941 | printf("Save jet QA histos: %s\n",decision[fSaveQAHistos]); | |
2942 | printf("Systematics: jet efficiency: %s with variation %1.0f%%\n",decision[fSysJetTrigEff],fVaryJetTrigEff*100); | |
2943 | printf("Systematics: tracking efficiency: %s with variation %1.0f%%\n",decision[fSysTrkEff],fVaryTrkEff*100); | |
2944 | printf("Systematics: track pt resolution: %s with variation %1.0f%%\n",decision[fSysTrkPtRes],fVaryTrkPtRes*100); | |
2945 | printf("Systematics: track-cluster matching: %s with |dEta|<%2.3f, |dPhi|<%2.3f\n",decision[fSysTrkClsMth],fCutdEta,fCutdPhi); | |
2946 | printf("Systematics: EMCal non-linearity: %s\n",decision[fSysNonLinearity]); | |
2947 | printf("Systematics: EMCal energy scale: %s with uncertainty of %1.0f%%\n",decision[fSysClusterEScale],fVaryClusterEScale*100); | |
2948 | printf("Systematics: EMCal energy resolution: %s with uncertainty of %1.0f%%\n",decision[fSysClusterERes],fVaryClusterERes*100); | |
2949 | printf("Smear lhc12a15a: %s\n",decision[fSmearMC]); | |
7c915297 | 2950 | printf("Run UE analysis: %s\n",decision[fRunUE]); |
2951 | printf("Run secondaries: %s\n",decision[fRunSecondaries]); | |
b6e50056 | 2952 | printf("=======================================\n\n"); |
2953 | } | |
2954 | ||
2955 | //________________________________________________________________________ | |
2956 | void AliAnalysisTaskFullppJet::CheckExoticEvent() | |
2957 | { | |
2958 | // | |
2959 | // Check if the event containts exotic clusters | |
2960 | // | |
7c915297 | 2961 | |
b6e50056 | 2962 | Double_t leadingE = 0; |
5543d2a2 | 2963 | Int_t nCls = fESD->GetNumberOfCaloClusters(); |
2964 | for(Int_t icl=0; icl<nCls; icl++) | |
b6e50056 | 2965 | { |
2966 | AliESDCaloCluster *cluster = fESD->GetCaloCluster(icl); | |
2967 | if(!cluster) continue; | |
2968 | if (!cluster->IsEMCAL()) continue; | |
2969 | if(fRecoUtil->IsExoticCluster(cluster, (AliVCaloCells*)fESD->GetEMCALCells()) && cluster->E()>leadingE) | |
2970 | { | |
2971 | leadingE = cluster->E(); | |
2972 | } | |
2973 | } | |
2974 | if(leadingE>3) fIsExoticEvent3GeV = kTRUE; | |
2975 | if(leadingE>5) fIsExoticEvent5GeV = kTRUE; | |
2976 | } | |
2977 | ||
7c915297 | 2978 | |
b6e50056 | 2979 | //________________________________________________________________________ |
7c915297 | 2980 | Bool_t AliAnalysisTaskFullppJet::HasPrimaryVertex() const |
b6e50056 | 2981 | { |
2982 | // | |
7c915297 | 2983 | // Check if the primary vertex exists |
b6e50056 | 2984 | // |
b6e50056 | 2985 | const AliESDVertex* vtx = fESD->GetPrimaryVertex(); |
2986 | if (!vtx || vtx->GetNContributors()<1) return kFALSE; | |
7c915297 | 2987 | else return kTRUE; |
2988 | } | |
b6e50056 | 2989 | |
b6e50056 | 2990 | |
7c915297 | 2991 | //________________________________________________________________________ |
2992 | Bool_t AliAnalysisTaskFullppJet::IsPrimaryVertexOk() const | |
2993 | { | |
2994 | // | |
2995 | // Check if the event vertex is good | |
2996 | // | |
b6e50056 | 2997 | return kTRUE; |
2998 | } | |
2999 | ||
3000 | //________________________________________________________________________ | |
3001 | Bool_t AliAnalysisTaskFullppJet::IsTPCOnlyVtx() const | |
3002 | { | |
3003 | // | |
3004 | // Check if the event only has valid TPC vertex | |
3005 | // | |
3006 | ||
7c915297 | 3007 | const Bool_t isPrimaryVtx = HasPrimaryVertex(); |
3008 | ||
3009 | Bool_t isGoodVtx = kTRUE; | |
3010 | AliESDVertex* goodvtx = const_cast<AliESDVertex*>(fESD->GetPrimaryVertexTracks()); | |
3011 | if(!goodvtx || goodvtx->GetNContributors()<1) | |
3012 | goodvtx = const_cast<AliESDVertex*>(fESD->GetPrimaryVertexSPD()); // SPD vertex | |
3013 | if(!goodvtx || !goodvtx->GetStatus()) isGoodVtx = kFALSE; // Event rejected | |
3014 | ||
3015 | if( isPrimaryVtx && !isGoodVtx ) | |
b6e50056 | 3016 | return kTRUE; |
3017 | else | |
3018 | return kFALSE; | |
3019 | } | |
3020 | ||
7c915297 | 3021 | //________________________________________________________________________ |
3022 | void AliAnalysisTaskFullppJet::CheckTPCOnlyVtx(const UInt_t trigger) | |
3023 | { | |
3024 | // | |
3025 | // Check the fraction of accepted events that have only TPC vertex | |
3026 | // | |
3027 | Int_t lTriggerType = -1; | |
3028 | if (trigger & AliVEvent::kMB) lTriggerType = 1; | |
3029 | if (trigger & AliVEvent::kFastOnly) lTriggerType = 0; | |
3030 | if (trigger & AliVEvent::kEMC1) lTriggerType = 2; | |
3031 | if(lTriggerType==-1) return; | |
3032 | fhEventStatTPCVtx->Fill(0.5+lTriggerType*3); | |
3033 | if(HasPrimaryVertex()) fhEventStatTPCVtx->Fill(1.5+lTriggerType*3); | |
3034 | if(IsTPCOnlyVtx()) fhEventStatTPCVtx->Fill(2.5+lTriggerType*3); | |
3035 | } | |
3036 | ||
b6e50056 | 3037 | //________________________________________________________________________ |
3038 | Bool_t AliAnalysisTaskFullppJet::IsLEDEvent() const | |
3039 | { | |
3040 | // | |
3041 | // Check if the event is contaminated by LED signal | |
3042 | // | |
3043 | ||
3044 | AliESDCaloCells *cells = fESD->GetEMCALCells(); | |
3045 | Short_t nCells = cells->GetNumberOfCells(); | |
3046 | Int_t nCellCount[2] = {0,0}; | |
3047 | for(Int_t iCell=0; iCell<nCells; iCell++) | |
3048 | { | |
3049 | Int_t cellId = cells->GetCellNumber(iCell); | |
3050 | Double_t cellE = cells->GetCellAmplitude(cellId); | |
3051 | Int_t sMod = fGeom->GetSuperModuleNumber(cellId); | |
3052 | ||
3053 | if(sMod==3 || sMod==4) | |
3054 | { | |
3055 | if(cellE>0.1) | |
3056 | nCellCount[sMod-3]++; | |
3057 | } | |
3058 | } | |
3059 | Bool_t isLED=kFALSE; | |
3060 | ||
3061 | if(fPeriod.CompareTo("lhc11a")==0) | |
3062 | { | |
3063 | if (nCellCount[1] > 100) | |
3064 | isLED = kTRUE; | |
3065 | Int_t runN = fESD->GetRunNumber(); | |
3066 | if (runN>=146858 && runN<=146860) | |
3067 | { | |
3068 | if(fTriggerType==0 && nCellCount[0]>=21) isLED=kTRUE; | |
3069 | if(fTriggerType==1 && nCellCount[0]>=35) isLED=kTRUE; | |
3070 | } | |
3071 | } | |
3072 | ||
3073 | return isLED; | |
3074 | } |