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