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