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