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