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