]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskFullppJet.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskFullppJet.cxx
... / ...
CommitLineData
1// **************************************
2// Full pp jet task - ESD input only
3// Extract the jet spectrum and all the
4// systematic uncertainties
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"
39#include "AliMCParticle.h"
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;
52using std::vector;
53
54ClassImp(AliAnalysisTaskFullppJet)
55
56const Float_t kRadius[3] = {0.4,0.2,0.3};
57const Double_t kPI = TMath::Pi();
58const Double_t kdRCut[3] = {0.25,0.1,0.15};
59
60//________________________________________________________________________
61AliAnalysisTaskFullppJet::AliAnalysisTaskFullppJet() :
62 AliAnalysisTaskSE("default"),
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.),
66 fCentrality(99), fZVtxMax(10),
67 fTriggerType(-1), fCheckTriggerMask(kTRUE), fIsTPCOnlyVtx(kFALSE),
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),
74 fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0),
75 fJetNEFMin(0.01), fJetNEFMax(0.98),
76 fSpotGoodJet(kFALSE), fFindChargedOnlyJet(kFALSE), fFindNeutralOnlyJet(kFALSE),
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.),
80 fSysTrkClsMth(kFALSE), fCutdEta(0.05), fCutdPhi(0.05),
81 fSysNonLinearity(kFALSE), fNonLinear(0x0), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0),
82 fSysClusterizer(0),
83 fNonStdBranch(""), fNonStdFile(""), fAlgorithm("aKt"), fRadius("0.4"), fRecombinationScheme(5),
84 fConstrainChInEMCal(kFALSE), fRejectNK(kFALSE), fRejectWD(kFALSE), fSmearMC(kFALSE), fTrkPtResData(0x0),
85 fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventCount(0x0), fhJetEventStat(0x0), fhEventStatTPCVtx(0x0), fhChunkQA(0x0)
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;
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 }
111 }
112
113 for(Int_t r=0; r<kNRadii; r++)
114 {
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 }
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)
140 {
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 }
170 }
171 }
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 }
185 }
186
187 for(Int_t s=0; s<3; s++)
188 {
189 for(Int_t a=0; a<2; a++)
190 {
191 for(Int_t r=0; r<kNRadii; r++)
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;
207 fhSecondaryResponse[j] = 0x0;
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),
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.),
223 fCentrality(99), fZVtxMax(10),
224 fTriggerType(-1), fCheckTriggerMask(kTRUE), fIsTPCOnlyVtx(kFALSE),
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),
231 fElectronRejection(kFALSE), fHadronicCorrection(kFALSE), fFractionHC(0.), fHCLowerPtCutMIP(1e4), fClusterEResolution(0x0),
232 fJetNEFMin(0.01), fJetNEFMax(0.98),
233 fSpotGoodJet(kFALSE), fFindChargedOnlyJet(kFALSE), fFindNeutralOnlyJet(kFALSE),
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.),
237 fSysTrkClsMth(kFALSE), fCutdEta(0.05), fCutdPhi(0.05),
238 fSysNonLinearity(kFALSE), fNonLinear(0x0), fSysClusterEScale(kFALSE), fVaryClusterEScale(0), fSysClusterERes(kFALSE), fVaryClusterERes(0),
239 fSysClusterizer(0),
240 fNonStdBranch(""), fNonStdFile(""), fAlgorithm("aKt"), fRadius("0.4"), fRecombinationScheme(5),
241 fConstrainChInEMCal(kFALSE), fRejectNK(kFALSE), fRejectWD(kFALSE), fSmearMC(kFALSE), fTrkPtResData(0x0),
242 fOutputList(0x0), fSaveQAHistos(kTRUE), fhJetEventCount(0x0), fhJetEventStat(0x0), fhEventStatTPCVtx(0x0), fhChunkQA(0x0)
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;
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 }
268 }
269
270 for(Int_t r=0; r<kNRadii; r++)
271 {
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 }
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)
297 {
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 }
327 }
328 }
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 }
342 }
343
344 for(Int_t s=0; s<3; s++)
345 {
346 for(Int_t a=0; a<2; a++)
347 {
348 for(Int_t r=0; r<kNRadii; r++)
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 }
360
361 for(Int_t j=0; j<3; j++)
362 {
363 fTrkEffFunc[j] = 0x0;
364 fhSecondaryResponse[j] = 0x0;
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 {
387 for(Int_t r=0; r<kNRadii; r++)
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];
403 if(fhSecondaryResponse[i]) delete fhSecondaryResponse[i];
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 }
410 for(Int_t r=0; r<kNRadii; r++)
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;
427 if(fNonLinear) delete fNonLinear;
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 //
450
451 if(fRunUE) fFindChargedOnlyJet = kTRUE;
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
465 fhJetEventStat = new TH1F("fhJetEventStat","Event statistics for jet analysis",12,0,12);
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");
475 fhJetEventStat->GetXaxis()->SetBinLabel(10,"MB+TVtx+10cm");
476 fhJetEventStat->GetXaxis()->SetBinLabel(11,"EMC+TVtx+10cm");
477 fhJetEventStat->GetXaxis()->SetBinLabel(12,"ALL-Pileup");
478 fOutputList->Add(fhJetEventStat);
479
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 }
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
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"};
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"};
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
556 for(Int_t r=0; r<kNRadii; r++)
557 {
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
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
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);
566 fOutputList->Add(fhNeutralPtInJet[i][r]);
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]);
570
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);
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
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 }
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 }
649
650 if(fRunUE)
651 {
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 }
666 }
667
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 }
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 }
685 }
686
687
688 if(fIsMC)
689 {
690 if(fStudyMcOverSubE)
691 {
692 for(Int_t r=0; r<kNRadii; r++)
693 {
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;
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 }
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 }
744 }
745
746 printf("\n=======================================\n");
747 printf("===== Jet task configuration ==========\n");
748
749 if(fNonStdBranch.Length()!=0)
750 {
751 const char* species[3] = {"in","ch","ne"};
752 const char* algorithm[2] = {"akt","kt"};
753 const char* radii[kNRadii] = {"04","02","03"};
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;
762 for(Int_t r=0; r<kNRadii; r++)
763 {
764 if(!fRadius.Contains("0.4") && r==0) continue;
765 if(!fRadius.Contains("0.2") && r==1) continue;
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 }
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 {
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 }
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);
809 if(fCheckTrkEffCorr && fAnaType==0)
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;
825 for(Int_t r=0; r<kNRadii; r++)
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
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 }
847
848 if(fCheckTriggerMask)
849 {
850 TString name = "TriggerCurve.root";
851 if(fAnaType==0) name = "/project/projectdirs/alice/marr/Analysis/2.76Jet/CorrectionFunctions/TriggerCurve.root";
852 TFile f (name.Data(),"read");
853 fTriggerMask = new TH2F(*((TH2F*)f.Get("lhc11a_TriggerMask")));
854 if(fOfflineTrigger)
855 {
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();
863 }
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
868 fGeom = AliEMCALGeometry::GetInstance("EMCAL_COMPLETEV1");
869 if(!fGeom)
870 {
871 AliError("No EMCal geometry is available!");
872 return;
873 }
874 fRecoUtil = new AliEMCALRecoUtils();
875 if(fRejectExoticCluster || fRejectExoticTrigger)
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 {
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);
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();
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
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 {
956 printf("ERROR: Could not retrieve MC event");
957 return;
958 }
959 fStack = fMC->Stack();
960 TParticle *particle = fStack->Particle(0);
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 }
966 }
967
968 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
969 if (!fESD)
970 {
971 AliError("fESD is not available");
972 return;
973 }
974
975 fhJetEventStat->Fill(0.5);
976 fhJetEventCount->Fill(0.5);
977 if(fIsMC)
978 {
979 GetMCInfo();
980 if(fVertexGenZ[0]) fVertexGenZ[0]->Fill(vtxTrueZ);
981 }
982
983 // Centrality, vertex, other event variables...
984 UInt_t trigger = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
985 if (trigger==0) return;
986 if(fCheckTPCOnlyVtx) CheckTPCOnlyVtx(trigger);
987 if(!fIsMC)
988 {
989 if (trigger & AliVEvent::kFastOnly) return; // Reject fast trigger cluster
990 if(fPeriod.CompareTo("lhc11a",TString::kIgnoreCase)==0)
991 {
992 if (trigger & AliVEvent::kMB) fTriggerType = 0;
993 else if(trigger & AliVEvent::kEMC1) fTriggerType = 1;
994 else fTriggerType = -1;
995 }
996 else if (fPeriod.CompareTo("lhc11c",TString::kIgnoreCase)==0 || fPeriod.CompareTo("lhc11d",TString::kIgnoreCase)==0)
997 {
998 if (trigger & AliVEvent::kINT7) fTriggerType = 0;
999 else if(trigger & AliVEvent::kEMC7) fTriggerType = 1;
1000 else fTriggerType = -1;
1001 }
1002 else
1003 {
1004 return;
1005 }
1006 }
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 }
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
1034 if (IsLEDEvent())
1035 {
1036 fhJetEventStat->Fill(8.5);
1037 return;
1038 }
1039
1040 fhJetEventStat->Fill(1.5+fTriggerType*2);
1041
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
1057 fIsTPCOnlyVtx = IsTPCOnlyVtx();
1058
1059 if(!fIsMC && fTriggerType==1)
1060 {
1061 // Check if event has valid trigger bit
1062 CheckEventTriggerBit();
1063 if(fIsEventTriggerBit)
1064 {
1065 fhJetEventStat->Fill(7.5);
1066 fhJetEventCount->Fill(9.5);
1067 }
1068 else return;
1069 }
1070 fhJetEventStat->Fill(2.5+fTriggerType*2);
1071 if(fIsTPCOnlyVtx) fhJetEventStat->Fill(9.5+fTriggerType);
1072
1073 // Check if event contains exotic clusters
1074 CheckExoticEvent();
1075
1076 // Write jet tree into AOD output in local analysis mode
1077 if(fAnaType==0) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
1078
1079 // Clean up arrays
1080 for(Int_t s=0; s<3; s++)
1081 {
1082 for(Int_t a=0; a<2; a++)
1083 {
1084 for(Int_t r=0; r<kNRadii; r++)
1085 {
1086 if(fJetTCA[s][a][r]) fJetTCA[s][a][r]->Delete();
1087 if(fDetJetFinder[s][a][r]) fDetJetFinder[s][a][r]->Clear();
1088
1089 if(s==0 && a==0)
1090 {
1091 if(fMcTruthAntikt[r]) fMcTruthAntikt[r]->Delete();
1092 if(fTrueJetFinder[r]) fTrueJetFinder[r]->Clear();
1093 }
1094 }
1095 }
1096 }
1097
1098 if(fVerbosity>5) printf("# of jets after clear: %d\n",fJetTCA[0][0][0]->GetEntries());
1099
1100 if(fTrackArray) fTrackArray->Delete();
1101 if(fClusterArray) fClusterArray->Delete();
1102 fMcPartArray->Reset();
1103
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
1131
1132 if(fIsMC)
1133 {
1134 for(Int_t r=0; r<kNRadii; r++)
1135 if(fTrueJetFinder[r]) ProcessMC(r); // find particle level jets
1136 }
1137
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 {
1159 // Feed in charged tracks
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);
1165 if(!t) continue;
1166 countTracks++;
1167 Double_t e = t->P();
1168 Double_t pt = t->Pt();
1169 if(s==1 && fRunUE && pt<1) continue;
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);
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());
1179 }
1180 if(fVerbosity>5) printf("[i] # of tracks filled: %d\n",countTracks);
1181 }
1182 if(isNe)
1183 {
1184 // Feed in EMCal clusters
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 }
1203 // Run jet finding
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
1214 //
1215
1216 Int_t radiusIndex = 0;
1217 TString arrayName = fJetArray->GetName();
1218 if(arrayName.Contains("_02_")) radiusIndex = 1;
1219 if(arrayName.Contains("_03_")) radiusIndex = 2;
1220 std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1221 if(fVerbosity>5 && radiusIndex==0) printf("[i] # of jets in %s : %d\n",fJetArray->GetName(),(Int_t)jetsIncl.size());
1222 AliAODJet *jet = 0;
1223 Int_t jetCount = 0;
1224 for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1225 {
1226 if(fVerbosity>10) printf("fastjet: eta=%f, phi=%f\n",jetsIncl[ij].eta(),jetsIncl[ij].phi()*TMath::RadToDeg());
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++;
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());
1233
1234 jet->GetRefTracks()->Clear();
1235 std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
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;
1239 Int_t leadIndex = -1;
1240 for(UInt_t ic=0; ic<constituents.size(); ic++)
1241 {
1242 if(fVerbosity>10 && radiusIndex==0) printf("ic=%d: user_index=%d, E=%f\n",ic,constituents[ic].user_index(),constituents[ic].E());
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);
1256 if(cluster->Chi2()>0.5) isTriggering = kTRUE;
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)
1282 {
1283 leadPt = constituents[ic].perp();
1284 leadIndex = ic;
1285 }
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);
1294 if(constituents[leadIndex].user_index()>0) jet->SetTrigger(AliAnalysisTaskFullppJet::kLeadCh);
1295
1296 if(jetsIncl[ij].E()>0) jet->SetNEF(neE/jetsIncl[ij].E());
1297 jet->SetPtLeading(leadPt);
1298 jet->SetPtSubtracted(leadTrkType,0);
1299
1300 Double_t effAErrCh = 0, effAErrNe = leadTrkType;
1301 Double_t chBgEnergy = 10, neBgEnergy = nNePart;
1302 if(!isTruth)
1303 {
1304 effAErrCh = GetMeasuredJetPtResolution(ij,jetFinder);
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 }
1320 }
1321
1322 jet->SetEffArea(leadPt, jetFinder->GetJetArea(ij),effAErrCh,effAErrNe);
1323 jet->SetBgEnergy(chBgEnergy,neBgEnergy);
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;
1326
1327 if(fVerbosity>5)
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
1357 }
1358 if(fVerbosity>5) printf("%s has %d jets\n",fJetArray->GetName(), fJetArray->GetEntries());
1359}
1360
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
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;
1409 if(vParticle->Charge()==0) { index=-1; }
1410 fMcPartArray->AddAt(ipart, countPart);
1411 countPart++;
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());
1416 }
1417 fMcPartArray->Set(countPart);
1418 fTrueJetFinder[r]->Run();
1419
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 }
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 {
1463 if(type<2 && fCheckTPCOnlyVtx && fIsTPCOnlyVtx)
1464 {
1465 fhJetInTPCOnlyVtx[type][r]->Fill(jetsIncl[ij].perp(),jetsIncl[ij].phi()*TMath::RadToDeg(),jetsIncl[ij].eta());
1466 }
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
1473 std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(ij);
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);
1512 if(type==1 && fIsExoticEvent3GeV) fhJetPtInExoticEvent[0][r]->Fill(jetPt);
1513 if(type==1 && fIsExoticEvent5GeV) fhJetPtInExoticEvent[1][r]->Fill(jetPt);
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;
1537 if(type<2)
1538 {
1539 AliESDCaloCluster *cluster = (AliESDCaloCluster *)fClusterArray->At(constituents[ic].user_index()*(-1)-1);
1540 Double_t clsE = cluster->E();
1541 if(cluster->Chi2()>0.5) fhTrigNeuPtInJet[type][r]->Fill(jetPt, partPt);
1542 if(fStudySubEInHC || fStudyMcOverSubE)
1543 {
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 }
1555 }
1556 }
1557 }
1558 else
1559 {
1560 if(partPt>leadChPt) {leadChPt = partPt;leadChIndex=ic;}
1561 fhChargedPtInJet[type][r]->Fill(jetPt, partPt);
1562 chPt += constituents[ic].perp();
1563 if(type<2)
1564 {
1565 AliESDtrack *track = (AliESDtrack*) fTrackArray->At(constituents[ic].user_index()-1);
1566 Int_t trkType = (Int_t)track->GetTRDQuality();
1567 totTrkPt[trkType] += partPt;
1568 Int_t clusterIndex = track->GetEMCALcluster();
1569 Int_t clusterPos = -1;
1570 Bool_t isExist = kFALSE;
1571 for(Int_t j=0; j<fClusterArray->GetEntriesFast(); j++)
1572 {
1573 AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(j);
1574 if( clusterIndex == cluster->GetID() )
1575 {
1576 isExist = kTRUE;
1577 clusterPos = j;
1578 break;
1579 }
1580 }
1581 if(isExist)
1582 {
1583 AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(clusterPos);
1584 Double_t subEtmp = cluster->GetDispersion();
1585 Double_t mipETmp = cluster->GetEmcCpvDistance();
1586 for(Int_t k=0; k<5; k++)
1587 {
1588 if(k==0)
1589 {
1590 if(mipETmp>cluster->E()) subTrkPtClean[k] += partPt;
1591 else subTrkPtAmbig[k] += partPt;
1592 }
1593 else
1594 {
1595 if(fraction[k]*subEtmp>cluster->E()) subTrkPtClean[k] += partPt;
1596 else subTrkPtAmbig[k] += partPt;
1597 }
1598 }
1599 }
1600 }
1601 }
1602
1603 Float_t dR = TMath::Sqrt( (partEta-jetEta)*(partEta-jetEta) + (partPhi-jetPhi)*(partPhi-jetPhi) );
1604 for(Int_t i=0; i<nBins; i++)
1605 {
1606 if( dR < radiusCut[i] )
1607 {
1608 eFraction[i] += partE;
1609 nPart[i]++;
1610 }
1611 }
1612 }
1613
1614 fhLeadNePtInJet[type][r]->Fill(jetPt, leadNePt);
1615 fhLeadChPtInJet[type][r]->Fill(jetPt, leadChPt);
1616 if(type<2 && leadChIndex>-1)
1617 {
1618 Double_t cz = GetZ(constituents[leadChIndex].px(),constituents[leadChIndex].py(),constituents[leadChIndex].pz(),jetsIncl[ij].px(),jetsIncl[ij].py(),jetsIncl[ij].pz());
1619 fhChLeadZVsJetPt[type][r]->Fill(jetPt, cz);
1620 }
1621
1622 if((fStudySubEInHC||fStudyMcOverSubE) && type<2)
1623 {
1624 for(Int_t i=0; i<5; i++)
1625 {
1626 if(fStudySubEInHC)
1627 {
1628 fhSubClsEVsJetPt[type][r][i]->Fill(jetPt-subClsE[4],subClsE[i]/(jetPt-subClsE[4]));
1629 fhHCTrkPtClean[type][r][i]->Fill(jetPt-subClsE[4],subTrkPtClean[i]/(jetPt-subClsE[4]));
1630 fhHCTrkPtAmbig[type][r][i]->Fill(jetPt-subClsE[4],subTrkPtAmbig[i]/(jetPt-subClsE[4]));
1631 }
1632 if(type==0 && fIsMC && fStudyMcOverSubE)
1633 {
1634 fHCOverSubE[r][i]->Fill(jetPt-subClsE[4],subClsE[i]-mcSubE);
1635 fHCOverSubEFrac[r][i]->Fill(jetPt-subClsE[4],(subClsE[i]-mcSubE)/(jetPt-subClsE[4]));
1636 }
1637 }
1638 }
1639 if(type<2 && chPt>0)
1640 {
1641 for(Int_t i=0; i<3; i++)
1642 {
1643 fRelTrkCon[type][r]->Fill(jetPt,totTrkPt[i]/chPt,i);
1644 }
1645 }
1646
1647
1648 for(Int_t ibin=0; ibin<nBins; ibin++)
1649 {
1650 Double_t fill1[3] = {jetPt,radiusCut[ibin]-0.005,eFraction[ibin]/jetE};
1651 fJetEnergyFraction[type][r]->Fill(fill1);
1652 Double_t fill2[3] = {jetPt,radiusCut[ibin]-0.005,static_cast<Double_t>(nPart[ibin])};
1653 fJetNPartFraction[type][r]->Fill(fill2);
1654 }
1655
1656 // Get the jet pt containing tracks or clusters above some threshold
1657 if(type<2)
1658 {
1659 Double_t thres[3] = {15,25,40};
1660 Int_t okCh[3] = {0,0,0};
1661 Int_t okNe[3] = {0,0,0};
1662 Double_t lowPt[2] = {0.,0.};
1663 for(UInt_t ic=0; ic<constituents.size(); ic++)
1664 {
1665 Float_t partPt = constituents[ic].perp();
1666
1667 if(partPt<0.3) lowPt[0] += partPt;
1668 else if(partPt<0.5) lowPt[1] += partPt;
1669
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 }
1692 for(Int_t k=0; k<2; k++) fhJetPtVsLowPtCons[type][r][k]->Fill(jetPt,lowPt[k]/jetPt);
1693 }
1694 }
1695}
1696
1697//________________________________________________________________________
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)
1826{
1827 //
1828 // Run analysis to estimate the underlying event
1829 // Adapted from Oliver
1830
1831 std::vector<fastjet::PseudoJet> jetsIncl = jetFinder->GetInclusiveJets();
1832 Double_t leadpt=0, leadPhi = -999, leadEta = -999;
1833 Int_t leadIndex = -1;
1834 for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1835 {
1836 Double_t jetEta = jetsIncl[ij].eta();
1837 Double_t jetPt = jetsIncl[ij].perp();
1838 Double_t jetPhi = jetsIncl[ij].phi();
1839 if(leadpt<jetPt)
1840 {
1841 leadpt = jetPt;
1842 leadPhi = jetPhi;
1843 leadEta = jetEta;
1844 leadIndex = ij;
1845 }
1846 }
1847 if(leadpt<1 || TMath::Abs(leadEta)>0.5 ) return;
1848
1849 Int_t backIndex = -1;
1850 Bool_t isBackToBack = kFALSE;
1851 for(UInt_t ij=0; ij<jetsIncl.size(); ij++)
1852 {
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)
1856 {
1857 backIndex = ij;
1858 isBackToBack = kTRUE;
1859 }
1860 }
1861
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}
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;
1999 std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
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;
2034 std::vector<fastjet::PseudoJet> constituents = jetFinder->GetJetConstituents(jetIndex);
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//________________________________________________________________________
2105Int_t AliAnalysisTaskFullppJet::FindSpatialMatchedJet(fastjet::PseudoJet jet, AliFJWrapper *jetFinder, Double_t &dEta, Double_t &dPhi, Double_t maxR)
2106{
2107 //
2108 // Find spatially matched detector and particle jets
2109 //
2110
2111 Int_t index=-1;
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++)
2115 {
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) );
2119 if(tmpR<maxR)
2120 {
2121 maxR=tmpR;
2122 index=ij;
2123 dEta = jet.eta()-jets_incl[ij].eta();
2124 dPhi = jet.phi()-jets_incl[ij].phi();
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();
2140 std::vector<fastjet::PseudoJet> constituents1 = jetFinder1->GetJetConstituents(index1);
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;
2148 std::vector<fastjet::PseudoJet> constituents2 = jetFinder2->GetJetConstituents(ij2);
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();
2203 if (pt < fTrkPtMin[fTriggerType] || pt > fTrkPtMax[fTriggerType] || TMath::Abs(eta) > fTrkEtaMax)
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
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);
2333 if (newCluster->E() < fClsEtMin[fTriggerType] || newCluster->E() > fClsEtMax[fTriggerType])
2334 {delete newCluster; continue;}
2335
2336 // Absolute scale
2337 if(fPeriod.Contains("lhc11a",TString::kIgnoreCase))
2338 {
2339 Double_t newE = newCluster->E() * 1.02;
2340 newCluster->SetE(newE);
2341 }
2342
2343 // Trigger efficiency systematic uncertainty
2344 if(fSysJetTrigEff)
2345 {
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 }
2354 }
2355
2356 // Cluster energy scale systematic uncertainty
2357 if(fSysClusterEScale)
2358 {
2359 Double_t newE = newCluster->E() * (1+fVaryClusterEScale);
2360 newCluster->SetE(newE);
2361 }
2362
2363 // Cluster energy resolution systematic uncertainty
2364 if(fSysClusterERes)
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)
2375 {
2376 Double_t oldE = newCluster->E();
2377 Double_t newE = oldE/fNonLinear->Eval(oldE) * 1.012;
2378 newCluster->SetE(newE);
2379 }
2380
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 }
2387
2388 // hadronic correction
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);
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 }
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));
2469 if(IsElectron(track,cluster->E())) //electrons
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 }
2478 else //hadrons
2479 {
2480 if(fHadronicCorrection)
2481 {
2482 MIPE += (trkP>0.27)?0.27:trkP; //MIP correction
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();
2523 if(labels)
2524 {
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;
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
2555 // Adpated from AliEMCalRecoUtils
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 //
2638 Info("Terminate","Terminate");
2639 AliAnalysisTaskSE::Terminate();
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 {
2657 // Check every cluster
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;
2665 fIsEventTriggerBit = 1;
2666 cluster->SetChi2(1);
2667 }
2668 else
2669 cluster->SetChi2(0);
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;
2689 fRecoUtil->GetMaxEnergyCell(fGeom, cells, cluster, absID, iSupMod, ieta, iphi, share); // Get the position of the most energetic cell in the cluster
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;
2698 if(fCheckTriggerMask && fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5 && iSupMod>-1) // check the trigger mask
2699 {
2700 Double_t clsE = cluster->E();
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
2703 else
2704 {
2705 Int_t bin = fTriggerCurve[iSupMod]->FindFixBin(clsE);
2706 pro = fTriggerCurve[iSupMod]->GetBinContent(bin)/fTriggerEfficiency[iSupMod]->Eval(10); // Read the probability from trigger turn-on curves
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
2775 Double_t resolution = track->Pt()*track->Pt()*TMath::Sqrt(track->GetSigma1Pt2());
2776 Double_t smear = resolution*TMath::Sqrt((1+fVaryTrkPtRes)*(1+fVaryTrkPtRes)-1);
2777 return fRandomGen->Gaus(0, smear);
2778
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
2824 if(fCheckTriggerMask && fTriggerMask->GetBinContent(globCol+1,globRow+1)<0.5) continue;
2825 trigInCut = 0;
2826 for( Int_t i = 0; i < ntimes; i++ )
2827 {
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);
2848 if(!cluster || !cluster->IsEMCAL()) continue;
2849 cluster->SetChi2(0);
2850 if(!IsGoodCluster(cluster)) continue;
2851
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
2861 if(fCheckTriggerMask && fTriggerMask->GetBinContent(feta+1,fphi+1)>0.5)
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"};
2902 const char *type[2] = {"local","grid"};
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());
2910 printf("Reject SPD pileup: %s\n",decision[fRejectPileup]);
2911 printf("Reject exotic triggered events: %s\n",decision[fRejectExoticTrigger]);
2912 printf("Is this MC data: %s\n",decision[fIsMC]);
2913 printf("Only find charged jets in MC: %s\n",decision[fChargedMC]);
2914 printf("Analyze on local or grid? %s\n",type[fAnaType]);
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]);
2950 printf("Run UE analysis: %s\n",decision[fRunUE]);
2951 printf("Run secondaries: %s\n",decision[fRunSecondaries]);
2952 printf("=======================================\n\n");
2953}
2954
2955//________________________________________________________________________
2956void AliAnalysisTaskFullppJet::CheckExoticEvent()
2957{
2958 //
2959 // Check if the event containts exotic clusters
2960 //
2961
2962 Double_t leadingE = 0;
2963 Int_t nCls = fESD->GetNumberOfCaloClusters();
2964 for(Int_t icl=0; icl<nCls; icl++)
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
2978
2979//________________________________________________________________________
2980Bool_t AliAnalysisTaskFullppJet::HasPrimaryVertex() const
2981{
2982 //
2983 // Check if the primary vertex exists
2984 //
2985 const AliESDVertex* vtx = fESD->GetPrimaryVertex();
2986 if (!vtx || vtx->GetNContributors()<1) return kFALSE;
2987 else return kTRUE;
2988}
2989
2990
2991//________________________________________________________________________
2992Bool_t AliAnalysisTaskFullppJet::IsPrimaryVertexOk() const
2993{
2994 //
2995 // Check if the event vertex is good
2996 //
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
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 )
3016 return kTRUE;
3017 else
3018 return kFALSE;
3019}
3020
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
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}