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