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