AliPhysics  b0c77bb (b0c77bb)
AliAnalysisTaskEmcalHfeTagging.cxx
Go to the documentation of this file.
1 //
2 // HFE tagged jet shape analysis task.
3 //
4 // Authors: D. Caffarri, L. Cunqueiro (jet); D. Godoy (HFE)
5 //
6 #include <TClonesArray.h>
7 #include <TH1F.h>
8 #include <TH2F.h>
9 #include <TH3F.h>
10 #include <THnSparse.h>
11 #include <TTree.h>
12 #include <TList.h>
13 #include <TLorentzVector.h>
14 #include <TProfile.h>
15 #include <TChain.h>
16 #include <TSystem.h>
17 #include <TFile.h>
18 #include <TKey.h>
19 #include <AliAnalysisDataSlot.h>
20 #include <AliAnalysisDataContainer.h>
21 #include "TMatrixD.h"
22 #include "TMatrixDSym.h"
23 #include "TMatrixDSymEigen.h"
24 #include "TVector3.h"
25 #include "TVector2.h"
26 #include "AliVCluster.h"
27 #include "AliVTrack.h"
28 #include "AliEmcalJet.h"
29 #include "AliRhoParameter.h"
30 #include "AliLog.h"
31 #include "AliEmcalParticle.h"
32 #include "AliMCEvent.h"
33 #include "AliGenPythiaEventHeader.h"
34 #include "AliAODMCHeader.h"
35 #include "AliAODHandler.h"
36 #include "AliAODEvent.h"
37 #include "AliMCEvent.h"
38 #include "AliAnalysisManager.h"
39 #include "AliJetContainer.h"
40 #include "AliParticleContainer.h"
41 #include "AliEmcalPythiaInfo.h"
42 #include "TRandom3.h"
43 #include "AliAODInputHandler.h"
44 #include "AliPIDResponse.h"
45 #include "AliCFContainer.h"
46 #include "AliCFManager.h"
47 #include "AliMultiEventInputHandler.h"
48 #include "AliKFParticle.h"
49 #include "AliKFVertex.h"
50 #include "AliMultSelection.h"
51 #include "AliAnalysisUtils.h"
52 
53 #include "AliAODEvent.h"
55 
56 using std::cout;
57 using std::endl;
58 
60 
61 //________________________________________________________________________
64 fAOD(0),
65 fVevent(0),
66 fpidResponse(0),
67 fTracksTender(0),
68 pVtx(0),
69 spdVtx(0),
70 fMC(0),
71 fStack(0),
72 fMCparticle(0),
73 fMCarray(0),
74 fMCheader(0),
75 fContainer(0),
76 fMinFractionShared(0),
77 fJetShapeType(kData),
78 fJetShapeSub(kNoSub),
79 fJetSelection(kInclusive),
80 fPtThreshold(-9999.),
81 fRMatching(0.2),
82 fSelectedShapes(0),
83 fminpTTrig(20.),
84 fmaxpTTrig(50.),
85 fangWindowRecoil(0.6),
86 fSemigoodCorrect(0),
87 fHolePos(0),
88 fHoleWidth(0),
89 fCentSelectOn(kTRUE),
90 fCentMin(0),
91 fCentMax(10),
92 fOneConstSelectOn(kFALSE),
93 fDerivSubtrOrder(0),
94 fMCweight(0),
95 fRunNumber(0),
96 fAssPtCut(0.1),
97 fITSncut(3),
98 fAssTPCnCut(60),
99 fTPCnCut(100),
100 fAssITSrefitCut(kTRUE),
101 fUseTender(kFALSE),
102 fSigmaTOFcut(3.),
103 fSigmaTPCcut(-1.),
104 fDcaXYcut(1.),
105 fDcaZcut(2.),
106 fIMcut(0.1),
107 fNeventV0(0),
108 fNeventT0(0),
109 fh2ResponseUW(0x0),
110 fh2ResponseW(0x0),
111 fPhiJetCorr6(0x0),
112 fPhiJetCorr7(0x0),
113 fEtaJetCorr6(0x0),
114 fEtaJetCorr7(0x0),
115 fPtJetCorr(0x0),
116 fPtJet(0x0),
117 fPtGenJet(0),
118 fPhiJet(0x0),
119 fEtaJet(0x0),
120 fEtaPhiJet(0x0),
121 fAreaJet(0x0),
122 fJetProbDensityDetPart(0x0),
123 fJetProbDensityPartDet(0x0),
124 fNbOfConstvspT(0x0),
125 fnTPCnTOFnocut(0),
126 fnTPCnocutP(0),
127 fnTOFnocutP(0),
128 fnTPCcutP(0),
129 fnTPCcutPt(0),
130 fnULSmLSpairsPerElectron(0),
131 fnPartPerJet(0),
132 fnElecOverPartPerJet(0),
133 fnInclElecPerJet(0),
134 fnPhotElecPerJet(0),
135 fnIncSubPhotElecPerJet(0),
136 fnTrueElecPerJet(0),
137 fnTrueHFElecPerJet(0),
138 fnTruePElecPerJet(0),
139 fPi0PtGen(0),
140 fPi0PtEnh(0),
141 fEtaPtGen(0),
142 fEtaPtEnh(0),
143 fGenHfePt(0),
144 fGenPePt(0),
145 fPtP(0),
146 fptJetIE(0),
147 fptJetPE(0),
148 fptJetHFE(0),
149 fptRecPE(0),
150 fptTruePE(0),
151 fptWrongPE(0),
152 fPtTrack(0),
153 fPhiTrack(0x0),
154 fEtaTrack(0x0),
155 fEtaPhiTrack(0x0),
156 fPhiRecElecTPC(0x0),
157 fEtaRecElecTPC(0x0),
158 fEtaPhiRecElecTPC(0x0),
159 fPhiRecElecEMCal(0x0),
160 fEtaRecElecEMCal(0x0),
161 fEtaPhiRecElecEMCal(0x0),
162 fPhiTrueElec(0x0),
163 fEtaTrueElec(0x0),
164 fEtaPhiTrueElec(0x0),
165 fnEovPelecTPCcut(0x0),
166 fnEovPelecTPCEMCalcut(0x0),
167 fnEovPbackg(0x0),
168 fnClsE(0x0),
169 fnM20(0x0),
170 fnM02(0x0),
171 fnClsTime(0x0),
172 fTreeObservableTagging(0)
173 
174 {
175  for(Int_t i=0;i<26;i++){
176  fShapesVar[i]=0;
177  }
178 
179  for(Int_t i=0;i<5;i++){
180  fptTrueHFEeffTPCTOF[i] = NULL;
181  fptTrueHFEeffEMCal[i] = NULL;
182  }
183 
184  for(Int_t i=0;i<5;i++){
185  fInvmassLS[i] = NULL;
186  fInvmassULS[i] = NULL;
187  fUlsLsElecPt[i] = NULL;
188  fTotElecPt[i] = NULL;
189  for(Int_t j=0;j<18;j++){
190  fnTPCTrueParticles[i][j] = NULL;
191  fnTPCSigma[i][j] = NULL;
192  }
193  }
194  SetMakeGeneralHistograms(kTRUE);
195  DefineOutput(1, TList::Class());
196  DefineOutput(2, TTree::Class());
197 
198 }
199 
200 //________________________________________________________________________
202 AliAnalysisTaskEmcalJet(name, kTRUE),
203 fAOD(0),
204 fVevent(0),
205 fpidResponse(0),
206 fTracksTender(0),
207 pVtx(0),
208 spdVtx(0),
209 fMC(0),
210 fStack(0),
211 fMCparticle(0),
212 fMCarray(0),
213 fMCheader(0),
214 fContainer(0),
215 fMinFractionShared(0),
216 fJetShapeType(kData),
217 fJetShapeSub(kNoSub),
218 fJetSelection(kInclusive),
219 fPtThreshold(-9999.),
220 fRMatching(0.2),
221 fSelectedShapes(0),
222 fminpTTrig(20.),
223 fmaxpTTrig(50.),
224 fangWindowRecoil(0.6),
225 fSemigoodCorrect(0),
226 fHolePos(0),
227 fHoleWidth(0),
228 fCentSelectOn(kTRUE),
229 fCentMin(0),
230 fCentMax(10),
231 fOneConstSelectOn(kFALSE),
232 fDerivSubtrOrder(0),
233 fMCweight(0),
234 fRunNumber(0),
235 fAssPtCut(0.1),
236 fITSncut(3),
237 fAssTPCnCut(60),
238 fTPCnCut(100),
239 fAssITSrefitCut(kTRUE),
240 fUseTender(kFALSE),
241 fSigmaTOFcut(3.),
242 fSigmaTPCcut(-1.),
243 fDcaXYcut(1.),
244 fDcaZcut(2.),
245 fIMcut(0.1),
246 fNeventV0(0),
247 fNeventT0(0),
248 fh2ResponseUW(0x0),
249 fh2ResponseW(0x0),
250 fPhiJetCorr6(0x0),
251 fPhiJetCorr7(0x0),
252 fEtaJetCorr6(0x0),
253 fEtaJetCorr7(0x0),
254 fPtJetCorr(0x0),
255 fPtJet(0x0),
256 fPtGenJet(0),
257 fPhiJet(0x0),
258 fEtaJet(0x0),
259 fEtaPhiJet(0x0),
260 fAreaJet(0x0),
261 fJetProbDensityDetPart(0x0),
262 fJetProbDensityPartDet(0x0),
263 fNbOfConstvspT(0x0),
264 fnTPCnTOFnocut(0),
265 fnTPCnocutP(0),
266 fnTOFnocutP(0),
267 fnTPCcutP(0),
268 fnTPCcutPt(0),
269 fnULSmLSpairsPerElectron(0),
270 fnPartPerJet(0),
271 fnElecOverPartPerJet(0),
272 fnInclElecPerJet(0),
273 fnPhotElecPerJet(0),
274 fnIncSubPhotElecPerJet(0),
275 fnTrueElecPerJet(0),
276 fnTrueHFElecPerJet(0),
277 fnTruePElecPerJet(0),
278 fPi0PtGen(0),
279 fPi0PtEnh(0),
280 fEtaPtGen(0),
281 fEtaPtEnh(0),
282 fGenHfePt(0),
283 fGenPePt(0),
284 fPtP(0),
285 fptJetIE(0),
286 fptJetPE(0),
287 fptJetHFE(0),
288 fptRecPE(0),
289 fptTruePE(0),
290 fptWrongPE(0),
291 fPtTrack(0),
292 fPhiTrack(0x0),
293 fEtaTrack(0x0),
294 fEtaPhiTrack(0x0),
295 fPhiRecElecTPC(0x0),
296 fEtaRecElecTPC(0x0),
297 fEtaPhiRecElecTPC(0x0),
298 fPhiRecElecEMCal(0x0),
299 fEtaRecElecEMCal(0x0),
300 fEtaPhiRecElecEMCal(0x0),
301 fPhiTrueElec(0x0),
302 fEtaTrueElec(0x0),
303 fEtaPhiTrueElec(0x0),
304 fnEovPelecTPCcut(0x0),
305 fnEovPelecTPCEMCalcut(0x0),
306 fnEovPbackg(0x0),
307 fnClsE(0x0),
308 fnM20(0x0),
309 fnM02(0x0),
310 fnClsTime(0x0),
311 fTreeObservableTagging(0)
312 
313 {
314  // Standard constructor.
315  for(Int_t i=0;i<26;i++){
316  fShapesVar[i]=0;
317  }
318 
319  for(Int_t i=0;i<5;i++){
320  fptTrueHFEeffTPCTOF[i] = NULL;
321  fptTrueHFEeffEMCal[i] = NULL;
322  }
323 
324  for(Int_t i=0;i<5;i++){
325  fInvmassLS[i] = NULL;
326  fInvmassULS[i] = NULL;
327  fUlsLsElecPt[i] = NULL;
328  fTotElecPt[i] = NULL;
329  for(Int_t j=0;j<18;j++){
330  fnTPCTrueParticles[i][j] = NULL;
331  fnTPCSigma[i][j] = NULL;
332  }
333  }
335 
336  DefineOutput(1, TList::Class());
337  DefineOutput(2, TTree::Class());
338 
339 }
340 
341 //________________________________________________________________________
343 {
344  // Destructor.
345 }
346 
347 //________________________________________________________________________
349 {
350  // Create user output.
351 
353 
354  Bool_t oldStatus = TH1::AddDirectoryStatus();
355  TH1::AddDirectory(kFALSE);
356 
357  fNeventV0 = new TH1F("fNeventV0","Number of Events (V0)",5,-0.5,4.5);
358  fOutput->Add(fNeventV0);
359 
360  fNeventT0 = new TH1F("fNeventT0","Number of Events (T0)",5,-0.5,4.5);
361  fOutput->Add(fNeventT0);
362 
363 
364  fh2ResponseUW= new TH2F("fh2ResponseUW", "fh2ResponseUW", 100, 0, 200, 100, 0, 200);
365  fOutput->Add(fh2ResponseUW);
366  fh2ResponseW= new TH2F("fh2ResponseW", "fh2ResponseW", 100, 0, 200, 100, 0, 200);
367  fOutput->Add(fh2ResponseW);
368  fPhiJetCorr6= new TH2F("fPhiJetCorr6", "fPhiJetCorr6", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
369  fOutput->Add(fPhiJetCorr6);
370  fEtaJetCorr6= new TH2F("fEtaJetCorr6", "fEtaJetCorr6", 50, -1.5, 1.5, 50, -1.5, 1.5);
371  fOutput->Add(fEtaJetCorr6);
372 
373  fPhiJetCorr7= new TH2F("fPhiJetCorr7", "fPhiJetCorr7", 50, 0, 2*TMath::Pi(), 50, 0, 2*TMath::Pi());
374  fOutput->Add(fPhiJetCorr7);
375  fEtaJetCorr7= new TH2F("fEtaJetCorr7", "fEtaJetCorr7", 50, -1.5, 1.5, 50, -1.5, 1.5);
376  fOutput->Add(fEtaJetCorr7);
377 
378  fPtJetCorr= new TH2F("fPtJetCorr", "fPtJetCorr", 100, 0, 200, 100, 0, 200);
379  fOutput->Add(fPtJetCorr);
380 
381  fPtJet= new TH1F("fPtJet", "fPtJet", 100, 0, 200);
382  fOutput->Add(fPtJet);
383 
384  fPtGenJet= new TH1F("fPtGenJet", "fPtGenJet", 100, 0, 200);
385  fOutput->Add(fPtGenJet);
386 
387  fPhiJet= new TH2F("fPhiJet", "fPhiJet", 100, 0, 200, 100, 0, TMath::TwoPi());
388  fOutput->Add(fPhiJet);
389 
390  fEtaJet= new TH2F("fEtaJet", "fEtaJet", 100, 0, 200, 100, -1,1);
391  fOutput->Add(fEtaJet);
392 
393  fEtaPhiJet= new TH2F("fEtaPhiJet", "fEtaPhiJet", 100, 0, TMath::TwoPi(), 100, -1,1);
394  fOutput->Add(fEtaPhiJet);
395 
396  fAreaJet= new TH2F("fAreaJet", "fAreaJet", 100, 0, 200, 100, 0,1.5);
397  fOutput->Add(fAreaJet);
398 
399  fNbOfConstvspT=new TH2F("fNbOfConstvspT", "fNbOfConstvspT", 100, 0, 100, 200, 0, 200);
400  fOutput->Add(fNbOfConstvspT);
401 
402  fnTPCnTOFnocut=new TH2F("fnTPCnTOFnocut", "fnTPCnTOFnocut", 200, -10, 10, 200, -10, 10);
403  fOutput->Add(fnTPCnTOFnocut);
404 
405  fnTPCnocutP=new TH2F("fnTPCnocutP", "fnTPCnocutP", 50, 0, 5, 200, -10, 10);
406  fOutput->Add(fnTPCnocutP);
407 
408  fnTOFnocutP=new TH2F("fnTOFnocutP", "fnTOFnocutP", 50, 0, 5, 200, -10, 10);
409  fOutput->Add(fnTOFnocutP);
410 
411  fnTPCcutP=new TH2F("fnTPCcutP", "fnTPCcutP", 50, 0, 5, 200, -10, 10);
412  fOutput->Add(fnTPCcutP);
413 
414  fnTPCcutPt=new TH2F("fnTPCcutPt", "fnTPCcutPt", 50, 0, 5, 200, -10, 10);
415  fOutput->Add(fnTPCcutPt);
416 
417  fnULSmLSpairsPerElectron =new TH2F("fnULSmLSpairsPerElectron", "fnULSmLSpairsPerElectron", 50, 0, 5, 100, 0, 100);
419 
420  fnPartPerJet=new TH1F("fnPartPerJet", "fnPartPerJet", 500,0,500);
421  fOutput->Add(fnPartPerJet);
422 
423  fnElecOverPartPerJet=new TH1F("fnElecOverPartPerJet", "fnElecOverPartPerJet", 100,0,0.1);
425 
426  fnInclElecPerJet=new TH1F("fnInclElecPerJet", "fnInclElecPerJet", 100,0,100);
428 
429  fnPhotElecPerJet=new TH1F("fnPhotElecPerJet", "fnPhotElecPerJet", 100,0,100);
431 
432  fnIncSubPhotElecPerJet=new TH1F("fnIncSubPhotElecPerJet", "fnIncSubPhotElecPerJet", 101,-1,100);
434 
435  fnTrueElecPerJet=new TH1F("fnTrueElecPerJet", "fnTrueElecPerJet", 100,0,100);
437 
438  fnTrueHFElecPerJet=new TH1F("fnTrueHFElecPerJet", "fnTrueHFElecPerJet", 100,0,100);
440 
441  fnTruePElecPerJet=new TH1F("fnTruePElecPerJet", "fnTruePElecPerJet", 100,0,100);
443 
444  fJetProbDensityDetPart=new TH2F("fJetProbDensityDetPart", "fJetProbDensityDetPart",200,-2,2,100, 0,200);
446 
447  fJetProbDensityPartDet=new TH2F("fJetProbDensityPartDet", "fJetProbDensityPartDet",200,-2,2,100, 0,200);
449 
450  int nbin = 59;
451  double xbins[60] = {0.01,0.1,0.12,0.14,0.16,0.18,0.2,0.25,0.3,0.35,
452  0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,
453  0.9,0.95,1,1.1,1.2,1.3,1.4,1.5,1.6,1.7,
454  1.8,1.9,2,2.2,2.4,2.6,2.8,3,3.2,3.4,
455  3.6,3.8,4,4.5,5,5.5,6,6.5,7,8,
456  9,10,11,12,13,14,15,16,18,20};
457 
458  fPi0PtGen = new TH1F("fPi0PtGen","fPi0PtGen",nbin,xbins);
459  fOutput->Add(fPi0PtGen);
460 
461  fPi0PtEnh = new TH1F("fPi0PtEnh","fPi0PtEnh",nbin,xbins);
462  fOutput->Add(fPi0PtEnh);
463 
464  fEtaPtGen = new TH1F("fEtaPtGen", "fEtaPtGen",nbin,xbins);
465  fOutput->Add(fEtaPtGen);
466 
467  fEtaPtEnh = new TH1F("fEtaPtEnh", "fEtaPtEnh",nbin,xbins);
468  fOutput->Add(fEtaPtEnh);
469 
470  fGenHfePt = new TH1F("fGenHfePt","fGenHfePt",nbin,xbins);
471  fOutput->Add(fGenHfePt);
472 
473  fGenPePt = new TH1F("fGenPePt", "fGenPePt",nbin,xbins);
474  fOutput->Add(fGenPePt);
475 
476  Double_t ptRange[34] = {0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4,
477  1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 5.,
478  6., 8., 10., 12., 14., 16., 19., 22., 26., 30.,
479  35., 40., 45., 50.};
480 
481  fPtP=new TH2F("fPtP", "fPtP", 33,ptRange,33,ptRange);
482  fOutput->Add(fPtP);
483 
484  for(Int_t i=0;i<5;i++){
485 
486  fInvmassLS[i] = new TH2F(Form("fInvmassLS%d",i), Form("fInvmassLS%d",i), 33,ptRange, 100, 0, 0.5);
487  fOutput->Add(fInvmassLS[i]);
488 
489  fInvmassULS[i] = new TH2F(Form("fInvmassULS%d",i), Form("fInvmassULS%d",i), 33,ptRange, 100, 0, 0.5);
490  fOutput->Add(fInvmassULS[i]);
491 
492  fUlsLsElecPt[i] = new TH1F(Form("fUlsLsElecPt%d",i), Form("fUlsLsElecPt%d",i), 33,ptRange);
493  fOutput->Add(fUlsLsElecPt[i]);
494 
495  fTotElecPt[i] = new TH1F(Form("fTotElecPt%d",i), Form("fTotElecPt%d",i), 33,ptRange);
496  fOutput->Add(fTotElecPt[i]);
497 
498 
499  for(Int_t j=0;j<18;j++){
500  fnTPCTrueParticles[i][j] = new TH2F(Form("fnTPCTrueParticles%d%d",i,j),Form("fnTPCTrueParticles%d%d",i,j), 7,0,7,100,-15,15);
501  fOutput->Add(fnTPCTrueParticles[i][j]);
502 
503  fnTPCSigma[i][j] = new TH1F(Form("fnTPCSigma%d%d",i,j),Form("fnTPCSigma%d%d",i,j), 100,-15,15);
504  fOutput->Add(fnTPCSigma[i][j]);
505  }
506  }
507 
508  for(Int_t i=0;i<5;i++){
509  fptTrueHFEeffTPCTOF[i] = new TH1F(Form("fptTrueHFEeffTPCTOF%d",i), Form("fptTrueHFEeffTPCTOF%d",i), 33,ptRange);
510  fOutput->Add(fptTrueHFEeffTPCTOF[i]);
511 
512  fptTrueHFEeffEMCal[i] = new TH1F(Form("fptTrueHFEeffEMCal%d",i), Form("fptTrueHFEeffEMCal%d",i), 33,ptRange);
513  fOutput->Add(fptTrueHFEeffEMCal[i]);
514  }
515 
516  fptJetIE= new TH1F("fptJetIE", "fptJetIE", 100, 0, 200);
517  fOutput->Add(fptJetIE);
518 
519  fptJetPE= new TH1F("fptJetPE", "fptJetPE", 100, 0, 200);
520  fOutput->Add(fptJetPE);
521 
522  fptJetHFE= new TH1F("fptJetHFE", "fptJetHFE", 100, 0, 200);
523  fOutput->Add(fptJetHFE);
524 
525  fptRecPE= new TH1F("fptRecPE", "fptRecPE", 33,ptRange);
526  fOutput->Add(fptRecPE);
527 
528  fptTruePE= new TH1F("fptTruePE", "fptTruePE", 33,ptRange);
529  fOutput->Add(fptTruePE);
530 
531  fptWrongPE= new TH1F("fptWrongPE", "fptWrongPE", 33,ptRange);
532  fOutput->Add(fptWrongPE);
533 
534  fPtTrack= new TH1F("fPtTrack", "fPtTrack", 100, 0, 200);
535  fOutput->Add(fPtTrack);
536 
537  fPhiTrack= new TH2F("fPhiTrack", "fPhiTrack", 100, 0, 200, 100, 0, TMath::TwoPi());
538  fOutput->Add(fPhiTrack);
539 
540  fEtaTrack= new TH2F("fEtaTrack", "fEtaTrack", 100, 0, 200, 100, -1,1);
541  fOutput->Add(fEtaTrack);
542 
543  fEtaPhiTrack= new TH2F("fEtaPhiTrack", "fEtaPhiTrack", 100, 0, TMath::TwoPi(), 100, -1,1);
544  fOutput->Add(fEtaPhiTrack);
545 
546  fPhiRecElecTPC= new TH2F("fPhiRecElecTPC", "fPhiRecElecTPC", 100, 0, 50, 100, 0, TMath::TwoPi());
547  fOutput->Add(fPhiRecElecTPC);
548 
549  fEtaRecElecTPC= new TH2F("fEtaRecElecTPC", "fEtaRecElecTPC", 100, 0, 50, 100, -1,1);
550  fOutput->Add(fEtaRecElecTPC);
551 
552  fEtaPhiRecElecTPC= new TH2F("fEtaPhiRecElecTPC", "fEtaPhiRecElecTPC", 100, 0, TMath::TwoPi(), 100, -1,1);
554 
555  fPhiRecElecEMCal= new TH2F("fPhiRecElecEMCal", "fPhiRecElecEMCal", 100, 0, 50, 100, 0, TMath::TwoPi());
557 
558  fEtaRecElecEMCal= new TH2F("fEtaRecElecEMCal", "fEtaRecElecEMCal", 100, 0, 50, 100, -1,1);
560 
561  fEtaPhiRecElecEMCal= new TH2F("fEtaPhiRecElecEMCal", "fEtaPhiRecElecEMCal", 100, 0, TMath::TwoPi(), 100, -1,1);
563 
564  fPhiTrueElec= new TH2F("fPhiTrueElec", "fPhiTrueElec", 100, 0, 50, 100, 0, TMath::TwoPi());
565  fOutput->Add(fPhiTrueElec);
566 
567  fEtaTrueElec= new TH2F("fEtaTrueElec", "fEtaTrueElec", 100, 0, 50, 100, -1,1);
568  fOutput->Add(fEtaTrueElec);
569 
570  fEtaPhiTrueElec= new TH2F("fEtaPhiTrueElec", "fEtaPhiTrueElec", 100, 0, TMath::TwoPi(), 100, -1,1);
571  fOutput->Add(fEtaPhiTrueElec);
572 
573  fnEovPelecTPCcut = new TH2F("fnEovPelecTPCcut", "fnEovPelecTPCcut", 100, 0, 100, 40,0,2);
575 
576  fnEovPelecTPCEMCalcut = new TH2F("fnEovPelecTPCEMCalcut", "fnEovPelecTPCEMCalcut", 100, 0, 100,40,0,2);
578 
579  fnEovPbackg = new TH2F("fnEovPbackg", "fnEovPbackg", 100, 0, 100,40,0,2);
580  fOutput->Add(fnEovPbackg);
581 
582  fnClsE = new TH2F("fnClsE", "fnClsE", 100, 0, 100, 100, 0,100);
583  fOutput->Add(fnClsE);
584 
585  fnM20 = new TH2F("fnM20", "fnM20", 100, 0, 100, 100, 0,2);
586  fOutput->Add(fnM20);
587 
588  fnM02 = new TH2F("fnM02", "fnM02", 100, 0, 100, 100, 0,2);
589  fOutput->Add(fnM02);
590 
591  fnClsTime = new TH2F("fnClsTime", "fnClsTime", 100, 0, 100, 100, -200,200);
592  fOutput->Add(fnClsTime);
593 
594 
595  // =========== Switch on Sumw2 for all histos ===========
596  for (Int_t i=0; i<fOutput->GetEntries(); ++i) {
597  TH1 *h1 = dynamic_cast<TH1*>(fOutput->At(i));
598  if (h1){
599  h1->Sumw2();
600  continue;
601  }
602  THnSparse *hn = dynamic_cast<THnSparse*>(fOutput->At(i));
603  if(hn)hn->Sumw2();
604  }
605 
606 
607  TH1::AddDirectory(oldStatus);
608  const Int_t nVar = 26;
609 
610  fTreeObservableTagging = new TTree(GetOutputSlot(2)->GetContainer()->GetName(), GetOutputSlot(2)->GetContainer()->GetName());
611 
612  TString *fShapesVarNames = new TString [nVar];
613 
614  fShapesVarNames[0] = "partonCode";
615  fShapesVarNames[1] = "ptJet";
616  fShapesVarNames[2] = "ptDJet";
617  fShapesVarNames[3] = "mJet";
618  // fShapesVarNames[4] = "nbOfConst";
619  fShapesVarNames[4] = "angularity";
620  fShapesVarNames[5] = "circularity";
621  fShapesVarNames[6] = "lesub";
622  fShapesVarNames[7] = "coronna";
623 
624  fShapesVarNames[8] = "ptJetMatch";
625  fShapesVarNames[9] = "ptDJetMatch";
626  fShapesVarNames[10] = "mJetMatch";
627  // fShapesVarNames[12] = "nbOfConstMatch";
628  fShapesVarNames[11] = "angularityMatch";
629  fShapesVarNames[12] = "circularityMatch";
630  fShapesVarNames[13] = "lesubMatch";
631  fShapesVarNames[14] = "coronnaMatch";
632  fShapesVarNames[15] = "weightPythia";
633  //fShapesVarNames[14]="ntrksEvt";
634  fShapesVarNames[16] = "rhoVal";
635  fShapesVarNames[17] = "rhoMassVal";
636  fShapesVarNames[18] = "ptUnsub";
637  fShapesVarNames[19] = "pElec";
638  fShapesVarNames[20] = "ptElec";
639  fShapesVarNames[21] = "nInclElec";
640  fShapesVarNames[22] = "nPhotElec";
641  fShapesVarNames[23] = "hasElec";
642  fShapesVarNames[24] = "nTrueElectronsMC";
643  fShapesVarNames[25] = "nTrueHFElecMC";
644 
645 
646  for(Int_t ivar=0; ivar < nVar; ivar++){
647  //cout<<"looping over variables"<<endl;
648  fTreeObservableTagging->Branch(fShapesVarNames[ivar].Data(), &fShapesVar[ivar], Form("%s/F", fShapesVarNames[ivar].Data()));
649  }
650 
651  PostData(1,fOutput);
652  PostData(2,fTreeObservableTagging);
653 
654  delete [] fShapesVarNames;
655 }
656 
657 //________________________________________________________________________
659 {
660  // // Run analysis code here, if needed. It will be executed before FillHistograms().
661 
662  fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
663  if (!fAOD){
664  printf("ERROR: fAOD not available\n");
665  return kFALSE;
666  }
667 
668  fVevent = dynamic_cast<AliVEvent*>(InputEvent());
669  if(!fVevent){
670  printf("ERROR: fVEvent not available\n");
671  return kFALSE;
672  }
673 
674  fRunNumber = fVevent->GetRunNumber();
675 
676  //remaining event selection
677  pVtx = fVevent->GetPrimaryVertex();
678  spdVtx = fVevent->GetPrimaryVertexSPD();
679 
680  fMCarray = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
681  fMCheader = dynamic_cast<AliAODMCHeader*>(fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
682 
683  fNeventV0->Fill(0);
684  if (!fMCheader){
685  TString firedTriggerClasses = static_cast<const AliAODEvent*>(InputEvent())->GetFiredTriggerClasses();
686  if ((firedTriggerClasses.Contains("C0TVX-B-NOPF-CENT"))){
687  fNeventT0->Fill(0);
688  }
689  }
690 
691  Int_t NpureMC = -1, NpureMCproc = -1;
692 
693  if (fMCheader){
694  TList *lh=fMCheader->GetCocktailHeaders();
695  if(lh){
696  for(int igene=0; igene<lh->GetEntries(); igene++){
697  AliGenEventHeader* gh=(AliGenEventHeader*)lh->At(igene);
698  if(gh){
699  // cout << "<------- imc = "<< gh->GetName() << endl;
700  if(igene==0)NpureMC = gh->NProduced(); // generate by PYTHIA or HIJING
701  NpureMCproc += gh->NProduced();
702  }
703  }
704  }
705  }
706  // pileup rejection with SPD multiple vertices and track multivertexer (using default configuration)
707 
708  Bool_t isPileupfromSPDmulbins = kFALSE, isPileupFromMV = kFALSE;
709  isPileupfromSPDmulbins = fAOD->IsPileupFromSPDInMultBins();
710  if (isPileupfromSPDmulbins) return kFALSE;
711  fNeventV0->Fill(1);
712 
713  AliAnalysisUtils utils;
714  utils.SetMinPlpContribMV(5);
715  utils.SetMaxPlpChi2MV(5.);
716  utils.SetMinWDistMV(15);
717  utils.SetCheckPlpFromDifferentBCMV(kFALSE);
718  isPileupFromMV = utils.IsPileUpMV(fAOD);
719  if (isPileupFromMV) return kFALSE;
720  fNeventV0->Fill(2);
721 
722  // Selection of pi0 and eta in MC to compute the weight
723 
724  Bool_t isPrimary = kFALSE, isFromLMdecay = kTRUE, isFromHFdecay=kTRUE, hasMother = kTRUE;
725 
726  Double_t MCweight = 1.;
727  Int_t iDecay = 0;
728 
729  if(fMCarray){
730  //Int_t nParticles = fMCarray->GetEntries();
731  for (Int_t iParticle = 0; iParticle < NpureMCproc; iParticle++) {
732  AliAODMCParticle* particle = (AliAODMCParticle*) fMCarray->At(iParticle);
733  int fPDG = particle->GetPdgCode();
734  Double_t pTMC = particle->Pt();
735  Double_t phiMC = particle->Phi();
736  Double_t etaMC = particle->Eta();
737 
738  Bool_t iEnhance = kFALSE;
739  if(iParticle>=NpureMC)iEnhance = kTRUE;
740 
741  //if (TMath::Abs(etaMC)>1.2)continue;
742 
743  isPrimary = IsPrimary(particle);
744  isFromLMdecay = IsFromLMdecay(particle);
745  isFromHFdecay = IsFromHFdecay(particle);
746  hasMother = HasMother(particle);
747 
748  MCweight = 1.;
749  iDecay = 0;
750 
751  GetWeightAndDecay(particle,iDecay,MCweight);
752 
753 
754  if (TMath::Abs(etaMC)<0.8 && TMath::Abs(fPDG)==11){
755  fPhiTrueElec->Fill(pTMC,phiMC);
756  fEtaTrueElec->Fill(pTMC,etaMC);
757  if (pTMC > 0.5) fEtaPhiTrueElec->Fill(phiMC,etaMC);
758 
759  if (isFromHFdecay) fGenHfePt->Fill(pTMC);
760  if (iDecay>0 && iDecay<7) fGenPePt->Fill(pTMC);
761  }
762 
763 
764  Double_t yMC = particle->Y();
765  if (TMath::Abs(yMC)>1.0)continue;
766 
767  if (isPrimary){
768  if (!hasMother || (!isFromLMdecay && !isFromHFdecay)){
769  if(fPDG==111 && iEnhance==kFALSE) fPi0PtGen->Fill(pTMC);
770  if(fPDG==111 && iEnhance==kTRUE) fPi0PtEnh->Fill(pTMC);
771 
772  if(fPDG==221 && iEnhance==kFALSE) fEtaPtGen->Fill(pTMC);
773  if(fPDG==221 && iEnhance==kTRUE) fEtaPtEnh->Fill(pTMC);
774  }
775  }
776  }
777  }//MC
778  return kTRUE;
779 }
780 //________________________________________________________________________
782 {
783  // Fill histograms.
784  //cout<<"base container"<<endl;
785  AliEmcalJet* jet1 = NULL;
786  AliJetContainer *jetCont = GetJetContainer(0);
787 
788  Float_t kWeight=1;
789  if (fCentSelectOn)
790  if ((fCent>fCentMax) || (fCent<fCentMin)) return 0;
791 
792  // list of kink mothers
793  Int_t nVertices = 1;
794  nVertices = fAOD->GetNumberOfVertices();
795  Double_t listofmotherkink[nVertices];
796  Int_t nMotherKink = 0;
797  for(Int_t ivertex=0; ivertex < nVertices; ivertex++) {
798  AliAODVertex *vertex = fAOD->GetVertex(ivertex);
799  if(!vertex) continue;
800  if(vertex->GetType()==AliAODVertex::kKink) {
801  AliAODTrack *mother = (AliAODTrack *) vertex->GetParent();
802  if(!mother) continue;
803  Int_t idmother = mother->GetID();
804  listofmotherkink[nMotherKink] = idmother;
805  nMotherKink++;
806  }
807  }
808 
809  Float_t rhoVal=0, rhoMassVal = 0.;
810 
811  if(jetCont) {
812  jetCont->ResetCurrentID();
814  //rho
815  AliRhoParameter* rhoParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject("RhoSparseR020"));
816  if (!rhoParam) {
817  Printf("%s: Could not retrieve rho %s (some histograms will be filled with zero)!", GetName(), jetCont->GetRhoName().Data());
818  } else rhoVal = rhoParam->GetVal();
819  //rhom
820  AliRhoParameter* rhomParam = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject("RhoMassSparseR020"));
821  if (!rhomParam) {
822  Printf("%s: Could not retrieve rho_m %s (some histograms will be filled with zero)!", GetName(), jetCont->GetRhoMassName().Data());
823  } else rhoMassVal = rhomParam->GetVal();
824  }
825 
826  while((jet1 = jetCont->GetNextAcceptJet())) {
827  if (!jet1) continue;
828  AliEmcalJet* jet2 = 0x0;
829  AliEmcalJet* jet3 = 0x0;
830  fPtJet->Fill(jet1->Pt());
831  fPhiJet->Fill(jet1->Pt(),jet1->Phi());
832  fEtaJet->Fill(jet1->Pt(),jet1->Eta());
833  if (jet1->Pt() > 5.) fEtaPhiJet->Fill(jet1->Phi(),jet1->Eta());
834  fAreaJet->Fill(jet1->Pt(),jet1->Area());
835  AliEmcalJet *jetUS = NULL;
836  Int_t ifound=0, ilab=-1;
837 
838  fShapesVar[0] = 0.;
840  //AliJetContainer *jetContTrue = GetJetContainer(1);
841  AliJetContainer *jetContUS = GetJetContainer(2);
842 
844  for(Int_t i = 0; i<jetContUS->GetNJets(); i++) {
845  jetUS = jetContUS->GetJet(i);
846  if(jetUS->GetLabel()==jet1->GetLabel()) {
847  ifound++;
848  if(ifound==1) ilab = i;
849  }
850  }
851  if(ilab==-1) continue;
852  jetUS=jetContUS->GetJet(ilab);
853  jet2=jetUS->ClosestJet();
854  }
855 
857  if (!jet2) {
858  Printf("jet2 does not exist, returning");
859  continue;
860  }
861 
862  //AliJetContainer *jetContPart=GetJetContainer(3);
863  jet3=jet2->ClosestJet();
864 
865  if(!jet3){
866  Printf("jet3 does not exist, returning");
867  continue;
868  }
869  //cout<<"jet 3 exists"<<jet3->Pt()<<endl;
870 
871 
872  fh2ResponseUW->Fill(jet1->Pt(),jet2->Pt());
873 
874  Double_t fraction=0;
875  if(!(fJetShapeSub==kConstSub)) fraction = jetCont->GetFractionSharedPt(jet1);
876  if(fJetShapeSub==kConstSub) fraction = jetContUS->GetFractionSharedPt(jetUS);
877  //if (fraction > 0.1) cout<<"***** hey a jet matched with fraction"<<fraction<<" "<<jet1->Pt()<<" "<<jet2->Pt()<<" "<<fCent<<endl;
878 
879  if(fraction<fMinFractionShared) continue;
880  //InputEvent()->Print();
881 
882  }
883 
885  jet3 = jet1->ClosestJet();
886  if (!jet3) continue;
887 
888  fh2ResponseUW->Fill(jet1->Pt(),jet3->Pt());
889 
890  Double_t probDensDetPart = -999., probDensPartDet = -999.;
891 
892  if (jet1->Pt()>0) probDensPartDet = (jet3->Pt()-jet1->Pt())/jet1->Pt();
893  if (jet3->Pt()>0) probDensDetPart = (jet1->Pt()-jet3->Pt())/jet3->Pt();
894 
895  fJetProbDensityDetPart->Fill(probDensDetPart,jet3->Pt());
896  fJetProbDensityPartDet->Fill(probDensPartDet,jet1->Pt());
897 
898  }
899 
900  if (fJetShapeType == kGenOnTheFly){
901  const AliEmcalPythiaInfo *partonsInfo = 0x0;
902  partonsInfo = GetPythiaInfo();
903  Double_t jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi6());
904  Double_t detap1=(jet1->Eta())-(partonsInfo->GetPartonEta6());
905  kWeight=partonsInfo->GetPythiaEventWeight();
906  fh2ResponseW->Fill(jet1->Pt(),jet1->Pt(),kWeight);
907 
908  Float_t dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
909  fEtaJetCorr6->Fill(jet1->Eta(), partonsInfo->GetPartonEta6());
910  fPhiJetCorr6->Fill(jet1->Phi(), partonsInfo->GetPartonPhi6());
911  if(dRp1 < fRMatching) {
912  fShapesVar[0] = partonsInfo->GetPartonFlag6();
913  fPtJetCorr ->Fill(partonsInfo->GetPartonPt6(), jet1->Pt());
914  }
915  else {
916  jp1=RelativePhi(jet1->Phi(),partonsInfo->GetPartonPhi7());
917  detap1=(jet1->Eta())-(partonsInfo->GetPartonEta7());
918  dRp1 = TMath::Sqrt(jp1 * jp1 + detap1 * detap1);
919  fEtaJetCorr7->Fill(jet1->Eta(), partonsInfo->GetPartonEta7());
920  fPhiJetCorr7->Fill(jet1->Phi(), partonsInfo->GetPartonPhi7());
921  if(dRp1 < fRMatching) {
922  fShapesVar[0] = partonsInfo->GetPartonFlag7();
923  fPtJetCorr->Fill(partonsInfo->GetPartonPt7(), jet1->Pt());
924  }
925  else fShapesVar[0]=0;
926  }
927  }
928 
929  Double_t ptSubtracted = 0;
930  if (fJetShapeSub==kConstSub) ptSubtracted= jet1->Pt();
931 
932  else if (fJetShapeSub==kDerivSub) {
933  ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
934  }
935 
936  else if (fJetShapeSub==kNoSub) {
937  if ((fJetShapeType==kData) || (fJetShapeType==kDetEmbPartPythia)) ptSubtracted=jet1->Pt()-GetRhoVal(0)*jet1->Area();
938  else if ((fJetShapeType==kPythiaDef) || (fJetShapeType==kMCTrue) || (fJetShapeType==kGenOnTheFly)) ptSubtracted= jet1->Pt();
939  }
940 
941  //Printf("ptSubtracted=%f,fPtThreshold =%f ", ptSubtracted, fPtThreshold);
942  if (ptSubtracted < fPtThreshold) continue;
943 
944  if (fOneConstSelectOn == kTRUE) fNbOfConstvspT->Fill(GetJetNumberOfConstituents(jet1,0), ptSubtracted);
945  if ((fCentSelectOn == kFALSE) && (jet1->GetNumberOfTracks() <= 1)) continue;
946 
947  fShapesVar[1] = ptSubtracted;
948  fShapesVar[2] = GetJetpTD(jet1,0);
949  fShapesVar[3] = GetJetMass(jet1,0);
950  fShapesVar[4] = GetJetAngularity(jet1,0);
951  fShapesVar[5] = GetJetCircularity(jet1,0);
952  fShapesVar[6] = GetJetLeSub(jet1,0);
953  fShapesVar[7] = GetJetCoronna(jet1,0);
954 
955  Float_t ptMatch=0., ptDMatch=0., massMatch=0., angulMatch=0.,circMatch=0., lesubMatch=0., coronnaMatch=0;
956  //Float constMatch=0., sigma2Match=0.;
957  Int_t kMatched = 0;
958 
959  if (fJetShapeType==kPythiaDef) {
960  kMatched =1;
961  if(fJetShapeSub==kConstSub) kMatched = 3;
962 
963  if (!jet3) {
964  Printf("jet3 does not exist, returning");
965  continue;
966  }
967 
968  ptMatch=jet3->Pt();
969  ptDMatch=GetJetpTD(jet3, kMatched);
970  massMatch=GetJetMass(jet3,kMatched);
971  //constMatch=1.*GetJetNumberOfConstituents(jet2,kMatched);
972  angulMatch=GetJetAngularity(jet3, kMatched);
973  circMatch=GetJetCircularity(jet3, kMatched);
974  lesubMatch=GetJetLeSub(jet3, kMatched);
975  coronnaMatch=GetJetCoronna(jet3,kMatched);
976  //sigma2Match = GetSigma2(jet2, kMatched);
977  }
978 
980  if(fJetShapeSub==kConstSub) kMatched = 3;
981  if(fJetShapeSub==kDerivSub) kMatched = 2;
982  ptMatch=jet3->Pt();
983  ptDMatch=GetJetpTD(jet3, kMatched);
984  massMatch=GetJetMass(jet3,kMatched);
985  // constMatch=1.*GetJetNumberOfConstituents(jet3,kMatched);
986  angulMatch=GetJetAngularity(jet3, kMatched);
987  circMatch=GetJetCircularity(jet3, kMatched);
988  lesubMatch=GetJetLeSub(jet3, kMatched);
989  coronnaMatch = GetJetCoronna(jet3, kMatched);
990  }
991 
992 
993 
995  kMatched = 0;
996  ptMatch=0.;
997  ptDMatch=0.;
998  massMatch=0.;
999  //constMatch=0.;
1000  angulMatch=0.;
1001  circMatch=0.;
1002  lesubMatch=0.;
1003  coronnaMatch =0.;
1004 
1005  }
1006 
1007  fShapesVar[8] = ptMatch;
1008  fShapesVar[9] = ptDMatch;
1009  fShapesVar[10] = massMatch;
1010  fShapesVar[11] = angulMatch;
1011  fShapesVar[12] = circMatch;
1012  fShapesVar[13] = lesubMatch;
1013  fShapesVar[14] = coronnaMatch;
1014  fShapesVar[15] = kWeight;
1015  fShapesVar[16] = rhoVal;
1016  fShapesVar[17] = rhoMassVal;
1017  fShapesVar[18] = jet1->Pt();
1018 
1019  Int_t nInclusiveElectrons = 0, nPhotonicElectrons = 0, nTrueElectronsMC= 0, nTrueHFElecMC= 0;
1020  Double_t pElec = 0., ptElec = 0.;
1021  Bool_t hasElectrons = kFALSE;
1022 
1023  GetNumberOfElectrons(jet1, 0,nMotherKink,listofmotherkink,nInclusiveElectrons,nPhotonicElectrons,pElec,ptElec,hasElectrons);
1024 
1025  // generated HFE jets
1026 
1027  AliVParticle *vp1 = 0x0;
1028  Int_t pdgCode = 0;
1029  Int_t elecCounter = 0;
1030 
1031  if(fMCarray){
1032 
1033  GetNumberOfTrueElectrons(jet1, 0,nMotherKink,listofmotherkink,nTrueElectronsMC,nTrueHFElecMC);
1034 
1035  for(UInt_t i = 0; i < jet1->GetNumberOfTracks(); i++) {
1036  vp1 = static_cast<AliVParticle*>(jet1->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1037 
1038  if (!vp1){
1039  Printf("AliVParticle associated to constituent not found");
1040  continue;
1041  }
1042 
1043  pdgCode = vp1->PdgCode();
1044  if (TMath::Abs(pdgCode)==11) elecCounter++;
1045  }
1046  if (elecCounter==1) fPtGenJet->Fill(jet1->Pt());
1047  }
1048 
1049  fShapesVar[19] = pElec;
1050  fShapesVar[20] = ptElec;
1051  fShapesVar[21] = nInclusiveElectrons;
1052  fShapesVar[22] = nPhotonicElectrons;
1053  fShapesVar[23] = hasElectrons;
1054  fShapesVar[24] = nTrueElectronsMC;
1055  fShapesVar[25] = nTrueHFElecMC;
1056 
1057  fTreeObservableTagging->Fill();
1058 
1059 
1060  }//jet loop
1061  }
1062  return kTRUE;
1063 }
1064 
1065 //________________________________________________________________________
1067  //calc subtracted jet mass
1068  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
1069  if (fDerivSubtrOrder == 1) return jet->GetShapeProperties()->GetFirstOrderSubtracted();
1070  else return jet->GetShapeProperties()->GetSecondOrderSubtracted();
1071  else
1072  return jet->M();
1073 }
1074 
1075 //________________________________________________________________________
1076 void AliAnalysisTaskEmcalHfeTagging::GetNumberOfElectrons(AliEmcalJet *jet, Int_t jetContNb , Int_t nMother, Double_t listMother[] , Int_t &nIncElec, Int_t &nPhotElec, Double_t &pElec, Double_t &ptElec, Bool_t &hasElec){
1077  // count the number of inclusive electrons per jet and per event
1078 
1079  AliVParticle *vp1 = 0x0;
1080  Int_t nIE=0, nPairs=0, nPE=0, sub = -1;
1081  Float_t ratioElec = -1.;
1082  Double_t pte=0., pe=0.;
1083  Bool_t hasElecCand = kFALSE;
1084 
1085  Double_t jetPt = jet->Pt();
1086 
1087  Double_t ptRange[34] = {0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4,
1088  1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 5.,
1089  6., 8., 10., 12., 14., 16., 19., 22., 26., 30.,
1090  35., 40., 45., 50.};
1091  Double_t ptJetRange[6] = {5,20,40,60,80,120};
1092 
1093  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1094  if (jet->GetNumberOfTracks()){
1095 
1096  fnPartPerJet->Fill(jet->GetNumberOfTracks());
1097  fpidResponse = fInputHandler->GetPIDResponse();
1098  if(!fpidResponse){
1099  printf("ERROR: pid response not available\n");
1100  }
1101 
1102  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
1103  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1104 
1105  if (!vp1){
1106  Printf("AliVParticle associated to constituent not found");
1107  continue;
1108  }
1109 
1110  AliVTrack *vtrack = dynamic_cast<AliVTrack*>(vp1);
1111  if (!vtrack) {
1112  printf("ERROR: Could not receive track%d\n", i);
1113  continue;
1114  }
1115 
1116  AliAODTrack *track = dynamic_cast<AliAODTrack*>(vtrack);
1117  if (!track) continue;
1118 
1119  // track cuts for electron identification
1120  Bool_t passTrackCut = kFALSE;
1121  passTrackCut = InclElecTrackCuts(pVtx,track,nMother,listMother);
1122  if (!passTrackCut) continue;
1123 
1124  Double_t p=-9., pt=-9., fTPCnSigma=-99., fTOFnSigma=-99., phi = -9., eta = -99.;
1125  p = track->P();
1126  pt = track->Pt();
1127  phi = track->Phi();
1128  eta = track->Eta();
1129 
1130  fPtTrack->Fill(pt);
1131  fPhiTrack->Fill(pt,phi);
1132  fEtaTrack->Fill(pt,eta);
1133  if (pt > 0.5) fEtaPhiTrack->Fill(phi,eta);
1134 
1135  nPairs=0;
1136 
1137  fTPCnSigma = fpidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
1138  fTOFnSigma = fpidResponse->NumberOfSigmasTOF(track, AliPID::kElectron);
1139 
1140  fnTPCnTOFnocut->Fill(fTPCnSigma,fTOFnSigma);
1141  fnTPCnocutP->Fill(p,fTPCnSigma);
1142  fnTOFnocutP->Fill(p,fTOFnSigma);
1143  if (TMath::Abs(fTOFnSigma)<fSigmaTOFcut){
1144  fnTPCcutP->Fill(p,fTPCnSigma);
1145  fnTPCcutPt->Fill(pt,fTPCnSigma);
1146 
1147  for (Int_t l=0;l<5;l++){// pt jet range
1148  for(Int_t k=0;k<18;k++){// pt electron range (up to 4 GeV/c)
1149  if (jetPt>=ptJetRange[l] && jetPt<ptJetRange[l+1] && p>=ptRange[k] && p<ptRange[k+1]) fnTPCSigma[l][k]->Fill(fTPCnSigma);
1150  }
1151  }
1152  }
1153 
1154 
1155  if(TMath::Abs(fTPCnSigma)<3.5) hasElecCand = kTRUE;
1156 
1157  if (TMath::Abs(fTOFnSigma)<fSigmaTOFcut) fPtP->Fill(pt,p);
1158 
1159  if (fTPCnSigma>fSigmaTPCcut && fTPCnSigma<3 && TMath::Abs(fTOFnSigma)<fSigmaTOFcut && pt > 0.5 && pt < 4){
1160  fPhiRecElecTPC->Fill(pt,phi);
1161  fEtaRecElecTPC->Fill(pt,eta);
1162  fEtaPhiRecElecTPC->Fill(phi,eta);
1163 
1164  nIE++;
1165  pte=pt;
1166  pe=p;
1167  nPairs = GetNumberOfPairs(jet,track,pVtx,nMother,listMother);
1168  if (nPairs>0) nPE++;
1169  }
1170 
1171  // Electron ID with EMCal
1172 
1173  Double_t clsE = -9., m02 = -9., m20 = -9., clsTime = -9, EovP = -9.;
1174 
1175  Double_t emcphimim = 1.39;
1176  Double_t emcphimax = 3.265;
1177 
1178  Int_t clsId = track->GetEMCALcluster();
1179  if (clsId>0){
1180  AliVCluster *cluster=0x0;
1181  cluster = (AliVCluster*)fVevent->GetCaloCluster(clsId);//tender is not used at the moment
1182 
1183  if(cluster && cluster->IsEMCAL() && phi>emcphimim && phi<emcphimax){
1184  clsE = cluster->E();
1185  m20 = cluster->GetM20();
1186  m02 = cluster->GetM02();
1187  clsTime = cluster->GetTOF()*1e+9; // ns
1188  }
1189  }
1190 
1191  EovP = clsE/p;
1192 
1193  if (fTPCnSigma>fSigmaTPCcut && fTPCnSigma<3) fnEovPelecTPCcut->Fill(pt,EovP);
1194  if (fTPCnSigma>fSigmaTPCcut && fTPCnSigma<3 && m20 > 0.03 && m20 < 0.3) fnEovPelecTPCEMCalcut->Fill(pt,EovP);
1195  if (fTPCnSigma < -3.5) fnEovPbackg->Fill(pt,EovP);
1196  fnClsE->Fill(pt,clsE);
1197  fnM20->Fill(pt,m20);
1198  fnM02->Fill(pt,m02);
1199  fnClsTime->Fill(pt,clsTime);
1200 
1201  if (clsE>0.9 && clsE<1.2 & m20 > 0.03 && m20 < 0.3 && pt > 4 && pt < 50){
1202  fPhiRecElecEMCal->Fill(pt,phi);
1203  fEtaRecElecEMCal->Fill(pt,eta);
1204  fEtaPhiRecElecEMCal->Fill(phi,eta);
1205 
1206  nIE++;
1207  pte=pt;
1208  pe=p;
1209  nPairs = GetNumberOfPairs(jet,track,pVtx,nMother,listMother);
1210  if (nPairs>0) nPE++;
1211  }
1212  }//tagged track
1213  }
1214 
1215  sub = nIE - nPE;
1216  if (jet->GetNumberOfTracks()>0) ratioElec = nIE/jet->GetNumberOfTracks();
1217  fnInclElecPerJet->Fill(nIE);
1218  fnPhotElecPerJet->Fill(nPE);
1219  fnIncSubPhotElecPerJet->Fill(sub);
1220  fnElecOverPartPerJet->Fill(ratioElec);
1221 
1222  nIncElec = nIE;
1223  nPhotElec = nPE;
1224  pElec = pe;//to be used for jets with only one IE
1225  ptElec = pte;//to be used for jets with only one IE
1226  hasElec = hasElecCand;
1227 }
1228 //________________________________________________________________________
1229 void AliAnalysisTaskEmcalHfeTagging::GetNumberOfTrueElectrons(AliEmcalJet *jet, Int_t jetContNb , Int_t nMother, Double_t listMother[] , Int_t &nTrueElec, Int_t &nTrueHFElec){
1230  // count the number of inclusive and HF electrons per jet and per event (true MC)
1231 
1232  AliVParticle *vp1 = 0x0;
1233  Int_t nIE=0, nHFE=0, nPE=0, PartId=0, nPairs=0, iDecay = 0;
1234  Double_t p=-9., pt=-9., fTPCnSigma=-99., fTOFnSigma=-99., MCweight = 1., eta = -99., phi = -99.;
1235  Double_t ptRange[34] = {0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4,
1236  1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.5, 4., 5.,
1237  6., 8., 10., 12., 14., 16., 19., 22., 26., 30.,
1238  35., 40., 45., 50.};
1239  Double_t ptJetRange[6] = {5,20,40,60,80,120};
1240 
1241  Bool_t isFromHFdecay=kFALSE;
1242  Bool_t isFromLMdecay=kFALSE;
1243  Bool_t passTrackCut = kFALSE;
1244 
1245  Double_t jetPt = jet->Pt();
1246 
1247  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1248  if (jet->GetNumberOfTracks()){
1249 
1250  fpidResponse = fInputHandler->GetPIDResponse();
1251  if(!fpidResponse){
1252  printf("ERROR: pid response not available\n");
1253  }
1254 
1255  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
1256  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1257 
1258  if (!vp1){
1259  Printf("AliVParticle associated to constituent not found");
1260  continue;
1261  }
1262 
1263  AliVTrack *vtrack = dynamic_cast<AliVTrack*>(vp1);
1264  if (!vtrack) {
1265  printf("ERROR: Could not receive track%d\n", i);
1266  continue;
1267  }
1268 
1269  AliAODTrack *track = dynamic_cast<AliAODTrack*>(vtrack);
1270  if (!track) continue;
1271 
1272  passTrackCut = kFALSE;
1273  isFromHFdecay=kFALSE;
1274  isFromLMdecay=kFALSE;
1275 
1276  p = track->P();
1277  pt = track->Pt();
1278  eta = track->Eta();
1279  phi = track->Phi();
1280 
1281  fTPCnSigma = fpidResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
1282  fTOFnSigma = fpidResponse->NumberOfSigmasTOF(track, AliPID::kElectron);
1283 
1284  nPairs=0;
1285  MCweight = 1.;
1286  iDecay = 0;
1287 
1288  // Electron ID with EMCal
1289 
1290  Double_t EovP = -9., clsE = -9., m20 = -9.;
1291 
1292  Double_t emcphimim = 1.39;
1293  Double_t emcphimax = 3.265;
1294 
1295  Int_t clsId = track->GetEMCALcluster();
1296  if (clsId>0){
1297  AliVCluster *cluster=0x0;
1298  cluster = (AliVCluster*)fVevent->GetCaloCluster(clsId);//tender is not used at the moment
1299 
1300  if(cluster && cluster->IsEMCAL() && phi > emcphimim && phi < emcphimax){
1301  clsE = cluster->E();
1302  m20 = cluster->GetM20();
1303  }
1304  }
1305 
1306  EovP = clsE/p;
1307 
1308  if(fMCarray){
1309  Int_t label = track->GetLabel();
1310  if(label!=0){
1311  fMCparticle = (AliAODMCParticle*) fMCarray->At(TMath::Abs(track->GetLabel()));
1312  if(fMCparticle){
1313  Int_t partPDG = fMCparticle->GetPdgCode();
1314 
1315  GetWeightAndDecay(fMCparticle,iDecay,MCweight);
1316  isFromHFdecay = IsFromHFdecay(fMCparticle);
1317  isFromLMdecay = IsFromLMdecay(fMCparticle);
1318 
1319  if (TMath::Abs(partPDG)==11 && isFromHFdecay){
1320  fptTrueHFEeffTPCTOF[0]->Fill(pt);
1321  fptTrueHFEeffEMCal[0]->Fill(pt);
1322  }
1323 
1324  // track cuts
1325  passTrackCut = InclElecTrackCuts(pVtx,track,nMother,listMother);
1326  if (!passTrackCut) continue;
1327 
1328  if (TMath::Abs(partPDG)==11 && isFromHFdecay) fptTrueHFEeffTPCTOF[1]->Fill(pt);
1329  if (TMath::Abs(partPDG)==11 && isFromHFdecay ) fptTrueHFEeffEMCal[1]->Fill(pt);
1330 
1331  //check sigma_TPC in MC for different particles
1332  if (TMath::Abs(fTOFnSigma)<fSigmaTOFcut){
1333  if (TMath::Abs(partPDG)==11) PartId = 1; // electrons
1334  if (TMath::Abs(partPDG)==13) PartId = 2; // muons
1335  if (TMath::Abs(partPDG)==321) PartId = 3; // kaons
1336  if (TMath::Abs(partPDG)==2212) PartId = 4; // protons
1337  if (TMath::Abs(partPDG)==211) PartId = 5; // pions
1338 
1339  for (Int_t l=0;l<5;l++){// pt jet range
1340  for(Int_t k=0;k<18;k++){// pt electron range (up to 4 GeV/c)
1341  if (jetPt>=ptJetRange[l] && jetPt<ptJetRange[l+1] && p>=ptRange[k] && p<ptRange[k+1]) fnTPCTrueParticles[l][k]->Fill(PartId,fTPCnSigma);
1342  }
1343  }
1344  }
1345 
1346  if (TMath::Abs(partPDG)==11 && isFromHFdecay && fTPCnSigma>fSigmaTPCcut && fTPCnSigma<3) fptTrueHFEeffTPCTOF[2]->Fill(pt);
1347  if (TMath::Abs(partPDG)==11 && isFromHFdecay && TMath::Abs(fTOFnSigma)<fSigmaTOFcut) fptTrueHFEeffTPCTOF[3]->Fill(pt);
1348  if (TMath::Abs(partPDG)==11 && isFromHFdecay && fTPCnSigma>fSigmaTPCcut && fTPCnSigma<3 && TMath::Abs(fTOFnSigma)<fSigmaTOFcut) fptTrueHFEeffTPCTOF[4]->Fill(pt);
1349 
1350  if (TMath::Abs(partPDG)==11 && isFromHFdecay && m20 > 0.03 && m20 < 0.3) fptTrueHFEeffEMCal[2]->Fill(pt);
1351  if (TMath::Abs(partPDG)==11 && isFromHFdecay && clsE>0.9 && clsE<1.2) fptTrueHFEeffEMCal[3]->Fill(pt);
1352  if (TMath::Abs(partPDG)==11 && isFromHFdecay && clsE>0.9 && clsE<1.2 & m20 > 0.03 && m20 < 0.3) fptTrueHFEeffEMCal[4]->Fill(pt);
1353 
1354  if (TMath::Abs(partPDG)==11) nIE++;
1355  if (TMath::Abs(partPDG)==11 && isFromHFdecay) nHFE++;
1356  if (TMath::Abs(partPDG)==11 && iDecay>0 && iDecay<7) nPE++;
1357 
1358  // TPC-TOF
1359  if (fTPCnSigma>fSigmaTPCcut && fTPCnSigma<3 && TMath::Abs(fTOFnSigma)<fSigmaTOFcut && pt > 0.5 && pt < 4){
1360  nPairs = GetNumberOfPairs(jet,track,pVtx,nMother,listMother);
1361  if (nPairs>0) fptRecPE->Fill(pt,MCweight);
1362  if (nPairs>0 && (iDecay<1 || iDecay>6)) fptWrongPE->Fill(pt,MCweight);
1363  if (iDecay>0 && iDecay<7) fptTruePE->Fill(pt,MCweight);
1364 
1365  for (Int_t l=0;l<5;l++){// pt jet range
1366  if (jetPt>=ptJetRange[l] && jetPt<ptJetRange[l+1]){
1367  if (iDecay>0 && iDecay<7) fTotElecPt[l]->Fill(pt,MCweight);
1368  if (nPairs>0) fUlsLsElecPt[l]->Fill(pt,MCweight);
1369  }
1370  }//jet pt
1371  }// PID cuts
1372 
1373  //EMCal
1374  if (clsE>0.9 && clsE<1.2 & m20 > 0.03 && m20 < 0.3 && pt > 4 && pt < 50 ){
1375  nPairs = GetNumberOfPairs(jet,track,pVtx,nMother,listMother);
1376  if (nPairs>0) fptRecPE->Fill(pt,MCweight);
1377  if (nPairs>0 && (iDecay<1 || iDecay>6)) fptWrongPE->Fill(pt,MCweight);
1378  if (iDecay>0 && iDecay<7) fptTruePE->Fill(pt,MCweight);
1379 
1380  for (Int_t l=0;l<5;l++){// pt jet range
1381  if (jetPt>=ptJetRange[l] && jetPt<ptJetRange[l+1]){
1382  if (iDecay>0 && iDecay<7) fTotElecPt[l]->Fill(pt,MCweight);
1383  if (nPairs>0) fUlsLsElecPt[l]->Fill(pt,MCweight);
1384  }
1385  }//jet pt
1386  }// PID cuts
1387  }
1388  }
1389  }
1390  }//track loop
1391  }
1392 
1393  if (nIE==1) fptJetIE->Fill(jet->Pt());
1394  if (nIE==1 && nPE==1) fptJetPE->Fill(jet->Pt());
1395  if (nIE==1 && nHFE==1) fptJetHFE->Fill(jet->Pt());
1396 
1397  fnTrueElecPerJet->Fill(nIE);
1398  fnTrueHFElecPerJet->Fill(nHFE);
1399  fnTruePElecPerJet->Fill(nPE);
1400 
1401  nTrueElec = nIE;
1402  nTrueHFElec = nHFE;
1403 }
1404 
1405 //_________________________________________
1406 Int_t AliAnalysisTaskEmcalHfeTagging::GetNumberOfPairs(AliEmcalJet *jet, AliAODTrack *track,const AliVVertex *pVtx, Int_t nMother, Double_t listMother[])
1407 {
1408  //Get number of ULS and LS pairs per event
1409 
1410  Int_t ntracks = 0;
1411  ntracks = fVevent->GetNumberOfTracks();
1412 
1413  Int_t nULSpairs = 0;
1414  Int_t nLSpairs = 0;
1415  Int_t sub = -1;
1416 
1417  Double_t ptJetRange[6] = {5,20,40,60,80,120};
1418  Double_t jetPt = jet->Pt();
1419 
1420  for(Int_t jTracks = 0; jTracks < ntracks; jTracks++){
1421 
1422  AliVParticle* Vassotrack = fVevent->GetTrack(jTracks);
1423 
1424  if (!Vassotrack) {
1425  printf("ERROR: Could not receive associated track %d\n", jTracks);
1426  continue;
1427  }
1428  AliAODTrack *trackAsso = dynamic_cast<AliAODTrack*>(Vassotrack);
1429  if(!trackAsso) continue;
1430 
1431  if((track->GetID())==(trackAsso->GetID())) continue;
1432 
1433  // track cuts
1434  Bool_t passAssoTrackCut = kFALSE;
1435  passAssoTrackCut = PhotElecTrackCuts(pVtx,trackAsso,nMother,listMother);
1436  if (!passAssoTrackCut) continue;
1437 
1438  // tagged particle
1439  Double_t p=-9.,pt=-9.,eta =-9.,phi=-9.;
1440  Int_t charge = 0;
1441 
1442  pt = track->Pt();
1443  p = track->P();
1444  phi = track->Phi();
1445  eta = track->Eta();
1446  charge = track->Charge();
1447 
1448 
1449  // associated particle variables
1450  Double_t pAsso=-9.,ptAsso=-9.,etaAsso =-9.,fITSnSigmaAsso=-9.,fTPCnSigmaAsso=-9.,phiAsso=-9.;
1451  Int_t chargeAsso = 0;
1452 
1453  ptAsso = trackAsso->Pt();
1454  pAsso = trackAsso->P();
1455  phiAsso = trackAsso->Phi();
1456  etaAsso = trackAsso->Eta();
1457  chargeAsso = trackAsso->Charge();
1458 
1459 
1460  // looser PID cuts
1461  fITSnSigmaAsso = fpidResponse->NumberOfSigmasITS(trackAsso, AliPID::kElectron);
1462  fTPCnSigmaAsso = fpidResponse->NumberOfSigmasTPC(trackAsso, AliPID::kElectron);
1463 
1464  if(TMath::Abs(fTPCnSigmaAsso)>3.5) continue;
1465 
1466  // invariant mass
1467  Bool_t fFlagLS=kFALSE, fFlagULS=kFALSE;
1468  Double_t openingAngle = -999., mass=999., width = -999;
1469 
1470  Int_t fPDGe1 = 11; Int_t fPDGe2 = 11;
1471  if(charge>0) fPDGe1 = -11;
1472  if(chargeAsso>0) fPDGe2 = -11;
1473 
1474  if(charge == chargeAsso) fFlagLS = kTRUE;
1475  if(charge != chargeAsso) fFlagULS = kTRUE;
1476 
1477  AliKFParticle::SetField(fVevent->GetMagneticField());
1478  AliKFParticle ge1(*track, fPDGe1);
1479  AliKFParticle ge2(*trackAsso, fPDGe2);
1480  AliKFParticle recg(ge1, ge2);
1481 
1482  if(recg.GetNDF()<1) continue;
1483  Double_t chi2recg = recg.GetChi2()/recg.GetNDF();
1484  if(TMath::Sqrt(TMath::Abs(chi2recg))>3.) continue;
1485 
1486  openingAngle = ge1.GetAngle(ge2);
1487  //if(openingAngle > fOpeningAngleCut) continue;
1488 
1489  Int_t MassCorrect=-9;
1490  MassCorrect = recg.GetMass(mass,width);
1491 
1492  for (Int_t l=0;l<5;l++){// pt jet range
1493  if (jetPt>=ptJetRange[l] && jetPt<ptJetRange[l+1] && fFlagULS) fInvmassULS[l]->Fill(pt,mass);
1494  if (jetPt>=ptJetRange[l] && jetPt<ptJetRange[l+1] && fFlagLS) fInvmassLS[l]->Fill(pt,mass);
1495  }
1496 
1497  if(mass<fIMcut && fFlagULS)
1498  nULSpairs++;
1499 
1500  if(mass<fIMcut && fFlagLS)
1501  nLSpairs++;
1502 
1503  }//track loop
1504 
1505  sub = nULSpairs-nLSpairs;
1506  fnULSmLSpairsPerElectron->Fill(track->Pt(),sub);
1507 
1508  return sub;
1509 }
1510 //_________________________________________
1512 {
1513  // Check if the mother comes from heavy-flavour decays
1514  Bool_t isHFdecay = kFALSE;
1515 
1516  Int_t idMother = particle->GetMother();
1517  if (idMother>0){
1518  AliAODMCParticle* mother = (AliAODMCParticle*) fMCarray->At(idMother);
1519  Int_t motherPDG = TMath::Abs(mother->GetPdgCode());
1520 
1521  // c decays
1522  if( (((motherPDG/1000)%10) == 4) || (((motherPDG/100)%10) == 4) ) isHFdecay = kTRUE;
1523 
1524  // b decays
1525  if( (((motherPDG/1000)%10) == 5) || (((motherPDG/100)%10) == 5) ) isHFdecay = kTRUE;
1526  }
1527 
1528  return isHFdecay;
1529 }
1530 //_________________________________________
1532 {
1533  // Check if mother comes from light-meson decays
1534  Bool_t isLMdecay = kFALSE;
1535 
1536  Int_t idMother = particle->GetMother();
1537  if (idMother>0){
1538  AliAODMCParticle* mother = (AliAODMCParticle*) fMCarray->At(idMother);
1539  Int_t motherPDG = TMath::Abs(mother->GetPdgCode());
1540 
1541  if(motherPDG == 111 || motherPDG == 221 || motherPDG==223 || motherPDG==333 || motherPDG==331 ||
1542  motherPDG==113 || motherPDG==213 || motherPDG==313 || motherPDG==323) isLMdecay = kTRUE;
1543  }
1544 
1545  return isLMdecay;
1546 }
1547 //_________________________________________
1549 {
1550  // Check if particle is primary
1551  Bool_t isprimary = kFALSE;
1552  if (particle->IsPrimary()) isprimary = kTRUE;
1553 
1554  return isprimary;
1555 }
1556 //_________________________________________
1558 {
1559  // Check if the particle has mother
1560  Bool_t hasMother = kTRUE;
1561 
1562  Int_t idMother = particle->GetMother();
1563  if (idMother < 0) hasMother = kFALSE;
1564 
1565  return hasMother;
1566 }
1567 //_________________________________________
1569 {
1570  //Get Pi0 weight
1571  double weight = 0;
1572  double parPi0_enh[5] = {0,0,0,0,0};
1573 
1574  if (fRunNumber >= 252235 && fRunNumber <= 264347){//LHC17h8b (o/p)
1575  parPi0_enh[0] = 11.4919;
1576  parPi0_enh[1] = 0.021972;
1577  parPi0_enh[2] = 0.133756;
1578  parPi0_enh[3] = 1.67114;
1579  parPi0_enh[4] = 5.80327;
1580  }
1581 
1582  if (fRunNumber >= 256504 && fRunNumber <= 259888){//LHC18f4b (k/l)
1583  parPi0_enh[0] = 17.2349;
1584  parPi0_enh[1] = 0.101117;
1585  parPi0_enh[2] = 0.165525;
1586  parPi0_enh[3] = 1.52289;
1587  parPi0_enh[4] = 5.72523;
1588  }
1589 
1590  if (fMCweight==1) weight = parPi0_enh[0] / TMath::Power((exp(parPi0_enh[1]*mcPi0pT - parPi0_enh[2]*mcPi0pT*mcPi0pT) + (mcPi0pT/parPi0_enh[3])) , parPi0_enh[4]);
1591 
1592  return weight;
1593 }
1594 //_________________________________________
1596 {
1597  //Get Pi0 weight
1598  double weight = 0;
1599  double parEta_enh[5] = {0,0,0,0,0};
1600 
1601  if (fRunNumber >= 252235 && fRunNumber <= 264347){//LHC17h8b (o/p)
1602  parEta_enh[0] = 3.89483;
1603  parEta_enh[1] = 2.8197;
1604  parEta_enh[2] = 0.0662309;
1605  parEta_enh[3] = 0.00563749;
1606  parEta_enh[4] = 0.435908;
1607  }
1608 
1609  if (fRunNumber >= 256504 && fRunNumber <= 259888){//LHC18f4b (k/l)
1610  parEta_enh[0] = 0.328435;
1611  parEta_enh[1] = -0.450847;
1612  parEta_enh[2] = 0.00972497;
1613  parEta_enh[3] = 3.13337;
1614  parEta_enh[4] = 5.91296;
1615  }
1616 
1617  if (fMCweight==1) weight = parEta_enh[0] / TMath::Power((exp(parEta_enh[1]*mcEtapT - parEta_enh[2]*mcEtapT*mcEtapT) + (mcEtapT/parEta_enh[3])), parEta_enh[4]);
1618 
1619  return weight;
1620 }
1621 //________________________________________________________________________
1622 void AliAnalysisTaskEmcalHfeTagging::GetWeightAndDecay(AliAODMCParticle *particle, Int_t &decay, Double_t &weight)
1623 {
1624  //Get decay channel and weight for pi0 and eta (MC with enhanced signal)
1625  Double_t w = 0.;
1626  Int_t d = 0;
1627  Int_t partPDG = particle->GetPdgCode();
1628 
1629  if (TMath::Abs(partPDG)==11){
1630  Int_t idMother = particle->GetMother();
1631 
1632  if (idMother>0){
1633  AliAODMCParticle* mother = (AliAODMCParticle*) fMCarray->At(idMother);
1634  Int_t motherPDG = mother->GetPdgCode();
1635  Double_t motherPt = mother->Pt();
1636 
1637  Bool_t isMotherPrimary = IsPrimary(mother);
1638  Bool_t isMotherFromHF = IsFromHFdecay(mother);
1639  Bool_t isMotherFromLM = IsFromLMdecay(mother);
1640  Bool_t hasSecondMother = HasMother(mother);
1641 
1642  if (motherPDG==111 && isMotherPrimary){// pi0 -> e
1643  if (!hasSecondMother || (!isMotherFromHF && !isMotherFromLM)){
1644  d = 1;
1645  w = GetPi0weight(motherPt);
1646  }
1647  }
1648 
1649  if (motherPDG==221 && isMotherPrimary){// eta -> e
1650  if (!hasSecondMother || (!isMotherFromHF && !isMotherFromLM)){
1651  d = 2;
1652  w = GetEtaweight(motherPt);
1653  }
1654  }
1655 
1656  Int_t idSecondMother = mother->GetMother();
1657 
1658  if (idSecondMother>0){
1659  AliAODMCParticle* secondMother = (AliAODMCParticle*) fMCarray->At(idSecondMother);
1660  Int_t secondMotherPDG = secondMother->GetPdgCode();
1661  Double_t secondMotherPt = secondMother->Pt();
1662 
1663  Bool_t isSecondMotherPrimary = IsPrimary(secondMother);
1664  Bool_t isSecondMotherFromHF = IsFromHFdecay(secondMother);
1665  Bool_t isSecondMotherFromLM = IsFromLMdecay(secondMother);
1666  Bool_t hasThirdMother = HasMother(secondMother);
1667 
1668  if (motherPDG==22 && secondMotherPDG==111 && isSecondMotherPrimary){ //pi0 -> g -> e
1669  if (!hasThirdMother || (!isSecondMotherFromHF && !isSecondMotherFromLM)){
1670  d = 3;
1671  w = GetPi0weight(secondMotherPt);
1672  }
1673  }
1674 
1675  if (motherPDG==22 && secondMotherPDG==221 && isSecondMotherPrimary){ //eta -> g -> e
1676  if (!hasThirdMother || (!isSecondMotherFromHF && !isSecondMotherFromLM)){
1677  d = 4;
1678  w = GetEtaweight(secondMotherPt);
1679  }
1680  }
1681 
1682  if (motherPDG==111 && secondMotherPDG==221 && isSecondMotherPrimary){ //eta -> pi0 -> e
1683  if (!hasThirdMother || (!isSecondMotherFromHF && !isSecondMotherFromLM)){
1684  d = 5;
1685  w = GetEtaweight(secondMotherPt);
1686  }
1687  }
1688 
1689  Int_t idThirdMother = secondMother->GetMother();
1690 
1691  if (idThirdMother>0){
1692  AliAODMCParticle* thirdMother = (AliAODMCParticle*) fMCarray->At(idThirdMother);
1693  Int_t thirdMotherPDG = thirdMother->GetPdgCode();
1694  Double_t thirdMotherPt = thirdMother->Pt();
1695 
1696  Bool_t isThirdMotherPrimary = IsPrimary(thirdMother);
1697  Bool_t isThirdMotherFromHF = IsFromHFdecay(thirdMother);
1698  Bool_t isThirdMotherFromLM = IsFromLMdecay(thirdMother);
1699  Bool_t hasFourthMother = HasMother(thirdMother);
1700 
1701  if (motherPDG==22 && secondMotherPDG==111 && thirdMotherPDG==221 && isThirdMotherPrimary){//eta->pi0->g-> e
1702  if (!hasFourthMother || (!isThirdMotherFromHF && !isThirdMotherFromLM)){
1703  d = 6;
1704  w = GetEtaweight(thirdMotherPt);
1705  }
1706  }
1707  }//third mother
1708  }//second mother
1709  }//mother
1710  }// if electron
1711  decay = d;
1712  weight = w;
1713 }
1714 //________________________________________________________________________
1716 
1717  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1718  if (!jet->GetNumberOfTracks())
1719  return 0;
1720  Double_t den=0.;
1721  Double_t num = 0.;
1722  AliVParticle *vp1 = 0x0;
1723  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
1724  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1725 
1726  if (!vp1){
1727  Printf("AliVParticle associated to constituent not found");
1728  continue;
1729  }
1730 
1731  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
1732  Double_t dr2 = (vp1->Eta()-jet->Eta())*(vp1->Eta()-jet->Eta()) + dphi*dphi;
1733  Double_t dr = TMath::Sqrt(dr2);
1734  num=num+vp1->Pt()*dr;
1735  den=den+vp1->Pt();
1736  }
1737  return num/den;
1738 }
1739 
1740 //________________________________________________________________________
1742 
1743  if((fJetShapeSub==kDerivSub) && (jetContNb==0))
1746  else
1747  return Angularity(jet, jetContNb);
1748 
1749 }
1750 
1751 //____________________________________________________________________________
1752 
1754 
1755  AliTrackContainer *PartCont = NULL;
1756  AliParticleContainer *PartContMC = NULL;
1757 
1758 
1759  if (fJetShapeSub==kConstSub){
1761  else PartCont = GetTrackContainer(1);
1762  }
1763  else{
1765  else PartCont = GetTrackContainer(0);
1766  }
1767 
1768  TClonesArray *TracksArray = NULL;
1769  TClonesArray *TracksArrayMC = NULL;
1770 
1771  if (fJetShapeType == AliAnalysisTaskEmcalHfeTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
1772  else TracksArray = PartCont->GetArray();
1773 
1775  if(!PartContMC || !TracksArrayMC) return -2;
1776  }
1777  else {
1778  if(!PartCont || !TracksArray) return -2;
1779  }
1780 
1781 
1782  AliAODTrack *Track = 0x0;
1783  Float_t sumpt=0;
1784  Int_t NTracks=0;
1785  if (fJetShapeType == AliAnalysisTaskEmcalHfeTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
1786  else NTracks = TracksArray->GetEntriesFast();
1787 
1788  for(Int_t i=0; i < NTracks; i++){
1790  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
1791  if (!Track) continue;
1792  if(TMath::Abs(Track->Eta())>0.9) continue;
1793  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
1794  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
1795  Double_t dr = TMath::Sqrt(dr2);
1796  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
1797  }
1798  }
1799  else{
1800  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
1801  if (!Track) continue;
1802  if(TMath::Abs(Track->Eta())>0.9) continue;
1803  if (Track->Pt()<0.15) continue;
1804  Double_t dphi = RelativePhi(Track->Phi(),jet->Phi());
1805  Double_t dr2 = (Track->Eta()-jet->Eta())*(Track->Eta()-jet->Eta()) + dphi*dphi;
1806  Double_t dr = TMath::Sqrt(dr2);
1807  if((dr>=0.8) && (dr<1)) sumpt=sumpt+Track->Pt();
1808 
1809  }
1810  }
1811  }
1812  return sumpt;
1813 }
1814 
1815 //________________________________________________________________________
1817 
1818  if((fJetShapeSub==kDerivSub) && (jetContNb==0)) return -2;
1819  else
1820  return Coronna(jet, jetContNb);
1821 
1822 }
1823 
1824 //________________________________________________________________________
1826 
1827  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1828  if (!jet->GetNumberOfTracks())
1829  return 0;
1830  Double_t den=0.;
1831  Double_t num = 0.;
1832  AliVParticle *vp1 = 0x0;
1833  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
1834  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1835 
1836  if (!vp1){
1837  Printf("AliVParticle associated to constituent not found");
1838  continue;
1839  }
1840 
1841  num=num+vp1->Pt()*vp1->Pt();
1842  den=den+vp1->Pt();
1843  }
1844  return TMath::Sqrt(num)/den;
1845 }
1846 
1847 //________________________________________________________________________
1849  //calc subtracted jet mass
1850  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
1852  else return jet->GetShapeProperties()->GetSecondOrderSubtractedpTD();
1853  else
1854  return PTD(jet, jetContNb);
1855 
1856 }
1857 
1858 //_____________________________________________________________________________
1860 
1861  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1862  if (!jet->GetNumberOfTracks())
1863  return 0;
1864  Double_t mxx = 0.;
1865  Double_t myy = 0.;
1866  Double_t mxy = 0.;
1867  int nc = 0;
1868  Double_t sump2 = 0.;
1869  Double_t pxjet=jet->Px();
1870  Double_t pyjet=jet->Py();
1871  Double_t pzjet=jet->Pz();
1872 
1873 
1874  //2 general normalized vectors perpendicular to the jet
1875  TVector3 ppJ1(pxjet, pyjet, pzjet);
1876  TVector3 ppJ3(- pxjet* pzjet, - pyjet * pzjet, pxjet * pxjet + pyjet * pyjet);
1877  ppJ3.SetMag(1.);
1878  TVector3 ppJ2(-pyjet, pxjet, 0);
1879  ppJ2.SetMag(1.);
1880  AliVParticle *vp1 = 0x0;
1881  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
1882  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
1883 
1884  if (!vp1){
1885  Printf("AliVParticle associated to constituent not found");
1886  continue;
1887  }
1888 
1889  TVector3 pp(vp1->Px(), vp1->Py(), vp1->Pz());
1890 
1891  //local frame
1892  TVector3 pLong = pp.Dot(ppJ1) / ppJ1.Mag2() * ppJ1;
1893  TVector3 pPerp = pp - pLong;
1894  //projection onto the two perpendicular vectors defined above
1895 
1896  Float_t ppjX = pPerp.Dot(ppJ2);
1897  Float_t ppjY = pPerp.Dot(ppJ3);
1898  Float_t ppjT = TMath::Sqrt(ppjX * ppjX + ppjY * ppjY);
1899  if(ppjT<=0) return 0;
1900 
1901  mxx += (ppjX * ppjX / ppjT);
1902  myy += (ppjY * ppjY / ppjT);
1903  mxy += (ppjX * ppjY / ppjT);
1904  nc++;
1905  sump2 += ppjT;}
1906 
1907  if(nc<2) return 0;
1908  if(sump2==0) return 0;
1909  // Sphericity Matrix
1910  Double_t ele[4] = {mxx / sump2, mxy / sump2, mxy / sump2, myy / sump2};
1911  TMatrixDSym m0(2,ele);
1912 
1913  // Find eigenvectors
1914  TMatrixDSymEigen m(m0);
1915  TVectorD eval(2);
1916  TMatrixD evecm = m.GetEigenVectors();
1917  eval = m.GetEigenValues();
1918  // Largest eigenvector
1919  int jev = 0;
1920  // cout<<eval[0]<<" "<<eval[1]<<endl;
1921  if (eval[0] < eval[1]) jev = 1;
1922  TVectorD evec0(2);
1923  // Principle axis
1924  evec0 = TMatrixDColumn(evecm, jev);
1925  Double_t compx=evec0[0];
1926  Double_t compy=evec0[1];
1927  TVector2 evec(compx, compy);
1928  Double_t circ=0;
1929  if(jev==1) circ=2*eval[0];
1930  if(jev==0) circ=2*eval[1];
1931 
1932  return circ;
1933 
1934 
1935 
1936 }
1937 
1938 //________________________________________________________________________
1940  //calc subtracted jet mass
1941 
1942  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
1945  else
1946  return Circularity(jet, jetContNb);
1947 
1948 }
1949 
1950 //________________________________________________________________________
1952 
1953  AliJetContainer *jetCont = GetJetContainer(jetContNb);
1954  if (!jet->GetNumberOfTracks())
1955  return 0;
1956  Double_t den=0.;
1957  Double_t num = 0.;
1958  AliVParticle *vp1 = 0x0;
1959  AliVParticle *vp2 = 0x0;
1960  std::vector<int> ordindex;
1961  ordindex=jet->GetPtSortedTrackConstituentIndexes(jetCont->GetParticleContainer()->GetArray());
1962  //Printf("Nbof const = %d", jet->GetNumberOfTracks());
1963  //Printf("ordindex[0] = %d, ordindex[1] = %d", ordindex[0], ordindex[1]);
1964 
1965  if(ordindex.size()<2) return -1;
1966 
1967  vp1 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[0], jetCont->GetParticleContainer()->GetArray()));
1968  if (!vp1){
1969  Printf("AliVParticle associated to Leading constituent not found");
1970  return -1;
1971  }
1972 
1973  vp2 = static_cast<AliVParticle*>(jet->TrackAt(ordindex[1], jetCont->GetParticleContainer()->GetArray()));
1974  if (!vp2){
1975  Printf("AliVParticle associated to Subleading constituent not found");
1976  return -1;
1977  }
1978 
1979 
1980  num=vp1->Pt();
1981  den=vp2->Pt();
1982  //Printf("vp1->Pt() =%f, vp2->Pt() =%f", vp1->Pt(), vp2->Pt());
1983 
1984  return num-den;
1985 }
1986 
1987 //________________________________________________________________________
1989  //calc subtracted jet mass
1990 
1991  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
1993  else return jet->GetShapeProperties()->GetSecondOrderSubtractedLeSub();
1994  else
1995  return LeSub(jet, jetContNb);
1996 
1997 }
1998 
1999 //________________________________________________________________________
2001  //calc subtracted jet mass
2002 
2003  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
2006  else
2007  return jet->GetNumberOfTracks();
2008 
2009 }
2010 
2011 
2012 //______________________________________________________________________________
2014 
2015  AliJetContainer *jetCont = GetJetContainer(jetContNb);
2016  if (!jet->GetNumberOfTracks())
2017  return 0;
2018  Double_t mxx = 0.;
2019  Double_t myy = 0.;
2020  Double_t mxy = 0.;
2021  int nc = 0;
2022  Double_t sump2 = 0.;
2023 
2024  AliVParticle *vp1 = 0x0;
2025  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {
2026  vp1 = static_cast<AliVParticle*>(jet->TrackAt(i, jetCont->GetParticleContainer()->GetArray()));
2027 
2028  if (!vp1){
2029  Printf("AliVParticle associated to constituent not found");
2030  continue;
2031  }
2032 
2033  Double_t ppt=vp1->Pt();
2034  Double_t dphi = RelativePhi(vp1->Phi(),jet->Phi());
2035 
2036  Double_t deta = vp1->Eta()-jet->Eta();
2037  mxx += ppt*ppt*deta*deta;
2038  myy += ppt*ppt*dphi*dphi;
2039  mxy -= ppt*ppt*deta*TMath::Abs(dphi);
2040  nc++;
2041  sump2 += ppt*ppt;
2042 
2043  }
2044  if(nc<2) return 0;
2045  if(sump2==0) return 0;
2046  // Sphericity Matrix
2047  Double_t ele[4] = {mxx , mxy , mxy , myy };
2048  TMatrixDSym m0(2,ele);
2049 
2050  // Find eigenvectors
2051  TMatrixDSymEigen m(m0);
2052  TVectorD eval(2);
2053  TMatrixD evecm = m.GetEigenVectors();
2054  eval = m.GetEigenValues();
2055  // Largest eigenvector
2056  int jev = 0;
2057  // cout<<eval[0]<<" "<<eval[1]<<endl;
2058  if (eval[0] < eval[1]) jev = 1;
2059  TVectorD evec0(2);
2060  // Principle axis
2061  evec0 = TMatrixDColumn(evecm, jev);
2062  Double_t compx=evec0[0];
2063  Double_t compy=evec0[1];
2064  TVector2 evec(compx, compy);
2065  Double_t sig=0;
2066  if(jev==1) sig=TMath::Sqrt(TMath::Abs(eval[0])/sump2);
2067  if(jev==0) sig=TMath::Sqrt(TMath::Abs(eval[1])/sump2);
2068 
2069  return sig;
2070 
2071 }
2072 
2073 //________________________________________________________________________
2075  //calc subtracted jet mass
2076 
2077  if((fJetShapeSub==kDerivSub)&&(jetContNb==0))
2079  else return jet->GetShapeProperties()->GetSecondOrderSubtractedSigma2();
2080  else
2081  return Sigma2(jet, jetContNb);
2082 
2083 }
2084 
2085 //________________________________________________________________________
2087 
2088  AliTrackContainer *PartCont = NULL;
2089  AliParticleContainer *PartContMC = NULL;
2090 
2091 
2092  if (fJetShapeSub==kConstSub){
2094  else PartCont = GetTrackContainer(1);
2095  }
2096  else{
2098  else PartCont = GetTrackContainer(0);
2099  }
2100 
2101  TClonesArray *TracksArray = NULL;
2102  TClonesArray *TracksArrayMC = NULL;
2103 
2104  if (fJetShapeType == AliAnalysisTaskEmcalHfeTagging::kGenOnTheFly) TracksArrayMC = PartContMC->GetArray();
2105  else TracksArray = PartCont->GetArray();
2106 
2108  if(!PartContMC || !TracksArrayMC) return -99999;
2109  }
2110  else {
2111  if(!PartCont || !TracksArray) return -99999;
2112  }
2113 
2114 
2115  AliAODTrack *Track = 0x0;
2116 
2117 
2118 
2119  //TList *trackList = new TList();
2120  Int_t triggers[100];
2121  for (Int_t iTrigger=0; iTrigger<100; iTrigger++) triggers[iTrigger] = 0;
2122  Int_t iTT = 0;
2123  Int_t NTracks=0;
2124  if (fJetShapeType == AliAnalysisTaskEmcalHfeTagging::kGenOnTheFly) NTracks = TracksArrayMC->GetEntriesFast();
2125  else NTracks = TracksArray->GetEntriesFast();
2126 
2127  for(Int_t i=0; i < NTracks; i++){
2129  if((Track = static_cast<AliAODTrack*>(PartContMC->GetAcceptParticle(i)))){
2130  if (!Track) continue;
2131  if(TMath::Abs(Track->Eta())>0.9) continue;
2132  if (Track->Pt()<0.15) continue;
2133  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
2134  triggers[iTT] = i;
2135  iTT++;
2136  }
2137  }
2138  }
2139  else{
2140  if((Track = static_cast<AliAODTrack*>(PartCont->GetAcceptTrack(i)))){
2141  if (!Track) continue;
2142  if(TMath::Abs(Track->Eta())>0.9) continue;
2143  if (Track->Pt()<0.15) continue;
2144  if ((Track->Pt() >= minpT) && (Track->Pt()< maxpT)) {
2145  triggers[iTT] = i;
2146  iTT++;
2147  }
2148  }
2149  }
2150  }
2151 
2152 
2153  if (iTT == 0) return -99999;
2154  Int_t nbRn = 0, index = 0 ;
2155  TRandom3* random = new TRandom3(0);
2156  nbRn = random->Integer(iTT);
2157  index = triggers[nbRn];
2158  //Printf("iTT Total= %d, nbRn = %d, Index = %d",iTT, nbRn, index );
2159  return index;
2160 
2161 }
2162 
2163 //__________________________________________________________________________________
2165 
2166  if (vphi < -1*TMath::Pi()) vphi += (2*TMath::Pi());
2167  else if (vphi > TMath::Pi()) vphi -= (2*TMath::Pi());
2168  if (mphi < -1*TMath::Pi()) mphi += (2*TMath::Pi());
2169  else if (mphi > TMath::Pi()) mphi -= (2*TMath::Pi());
2170  double dphi = mphi-vphi;
2171  if (dphi < -1*TMath::Pi()) dphi += (2*TMath::Pi());
2172  else if (dphi > TMath::Pi()) dphi -= (2*TMath::Pi());
2173  return dphi;//dphi in [-Pi, Pi]
2174 }
2175 
2176 
2177 //________________________________________________________________________
2179  //
2180  // retrieve event objects
2181  //
2183  return kFALSE;
2184 
2185  return kTRUE;
2186 }
2187 
2188 //_________________________________________
2189 Bool_t AliAnalysisTaskEmcalHfeTagging::InclElecTrackCuts(const AliVVertex *pVietx,AliAODTrack *ietrack,Int_t nMother, Double_t listMother[])
2190 {
2191  // track cuts for inclusive electrons
2192 
2193  if(!ietrack->TestFilterMask(AliAODTrack::kTrkGlobalNoDCA)) return kFALSE;
2194  if(TMath::Abs(ietrack->Eta())>0.7) return kFALSE;
2195  if(ietrack->GetTPCNcls() < fTPCnCut) return kFALSE;
2196  if (ietrack->GetITSNcls() < fITSncut) return kFALSE;
2197  if(!ietrack->IsOn(AliAODTrack::kITSrefit)) return kFALSE;
2198  if(!ietrack->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
2199  if(!(ietrack->HasPointOnITSLayer(0) && ietrack->HasPointOnITSLayer(1))) return kFALSE;
2200  if(ietrack->GetTPCFoundFraction() < 0.6) return kFALSE;
2201  if(ietrack->Chi2perNDF() > 4) return kFALSE;
2202 
2203  Bool_t kinkmotherpass = kTRUE;
2204  for(Int_t kinkmother = 0; kinkmother < nMother; kinkmother++) {
2205  if(ietrack->GetID() == listMother[kinkmother]) {
2206  kinkmotherpass = kFALSE;
2207  continue;
2208  }
2209  }
2210  if(!kinkmotherpass) return kFALSE;
2211 
2212  Double_t d0z0[2]={-999,-999}, cov[3];
2213 
2214  if(ietrack->PropagateToDCA(pVietx, fVevent->GetMagneticField(), 20., d0z0, cov))
2215  if(TMath::Abs(d0z0[0]) > fDcaXYcut || TMath::Abs(d0z0[1]) > fDcaZcut) return kFALSE;
2216 
2217  return kTRUE;
2218 
2219 }
2220 //_________________________________________
2221 Bool_t AliAnalysisTaskEmcalHfeTagging::PhotElecTrackCuts(const AliVVertex *pVaetx,AliAODTrack *aetrack,Int_t nMother, Double_t listMother[])
2222 {
2223  // track cuts for associate tracks of photonic electrons
2224 
2225  if(!aetrack->TestFilterMask(AliAODTrack::kTrkTPCOnly)) return kFALSE;
2226  if(aetrack->Pt() < fAssPtCut) return kFALSE;
2227  if(TMath::Abs(aetrack->Eta())>0.9) return kFALSE;
2228  if(aetrack->GetTPCNcls() < fAssTPCnCut) return kFALSE;
2229  if (fAssITSrefitCut && !(aetrack->IsOn(AliAODTrack::kITSrefit))) return kFALSE;
2230  if(!aetrack->IsOn(AliAODTrack::kTPCrefit)) return kFALSE;
2231 
2232  return kTRUE;
2233 
2234 }
2235 
2236 //_______________________________________________________________________
2238 {
2239  // Called once at the end of the analysis.
2240 
2241  // fTreeObservableTagging = dynamic_cast<TTree*>(GetOutputData(1));
2242  // if (!fTreeObservableTagging){
2243  // Printf("ERROR: fTreeObservableTagging not available");
2244  // return;
2245  // }
2246 
2247 }
2248 
Int_t charge
Float_t Angularity(AliEmcalJet *jet, Int_t jetContNb)
Double_t GetFirstOrderSubtractedAngularity() const
Float_t GetJetLeSub(AliEmcalJet *jet, Int_t jetContNb)
Double_t Area() const
Definition: AliEmcalJet.h:130
Bool_t IsFromHFdecay(AliAODMCParticle *particle)
const TString & GetRhoName() const
double Double_t
Definition: External.C:58
Double_t GetSecondOrderSubtractedSigma2() const
Definition: External.C:236
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
Bool_t FillHistograms()
Function filling histograms.
Float_t GetPartonEta6() const
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:341
AliJetContainer * GetJetContainer(Int_t i=0) const
Float_t GetPartonEta7() const
Double_t GetSecondOrderSubtractedConstituent() const
Double_t Eta() const
Definition: AliEmcalJet.h:121
Double_t Py() const
Definition: AliEmcalJet.h:107
Double_t Phi() const
Definition: AliEmcalJet.h:117
Double_t mass
Float_t GetPythiaEventWeight() const
Container with name, TClonesArray and cuts for particles.
Double_t GetPi0weight(Double_t mcPi0pT) const
Int_t GetLabel() const
Definition: AliEmcalJet.h:124
Float_t GetPartonPhi7() const
Bool_t RetrieveEventObjects()
Retrieve common objects from event.
Bool_t PhotElecTrackCuts(const AliVVertex *pVtx, AliAODTrack *aetrack, Int_t nMother, Double_t listMother[])
Float_t Circularity(AliEmcalJet *jet, Int_t jetContNb)
Bool_t HasMother(AliAODMCParticle *particle)
Int_t GetNumberOfPairs(AliEmcalJet *jet, AliAODTrack *track, const AliVVertex *pVtx, Int_t nMother, Double_t listMother[])
Float_t GetJetCircularity(AliEmcalJet *jet, Int_t jetContNb)
Int_t GetPartonFlag7() const
Container for particles within the EMCAL framework.
Float_t Sigma2(AliEmcalJet *jet, Int_t jetContNb)
Float_t GetPartonPhi6() const
Int_t TrackAt(Int_t idx) const
Definition: AliEmcalJet.h:160
UShort_t GetNumberOfTracks() const
Definition: AliEmcalJet.h:139
TString kData
Declare data MC or deltaAOD.
Double_t Px() const
Definition: AliEmcalJet.h:106
Float_t PTD(AliEmcalJet *jet, Int_t jetContNb)
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
const TString & GetRhoMassName() const
Int_t GetPartonFlag6() const
Double_t GetEtaweight(Double_t mcEtapT) const
AliParticleContainer * GetParticleContainer() const
Double_t GetFirstOrderSubtractedConstituent() const
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
Bool_t IsPrimary(AliAODMCParticle *particle)
Int_t GetNJets() const
Float_t GetJetMass(AliEmcalJet *jet, Int_t jetContNb)
Double_t GetSecondOrderSubtractedAngularity() const
Bool_t InclElecTrackCuts(const AliVVertex *pVtx, AliAODTrack *ietrack, Int_t nMother, Double_t listMother[])
Double_t RelativePhi(Double_t mphi, Double_t vphi)
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
Double_t fCent
!event centrality
AliEmcalJet * GetNextAcceptJet()
void GetNumberOfElectrons(AliEmcalJet *jet, Int_t jetContNb, Int_t nMother, Double_t listMother[], Int_t &nIncElec, Int_t &nPhotElec, Double_t &pElec, Double_t &ptElec, Bool_t &hasElec)
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Float_t LeSub(AliEmcalJet *jet, Int_t jetContNb)
Double_t Pt() const
Definition: AliEmcalJet.h:109
Double_t GetRhoVal(Int_t i=0) const
void GetWeightAndDecay(AliAODMCParticle *particle, Int_t &decay, Double_t &weight)
Double_t GetFirstOrderSubtractedCircularity() const
void Track()
decay
Definition: HFPtSpectrum.C:41
AliEmcalList * fOutput
!output list
Float_t GetJetCoronna(AliEmcalJet *jet, Int_t jetContNb)
Double_t GetSecondOrderSubtractedCircularity() const
AliTrackContainer * GetTrackContainer(Int_t i=0) const
Store some informaion about a Pythia eventThis class is used to store some information about a Pythia...
void SetMakeGeneralHistograms(Bool_t g)
Float_t GetJetAngularity(AliEmcalJet *jet, Int_t jetContNb)
Float_t GetSigma2(AliEmcalJet *jet, Int_t jetContNb)
virtual AliVTrack * GetAcceptTrack(Int_t i=-1) const
Base task in the EMCAL jet framework.
Represent a jet reconstructed using the EMCal jet framework.
Definition: AliEmcalJet.h:51
void GetNumberOfTrueElectrons(AliEmcalJet *jet, Int_t jetContNb, Int_t nMother, Double_t listMother[], Int_t &nTrueElec, Int_t &nTrueHFElec)
const AliEmcalPythiaInfo * GetPythiaInfo() const
Double_t Pz() const
Definition: AliEmcalJet.h:108
Double_t GetFractionSharedPt(const AliEmcalJet *jet, AliParticleContainer *cont2=0x0) const
Declaration of class AliEmcalPythiaInfo.
const char Option_t
Definition: External.C:48
void UserCreateOutputObjects()
Main initialization function on the worker.
Double_t GetFirstOrderSubtractedSigma2() const
bool Bool_t
Definition: External.C:53
Bool_t IsFromLMdecay(AliAODMCParticle *particle)
AliEmcalJetShapeProperties * GetShapeProperties() const
Definition: AliEmcalJet.h:375
Int_t SelectTrigger(Float_t minpT, Float_t maxpT)
Float_t GetJetpTD(AliEmcalJet *jet, Int_t jetContNb)
Float_t Coronna(AliEmcalJet *jet, Int_t jetContNb)
Float_t GetPartonPt7() const
Double_t M() const
Definition: AliEmcalJet.h:120
Container for jet within the EMCAL jet framework.
Definition: External.C:196
std::vector< int > GetPtSortedTrackConstituentIndexes(TClonesArray *tracks) const
Float_t GetPartonPt6() const
Float_t GetJetNumberOfConstituents(AliEmcalJet *jet, Int_t jetContNb)
AliEmcalJet * GetJet(Int_t i) const