AliPhysics  master (3d17d9d)
AliAnalysisTaskHFJetIPQA.cxx
Go to the documentation of this file.
1 #include "TList.h"
2 #include "TMatrixD.h"
3 #include "TParticle.h"
4 #include "TMath.h"
5 #include "TRandom3.h"
6 #include "TVector3.h"
7 #include "TGraph.h"
8 #include "TFile.h"
9 #include "TVector3.h"
10 #include "TLorentzVector.h"
11 #include "Math/SVector.h"
12 #include "Math/SMatrix.h"
13 #include "TCanvas.h"
14 #include "TPaveText.h"
15 #include "AliEmcalJet.h"
16 #include "AliParticleContainer.h"
17 #include "AliAnalysisUtils.h"
18 #include "AliExternalTrackParam.h"
19 #include "AliAODEvent.h"
20 #include "AliAODTrack.h"
21 #include "AliAODVertex.h"
22 #include "AliAODPidHF.h"
23 #include "AliAODMCParticle.h"
24 #include "AliAODTracklets.h"
25 #include "AliMCEvent.h"
26 #include "AliESDEvent.h"
27 #include "AliESDUtils.h"
28 #include "AliMCEventHandler.h"
29 #include "AliStack.h"
30 #include "AliPIDResponse.h"
31 #include "AliLog.h"
32 #include "AliAnalysisManager.h"
33 #include "AliInputEventHandler.h"
34 #include "AliVEventHandler.h"
35 #include "AliVVertex.h"
36 #include "AliVParticle.h"
37 #include "AliAODMCHeader.h"
38 #include "AliJetContainer.h"
39 #include "AliGenEventHeader.h"
40 #include "AliVertexerTracks.h"
41 #include "AliEmcalList.h"
42 #include "THnSparse.h"
43 #include "TObjectTable.h"
44 #include "AliAnalysisTaskEmcalJet.h"
45 
46 
47 //***********************************//
48 #include "AliAnalysisTaskHFJetIPQA.h"
49 //***********************************//
50 #include <vector>
51 #include <algorithm>
52 #include <stdio.h>
53 #include <stdlib.h>
54 #include "AliAnalysisHelperJetTasks.h"
55 #include "AliGenPythiaEventHeader.h"
56 #include "TChain.h"
57 #include <map>
58 using std::min;
59 using std::cout;
60 using std::endl;
61 using std::vector;
62 using std::pair;
63 using std::map;
64 ClassImp(AliAnalysisTaskHFJetIPQA)
65 
66 AliAnalysisTaskHFJetIPQA::AliAnalysisTaskHFJetIPQA():
67 AliAnalysisTaskEmcalJet(),
68 fEventCuts(0),
69 fHistManager(),
70 fEventVertex(nullptr),
71 fPidResponse(nullptr),
72 jetconrec(nullptr),
73 jetcongen(nullptr),
74 fRunSmearing(kFALSE),
75 fUsePIDJetProb(kFALSE),
76 fDoMCCorrection(kFALSE),
77 fDoUnderlyingEventSub(kFALSE),
78 fApplyV0Rej(kFALSE),
79 fDoFlavourMatching(kFALSE),
80 fParam_Smear_Sigma(1.),
81 fParam_Smear_Mean(0.),
82 fGlobalVertex(kFALSE),
83 fDoNotCheckIsPhysicalPrimary(kFALSE),
84 fDoJetProb(kFALSE),
85 fFillCorrelations(kFALSE),
86 fDoLundPlane(kFALSE),
87 fDoTCTagging(0),
88 fDoProbTagging(0),
89 fDoMCEffs(0),
90 fUseSignificance(kTRUE),
91 kTagLevel(3),
92 fFracs(0),
93 fXsectionWeightingFactor(1),
94 fProductionNumberPtHard(-1),
95 fNThresholds(1),
96 sTemplateFlavour(0),
97 fJetRadius(0.4),
98 fDaughtersRadius(1),
99 fNoJetConstituents(0),
100 fTCThresholdPtFixed(0.008),
101 fGraphMean(nullptr),
102 fGraphSigmaData(nullptr),
103 fGraphSigmaMC(nullptr),
104 fGraphXi(nullptr),
105 fGraphOmega(nullptr),
106 fK0Star(nullptr),
107 fPhi(nullptr),
108 fGeant3FlukaProton(nullptr),
109 fGeant3FlukaAntiProton(nullptr),
110 fGeant3FlukaLambda(nullptr),
111 fGeant3FlukaAntiLambda(nullptr),
112 fGeant3FlukaKMinus(nullptr),
113 h1DThresholdsFirst(0),
114 h1DThresholdsSecond(0),
115 h1DThresholdsThird(0),
116 h2DProbLookup(0),
117 h2DProbDistsUnid(0),
118 h2DProbDistsudsg(0),
119 h2DProbDistsc(0),
120 h2DProbDistsb(0),
121 h2DProbDistsudsgV0(0),
122 h2DProbDistscV0(0),
123 h2DLNProbDistsUnid(0),
124 h2DLNProbDistsudsg(0),
125 h2DLNProbDistsc(0),
126 h2DLNProbDistsb(0),
127 h2DLNProbDistsudsgV0(0),
128 h2DLNProbDistscV0(0),
129 h1DProbThresholds(0),
130 cCuts(0),
131 fh1DCutInclusive(0),
132 fh1dCutudg(0),
133 fh1dCutc(0),
134 fh1dCutb(0),
135 fh1dCuts(0),
136 fh1dTracksAccepeted(0),
137 fh1dCutsPrinted(0),
138 fHLundIterative(nullptr),
139 fhnV0InJetK0s(nullptr),
140 fhnV0InJetLambda(nullptr),
141 fhnV0InJetALambda(nullptr),
142 fh1V0CounterCentK0s(nullptr),
143 fh1V0CounterCentLambda(nullptr),
144 fh1V0CounterCentALambda(nullptr),
145 fh2dKshortMassVsPt(nullptr),
146 fh2dLamdaMassVsPt(nullptr),
147 fh2dAnLamdaMassVsPt(nullptr),
148 h1DV0FalseRec(nullptr),
149 h1DV0TrueRec(nullptr),
150 h1DV0TrueDataDef(nullptr),
151 h1DV0TrueMCDef(nullptr),
152 fh1dKshortPtMC(nullptr),
153 fh1dLamdaPtMC(nullptr),
154 fh1dAnLamdaPtMC(nullptr),
155 fh2dKshortPtVsJetPtMC(nullptr),
156 fh2dLamdaPtVsJetPtMC(nullptr),
157 fh2dAnLamdaPtVsJetPtMC(nullptr),
158 fMCArray(nullptr),
159 fMCEvent(nullptr),
160 fESDTrackCut(nullptr),
161 fVertexer(nullptr),
162 fV0CandidateArray(nullptr),
163 fMcEvtSampled(kFALSE),
164 fBackgroundFactorLinus{0},
165 fPUdsgJet(100),fPSJet(100),fPCJet(100),fPBJet(100),
166 fJetCont(10),
167 fAnalysisCuts{0},
168 fCombined(nullptr),
169 fMCglobalDCAxyShift(0.0008),
170 fMCglobalDCASmear(1),
171 fVertexRecalcMinPt(1.0),
172 fHardCutOff(0),
173 fn1_mix(-999.),
174 fn2_mix(-999.),
175 fn3_mix(-999.),
176 fIsMixSignalReady_n1(kFALSE),
177 fIsMixSignalReady_n2(kFALSE),
178 fIsMixSignalReady_n3(kFALSE),
179 fIsSameEvent_n1(kFALSE),
180 fIsSameEvent_n2(kFALSE),
181 fIsSameEvent_n3(kFALSE),
182 fUseTreeForCorrelations(kFALSE),
183 fCorrelationCrossCheck(nullptr),
184 fTREE_n1(-99.),
185 fTREE_n2(-99.),
186 fTREE_n3(-99.),
187 fTREE_pt(-1.)
188 
189 {
190  SetMakeGeneralHistograms(kTRUE);
191  SetDefaultAnalysisCuts();
192  SetDefaultV0Cuts();
193  SetNeedEmcalGeom(kFALSE);
194  SetOffTrigger(AliVEvent::kINT7);
195  SetVzRange(-10,10);
196  DefineOutput(1, AliEmcalList::Class()) ;
197 }
198 AliAnalysisTaskHFJetIPQA::AliAnalysisTaskHFJetIPQA(const char *name):
199 AliAnalysisTaskEmcalJet(name, kTRUE),
200 fEventCuts(0),
201 fHistManager(name),
202 fEventVertex(nullptr),
203 fPidResponse(nullptr),
204 jetconrec(nullptr),
205 jetcongen(nullptr),
206 fRunSmearing(kFALSE),
207 fUsePIDJetProb(kFALSE),
208 fDoMCCorrection(kFALSE),
209 fDoUnderlyingEventSub(kFALSE),
210 fApplyV0Rej(kFALSE),
211 fDoFlavourMatching(kFALSE),
212 fParam_Smear_Sigma(1.),
213 fParam_Smear_Mean(0.),
214 fGlobalVertex(kFALSE),
215 fDoNotCheckIsPhysicalPrimary(kFALSE),
216 fDoJetProb(kFALSE),
217 fFillCorrelations(kFALSE),
218 fDoLundPlane(kFALSE),
219 fDoTCTagging(0),
220 fDoProbTagging(0),
221 fDoMCEffs(0),
222 fUseSignificance(kTRUE),
223 kTagLevel(3),
224 fFracs(0),
225 fXsectionWeightingFactor(1.),
226 fProductionNumberPtHard(-1),
227 fNThresholds(1),
228 sTemplateFlavour(0),
229 fJetRadius(0.4),
230 fDaughtersRadius(1),
231 fNoJetConstituents(0),
232 fTCThresholdPtFixed(0.008),
233 fGraphMean(nullptr),
234 fGraphSigmaData(nullptr),
235 fGraphSigmaMC(nullptr),
236 fGraphXi(nullptr),
237 fGraphOmega(nullptr),
238 fK0Star(nullptr),
239 fPhi(nullptr),
240 fGeant3FlukaProton(nullptr),
241 fGeant3FlukaAntiProton(nullptr),
242 fGeant3FlukaLambda(nullptr),
243 fGeant3FlukaAntiLambda(nullptr),
244 fGeant3FlukaKMinus(nullptr),
245 h1DThresholdsFirst(0),
246 h1DThresholdsSecond(0),
247 h1DThresholdsThird(0),
248 h2DProbLookup(0),
249 h2DProbDistsUnid(0),
250 h2DProbDistsudsg(0),
251 h2DProbDistsc(0),
252 h2DProbDistsb(0),
253 h2DProbDistsudsgV0(0),
254 h2DProbDistscV0(0),
255 h2DLNProbDistsUnid(0),
256 h2DLNProbDistsudsg(0),
257 h2DLNProbDistsc(0),
258 h2DLNProbDistsb(0),
259 h2DLNProbDistsudsgV0(0),
260 h2DLNProbDistscV0(0),
261 h1DProbThresholds(0),
262 cCuts(0),
263 fh1DCutInclusive(0),
264 fh1dCutudg(0),
265 fh1dCutc(0),
266 fh1dCutb(0),
267 fh1dCuts(0),
268 fh1dTracksAccepeted(0),
269 fh1dCutsPrinted(0),
270 fHLundIterative(nullptr),
271 fhnV0InJetK0s(nullptr),
272 fhnV0InJetLambda(nullptr),
273 fhnV0InJetALambda(nullptr),
274 fh1V0CounterCentK0s(nullptr),
275 fh1V0CounterCentLambda(nullptr),
276 fh1V0CounterCentALambda(nullptr),
277 fh2dKshortMassVsPt(nullptr),
278 fh2dLamdaMassVsPt(nullptr),
279 fh2dAnLamdaMassVsPt(nullptr),
280 h1DV0FalseRec(nullptr),
281 h1DV0TrueRec(nullptr),
282 h1DV0TrueDataDef(nullptr),
283 h1DV0TrueMCDef(nullptr),
284 fh1dKshortPtMC(nullptr),
285 fh1dLamdaPtMC(nullptr),
286 fh1dAnLamdaPtMC(nullptr),
287 fh2dKshortPtVsJetPtMC(nullptr),
288 fh2dLamdaPtVsJetPtMC(nullptr),
289 fh2dAnLamdaPtVsJetPtMC(nullptr),
290 fMCArray(nullptr),
291 fMCEvent(nullptr),
292 fESDTrackCut(nullptr),
293 fVertexer(nullptr),
294 fV0CandidateArray(nullptr),
295 fMcEvtSampled(kFALSE),
296 fBackgroundFactorLinus{0},
297 fPUdsgJet(100),fPSJet(100),fPCJet(100),fPBJet(100),
298 fJetCont(10),
299 fAnalysisCuts{0},
300 fCombined(nullptr),
301 fMCglobalDCAxyShift(0.000668),
302 fMCglobalDCASmear(1),
303 fVertexRecalcMinPt(1.0),
304 fHardCutOff(0),
305 fn1_mix(-999.),
306 fn2_mix(-999.),
307 fn3_mix(-999.),
308 fIsMixSignalReady_n1(kFALSE),
309 fIsMixSignalReady_n2(kFALSE),
310 fIsMixSignalReady_n3(kFALSE),
311 fIsSameEvent_n1(kFALSE),
312 fIsSameEvent_n2(kFALSE),
313 fIsSameEvent_n3(kFALSE),
314 fUseTreeForCorrelations(kFALSE),
315 fCorrelationCrossCheck(nullptr),
316 fTREE_n1(-99.),
317 fTREE_n2(-99.),
318 fTREE_n3(-99.),
319 fTREE_pt(-1.)
320 {
321  SetNeedEmcalGeom(kFALSE);
322  SetOffTrigger(AliVEvent::kINT7);
323  SetVzRange(-10,10);
324  SetMakeGeneralHistograms(kTRUE);
325  SetDefaultAnalysisCuts();
326  SetDefaultV0Cuts();
327  DefineOutput(1, AliEmcalList::Class()) ;
328 }
329 
335 void AliAnalysisTaskHFJetIPQA::ChangeDefaultCutTo(AliAnalysisTaskHFJetIPQA::bCuts cutname, Double_t newcutvalue){
336  fAnalysisCuts[cutname] =newcutvalue;
337 }
343 void AliAnalysisTaskHFJetIPQA::SetDefaultAnalysisCuts(){
344  //DCA
345  fAnalysisCuts[bAnalysisCut_DCAJetTrack] = 0.07;
346  fAnalysisCuts[bAnalysisCut_MaxDecayLength] = 5.;
347  fAnalysisCuts[bAnalysisCut_MaxDCA_XY] = 1.;
348  fAnalysisCuts[bAnalysisCut_MaxDCA_Z] = 2.;
349 
350  //Vertex
351  fAnalysisCuts[bAnalysisCut_NContibutors] = 3 ;
352  //fAnalysisCuts[bAnalysisCut_RelError_Y] = 0.2;
353  //fAnalysisCuts[bAnalysisCut_RelError_Z] = 0.2;
354  //fAnalysisCuts[bAnalysisCut_Sigma_Y] = 0.3;
355  //fAnalysisCuts[bAnalysisCut_Sigma_Z] = 0.3;
356  //fAnalysisCuts[bAnalysisCut_SigmaDiamond] = 2.;
357  fAnalysisCuts[bAnalysisCut_MaxVtxZ] = 10.;
358  fAnalysisCuts[bAnalysisCut_Z_Chi2perNDF] =3.5*3.5;
359  fAnalysisCuts[bAnalysisCut_MinNewVertexContrib] = 1;
360 
361  //Tracks
362  fAnalysisCuts[bAnalysisCut_MinTrackPt] =.5;
363  //fAnalysisCuts[bAnalysisCut_MinTrackPtMC] =.5;
364  fAnalysisCuts[bAnalysisCut_MinTPCClus] =100;
365  fAnalysisCuts[bAnalysisCut_MinITSLayersHit] =4;
366  fAnalysisCuts[bAnalysisCut_MinTrackChi2] =5;
367  fAnalysisCuts[bAnalysisCut_HasSDD]=1;
368  fAnalysisCuts[bAnalysisCut_HasTPCrefit]=1;
369  fAnalysisCuts[bAnalysisCut_HasITSrefit]=1;
370  fAnalysisCuts[bAnalysisCut_KinkCand]=1;
371 
372  //Jet Cuts
373  fAnalysisCuts[bAnalysisCut_MinJetPt] =0; //only for settings output. Not really used as cuts are done in .C file
374  fAnalysisCuts[bAnalysisCut_MaxJetPt] =1000;
375  fAnalysisCuts[bAnalysisCut_MinJetEta] =-0.9;
376  fAnalysisCuts[bAnalysisCut_MaxJetEta] =0.9;
377 
378  //Events
379  fAnalysisCuts[bAnalysisCut_PtHardAndJetPtFactor] =3;
380 }
381 
382 void AliAnalysisTaskHFJetIPQA::SetDefaultV0Cuts(){
383  fV0Cuts[DaughMaxEta]=0.8;
384  fV0Cuts[DaughMinPt]=0.15;
385  fV0Cuts[MinDCADaughWrtPV]=0.06;
386  fV0Cuts[MaxDCADaughvsDaugh]=1;
387  fV0Cuts[IsTPCRefitOn]=1;
388  fV0Cuts[DoPosNoTPCClusters]=1;
389  fV0Cuts[MinNoCrossedTPCRows]=70;
390  fV0Cuts[NoCrossedOverNoTPCClustersMin]=0.8;
391  fV0Cuts[NoCrossedOverNoTPCClustersMax]=1000;
392  fV0Cuts[IsKinkCand]=1;
393  fV0Cuts[MaxSigmadEdxTPC]=3;
394 
395  fV0Cuts[MaxV0Eta]=0;
396  fV0Cuts[MaxV0Rap]=0.5;
397  fV0Cuts[MinDecayRadius]=0.5;
398  fV0Cuts[MaxDecayRadius]=1000;
399  fV0Cuts[MaxCosPALambda]=0.995;
400  fV0Cuts[MinCosPAK0]=0.97;
401  fV0Cuts[MaxLifeTime]=5;
402  fV0Cuts[DoArmenteros]=1;
403  fV0Cuts[DoMassWindow]=1;
404  fV0Cuts[InvarMassWindowK0]=0.01;
405  fV0Cuts[InvarMassWindowLambda]=0.005;
406 
407  fV0Cuts[fAV0Cut]=0;
408  fV0Cuts[fBV0Cut]=0;
409  fV0Cuts[fCV0Cut]=9999;
410 }
411 
412 void AliAnalysisTaskHFJetIPQA::SmearTrack(AliAODTrack *track) {
413  if(!fIsPythia) return;
414  printf("Run Track Smearing.\n");
415  // Get reconstructed track parameters
416  AliExternalTrackParam et; et.CopyFromVTrack(track);
417  Double_t *param=const_cast<Double_t*>(et.GetParameter());
418  // Double_t *covar=const_cast<Double_t*>(et.GetCovariance());
419  // Get MC info
420  Int_t imc=track->GetLabel();
421  if (imc<=0) return;
422  const AliAODMCParticle *mc=static_cast<AliAODMCParticle*>(fMCArray->At(imc));
423  Double_t mcx[3];
424  Double_t mcp[3];
425  Double_t mccv[36]={0.};
426  Short_t mcc;
427  mc->XvYvZv(mcx);
428  mc->PxPyPz(mcp);
429  mcc=mc->Charge();
430  AliExternalTrackParam mct(mcx,mcp,mccv,mcc);
431  const Double_t *parammc=mct.GetParameter();
432  AliVertex vtx(mcx,1.,1);
433  // Correct reference points and frames according to MC
434  // TODO: B-Field correct?
435  // TODO: failing propagation....
436  et.PropagateToDCA(&vtx,track->GetBz(),10.);
437  et.Rotate(mct.GetAlpha());
438  // Select appropriate smearing functions
439  Double_t sd0rpn=fParam_Smear_Sigma;
440  Double_t sd0mrpn=fParam_Smear_Mean;//mu m
441  Double_t sd0zn =1.;
442  Double_t spt1n =1.;
443  Double_t sd0rpo=1;
444  Double_t sd0mrpo=0;
445  Double_t sd0zo =1.;
446  Double_t spt1o =1.;
447  // Use the same units (i.e. cm and GeV/c)! TODO: pt!
448  sd0rpo*=1.e-4;
449  sd0zo *=1.e-4;
450  sd0rpn*=1.e-4;
451  sd0zn *=1.e-4;
452  sd0mrpo*=1.e-4;
453  sd0mrpn*=1.e-4;
454  // Apply the smearing
455  Double_t d0zo =param [1];
456  Double_t d0zmc =parammc[1];
457  Double_t d0rpo =param [0];
458  Double_t d0rpmc=parammc[0];
459  Double_t pt1o =param [4];
460  Double_t pt1mc =parammc[4];
461  Double_t dd0zo =d0zo-d0zmc;
462  Double_t dd0zn =dd0zo *(sd0zo >0. ? (sd0zn /sd0zo ) : 1.);
463  Double_t d0zn =d0zmc+dd0zn;
464  Double_t dd0rpo=d0rpo-d0rpmc;
465  Double_t dd0rpn=dd0rpo*(sd0rpo>0. ? (sd0rpn/sd0rpo) : 1.);
466  Double_t dd0mrpn=sd0mrpn-sd0mrpo;
467  Double_t d0rpn =d0rpmc+dd0rpn-dd0mrpn;
468  Double_t dpt1o =pt1o-pt1mc;
469  Double_t dpt1n =dpt1o *(spt1o >0. ? (spt1n /spt1o ) : 1.);
470  Double_t pt1n =pt1mc+dpt1n;
471  param[0]=d0rpn;
472  param[1]=d0zn ;
473  param[4]=pt1n ;
474  // Copy the smeared parameters to the AOD track
475  Double_t x[3];
476  Double_t p[3];
477  et.GetXYZ(x);
478  et.GetPxPyPz(p);
479  Double_t cv[21];
480  et.GetCovarianceXYZPxPyPz(cv);
481  track->SetPosition(x,kFALSE);
482  track->SetP(p,kTRUE);
483  track->SetCovMatrix(cv);
484  // Mark the track as "improved" with a trick (this is done with a trick using layer 7 (ie the 8th))
485  UChar_t itsClusterMap = track->GetITSClusterMap();
486  SETBIT(itsClusterMap,7);
487  track->SetITSClusterMap(itsClusterMap);
488 }
489 
490 
491 int AliAnalysisTaskHFJetIPQA::GetMCTruth(AliAODTrack * track, int &motherpdg){
492  if(!fIsPythia) return 0;
493  int pdg = 0;
494  AliAODMCParticle *pMCAOD = nullptr;
495  if(track->GetLabel()< 0) return pdg;
496  pMCAOD = static_cast<AliAODMCParticle*>(fMCArray->At(track->GetLabel()));
497  if(!(pMCAOD)) return pdg;
498  pdg = pMCAOD->PdgCode();
499  motherpdg=0;
500  AliAODMCParticle *pMCAODmother = nullptr;
501  pMCAODmother = static_cast<AliAODMCParticle*>(fMCArray->At(pMCAOD->GetMother()));
502  if(!(pMCAOD)) return pdg;
503  motherpdg =pMCAODmother->PdgCode();
504  return pdg;
505 }
506 
507 
508 
509 Bool_t AliAnalysisTaskHFJetIPQA::FillTrackHistograms(AliVTrack *track, double *dca, double *cov, double weight)
510 {
511  FillHist("fh2dTracksImpParXY",GetValImpactParameter(kXY,dca,cov),track->Pt(),1); //*this->fXsectionWeightingFactor );
512  FillHist("fh2dTracksImpParXYZ",GetValImpactParameter(kXYZ,dca,cov),track->Pt(),1); //*this->fXsectionWeightingFactor );
513  FillHist("fh2dTracksImpParZ",dca[1],track->Pt(),1); //*this->fXsectionWeightingFactor );
514  FillHist("fh2dTracksImpParXYSignificance",GetValImpactParameter(kXYSig,dca,cov),track->Pt(),1); //*this->fXsectionWeightingFactor );
515  //FillHist("fh2dTracksImpParXYZSignificance",GetValImpactParameter(kXYZSig,dca,cov),track->Pt(),1); //*this->fXsectionWeightingFactor );
516  FillHist("fh2dTracksImpParZSignificance",GetValImpactParameter(kZSig,dca,cov),track->Pt(),1); //*this->fXsectionWeightingFactor );
517  FillHist("fh1dTracksImpParXY",GetValImpactParameter(kXY,dca,cov),1); //*this->fXsectionWeightingFactor );
518  FillHist("fh1dTracksImpParXYZ",GetValImpactParameter(kXYZ,dca,cov),1); //*this->fXsectionWeightingFactor );
519  FillHist("fh1dTracksImpParXYSignificance",GetValImpactParameter(kXYSig,dca,cov),1); //*this->fXsectionWeightingFactor );
520  //FillHist("fh1dTracksImpParXYZSignificance",GetValImpactParameter(kXYZSig,dca,cov),1); //*this->fXsectionWeightingFactor );
521  if(fIsPythia){
522  FillHist("fh1dTracksImpParXY_McCorr",GetValImpactParameter(kXY,dca,cov),weight); //*this->fXsectionWeightingFactor );
523  FillHist("fh1dTracksImpParXYZ_McCorr",GetValImpactParameter(kXYZ,dca,cov),weight); //*this->fXsectionWeightingFactor );
524  FillHist("fh1dTracksImpParXYSignificance_McCorr",GetValImpactParameter(kXYSig,dca,cov),weight); //*this->fXsectionWeightingFactor );
525  //FillHist("fh1dTracksImpParXYZSignificance_McCorr",GetValImpactParameter(kXYZSig,dca,cov),weight); //*this->fXsectionWeightingFactor );
526  FillHist("fh2dTracksImpParXY_McCorr",GetValImpactParameter(kXY,dca,cov),track->Pt(),weight); //*this->fXsectionWeightingFactor );
527  FillHist("fh2dTracksImpParXYZ_McCorr",GetValImpactParameter(kXYZ,dca,cov),track->Pt(),weight); //*this->fXsectionWeightingFactor );
528  //FillHist("fh2dTracksImpParXYZSignificance_McCorr",GetValImpactParameter(kXYZSig,dca,cov),track->Pt(),weight); //*this->fXsectionWeightingFactor );
529  }
530  return kTRUE;
531 }
532 
538 void AliAnalysisTaskHFJetIPQA::localtoglobal(Double_t alpha ,Double_t* local,Double_t* global)
539 {
540  global[0] = local[0]* TMath::Sin(alpha) + local[1] * TMath::Cos(alpha);
541  global[1] = local[0]* TMath::Cos(alpha) + local[1] * TMath::Sin(-alpha);
542  global[2] = local[2];
543  return;
544 }
550 /*void AliAnalysisTaskHFJetIPQA::EventwiseCleanup(){
551  fEtaBEvt.clear();
552  fPhiBEvt.clear();
553  fEtaCEvt.clear();
554  fPhiCEvt.clear();
555  fEtaUdsgEvt.clear();
556  fPhiUdsgEvt.clear();
557  fEtaSEvt.clear();
558  fPhiSEvt.clear();
559  fMcEvtSampled = kFALSE;
560 }
561 */
562 
563 void AliAnalysisTaskHFJetIPQA::FillRecHistograms(int jetflavour, double jetpt, double eta, double phi){
564  FillHist("fh1dJetRecPt",jetpt, 1); //this->fXsectionWeightingFactor );
565  FillHist("fh1dJetRecEtaPhiAccepted",eta,phi, 1); //this->fXsectionWeightingFactor );
566  FillHist("fh1dJetRecPtAccepted",jetpt, 1); //this->fXsectionWeightingFactor );
567 
568  if(fIsPythia){
569  if(jetflavour==0) FillHist("fh1dJetRecPtUnidentified",jetpt, 1); //this->fXsectionWeightingFactor );
570  else if(jetflavour==1)FillHist("fh1dJetRecPtudsg", jetpt, 1); //this->fXsectionWeightingFactor );
571  else if(jetflavour==2)FillHist("fh1dJetRecPtc", jetpt, 1); //this->fXsectionWeightingFactor );
572  else if(jetflavour==3)FillHist("fh1dJetRecPtb", jetpt, 1); //this->fXsectionWeightingFactor );
573  else if(jetflavour==4)FillHist("fh1dJetRecPts", jetpt, 1); //this->fXsectionWeightingFactor );
574  }
575 }
576 
577 void AliAnalysisTaskHFJetIPQA::FillGenHistograms(int jetflavour, AliEmcalJet* jetgen){
578  FillHist("fh1dJetGenPt",GetPtCorrectedMC(jetgen), 1); //this->fXsectionWeightingFactor);
579  if(jetflavour ==0) FillHist("fh1dJetGenPtUnidentified",GetPtCorrectedMC(jetgen), 1); // this->fXsectionWeightingFactor );
580  else if(jetflavour ==1) FillHist("fh1dJetGenPtudsg",GetPtCorrectedMC(jetgen), 1); //this->fXsectionWeightingFactor );
581  else if(jetflavour ==2) FillHist("fh1dJetGenPtc",GetPtCorrectedMC(jetgen), 1); //this->fXsectionWeightingFactor );
582  else if(jetflavour ==3) FillHist("fh1dJetGenPtb",GetPtCorrectedMC(jetgen), 1); //this->fXsectionWeightingFactor );
583  else if(jetflavour ==4) FillHist("fh1dJetGenPts",GetPtCorrectedMC(jetgen), 1); //this->fXsectionWeightingFactor );*/
584 }
585 
586 
587 void AliAnalysisTaskHFJetIPQA::FillIPTypePtHists(int jetflavour, double jetpt, bool* nTracks){
588  //Fill histograms for jets which have largest, second largest and third largest impact parameter
589  //tracks passing the selection criterion
590 
591  for (Int_t iN = 1 ; iN <=3 ;++iN){
592  if(!nTracks[iN]) continue;
593  FillHist(Form("fh1dJetRecPt_n_%i_%s_Accepted",iN,"all"),jetpt,1); //*this->fXsectionWeightingFactor );
594 
595  if(jetflavour==0) continue;
596  FillHist(Form("fh1dJetRecPt_n_%i_%s_Accepted",iN,sTemplateFlavour[jetflavour].Data()),jetpt,1); //*this->fXsectionWeightingFactor );
597  }
598 }
599 
600 void AliAnalysisTaskHFJetIPQA::FillIPTemplateHists(double jetpt, int iN,int jetflavour, double* params){
601  const char * stype [4] = {"fh2dJetSignedImpParXY","fh2dJetSignedImpParXYSignificance","fh2dJetSignedImpParXYZ","fh2dJetSignedImpParXYZSignificance"};
602  const char * subord [3] = {"First","Second","Third"};
603 
604  for (Int_t iType = 0 ;iType <2 ;++iType){
605  TString hname = Form("%s%s",stype[iType],subord[iN]);
606  FillHist(hname.Data(),jetpt,params[iType],1);
607  if(fIsPythia){
608  TString hnameflav = Form("%s%s%s",stype[iType],sTemplateFlavour[jetflavour].Data(),subord[iN]);
609  FillHist(hnameflav.Data(),jetpt,params[iType],1);
610  }
611  }
612 }
613 
614 void AliAnalysisTaskHFJetIPQA::FillTrackIPvsPt(int isV0, double pt, double IP, int jetflavour){
615  if(jetflavour==CV0||jetflavour==UDSGV0){
616  FillHist("fh1dTracksIPvsPt_V0JetTracks", pt, IP,1);
617  //printf("Filling fh1dTracksIPvsPt_V0JetTracks: isV0=%i, pt=%f, IP=%f, jetflavour=%i",isV0, pt, IP, jetflavour);
618  if((isV0==V0MC)||(isV0==V0TrueRec)){
619  FillHist("fh1dTracksIPvsPt_V0inV0Jet", pt, IP,1);
620  //printf("Filling fh1dTracksIPvsPt_V0inV0Jet: isV0=%i, pt=%f, IP=%f, jetflavour=%i",isV0, pt, IP, jetflavour);
621  }
622  }
623  if(jetflavour==B){
624  FillHist("fh1dTracksIPvsPt_B", pt, IP,1);
625  //printf("Filling fh1dTracksIPvsPt_B: isV0=%i, pt=%f, IP=%f, jetflavour=%i",isV0, pt, IP, jetflavour);
626  if((isV0==V0MC)||(isV0==V0TrueRec)){
627  FillHist("fh1dTracksIPvsPt_V0inBJet", pt, IP,1);
628  //printf("Filling fh1dTracksIPvsPt_V0inBJet: isV0=%i, pt=%f, IP=%f, jetflavour=%i",isV0, pt, IP, jetflavour);
629  }
630  }
631 }
632 
633 
634 void AliAnalysisTaskHFJetIPQA::FillTrackTypeResHists(){
635  printf("Filling track type resolution hists");
636 
637  /* if(GetImpactParameterWrtToJet((AliAODTrack*)trackV,(AliAODEvent*)InputEvent(),jetrec,dca,cov,xyzatcda,sign)){
638  if(fEventVertex) {
639  delete fEventVertex;
640  fEventVertex =nullptr;
641  }
642  dca[0]=fabs(dca[0]);
643  Double_t cursImParXYSig =TMath::Abs(GetValImpactParameter(kXYSig,dca,cov))*sign;
644  Double_t cursImParXYZSig =TMath::Abs(GetValImpactParameter(kXYZSig,dca,cov))*sign;
645 
646  Int_t corridx=-1;double ppt;
647  (fIsPythia&&fDoMCCorrection) ? TrackWeight = GetMonteCarloCorrectionFactor(trackV,corridx,ppt) : TrackWeight =1;
648  Double_t cursImParXY =TMath::Abs(GetValImpactParameter( kXY,dca,cov))*sign;
649  Double_t cursImParXYZ =TMath::Abs(GetValImpactParameter( kXYZ,dca,cov))*sign;
650 
651  if(fIsPythia){
652  if(is_udgjet){
653  if (IsTrackAcceptedJP((AliAODTrack*)trackV,6)){
654  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_6ITShits",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
655  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_6ITShits",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
656  }
657  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,5)){
658  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_5ITShits",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
659  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_5ITShits",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
660  }
661  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,4)){
662  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_4ITShits",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
663  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_4ITShits",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
664  }
665  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,3)){
666  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_3ITShits",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
667  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_3ITShits",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
668  }
669  }
670  if(IsFromElectron((AliAODTrack*)trackV)){
671  if(is_udgjet){
672  if (IsTrackAcceptedJP((AliAODTrack*)trackV,6)){
673  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_6ITShitsElectrons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
674  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_6ITShitsElectrons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
675  }
676  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,5)){
677  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_5ITShitsElectrons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
678  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_5ITShitsElectrons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
679  }
680  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,4)){
681  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_4ITShitsElectrons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
682  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_4ITShitsElectrons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
683  }
684  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,3)){
685  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_3ITShitsElectrons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
686  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_3ITShitsElectrons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
687  }
688  }
689  }
690  else if(IsFromPion((AliAODTrack*)trackV)){
691  if(is_udgjet){
692  if (IsTrackAcceptedJP((AliAODTrack*)trackV,6)){
693  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_6ITShitsPions",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
694  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_6ITShitsPions",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
695  }
696  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,5)){
697  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_5ITShitsPions",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
698  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_5ITShitsPions",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
699  }
700  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,4)){
701  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_4ITShitsPions",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
702  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_4ITShitsPions",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
703  }
704  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,3)){
705  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_3ITShitsPions",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
706  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_3ITShitsPions",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
707  }
708  }
709  }
710  else if(IsFromKaon((AliAODTrack*)trackV)){
711  if(is_udgjet){
712  if (IsTrackAcceptedJP((AliAODTrack*)trackV,6)){
713  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_6ITShitsKaons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
714  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_6ITShitsKaons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
715  }
716  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,5)){
717  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_5ITShitsKaons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
718  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_5ITShitsKaons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
719  }
720  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,4)){
721  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_4ITShitsKaons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
722  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_4ITShitsKaons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
723  }
724  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,3)){
725  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_3ITShitsKaons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
726  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_3ITShitsKaons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
727  }
728  }
729  }
730  else if(IsFromProton((AliAODTrack*)trackV)){
731  if(is_udgjet){
732  if (IsTrackAcceptedJP((AliAODTrack*)trackV,6)){
733  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_6ITShitsProtons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
734  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_6ITShitsProtons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
735  }
736  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,5)){
737  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_5ITShitsProtons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
738  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_5ITShitsProtons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
739  }
740  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,4)){
741  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_4ITShitsProtons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
742  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_4ITShitsProtons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
743  }
744  else if (IsTrackAcceptedJP((AliAODTrack*)trackV,3)){
745  FillHist("fh2dJetSignedImpParXYSignificanceudg_light_resfunction_3ITShitsProtons",trackV->Pt(),cursImParXYSig,TrackWeight); //this->fXsectionWeightingFactor );
746  //FillHist("fh2dJetSignedImpParXYZSignificanceudg_light_resfunction_3ITShitsProtons",trackV->Pt(),cursImParXYZSig,TrackWeight); //this->fXsectionWeightingFactor );
747  }
748  }
749  }
750  //Fill jet probability ipsig histograms for template fitting
751  const char * subtype_jp [4] = {"","udsg","c","b"};*/
752 }
753 
754 double AliAnalysisTaskHFJetIPQA::DoUESubtraction(AliJetContainer* &jetcongen, AliJetContainer* &jetconrec, AliEmcalJet* &jetrec, double jetpt){
755  //_________________________
756  //Underlying Event Subtraction
757  if((!(jetconrec->GetRhoParameter() == nullptr)))
758  {
759  printf("Correct for Underlying Event.\n");
760  jetpt = jetpt - jetconrec->GetRhoVal() * jetrec->Area();
761  }
762  if(fIsPythia){
763  if (jetrec->MatchedJet()) {
764  Double_t genpt = jetrec->MatchedJet()->Pt();
765  if((!(jetcongen->GetRhoParameter() == nullptr)))
766  {
767  printf("Correct for Underlying Event.\n");
768  genpt = genpt - jetcongen->GetRhoVal() * jetrec->MatchedJet()->Area();
769  }
770  FillHist("fh2dJetGenPtVsJetRecPt",genpt,jetpt,1); // this->fXsectionWeightingFactor );
771  }
772  }
773  return jetpt;
774 }
775 
776 Bool_t AliAnalysisTaskHFJetIPQA::IsEventAccepted(AliAODEvent *ev){
777  if(!fEventCuts.AcceptEvent(ev)){
778  return kFALSE;
779  }
780 
781  //if(!fMCRejectFilter) return true;
782  if(!(fIsPythia)) return true; // Only relevant for pt-hard production
783  AliDebugStream(1) << "Using custom MC outlier rejection" << std::endl;
784  //auto partjets =GetJetContainer("mcparticles");
785  AliJetContainer * partjets=static_cast<AliJetContainer*>(fJetCollArray.At(1));
786  if(!partjets){
787  printf("No particle container found\n");
788  return true;
789  }
790 
791  // Check whether there is at least one particle level jet with pt above n * event pt-hard
792  auto jetiter = partjets->accepted();
793  auto max = std::max_element(jetiter.begin(), jetiter.end(), [](const AliEmcalJet *lhs, const AliEmcalJet *rhs ) { return lhs->Pt() < rhs->Pt(); });
794  if(max != jetiter.end()) {
795  // At least one jet found with pt > n * pt-hard
796  AliDebugStream(1) << "Found max jet with pt " << (*max)->Pt() << " GeV/c" << std::endl;
797  if((*max)->Pt() > fAnalysisCuts[bAnalysisCut_PtHardAndJetPtFactor] * fPtHard){
798  //printf("Refuse jet with jetpt=%f, fPtHard=%f, fPtHardAndJetPtFactor=%f\n",(*max)->Pt(), fPtHard,fAnalysisCuts[bAnalysisCut_PtHardAndJetPtFactor]);
799  return false;
800  }
801  }
802  return true;
803 }
804 
805 /*void AliAnalysisTaskHFJetIPQA::SetIPVals(vector <SJetIpPati > sImpPar, bool* hasIPs, double* ipval){
806  if((int)sImpPar.size()>0) hasIPs[0]=kTRUE;
807  if((int)sImpPar.size()>1) hasIPs[1]=kTRUE;
808  if((int)sImpPar.size()>2) hasIPs[2]=kTRUE;
809 
810  if(hasIPs[0]){
811  ipval[0] =sImpPar.at(0).first;
812  //printf("HasIP0, ipval[0]=%f\n", ipval[0]);
813  }
814  if(hasIPs[1]){
815  ipval[1] =sImpPar.at(1).first;
816  //printf("HasIP1, ipval[1]=%f\n",ipval[1]);
817  }
818  if(hasIPs[2]){
819  ipval[2] =sImpPar.at(2).first;
820  //printf("HasIP2, ipval[2]=%f\n", ipval[2]);
821  }
822 }*/
823 
824 void AliAnalysisTaskHFJetIPQA::GetV0Properties(SV0Cand* & sV0, AliAODv0* &v0){
825  // particle masses from PDG
826  Double_t dMassPDGK0s = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
827  Double_t dMassPDGLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
828  sV0->fPt = TMath::Sqrt(v0->Pt2V0()); // transverse momentum of V0
829 
830  Double_t dSecVtxPos[3]; // V0 vertex position {x,y,z}
831  v0->GetSecondaryVtx(dSecVtxPos);
832  Double_t dPrimVtxPos[3]; // primary vertex position {x,y,z}
833  fEventVertex->GetXYZ(dPrimVtxPos);
834  Double_t dDecayPath[3];
835  for(Int_t iPos = 0; iPos < 3; iPos++)
836  dDecayPath[iPos] = dSecVtxPos[iPos] - dPrimVtxPos[iPos]; // vector of the V0 path
837  //Double_t dDecLen = TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1] + dDecayPath[2] * dDecayPath[2]); // path length L
838  Double_t dDecLen2D = TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1]); // transverse path length R
839  Double_t dROverPt = dDecLen2D / sV0->fPt; // R/pT
840 
841  const AliAODTrack* trackPos = (AliAODTrack*)v0->GetDaughter(0); // positive daughter track
842  const AliAODTrack* trackNeg = (AliAODTrack*)v0->GetDaughter(1); // negative daughter track
843  if(!trackPos || !trackPos){
844  sV0->bDaughsMissing=kTRUE;
845  }
846  //printf("GetV0Properties: Pointers DaughPos=%p, DaughNeg=%p\n", trackPos, trackNeg);
847  //set v0 parameters
848  sV0->bOnFly=v0->GetOnFlyStatus();
849  sV0->fDCAV0DaughvsDaugh=v0->DcaV0Daughters();
850  sV0->fPA=v0->CosPointingAngle(fEventVertex);
851  sV0->fDecayRadius=TMath::Sqrt(dSecVtxPos[0] * dSecVtxPos[0] + dSecVtxPos[1] * dSecVtxPos[1]);
852  sV0->fLifetimeK0 = dMassPDGK0s * dROverPt; // m*R/pT
853  sV0->fLifetimeLambda = dMassPDGLambda * dROverPt; // m*R/pT
854  sV0->fEta=v0->Eta();
855  sV0->fRapK0=v0->RapK0Short();
856  sV0->fRapLambda=v0->RapLambda();
857  sV0->fDecayLength3D=TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1] + dDecayPath[2] * dDecayPath[2]);
858  sV0->fDecayLength2D=TMath::Sqrt(dDecayPath[0] * dDecayPath[0] + dDecayPath[1] * dDecayPath[1]);
859  sV0->fArmenterosAlpha=v0->AlphaV0();
860  sV0->fArmenterosPt=v0->PtArmV0();
861  sV0->fMassK0=v0->MassK0Short();
862  sV0->fMassLambda=v0->MassLambda();
863  sV0->fMassAntilambda=v0->MassAntiLambda();
864  sV0->fSigmaPosPion = (fPidResponse ? TMath::Abs(fPidResponse->NumberOfSigmasTPC(trackPos, AliPID::kPion)) : 0.); // difference between measured and expected signal of the dE/dx in the TPC
865  sV0->fSigmaPosProton = (fPidResponse ? TMath::Abs(fPidResponse->NumberOfSigmasTPC(trackPos, AliPID::kProton)) : 0.);
866  sV0->fSigmaNegPion = (fPidResponse ? TMath::Abs(fPidResponse->NumberOfSigmasTPC(trackNeg, AliPID::kPion)) : 0.);
867  sV0->fSigmaNegProton = (fPidResponse ? TMath::Abs(fPidResponse->NumberOfSigmasTPC(trackNeg, AliPID::kProton)) : 0.);
868 }
869 
870 void AliAnalysisTaskHFJetIPQA::GetV0DaughProperties(SV0Daugh* & sTrack,AliAODv0* &v0, bool isPos){
871  const AliAODTrack* vTrack =0x0;
872  const AliAODVertex* prodVtxDaughter =0x0;
873  if(isPos){
874  vTrack =(AliAODTrack*)v0->GetDaughter(0); // positive daughter track
875  }
876  else{
877  vTrack =(AliAODTrack*)v0->GetDaughter(1); // positive daughter track
878  }
879  if(!vTrack) AliError(" There should be daughter tracks but GetV0DaughProperties does not get them!\n");
880  prodVtxDaughter = (AliAODVertex*)(vTrack->GetProdVertex()); // production vertex of the positive daughter track
881  Int_t cTypeVtxProd = prodVtxDaughter->GetType(); // type of the production vertex
882 
883  sTrack->fPt= vTrack->Pt();
884  sTrack->fEta=vTrack->Eta();
885  sTrack->iCharge=vTrack->Charge();
886  sTrack->iCrossedTPC=vTrack->GetTPCClusterInfo(2, 1);
887  sTrack->iNoTPCCluster=Double_t(vTrack->GetTPCNclsF());
888  if(isPos){
889  sTrack->fDCAtoPV=TMath::Abs(v0->DcaPosToPrimVertex());
890  //printf("v0: dcapos=%f\n",v0->DcaPosToPrimVertex());
891  }
892  else{
893  sTrack->fDCAtoPV=TMath::Abs(v0->DcaNegToPrimVertex());
894  //printf("v0: dcaneg=%f\n",v0->DcaNegToPrimVertex());
895  }
896  sTrack->bTPCRefitOn=vTrack->IsOn(AliAODTrack::kTPCrefit);
897  sTrack->bIsKink=(cTypeVtxProd == AliAODVertex::kKink);
898 }
899 
900 Bool_t AliAnalysisTaskHFJetIPQA::IsParticleInCone(const AliVParticle* part, const AliEmcalJet* jet, Double_t dRMax) {
901 // decides whether a particle is inside a jet cone
902  if(!part || !jet) AliError(Form("Particle or Jet missing: part=%p, jet=%p\n", part,jet));
903 
904  TVector3 vecMom2(jet->Px(), jet->Py(), jet->Pz());
905  TVector3 vecMom1(part->Px(), part->Py(), part->Pz());
906  Double_t dR = vecMom2.DeltaR(vecMom1); // = sqrt(dEta*dEta+dPhi*dPhi)
907  if(dR <= dRMax) return kTRUE;// momentum vectors of part1 and part2 are closer than dRMax
908  return kFALSE;
909 }
910 
911 //==========================================
912 void AliAnalysisTaskHFJetIPQA::FillV0Candidates(Bool_t isK, Bool_t isL, Bool_t isAL, Int_t iCut/*cut index*/){
913  if(isK){
914  fh1V0CounterCentK0s->Fill(iCut);
915  }
916  if(isL){
917  fh1V0CounterCentLambda->Fill(iCut);
918  }
919  if(isAL){
920  fh1V0CounterCentALambda->Fill(iCut);
921  }
922 }
923 
924 
925 void AliAnalysisTaskHFJetIPQA::GetV0MCTrueCandidates(AliAODEvent *fAODIn){
926 
927  AliAODMCHeader* headerMC = (AliAODMCHeader*)fAODIn->FindListObject(AliAODMCHeader::StdBranchName());
928 
929  if(!headerMC){
930  AliError("No MC header found!");
931  }
932  // get position of the MC primary vertex
933  Double_t dPrimVtxMCX = headerMC->GetVtxX();
934  Double_t dPrimVtxMCY = headerMC->GetVtxY();
935  Double_t dPrimVtxMCZ = headerMC->GetVtxZ();
936 
937  AliAODMCParticle *pAOD = 0;
938  AliEmcalJet * jetMC = 0x0;
939  double fJetPt=0;
940 
941  for (Int_t i=0; i<fMCArray->GetEntriesFast(); i++) {
942  pAOD = dynamic_cast<AliAODMCParticle*>(fMCArray->At(i));
943  if (!pAOD) continue;
944 
945  // Get the distance between the production point of the MC V0 particle and the primary vertex
946  Double_t dx = dPrimVtxMCX - pAOD->Xv();
947  Double_t dy = dPrimVtxMCY - pAOD->Yv();
948  Double_t dz = dPrimVtxMCZ - pAOD->Zv();
949  Double_t dDistPrimary = TMath::Sqrt(dx * dx + dy * dy + dz * dz);
950  Bool_t bV0MCIsPrimaryDist = (dDistPrimary < 0.01); // Is close enough to be considered primary-like?
951  if(!bV0MCIsPrimaryDist) continue;
952 
953  //Ask for PDG
954  Int_t id = pAOD->GetPdgCode();
955  Bool_t bV0 = ((id==3122) || (id==-3122) || (id==310));
956  if (!bV0) { pAOD=0; continue; }
957 
958  //Daughter Properties
959  int itrackPos = pAOD->GetDaughterLabel(0); // positive daughter track
960  int itrackNeg = pAOD->GetDaughterLabel(1); // negative daughter track
961 
962  AliAODMCParticle* trackPos=(AliAODMCParticle*) GetMCTrack(itrackPos);
963  AliAODMCParticle* trackNeg=(AliAODMCParticle*) GetMCTrack(itrackNeg);
964 
965  if((!trackPos)||(!trackNeg)) continue;
966  Double_t fPosPt=trackPos->Pt();
967  Double_t fNegPt=trackNeg->Pt();
968  Double_t fPosEta=trackPos->Eta();
969  Double_t fNegEta=trackNeg->Eta();
970 
971  //acceptance cuts
972  if (TMath::Abs(pAOD->Y()) > fV0Cuts[MaxV0Rap]) { pAOD=0; continue; }
973  if (pAOD->Pt() <0.15) continue;
974  if ((TMath::Abs(fPosEta) > fV0Cuts[DaughMaxEta])||(TMath::Abs(fNegEta) > fV0Cuts[DaughMaxEta])) continue;
975  if ((fPosPt < fV0Cuts[DaughMinPt])||(fNegPt < fV0Cuts[DaughMinPt])) continue;
976 
977 
978  if(id==310) {fh1dKshortPtMC->Fill(pAOD->Pt());}
979  if(id==3122) { fh1dLamdaPtMC->Fill(pAOD->Pt());}
980  if(id==-3122) { fh1dAnLamdaPtMC->Fill(pAOD->Pt());}
981 
982  if(!jetcongen) AliError("No MC jet container available!\n");
983  jetcongen->ResetCurrentID();
984  while ((jetMC = jetcongen->GetNextAcceptJet())){
985  fJetPt= jetMC->Pt();
986  //if(!(jetcongen->GetRhoParameter() == 0x0)){
987  // fJetPt = fJetPt - jetcongen->GetRhoVal() * jetMC->Area();
988  //}
989  if(fJetPt < 5.) continue;
990  bool bIsInCone=IsParticleInCone(pAOD, jetMC, fJetRadius);
991  if(!bIsInCone) continue;
992 
993  if(id==310) {
994  fh2dKshortPtVsJetPtMC->Fill(pAOD->Pt(), fJetPt);
995  //printf("Fount MCTrue K0s: id=%i, y=%f, pt=%f\n",id,pAOD->Y(),pAOD->Pt());
996  }
997  if(id==3122) {
998  fh2dLamdaPtVsJetPtMC->Fill(pAOD->Pt(), fJetPt);
999  //printf("Fount MCTrue Lambda: id=%i, y=%f, pt=%f\n",id,pAOD->Y(),pAOD->Pt());
1000  }
1001  if(id==-3122) {
1002  fh2dAnLamdaPtVsJetPtMC->Fill(pAOD->Pt(), fJetPt);
1003  //printf("Fount MCTrue ALambda: id=%i, y=%f, pt=%f\n",id,pAOD->Y(),pAOD->Pt());
1004  }
1005  }
1006  }
1007  //delete jetMC;
1008  jetMC=NULL;
1009  pAOD=NULL;
1010  //delete pAOD;
1011  headerMC=NULL;
1012  //delete headerMC;
1013 }
1014 
1016 Int_t AliAnalysisTaskHFJetIPQA::IsV0Daughter(const AliAODTrack* track){
1017  if(!track)return kFALSE;
1018  AliAODv0* v0aod = 0x0;
1019  int posid = -1;
1020  int negid = -1;
1021  int trackid = -1;
1022 
1023  if(!fV0CandidateArray) {return kFALSE;}
1024 
1025  const AliAODTrack* trackPos =0x0;
1026  const AliAODTrack* trackNeg =0x0; // negative daughter track
1027  AliAODMCParticle *v0MC=0x0;
1028  AliAODMCParticle *v0Daugh=0x0;
1029  int tracklabel=track->GetLabel();
1030  int posdaughlabel=-99;
1031  int negdaughlabel=-99;
1032  int iMCLabelMother=-99;
1033  int iMCPdgMother=-99;
1034  int iMCPdgDaughter=-99;
1035  int iV0Tag=V0No;
1036 
1037  //printf("--------------------------Start with %i candidated------------------------\n",fV0CandidateArray->GetEntriesFast());
1038  for(int i = 0; i < fV0CandidateArray->GetEntriesFast(); ++i) {
1039  v0aod = dynamic_cast<AliAODv0*>(fV0CandidateArray->At(i));
1040  if(!v0aod){
1041  continue;
1042  }
1043  trackPos=(AliAODTrack*) (AliAODTrack*)v0aod->GetDaughter(0); // positive daughter track
1044  trackNeg=(AliAODTrack*)(AliAODTrack*)v0aod->GetDaughter(1); // positive daughter track
1045  posid = trackPos->GetID();
1046  negid = trackNeg->GetID();
1047  trackid = track->GetID();
1048 
1049  //test the matching of IDs and compare to MC truth
1050  if(fIsPythia){
1051  posdaughlabel=trackPos->GetLabel();
1052  negdaughlabel=trackNeg->GetLabel();
1053 
1054  if((tracklabel==posdaughlabel)||(tracklabel==negdaughlabel)){
1055  //printf("Matched track to daughter: tracklabel=%i, posdaughlabel=%i, negdaughlabel=%i, trackID=%i, posID=%i, negIP=%i\n", tracklabel, posdaughlabel, negdaughlabel, trackid, posid, negid);
1056  }
1057  if((tracklabel == posdaughlabel)&&(trackid!=posid) ) {
1058  AliError(Form("Mismatch of trackid=%i and posdaughter id=%i! (tracklabel=%i, daughlabel=%i) Are you using hybrid tracks?",trackid, posid, tracklabel, posdaughlabel));
1059  }
1060  if((tracklabel == negdaughlabel)&&(trackid!=negid)){
1061  AliError(Form("Mismatch of trackid=%i and negdaughter id=%i! (tracklabel=%i, daughlabel=%i) Are you using hybrid tracks?",trackid, negid, tracklabel, negdaughlabel));
1062  }
1063  }
1064 
1065  if(posid == trackid || negid == trackid) {
1066  //printf("Reject V0 candidate: posid=%i, negid=%i, trackid=%i\n", posid, negid, trackid);
1067  iV0Tag=V0Rec;
1068  }
1069  }
1070  if(fIsPythia){
1071  if(tracklabel>(fMCArray->GetEntriesFast())||(tracklabel<0)) return 0;
1072  v0Daugh=dynamic_cast<AliAODMCParticle*>(fMCArray->At(tracklabel));
1073  if(!v0Daugh) return 0;
1074  iMCLabelMother=v0Daugh->GetMother();
1075  v0MC=dynamic_cast<AliAODMCParticle*>(fMCArray->At(iMCLabelMother));
1076  if(!v0MC){
1077  return 0;
1078  }
1079  iMCPdgMother=v0MC->GetPdgCode();
1080  iMCPdgDaughter=v0Daugh->GetPdgCode();
1081  //printf("TrackLabel=%i, iMCLabelMother=%i, pdgduaghter=%i,pdgmother=%i\n",tracklabel,iMCLabelMother, iMCPdgDaughter,iMCPdgMother);
1082  //printf("Daughter print!\n");
1083  //v0Daugh->Print();
1084  //printf("Mother print!\n");
1085  //v0MC->Print();
1086  if(((iMCPdgMother==310)||(iMCPdgMother==3122)||(iMCPdgMother==-3122))&&(iV0Tag!=V0Rec)){
1087  iV0Tag=V0MC;
1088  //printf("V0MC\n");
1089  }
1090  if(((iMCPdgMother==310)||(iMCPdgMother==3122)||(iMCPdgMother==-3122))&&(iV0Tag==V0Rec)){
1091  iV0Tag=V0TrueRec;
1092  //printf("TrueRec!\n");
1093  }
1094  }
1095  //printf("--------------------------End candidates------------------------\n");
1096  return iV0Tag;
1097 }
1098 
1099 void AliAnalysisTaskHFJetIPQA::SV0Cand::Print() const{
1100  printf("----------------- V0 -----------------\n");
1101  printf("bOnFly=%i, bDaughsMissing=%i\n", bOnFly, bDaughsMissing);
1102  printf("fDCAV0DaughvsDaugh=%f, fPA=%f, fDecayRadius=%f, fLifetimeK0=%f, fLifetimeLambda=%f\n", fDCAV0DaughvsDaugh,fPA,fDecayRadius,fLifetimeK0,fLifetimeLambda);
1103  printf("fEta=%f, fPt=%f, fRapK0=%f, fRapLambda=%f, fDecayLength3D=%f, fDecayLength2D=%f\n",fEta, fPt, fRapK0, fRapLambda, fDecayLength3D, fDecayLength3D);
1104  printf("fArmenterosAlpha=%f, fArmenterosPt=%f, fMassK0=%f, fMassLambda=%f, fMassAntilambda=%f\n",fArmenterosAlpha,fArmenterosPt,fMassK0,fMassLambda,fMassAntilambda);
1105  printf("fSigmaPosPion=%f fSigmaPosProton=%f fSigmaNegPion=%f fSigmaNegProton=%f\n", fSigmaPosPion,fSigmaPosProton,fSigmaNegPion,fSigmaNegProton);
1106  printf("bIsCandidateK0s=%i ,bIsCandidateLambda=%i, bIsCandidateALambda=%i, bIsInPeakK0s=%i bIsInPeakLambda=%i bIsInPeakALambda=%i\n", bIsCandidateK0s,bIsCandidateLambda,bIsCandidateALambda,bIsInPeakK0s,bIsInPeakLambda,bIsInPeakALambda);
1107  printf("bIsInConeJet=%i, bIsInConePerp=%i, bIsInConeRnd=%i, bIsInConeMed=%i, bIsOutsideCones=%i\n",bIsInConeJet,bIsInConePerp,bIsInConeRnd,bIsInConeMed,bIsOutsideCones);
1108  printf("-----------------------------------------\n");
1109 }
1110 
1111 void AliAnalysisTaskHFJetIPQA::SV0Daugh::Print() const{
1112  printf("----------------- V0 Daugh -----------------\n");
1113  printf("fPT=%f, fEta=%f, iCharge=%i\n", fPt, fEta, iCharge);
1114  printf("iCrossedTPC=%i iNoTPCCluster=%i fDCAtoPV=%f bTPCRefitOn=%i bIsKink=%i\n", iCrossedTPC, iNoTPCCluster, fDCAtoPV, bTPCRefitOn, bIsKink);
1115  printf("-----------------------------------------\n");
1116 }
1117 
1118 /*AliAODMCParticle* AliAnalysisTaskHFJetIPQA::GetMCTrack( const AliAODTrack* track){
1119  //
1120  // return MC track
1121  //
1122  if(!fIsPythia) return NULL;
1123  if(!fMCArray) { AliError("No fMCArray"); return NULL;}
1124  Int_t nStack = fMCArray->GetEntriesFast();
1125  Int_t iLabel = track->GetLabel(); // negative label indicate poor matching quality
1126  if(iLabel<0||(iLabel > nStack) ){
1127  AliError(Form("Stupid Label given iLabel=%i\n",iLabel));
1128  return NULL;
1129  }
1130  printf("MCTrack: iLabel=%i\n",iLabel);
1131  AliAODMCParticle *mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(iLabel));
1132  return mctrack;*/
1133 //}
1134 
1135 AliAODMCParticle* AliAnalysisTaskHFJetIPQA::GetMCTrack(int iLabel){
1136  //
1137  // return MC track
1138  //
1139  if(!fIsPythia) return NULL;
1140  if(!fMCArray) { AliError("No fMCArray"); return NULL;}
1141  Int_t nStack = fMCArray->GetEntriesFast();
1142 
1143  if((iLabel < 0) || (iLabel >= nStack)){
1144  //printf("Daugh not in array range: iLabel=%i\n", iLabel);
1145  return NULL;
1146  }
1147 
1148  AliAODMCParticle *mctrack = dynamic_cast<AliAODMCParticle *>(fMCArray->At(iLabel));
1149  return mctrack;
1150 }
1151 
1152 int AliAnalysisTaskHFJetIPQA::GetV0MCVeto(AliAODEvent* fAODIn, AliAODv0* v0, bool bIsCandidateK0s,bool bIsCandidateLambda, bool bIsCandidateALambda){
1153  // PDG codes of used particles
1154  Int_t iPdgCodePion = 211;
1155  Int_t iPdgCodeProton = 2212;
1156  Int_t iPdgCodeK0s = 310;
1157  Int_t iPdgCodeLambda = 3122;
1158 
1159  AliAODMCHeader* headerMC = (AliAODMCHeader*)fAODIn->FindListObject(AliAODMCHeader::StdBranchName());
1160  Double_t dPrimVtxMCX = 0., dPrimVtxMCY = 0., dPrimVtxMCZ = 0.; // position of the MC primary vertex
1161 
1162  if(!headerMC){
1163  AliError("No MC header found!");
1164  }
1165  // get position of the MC primary vertex
1166  dPrimVtxMCX = headerMC->GetVtxX();
1167  dPrimVtxMCY = headerMC->GetVtxY();
1168  dPrimVtxMCZ = headerMC->GetVtxZ();
1169 
1170  Int_t iNTracksMC = fMCArray->GetEntriesFast();
1171  if(!(bIsCandidateK0s) && !(bIsCandidateLambda) && !(bIsCandidateALambda)) return 0; // chosen candidates with any mass
1172 
1173  const AliAODTrack* postrack = (AliAODTrack*)v0->GetDaughter(0); // positive daughter track
1174  const AliAODTrack* negtrack = (AliAODTrack*)v0->GetDaughter(1); // positive daughter track
1175  Int_t iposLabel = postrack->GetLabel();
1176  Int_t inegLabel = negtrack->GetLabel();
1177  AliAODMCParticle* particleMCDaughterPos=(AliAODMCParticle*) GetMCTrack(iposLabel);
1178  AliAODMCParticle* particleMCDaughterNeg=(AliAODMCParticle*) GetMCTrack(inegLabel);
1179  if(!particleMCDaughterNeg || !particleMCDaughterPos) return 0;
1180 
1181  Int_t iPdgCodeDaughterPos = particleMCDaughterPos->GetPdgCode();
1182  Int_t iPdgCodeDaughterNeg = particleMCDaughterNeg->GetPdgCode();
1183  Int_t iIndexMotherPos = particleMCDaughterPos->GetMother();
1184  Int_t iIndexMotherNeg = particleMCDaughterNeg->GetMother();
1185  Double_t fPosEta=particleMCDaughterPos->Eta();
1186  Double_t fPosPt=particleMCDaughterPos->Pt();
1187  Double_t fNegEta=particleMCDaughterNeg->Eta();
1188  Double_t fNegPt=particleMCDaughterNeg->Pt();
1189 
1190  if((iIndexMotherNeg < 0) || (iIndexMotherNeg >= iNTracksMC) || (iIndexMotherPos < 0) || (iIndexMotherPos >= iNTracksMC)){
1191  //printf("Mother not in array range: iIndexMotherNeg=%i, iIndexMotherPos=%i,iNTracksMC=%i\n", iIndexMotherNeg, iIndexMotherPos, iNTracksMC);
1192  return 0;
1193  }
1194  AliAODMCParticle* particleMCMotherPos = (AliAODMCParticle*)fMCArray->At(iIndexMotherNeg);
1195  AliAODMCParticle* particleMCMotherNeg = (AliAODMCParticle*)fMCArray->At(iIndexMotherPos);
1196  int posidmother=particleMCMotherPos->GetPdgCode();
1197  int negidmother=particleMCMotherNeg->GetPdgCode();
1198  /*int posndaughs=particleMCMotherPos->GetNDaughters();
1199  int negndaughs=particleMCMotherNeg->GetNDaughters();
1200  int pos0daughs=((AliAODMCParticle*)fMCArray->At(particleMCMotherPos->GetDaughterLabel(0)))->GetPdgCode();
1201  int pos1daughs=((AliAODMCParticle*)fMCArray->At(particleMCMotherPos->GetDaughterLabel(1)))->GetPdgCode();
1202  int neg0daughs=((AliAODMCParticle*)fMCArray->At(particleMCMotherNeg->GetDaughterLabel(0)))->GetPdgCode();
1203  int neg1daughs=((AliAODMCParticle*)fMCArray->At(particleMCMotherNeg->GetDaughterLabel(1)))->GetPdgCode();
1204 
1205  int pos0daughslabel=((AliAODMCParticle*)fMCArray->At(particleMCMotherPos->GetDaughterLabel(0)))->GetLabel();
1206  int pos1daughslabel=((AliAODMCParticle*)fMCArray->At(particleMCMotherPos->GetDaughterLabel(1)))->GetLabel();
1207  int neg0daughslabel=((AliAODMCParticle*)fMCArray->At(particleMCMotherNeg->GetDaughterLabel(0)))->GetLabel();
1208  int neg1daughslabel=((AliAODMCParticle*)fMCArray->At(particleMCMotherNeg->GetDaughterLabel(1)))->GetLabel();*/
1209 
1210  if((posidmother != iPdgCodeK0s) && (TMath::Abs(posidmother) != iPdgCodeLambda)&&(negidmother != iPdgCodeK0s) && (TMath::Abs(negidmother) != iPdgCodeLambda)) return 0;
1211 
1212  if(iIndexMotherNeg != iIndexMotherPos){
1213  //printf("Mothers are different: iIndexMotherNeg=%i, iIndexMotherPos=%i, neglab=%i, poslab=%i, negid=%i, posid=%i, v0idneg=%i, v0idpos=%i\n", iIndexMotherNeg,iIndexMotherPos,inegLabel, iposLabel ,iPdgCodeDaughterNeg,iPdgCodeDaughterPos,negidmother,posidmother);
1214  //printf("posndaughs=%i, pos0daughs=%i, pos1daughs=%i, negndaughs=%i, neg0daughs=%i, neg1daughs=%i\n",posndaughs, pos0daughs, pos1daughs, negndaughs, neg0daughs, neg1daughs);
1215  //printf("pos0daughslabel=%i, pos1daughslabel=%i, neg0daughslabel=%i, neg1daughslabel=%i\n",pos0daughslabel, pos1daughslabel, neg0daughslabel, neg1daughslabel);
1216  return 0;
1217  }
1218  //_____________________
1219  // Mother Properties
1220  AliAODMCParticle* particleMCMother = (AliAODMCParticle*)fMCArray->At(iIndexMotherPos);
1221  if(!particleMCMother)return 0;
1222 
1223  Int_t iPdgCodeMother = particleMCMother->GetPdgCode();
1224  Double_t dPtV0Gen = particleMCMother->Pt();
1225  Double_t dRapV0Gen = particleMCMother->Y();
1226 
1227  // Acceptance Cuts
1228  if((TMath::Abs(dRapV0Gen) > fV0Cuts[MaxV0Rap])) return 0;
1229  if(dPtV0Gen <0.15) return 0;
1230  if((TMath::Abs(fPosEta) > fV0Cuts[DaughMaxEta])||(TMath::Abs(fNegEta) > fV0Cuts[DaughMaxEta])) return 0;
1231  if((fPosPt<fV0Cuts[DaughMinPt])||(fNegPt<fV0Cuts[DaughMinPt])) return 0;
1232 
1233  // Skip not interesting particles
1234  if((iPdgCodeMother != iPdgCodeK0s) && (TMath::Abs(iPdgCodeMother) != iPdgCodeLambda)) return 0;
1235 
1236  // Check identity of the MC mother particle and the decay channel
1237  Bool_t bV0MCIsK0s = ((iPdgCodeMother == iPdgCodeK0s) && (iPdgCodeDaughterPos == +iPdgCodePion) && (iPdgCodeDaughterNeg == -iPdgCodePion));
1238  Bool_t bV0MCIsLambda = ((iPdgCodeMother == +iPdgCodeLambda) && (iPdgCodeDaughterPos == +iPdgCodeProton) && (iPdgCodeDaughterNeg == -iPdgCodePion));
1239  Bool_t bV0MCIsALambda = ((iPdgCodeMother == -iPdgCodeLambda) && (iPdgCodeDaughterPos == +iPdgCodePion) && (iPdgCodeDaughterNeg == -iPdgCodeProton));
1240 
1241  // Get the distance between production point of the MC mother particle and the primary vertex
1242  Double_t dx = dPrimVtxMCX - particleMCMother->Xv();
1243  Double_t dy = dPrimVtxMCY - particleMCMother->Yv();
1244  Double_t dz = dPrimVtxMCZ - particleMCMother->Zv();
1245  Double_t dDistPrimary = TMath::Sqrt(dx * dx + dy * dy + dz * dz);
1246  Bool_t bV0MCIsPrimaryDist = (dDistPrimary < 0.01); // Is close enough to be considered primary-like?
1247 
1248  // K0s
1249  if(bIsCandidateK0s){ // selected candidates with any mass
1250  if(bV0MCIsK0s && bV0MCIsPrimaryDist){ // well reconstructed candidates
1251  //printf("Found K0s iPdgCodeMother=%i, iPdgCodeDaughterPos=%i, iPdgCodeDaughterNeg=%i, y=%f, ptv0=%f, etadaughpos=%f, ptdaughpos=%f,etadaughneg=%f, ptdaughneg=%f\n",
1252  //iPdgCodeMother,iPdgCodeDaughterPos,iPdgCodeDaughterNeg, dRapV0Gen,dPtV0Gen,fPosEta, fPosPt, fNegEta,fNegPt);
1253  return 1;
1254  }
1255  }
1256  // Lambda
1257  if(bIsCandidateLambda){ // selected candidates with any mass
1258  if(bV0MCIsLambda && bV0MCIsPrimaryDist){ // well reconstructed candidates
1259  //printf("Found Lambda iPdgCodeMother=%i, iPdgCodeDaughterPos=%i, iPdgCodeDaughterNeg=%i, y=%f, ptv0=%f, etadaughpos=%f, ptdaughpos=%f,etadaughneg=%f, ptdaughneg=%f\n",
1260  //iPdgCodeMother,iPdgCodeDaughterPos,iPdgCodeDaughterNeg, dRapV0Gen,dPtV0Gen,fPosEta, fPosPt, fNegEta,fNegPt);
1261  return 2;
1262  }
1263  }
1264 
1265  // anti-Lambda
1266  if(bIsCandidateALambda){ // selected candidates with any mass
1267  if(bV0MCIsALambda && bV0MCIsPrimaryDist){ // well reconstructed candidates
1268  //printf("Found ALambda iPdgCodeMother=%i, iPdgCodeDaughterPos=%i, iPdgCodeDaughterNeg=%i, y=%f, ptv0=%f, etadaughpos=%f, ptdaughpos=%f,etadaughneg=%f, ptdaughneg=%f\n",
1269  //iPdgCodeMother,iPdgCodeDaughterPos,iPdgCodeDaughterNeg, dRapV0Gen,dPtV0Gen,fPosEta, fPosPt, fNegEta,fNegPt);
1270  return 3;
1271  }
1272  }
1273 
1274  return 0;
1275 }
1276 
1277 void AliAnalysisTaskHFJetIPQA::SelectV0Candidates(AliAODEvent *fAODIn){
1278  AliAODv0* v0 = 0; // pointer to V0 candidates
1279  SV0Cand* sV0=new SV0Cand();
1280  SV0Daugh* sPosDaugh=new SV0Daugh();
1281  SV0Daugh* sNegDaugh=new SV0Daugh();
1282 
1283 
1284  fV0CandidateArray->Delete();//Reset the TClonesArray
1285 
1286  // Mean lifetime
1287  Int_t iNV0s = fAODIn->GetNumberOfV0s(); // get the number of V0 candidates
1288 
1289  //printf("################## Select candidates ########################\n");
1290 
1291  for(Int_t iV0 = 0; iV0 < iNV0s; iV0++){
1292  v0 = fAODIn->GetV0(iV0); // get next candidate from the list in AOD
1293  if(!v0) continue;
1294 
1295  sV0->Reset();
1296  sPosDaugh->Reset();
1297  sNegDaugh->Reset();
1298 
1299  //sV0->Print();
1300  //sPosDaugh->Print();
1301 
1302  GetV0Properties(sV0, v0);
1303  bool hasNoDaughters=sV0->bDaughsMissing;
1304  if(hasNoDaughters) continue;
1305 
1306  GetV0DaughProperties(sPosDaugh,v0, kTRUE);
1307  GetV0DaughProperties(sNegDaugh, v0, kFALSE);
1308  //v0->Print();
1309 
1310  Int_t iCutIndex = 0; // indicator of current selection step
1311  // 1
1312  // All V0 candidates
1313  //printf("Found V0 cand\n");
1314  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1315  iCutIndex++;
1316 
1317  // Start of global cuts
1318  // 2
1319  // Reconstruction method
1320  if(sV0->bOnFly) continue;
1321  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1322  iCutIndex++;
1323 
1324  // 3
1325  // Tracks TPC OK
1326  if(sPosDaugh->iCharge == sNegDaugh->iCharge) continue;// daughters have different charge?
1327  if(sNegDaugh->iCharge != -1) continue;// daughters have expected charge?
1328  if(sPosDaugh->iCharge != 1) continue;// daughters have expected charge?
1329 
1330  if(fV0Cuts[IsTPCRefitOn]){
1331  if(!sPosDaugh->bTPCRefitOn) continue;// TPC refit is ON?
1332  if(!sNegDaugh->bTPCRefitOn) continue;
1333  }
1334 
1335  if(fV0Cuts[IsKinkCand]){
1336  if(sPosDaugh->bIsKink) continue;// kink daughter rejection
1337  if(sNegDaugh->bIsKink) continue;
1338  }
1339 
1340  if(fV0Cuts[DoPosNoTPCClusters]){
1341  if(sPosDaugh->iNoTPCCluster <= 0.) continue;
1342  if(sNegDaugh->iNoTPCCluster <= 0.) continue;
1343  }
1344 
1345  if(fV0Cuts[MinNoCrossedTPCRows] > 0.){
1346  if(sPosDaugh->iCrossedTPC < fV0Cuts[MinNoCrossedTPCRows]) continue;// Crossed TPC padrows
1347  if(sNegDaugh->iCrossedTPC < fV0Cuts[MinNoCrossedTPCRows]) continue;
1348  }
1349 
1350  if(fV0Cuts[NoCrossedOverNoTPCClustersMin] > 0.){
1351  if(sPosDaugh->iCrossedTPC / sPosDaugh->iNoTPCCluster < fV0Cuts[NoCrossedOverNoTPCClustersMin]) continue;
1352  if(sNegDaugh->iCrossedTPC / sNegDaugh->iNoTPCCluster < fV0Cuts[NoCrossedOverNoTPCClustersMin]) continue;
1353  }
1354 
1355  if(fV0Cuts[NoCrossedOverNoTPCClustersMax] > 0.){
1356  if(sPosDaugh->iCrossedTPC / sPosDaugh->iNoTPCCluster > fV0Cuts[NoCrossedOverNoTPCClustersMax]) continue;
1357  if(sNegDaugh->iCrossedTPC / sNegDaugh->iNoTPCCluster > fV0Cuts[NoCrossedOverNoTPCClustersMax]) continue;
1358  }
1359 
1360  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1361  iCutIndex++;
1362 
1363  // 4
1364  // Daughters: transverse momentum cut
1365  if(fV0Cuts[DaughMinPt] > 0.){
1366  if((sPosDaugh->fPt < fV0Cuts[DaughMinPt]) || (sNegDaugh->fPt < fV0Cuts[DaughMinPt])) continue;
1367  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1368  }
1369  iCutIndex++;
1370 
1371  // 5
1372  // Daughters: Impact parameter of daughters to prim vtx
1373  if(fV0Cuts[MinDCADaughWrtPV] > 0.){
1374  if((sPosDaugh->fDCAtoPV < fV0Cuts[MinDCADaughWrtPV] ) || (sNegDaugh->fDCAtoPV < fV0Cuts[MinDCADaughWrtPV] )) continue;
1375  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1376  }
1377  iCutIndex++;
1378 
1379  // 6
1380  // Daughters: DCA
1381  if(fV0Cuts[MaxDCADaughvsDaugh] > 0.){
1382  if(sV0->fDCAV0DaughvsDaugh > fV0Cuts[MaxDCADaughvsDaugh]) continue;
1383  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1384  }
1385  iCutIndex++;
1386 
1387  // 7
1388  // V0: Cosine of the pointing angle
1389  if(fV0Cuts[MinCosPAK0] > 0.){
1390  if(sV0->fPA < fV0Cuts[MinCosPAK0]){
1391  sV0->bIsCandidateK0s = kFALSE;
1392  }
1393  }
1394  if(fV0Cuts[MaxCosPALambda] > 0.){
1395  if(sV0->fPA < fV0Cuts[MaxCosPALambda]){
1396  sV0->bIsCandidateLambda = kFALSE;
1397  sV0->bIsCandidateALambda = kFALSE;
1398  }
1399  }
1400  if(fV0Cuts[MinCosPAK0] > 0. || fV0Cuts[MaxCosPALambda] > 0.){
1401  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1402  }
1403  iCutIndex++;
1404 
1405 
1406  // 8
1407  // V0: Fiducial volume
1408  /*if(fdCutRadiusDecayMin > 0. && fdCutRadiusDecayMax > 0.)
1409  {
1410  if((dRadiusDecay < fdCutRadiusDecayMin) || (dRadiusDecay > fdCutRadiusDecayMax))
1411  continue;
1412  FillV0Candidates(bIsCandidateK0s, bIsCandidateLambda, bIsCandidateALambda, iCutIndex);
1413  }*/
1414 
1415  if(fV0Cuts[MinDecayRadius]>0){
1416  if(sV0->fDecayRadius<fV0Cuts[MinDecayRadius]) {
1417  continue;
1418  }
1419  }
1420  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1421  iCutIndex++;
1422 
1423  // 9
1424  // Daughters: pseudorapidity cut
1425  if(fV0Cuts[DaughMaxEta] > 0.){
1426  if((TMath::Abs(sPosDaugh->fEta) > fV0Cuts[DaughMaxEta]) || (TMath::Abs(sNegDaugh->fEta) > fV0Cuts[DaughMaxEta]))
1427  continue;
1428  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1429  }
1430  iCutIndex++;
1431 
1432  // 10
1433  // V0: rapidity cut
1434  if(fV0Cuts[MaxV0Rap] > 0.){
1435  if(TMath::Abs(sV0->fRapK0) > fV0Cuts[MaxV0Rap])
1436  sV0->bIsCandidateK0s = kFALSE;
1437  if(TMath::Abs(sV0->fRapLambda) > fV0Cuts[MaxV0Rap]){
1438  sV0->bIsCandidateLambda = kFALSE;
1439  sV0->bIsCandidateALambda = kFALSE;
1440  }
1441  }
1442  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1443  iCutIndex++;
1444 
1445  // 11
1446  // Lifetime cut
1447 
1448  // Mean lifetime
1449  Double_t dCTauK0s = 2.6844; // [cm] c*tau of K0S
1450  Double_t dCTauLambda = 7.89; // [cm] c*tau of Lambda
1451  if(fV0Cuts[MaxLifeTime] > 0.){
1452  if(sV0->fLifetimeK0 > fV0Cuts[MaxLifeTime] * dCTauK0s)
1453  sV0->bIsCandidateK0s = kFALSE;
1454  }
1455  if(fV0Cuts[MaxLifeTime] > 0.){
1456  if(sV0->fLifetimeLambda > fV0Cuts[MaxLifeTime] * dCTauLambda)
1457  {
1458  sV0->bIsCandidateLambda = kFALSE;
1459  sV0->bIsCandidateALambda = kFALSE;
1460  }
1461  }
1462  if(fV0Cuts[MaxLifeTime] > 0.){
1463  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1464  }
1465  iCutIndex++;
1466 
1467  // 12
1468  // Daughter PID
1469  if(fV0Cuts[MaxSigmadEdxTPC] > 0.){
1470  if(sV0->fSigmaPosPion > fV0Cuts[MaxSigmadEdxTPC] || sV0->fSigmaNegPion > fV0Cuts[MaxSigmadEdxTPC]){ // pi+, pi-
1471  sV0->bIsCandidateK0s = kFALSE;
1472  }
1473  if(sV0->fSigmaPosProton > fV0Cuts[MaxSigmadEdxTPC] || sV0->fSigmaNegPion > fV0Cuts[MaxSigmadEdxTPC]){ // p+, pi-
1474  sV0->bIsCandidateLambda = kFALSE;
1475  }
1476  if(sV0->fSigmaNegProton > fV0Cuts[MaxSigmadEdxTPC] || sV0->fSigmaPosPion > fV0Cuts[MaxSigmadEdxTPC]){ // p-, pi+
1477  sV0->bIsCandidateALambda = kFALSE;
1478  }
1479  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1480  }
1481  iCutIndex++;
1482 
1483  // 13
1484  // Armenteros-Podolanski cut
1485  if(fV0Cuts[DoArmenteros]){
1486  if(sV0->fArmenterosPt < TMath::Abs(0.2 * sV0->fArmenterosAlpha)){
1487  sV0->bIsCandidateK0s = kFALSE;
1488  }
1489  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1490  }
1491  iCutIndex++;
1492 
1493  // 14
1494  // Invariant mass peak selection
1495  Double_t dMassPeakWindowK0s = fV0Cuts[InvarMassWindowK0]; //0.010; // LF p-p
1496  Double_t dMassPeakWindowLambda = fV0Cuts[InvarMassWindowLambda]; //0.005; // LF p-p
1497  Double_t dMassPDGK0s = TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
1498  Double_t dMassPDGLambda = TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
1499 
1500  if(TMath::Abs(sV0->fMassK0 - dMassPDGK0s) < dMassPeakWindowK0s)
1501  sV0->bIsInPeakK0s = kTRUE;
1502  if(TMath::Abs(sV0->fMassLambda - dMassPDGLambda) < dMassPeakWindowLambda)
1503  sV0->bIsInPeakLambda = kTRUE;
1504  if(TMath::Abs(sV0->fMassAntilambda - dMassPDGLambda) < dMassPeakWindowLambda)
1505  sV0->bIsInPeakALambda = kTRUE;
1506 
1507  if(fV0Cuts[DoMassWindow]){
1508  if(sV0->bIsInPeakK0s){
1509  sV0->bIsCandidateLambda = kFALSE;
1510  sV0->bIsCandidateALambda = kFALSE;
1511  }
1512  if(sV0->bIsInPeakLambda){
1513  sV0->bIsCandidateK0s = kFALSE;
1514  }
1515  if(sV0->bIsInPeakALambda){
1516  sV0->bIsCandidateK0s = kFALSE;
1517  }
1518  FillV0Candidates(sV0->bIsCandidateK0s, sV0->bIsCandidateLambda, sV0->bIsCandidateALambda, iCutIndex);
1519  }
1520  iCutIndex++;
1521 
1522  if(sV0->bIsCandidateK0s) {fh2dKshortMassVsPt->Fill(sV0->fPt, sV0->fMassK0,1); }
1523  if(sV0->bIsCandidateLambda) {fh2dLamdaMassVsPt->Fill(sV0->fPt, sV0->fMassLambda,1); }
1524  if(sV0->bIsCandidateALambda) {fh2dAnLamdaMassVsPt->Fill(sV0->fPt, sV0->fMassAntilambda,1); }
1525 
1526  //Is true MC V0?
1527  if(fIsPythia){
1528  int iVetoDec=0;
1529  iVetoDec=(int) GetV0MCVeto(fAODIn,v0,sV0->bIsCandidateK0s,sV0->bIsCandidateLambda,sV0->bIsCandidateALambda);
1530  Double_t fIsMCTrueK0=0.;
1531  Double_t fIsMCTrueLambda=0.;
1532  Double_t fIsMCTrueALambda=0.;
1533  if(iVetoDec==1) fIsMCTrueK0=1.;
1534  if(iVetoDec==2) fIsMCTrueLambda=1.;
1535  if(iVetoDec==3) fIsMCTrueALambda=1.;
1536 
1537  //printf("VetoDecision=%i, fIsMCTrueK0=%f, fIsMCTrueLambda=%f, fIsMCTrueALamda=%f\n",iVetoDec,fIsMCTrueK0,fIsMCTrueLambda,fIsMCTrueALambda);
1538 
1539  AliEmcalJet * jetrec = 0x0;
1540  double fJetPt=0;
1541  if(!jetconrec)printf("No jet container with reconstructed jets!\n");
1542 
1543  while ((jetrec = jetconrec->GetNextAcceptJet())){
1544  fJetPt= jetrec->Pt();
1545  //if(!(fJetContainerData->GetRhoParameter() == 0x0)){
1546  // fJetPt = fJetPt - fJetContainerData->GetRhoVal() * jetrec->Area();
1547  //}
1548  if(fJetPt < 5.) continue;
1549  bool bIsInCone=IsParticleInCone(v0, jetrec, fJetRadius);
1550  //printf("VetoDecision=%i, fIsMCTrueK0=%f, fIsMCTrueLambda=%f, fIsMCTrueALamda=%f, bIsINCone=%i\n",iVetoDec,fIsMCTrueK0,fIsMCTrueLambda,fIsMCTrueALambda,bIsInCone);
1551 
1552  // make inclusive signed imp. parameter constituent histograms
1553  if(sV0->bIsCandidateK0s &&bIsInCone) {
1554  Double_t valueKInJC[5] = {sV0->fMassK0, sV0->fPt, sV0->fEta, fJetPt, fIsMCTrueK0};
1555  fhnV0InJetK0s->Fill(valueKInJC);
1556  //printf("massk0=%f, pt=%f, eta=%f, jetpt=%f, istrue=%f\n",sV0->fMassK0, sV0->fPt, sV0->fEta,fJetPt,fIsMCTrueK0);
1557  }
1558  if(sV0->bIsCandidateLambda && bIsInCone) {
1559  Double_t valueLInJC[5] = {sV0->fMassLambda, sV0->fPt, sV0->fEta, fJetPt,fIsMCTrueLambda};
1560  fhnV0InJetLambda->Fill(valueLInJC);
1561  //printf("masslambda=%f, pt=%f, eta=%f, jetpt=%f, istrue=%f\n",sV0->fMassK0, sV0->fPt, sV0->fEta,fJetPt,fIsMCTrueLambda);
1562  }
1563  if(sV0->bIsCandidateALambda && bIsInCone) {
1564  Double_t valueLInJC[5] = {sV0->fMassAntilambda, sV0->fPt, sV0->fEta, fJetPt,fIsMCTrueALambda};
1565  fhnV0InJetALambda->Fill(valueLInJC);
1566  //printf("massalambda=%f, pt=%f, eta=%f, jetpt=%f, istrue=%f\n",sV0->fMassK0, sV0->fPt, sV0->fEta,fJetPt,fIsMCTrueALambda);
1567  }
1568 
1569  }
1570  jetconrec->ResetCurrentID();
1571  jetrec=NULL;
1572  //delete jetrec;
1573 
1574  }
1575  //Store V0 candidate for later rejection of daughter tracks
1576  if(sV0->bIsCandidateK0s || sV0->bIsCandidateLambda || sV0->bIsCandidateALambda){
1577  int nBins=fV0CandidateArray->GetEntriesFast();
1578  //v0->Print();
1579  new((*fV0CandidateArray)[nBins]) AliAODv0(*v0);
1580  }
1581  }
1582  delete sV0;
1583  delete sPosDaugh;
1584  delete sNegDaugh;
1585  //printf("################## Select candidates End########################\n");
1586 
1587 }
1588 
1589 Bool_t AliAnalysisTaskHFJetIPQA::Run(){
1590 
1591  //*******************************
1592  //Selection
1593  FillGeneralHistograms();
1594  /*Vertex Pos Selection*/
1595 
1596  AliAODEvent *ev = dynamic_cast<AliAODEvent*>(InputEvent());
1597  fEventVertex = dynamic_cast<AliAODVertex*>(ev->GetPrimaryVertex());
1598  fIsSameEvent_n1 = kFALSE;
1599  fIsSameEvent_n2 = kFALSE;
1600  fIsSameEvent_n3 = kFALSE;
1601 
1602  IsEventAccepted(ev);
1603 
1604 
1605  Bool_t HasImpactParameter = kFALSE;
1606  Double_t dca[2] = {-99999,-99999};
1607  Double_t cov[3] = {-99999,-99999,-99999};
1608  Double_t TrackWeight = 1;
1609  AliVTrack* trackV = NULL;
1610  fIsEsd = (InputEvent()->IsA()==AliESDEvent::Class())? kTRUE : kFALSE;
1611  // EventwiseCleanup();
1612  if(fIsPythia){
1613  if(fIsEsd){
1614  fMCEvent = dynamic_cast<AliMCEvent*>(MCEvent()) ;
1615  if (!fMCEvent){
1616  AliError("Could not retrieve MC particles! Returning");
1617  return kFALSE;
1618  }
1619  }
1620  else{
1621  fMCArray = static_cast<TClonesArray*>(InputEvent()->FindListObject(AliAODMCParticle::StdBranchName()));
1622  if (!fMCArray){
1623  AliError("Could not retrieve AOD MC particles! Returning");
1624  return kFALSE;
1625  }
1626  }
1627  jetcongen = static_cast<AliJetContainer*>(fJetCollArray.At(1));
1628  jetcongen->ResetCurrentID();
1629  }
1630 
1631  jetconrec = static_cast<AliJetContainer*>(fJetCollArray.At(0));
1632  jetconrec->ResetCurrentID();
1633 
1634 
1635  if(fApplyV0Rej){
1636  SelectV0Candidates(ev);
1637  if(fIsPythia) GetV0MCTrueCandidates(ev);
1638  }
1639  IncHist("fh1dEventsAcceptedInRun",1);
1640 
1641  //if(fIsPythia)FillParticleCompositionEvent(); //Added for in cone vs inclusive comparison
1642 
1643  FillHist("fh1dNoParticlesPerEvent",InputEvent()->GetNumberOfTracks(),1);
1644 
1645  //************************
1646  //Investigate Track Corrections: Smearing & particle weighting
1647  for(long itrack= 0; itrack<InputEvent()->GetNumberOfTracks();++itrack)
1648  {
1649  trackV = static_cast<AliVTrack*>(InputEvent()->GetTrack(itrack));
1650  if(!trackV) {
1651  AliInfo("Could not retrieve Track");
1652  continue;
1653  }
1654  IncHist("fh1dTracksAccepeted",1);
1655  if(!IsTrackAccepted(trackV,-1)) {
1656  IncHist("fh1dTracksAccepeted",3);
1657  continue;
1658  }
1659  if(fRunSmearing)SmearTrack((AliAODTrack*)trackV);
1660  IncHist("fh1dTracksAccepeted",2);
1661  FillHist("fh2dAcceptedTracksEtaPhi",trackV->Eta(),trackV->Phi(),1); //this->fXsectionWeightingFactor );
1662  TrackWeight =1;dca[0]=-9999;dca[1]=-9999;cov[0]=-9999;cov[1]=-9999;cov[2]=-9999;
1663  HasImpactParameter =kFALSE;
1664  Double_t xyzatdca[3];
1665  if (GetImpactParameter(((AliAODTrack*)trackV),(AliAODEvent *)InputEvent(), dca, cov,xyzatdca))HasImpactParameter =kTRUE;
1666  if(fEventVertex) {
1667  delete fEventVertex;
1668  fEventVertex =nullptr;
1669  }
1670  if(!HasImpactParameter) continue;
1671  /*if(fIsPythia){
1672  Int_t corrpartidx =-1;
1673  double ppt;
1674  //if(fDoMCCorrection) TrackWeight *= GetMonteCarloCorrectionFactor(trackV,corrpartidx,ppt);
1675  }*/
1676  FillTrackHistograms(trackV,dca,cov,TrackWeight);
1677  }
1678  //**********************************
1679  //JetMatching between generated and reconstructed Jets
1680 
1681  //printf("In Program %f < jetpt <%f, %f < jeteta < %f\n",fAnalysisCuts[bAnalysisCut_MinJetPt],fAnalysisCuts[bAnalysisCut_MaxJetPt],fAnalysisCuts[bAnalysisCut_MinJetEta],fAnalysisCuts[bAnalysisCut_MaxJetEta] );
1682  FillHist("fh1dNoJetsPerEvent",jetconrec->GetNJets(),1);
1683  if (!jetconrec) return kFALSE;
1684  AliEmcalJet * jetgen = nullptr;
1685 
1686  if(fIsPythia){
1687  if(!MatchJetsGeometricDefault()) AliInfo("Jet matching did not succeed!");
1688  jetcongen->ResetCurrentID();
1689  while ((jetgen = jetcongen->GetNextAcceptJet()))
1690  {
1691  if (!jetgen) continue;
1692  //Int_t jetflavour =0;
1693  //Bool_t is_udgjet = kFALSE;
1694  //jetflavour =IsMCJetPartonFast(jetgen,fJetRadius,is_udgjet);
1695  //FillGenHistograms(jetflavour, jetgen);
1696  }
1697  jetcongen->ResetCurrentID();
1698  jetconrec->ResetCurrentID();
1699  }
1700 
1701  //*************************************************
1702  // Loop over reconstructed/matched jets for template creation and analysis
1703  AliEmcalJet * jetrec = nullptr;
1704  AliEmcalJet * jetmatched = nullptr;
1705  jetconrec->ResetCurrentID();
1706  Double_t jetpt=0;
1707  while ((jetrec = jetconrec->GetNextAcceptJet()))
1708  {//start jetloop
1709  if(!jetrec) continue;
1710  jetpt = jetrec->Pt();
1711  if(fDoUnderlyingEventSub)jetpt=DoUESubtraction(jetcongen, jetconrec,jetrec, jetpt);
1712 
1713  FillHist("fh1dJetArea",jetrec->Area(),1);
1714 
1715  //________________________
1716  //Determination of Jet Flavour
1717  Int_t jetflavour=0;
1718  Bool_t is_udgjet = kFALSE;
1719  if(fIsPythia){
1720  jetmatched = nullptr;
1721  jetmatched =jetrec->MatchedJet();
1722  if(jetmatched){
1723  jetflavour = IsMCJetPartonFast(jetmatched,fJetRadius,is_udgjet); //Event based association to save memory
1724  }
1725  else{
1726  jetflavour=Unid;
1727  }
1728  }
1729  FillRecHistograms( jetflavour, jetpt, jetrec->Eta(),jetrec->Phi());
1730  if(fDoLundPlane)RecursiveParents(jetrec, jetconrec);
1731 
1732  //_____________________________
1733  //Determination of impact parameters
1734  std::vector<SJetIpPati> sImpParXY,sImpParXYZ,sImpParXYSig,sImpParXYZSig;
1735  AliVParticle *vp=0x0;
1736  int NJetParticles=0; //Used for counting particles per jet
1737  int isV0=kFALSE;
1738  Double_t dca[2] = {-99999,-99999};
1739  Double_t cov[3] = {-99999,-99999,-99999};
1740  Double_t sign=0;
1741 
1742  for(UInt_t i = 0; i < jetrec->GetNumberOfTracks(); i++) {//start trackloop
1743  TrackWeight=1;
1744  isV0=kFALSE;
1745  Double_t xyzatcda[3];
1746 
1747  vp = static_cast<AliVParticle*>(jetrec->TrackAt(i, jetconrec->GetParticleContainer()->GetArray()));
1748  if (!vp){
1749  AliError("AliVParticle associated to constituent not found");
1750  continue;
1751  }
1752  AliVTrack *vtrack = dynamic_cast<AliVTrack*>(vp);
1753  if (!vtrack) {
1754  AliError(Form("Could not receive track%d\n", i));
1755  continue;
1756  }
1757  AliAODTrack *trackV = dynamic_cast<AliAODTrack*>(vtrack);
1758 
1759  if (!trackV || !jetrec) continue;
1760  if (fIsPythia&&!IsTrackAccepted((AliAODTrack*)trackV,jetflavour)) continue;
1761 
1762  isV0=IsV0Daughter(trackV);
1763  ++NJetParticles;
1764 
1765  //FillTrackTypeResHists();
1766 
1767  if(GetImpactParameterWrtToJet((AliAODTrack*)trackV,(AliAODEvent*)InputEvent(),jetrec,dca,cov,xyzatcda,sign, jetflavour)){
1768  if(fEventVertex) {
1769  delete fEventVertex;
1770  fEventVertex =nullptr;
1771  }
1772  Int_t corridx=-1;double ppt;
1773  //(fIsPythia&&fDoMCCorrection) ? TrackWeight = GetMonteCarloCorrectionFactor(trackV,corridx,ppt) : TrackWeight =1;
1774  dca[0]=fabs(dca[0]);
1775  Double_t cursImParXY =TMath::Abs(GetValImpactParameter( kXY,dca,cov))*sign;
1776  Double_t cursImParXYSig =TMath::Abs(GetValImpactParameter(kXYSig,dca,cov))*sign;
1777  Double_t cursImParXYZ =TMath::Abs(GetValImpactParameter( kXYZ,dca,cov))*sign;
1778  Double_t cursImParXYZSig =TMath::Abs(GetValImpactParameter(kXYZSig,dca,cov))*sign;
1779  FillHist("fh2dJetSignedImpParXY" ,jetpt,cursImParXY,TrackWeight); //*this->fXsectionWeightingFactor );
1780  FillHist("fh2dJetSignedImpParXYSignificance",jetpt,cursImParXYSig,TrackWeight); //*this->fXsectionWeightingFactor );
1781 
1782  const char * subtype [5] = {"Unidentified","udsg","c","b","s"};
1783  if(fIsPythia){
1784  FillHist(Form("fh2dJetSignedImpParXY%s",subtype[jetflavour]),jetpt,cursImParXY,TrackWeight); //*this->fXsectionWeightingFactor );
1785  FillHist(Form("fh2dJetSignedImpParXYSignificance%s",subtype[jetflavour]),jetpt,cursImParXYSig,TrackWeight); //*this->fXsectionWeightingFactor );
1786  }
1787  double fTrackPt=trackV->Pt();
1788  double fIPValue=fV0Cuts[fAV0Cut]*TMath::Exp(fV0Cuts[fBV0Cut]*fTrackPt)+fV0Cuts[fCV0Cut];
1789  //printf("trackpt=%f, IPValue=%f, TrueIP=%f, a=%f, b=%f, c=%f\n", fTrackPt, fIPValue,cursImParXYSig, fV0Cuts[fAV0Cut], fV0Cuts[fBV0Cut], fV0Cuts[fCV0Cut]);
1790  if(cursImParXYSig>fIPValue){
1791  //printf("Going into switch!\n");
1792  switch (isV0){
1793  case V0No:
1794  isV0=V0Rec;
1795  //printf("V0No set to V0Rec!\n");
1796  break;
1797  case V0MC:
1798  isV0=V0TrueRec;
1799  //printf("V0MC set to V0TrueRec!\n");
1800  break;
1801  }
1802  }
1803 
1804  SJetIpPati a(cursImParXY, TrackWeight,isV0,kFALSE,corridx,trackV->Pt()); sImpParXY.push_back(a);
1805  SJetIpPati b(cursImParXYZ, TrackWeight,isV0,kFALSE,corridx,trackV->Pt()); sImpParXYZ.push_back(b);
1806  SJetIpPati c(cursImParXYSig, TrackWeight,isV0,kFALSE,corridx,trackV->Pt());sImpParXYSig.push_back(c);
1807  SJetIpPati d(cursImParXYZSig, TrackWeight,isV0,kFALSE,corridx,trackV->Pt());sImpParXYZSig.push_back(d);
1808  //printf("curImParXY=%f, isV0=%i, pt=%f\n",sImpParXYSig.back().first,sImpParXYSig.back().is_V0, sImpParXYSig.back().trackpt);
1809 
1810  }
1811  }//end trackloop
1812 
1813  FillHist("fh1dParticlesPerJet",NJetParticles,1);
1814  //_________________________
1815  //Sorting of Impact Parameters
1816  bool hasIPs[4] ={kFALSE,kFALSE,kFALSE, kFALSE};
1817  Double_t ipval[4] = {-9999.,-9999.,-9999., -9999.};
1818  Double_t ipvalsig[4] = {-9999.,-9999.,-9999., -9999.};
1819 
1820  std::sort(sImpParXY.begin(),sImpParXY.end(), AliAnalysisTaskHFJetIPQA::mysort);
1821  std::sort(sImpParXYSig.begin(),sImpParXYSig.end(), AliAnalysisTaskHFJetIPQA::mysort);
1822  std::sort(sImpParXYZ.begin(),sImpParXYZ.end(), AliAnalysisTaskHFJetIPQA::mysort);
1823  std::sort(sImpParXYZSig.begin(),sImpParXYZSig.end(),AliAnalysisTaskHFJetIPQA::mysort);
1824 
1825  if((int)sImpParXYSig.size()>0) hasIPs[0]=kTRUE;
1826  if((int)sImpParXYSig.size()>1) hasIPs[1]=kTRUE;
1827  if((int)sImpParXYSig.size()>2) hasIPs[2]=kTRUE;
1828  if((int)sImpParXYSig.size()>3) hasIPs[3]=kTRUE;
1829 
1830  if(hasIPs[0]){
1831  ipvalsig[0] =sImpParXYSig.at(0).first;
1832  ipval[0] =sImpParXY.at(0).first;
1833  //printf("HasIP0, ipval[0]=%f, ipvalsig[0]=%f\n", ipval[0], ipvalsig[0]);
1834  }
1835  if(hasIPs[1]){
1836  ipvalsig[1] =sImpParXYSig.at(1).first;
1837  ipval[1] =sImpParXY.at(1).first;
1838  //printf("HasIP1, ipval[1]=%f, ipvalsig[1]=%f\n",ipval[1], ipvalsig[1]);
1839  }
1840  if(hasIPs[2]){
1841  ipvalsig[2] =sImpParXYSig.at(2).first;
1842  ipval[2] =sImpParXY.at(2).first;
1843  //printf("HasIP2, ipval[2]=%f, ipvalsig[2]=%f\n", ipval[2], ipvalsig[2]);
1844  }
1845  if(hasIPs[3]){
1846  ipvalsig[3] =sImpParXYSig.at(3).first;
1847  ipval[3] =sImpParXY.at(3).first;
1848  //printf("HasIP3, ipval[3]=%f, ipvalsig[3]=%f\n", ipval[3], ipvalsig[3]);
1849  }
1850 
1851  //if(hasIPs[0])printf("N=1: cursImParXY=%f, TrackWeight=%f,corridx=%i, pt=%f\n",sImpParXYSig.at(0).first, sImpParXYSig.at(0).second, sImpParXYSig.at(0).trackLabel, sImpParXYSig.at(0).trackpt);
1852  //if(hasIPs[1])printf("N=2: cursImParXY=%f, TrackWeight=%f, corridx=%i, pt=%f\n",sImpParXYSig.at(1).first, sImpParXYSig.at(1).second, sImpParXYSig.at(1).trackLabel, sImpParXYSig.at(1).trackpt);
1853  //if(hasIPs[2])printf("N=3: cursImParXY=%f, TrackWeight=%f, corridx=%i, pt=%f\n",sImpParXYSig.at(2).first, sImpParXYSig.at(2).second, sImpParXYSig.at(2).trackLabel, sImpParXYSig.at(2).trackpt);
1854  //printf("*********************************************************\n");
1855 
1856  //_____________________________
1857  //V0 tag decisions
1858  //printf("isV0=%i, jetflaovur=%i, jetpt=%f", );
1859  bool isV0Jet=kFALSE;
1860  if((hasIPs[0])&&(!fIsPythia)&&(sImpParXYSig[0].is_V0==V0Rec))isV0Jet=kTRUE;
1861  //printf("New jetflavour=%i, isV0Jet=%i\n",jetflavour, isV0Jet);
1862 
1863  if(fIsPythia){
1864  if(hasIPs[0])FillV0EfficiencyHists(sImpParXYSig[0].is_V0, jetflavour, jetpt, isV0Jet);
1865  for(long unsigned iTrack=0;iTrack<sImpParXYSig.size();iTrack++){
1866  FillTrackIPvsPt(sImpParXYSig[iTrack].is_V0,sImpParXYSig[iTrack].trackpt,sImpParXYSig[iTrack].first,jetflavour);
1867  }
1868 
1869  //_______________________________
1870  //IP Template Generation
1871  FillIPTypePtHists(jetflavour, jetpt, hasIPs);
1872  }
1873  for(int iN=0;iN<3;iN++){
1874  if(!hasIPs[iN]) continue;
1875  //printf("iN=%i, jetflavour=%i xy=%f, xysig=%f\n",iN,jetflavour,sImpParXY.at(iN).first,sImpParXYSig.at(iN).first);
1876  Double_t params [4] ={sImpParXY.at(iN).first,sImpParXYSig.at(iN).first,sImpParXYZ.at(iN).first,sImpParXYZSig.at(iN).first};
1877  FillIPTemplateHists(jetpt,iN,jetflavour, params);
1878  }
1879 
1880  //____________________________________________
1881  //TAGGING
1882  if(fDoTCTagging!=TCNo&&fDoProbTagging!=ProbNo) AliFatal("Don't do track counting and probability tagging simultaneously!");
1883 
1884  bool ** kTagDec=new bool*[fNThresholds];
1885  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
1886  kTagDec[iThresh]=new bool[6];
1887  for(int iType=0;iType<6;iType++){
1888  kTagDec[iThresh][iType]=0;
1889  }
1890  }
1891  //**************
1892  //MC Track Counting
1893  if(fDoTCTagging!=TCNo){
1894  if(fIsPythia||((!fIsPythia)&&(!isV0Jet))){
1895  //printf("isV0Jet=%i\n", isV0Jet);
1896  if(fUseSignificance){DoTCTagging(jetpt, hasIPs,ipvalsig, kTagDec);}
1897  else{DoTCTagging(jetpt, hasIPs,ipval, kTagDec);}
1898  }
1899  if(fDoMCEffs){
1900  //printf("Filling Efficiency hists\n");
1901  FillEfficiencyHists(kTagDec, jetflavour, jetpt,hasIPs[0]);
1902  }
1903  FillTaggedJetPtDistribution(kTagDec,jetpt);
1904  }
1905  //**************
1906  //Probability Dists
1907  double probval=0;
1908  probval=GetTrackProbability(jetpt,hasIPs, ipvalsig);
1909  //Generation of Track Probability Hists
1910  if(fDoJetProb&&(probval>0)){
1911  //printf("Doing Jet Probability!\n");
1912  FillProbabilityHists(jetpt, probval, jetflavour, kTagDec);
1913  FillProbThreshHists(probval, ipvalsig, jetpt, jetflavour, hasIPs, kTagDec);
1914  }
1915  //Probability Tagging
1916  /*if(fDoProbTagging!=ProbNo){
1917  DoProbTagging(probval, jetpt,kTagDec);
1918  FillEfficiencyHists(kTagDec, jetflavour, jetpt,hasIPs[0]);
1919  }*/
1920 
1921  if(sImpParXY.size()!=0){
1922  FillHist("fh2dNoAcceptedTracksvsJetArea",(int)sImpParXY.size(),jetrec->Area(),1);
1923  }
1924  sImpParXY.clear();
1925  sImpParXYSig.clear();
1926  sImpParXYZ.clear();
1927  sImpParXYZSig.clear();
1928 
1929  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
1930  delete kTagDec[iThresh];
1931  }
1932  delete [] kTagDec;
1933  }//end jetloop*/
1934  return kTRUE;
1935  }
1936 
1937 
1938 
1939  Double_t AliAnalysisTaskHFJetIPQA::GetLocalAlphaAOD(AliAODTrack * track)
1940  {
1941  AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1942  return etp.GetAlpha();
1943  }
1944  Double_t AliAnalysisTaskHFJetIPQA::GetLocalThetaAOD(AliAODTrack * track)
1945  {
1946  // convert to AliExternalTrackParam
1947  AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1948  // propagate
1949 
1950  Double_t dv[2],dcov[3];
1951  const Double_t kBeampiperadius=3;
1952  AliVEvent * eev = (AliVEvent*)InputEvent();
1953  const AliVVertex *vtxESDSkip =(const AliVVertex *) (InputEvent()->GetPrimaryVertex()) ;
1954  if(!vtxESDSkip) return -9999;
1955  if(!(etp.PropagateToDCA(vtxESDSkip, eev->GetMagneticField(), kBeampiperadius, dv, dcov)))return -9999.;
1956  // update track position and momentum
1957  return etp.Theta();
1958  }
1959 
1960  Double_t AliAnalysisTaskHFJetIPQA::GetTrackCurvature(AliAODTrack * track)
1961  {
1962  // convert to AliExternalTrackParam
1963  AliVEvent * eev = (AliVEvent*)InputEvent();
1964 
1965  AliExternalTrackParam etp; etp.CopyFromVTrack(track);
1966 
1967  return etp.GetC(eev->GetMagneticField());
1968  }
1969 
1970 
1985 void AliAnalysisTaskHFJetIPQA::SetUseMonteCarloWeighingLinus(TH1F *Pi0, TH1F *Eta, TH1F *EtaP, TH1F *Rho, TH1F *Phi, TH1F *Omega, TH1F *K0s, TH1F *Lambda, TH1F *ChargedPi, TH1F *ChargedKaon, TH1F *Proton, TH1F *D0, TH1F *DPlus, TH1F *DStarPlus, TH1F *DSPlus, TH1F *LambdaC, TH1F *BPlus, TH1F *B0, TH1F *LambdaB, TH1F *BStarPlus)
1986  {
1987  for(Int_t i =1 ; i< Pi0->GetNbinsX()+1;++i){
1988  fBackgroundFactorLinus[bIdxPi0][i-1] =Pi0->GetBinContent(i);
1989  fBackgroundFactorLinus[bIdxEta][i-1] =Eta->GetBinContent(i);
1990  fBackgroundFactorLinus[bIdxEtaPrime][i-1] =EtaP->GetBinContent(i);
1991  fBackgroundFactorLinus[bIdxRho][i-1] =Rho->GetBinContent(i);
1992  fBackgroundFactorLinus[bIdxPhi][i-1] =Phi->GetBinContent(i);
1993  fBackgroundFactorLinus[bIdxOmega][i-1] =Omega->GetBinContent(i);
1994  fBackgroundFactorLinus[bIdxK0s][i-1] =K0s->GetBinContent(i);
1995  fBackgroundFactorLinus[bIdxLambda][i-1] =Lambda->GetBinContent(i);
1996  fBackgroundFactorLinus[bIdxPi][i-1] =ChargedPi->GetBinContent(i);
1997  fBackgroundFactorLinus[bIdxKaon][i-1] =ChargedKaon->GetBinContent(i);
1998  fBackgroundFactorLinus[bIdxProton][i-1] =Proton->GetBinContent(i);
1999  fBackgroundFactorLinus[bIdxD0][i-1] =D0->GetBinContent(i);
2000  fBackgroundFactorLinus[bIdxDPlus][i-1] =DPlus->GetBinContent(i);
2001  fBackgroundFactorLinus[bIdxDStarPlus][i-1] =DStarPlus->GetBinContent(i);
2002  fBackgroundFactorLinus[bIdxDSPlus][i-1] =DSPlus->GetBinContent(i);
2003  fBackgroundFactorLinus[bIdxLambdaC][i-1] =LambdaC->GetBinContent(i);
2004  fBackgroundFactorLinus[bIdxBPlus][i-1] =BPlus->GetBinContent(i);
2005  fBackgroundFactorLinus[bIdxB0][i-1] =B0->GetBinContent(i);
2006  fBackgroundFactorLinus[bIdxLambdaB][i-1] =LambdaB->GetBinContent(i);
2007  fBackgroundFactorLinus[bIdxBStarPlus][i-1] =BStarPlus->GetBinContent(i);
2008  }
2009  return;
2010 }
2011 
2012 void AliAnalysisTaskHFJetIPQA::SetFlukaFactor(TGraph* GraphOmega, TGraph* GraphXi, TGraph* K0Star, TGraph* Phi)
2013  {
2014  fGraphOmega=(TGraph*)GraphOmega;
2015  fGraphXi=(TGraph*)GraphXi;
2016  fK0Star=(TGraph*)K0Star;
2017  fPhi=(TGraph*)Phi;
2018 
2019  return;
2020  }
2021 
2022 
2023 
2030  Bool_t AliAnalysisTaskHFJetIPQA::SetResFunctionPID(const char * filename){
2031  TFile * jetProbfile=TFile::Open(filename,"READ");
2032 
2033  int bins_low[5] = {0,1,2,4,6};
2034  int bins_high[5] = {1,2,4,6,255};
2035  int its_hits[4] = {6,5,4,3};
2036 
2037  const char * type[5] ={"Electron","Pion","Kaon","Proton",""};
2038  for (int k=0;k<5;++k){
2039  for (int i=0;i<4;++i){
2040  for (int j=0;j<5;++j){
2041  TGraph *fResulFkt = 0x0;
2042  if(k==4) fResulFkt = (TGraph*)jetProbfile->Get(Form("fResulFkt_ITS_%i_PT_%i_to_%i",its_hits[i],bins_low[j],bins_high[j]));
2043  else fResulFkt = (TGraph*)jetProbfile->Get(Form("fResulFkt_ITS_%i_PT_%i_to_%i_%s",its_hits[i],bins_low[j],bins_high[j],type[k]));
2044  if(!fResulFkt){
2045  return kFALSE;
2046  }
2047  else {
2048  fResolutionFunction [20*k + 4*j +i] = *fResulFkt;
2049  Printf("Added %i %i %i -> [%i][%i][%i]->[%i]",j,k,i,j,i,k,20*k + 4*j +i );
2050  delete fResulFkt;
2051  }
2052  }
2053  }
2054  }
2055  if(jetProbfile)jetProbfile->Close();
2056 
2057  fUsePIDJetProb =kTRUE;
2058  return kTRUE;
2059 
2060  }
2061 
2062  void AliAnalysisTaskHFJetIPQA::FillCorrelations(bool bn[3],double v[3], double jetpt ){
2063  //Fill all possible same jet distributions
2064  fTREE_n1 = -999;
2065  fTREE_n2 = -999;
2066  fTREE_n3 = -999;
2067  fTREE_pt = -999;
2068 
2069  if(fUseTreeForCorrelations){
2070  if(bn[0] && bn[1]){
2071  fTREE_n1 = v[0];
2072  fTREE_n2 = v[1];
2073  fTREE_n3 = v[2];
2074  fTREE_pt = jetpt;
2075  fCorrelationCrossCheck->Fill();
2076  }
2077  }
2078  else if (fFillCorrelations && !fUseTreeForCorrelations )
2079  {
2080  if (bn[0] && bn[1]) {
2081  FillHist("fh2dInclusiveCorrelationN1N2",v[0],v[1]);
2082  if(jetpt>10 && jetpt <20) FillHist("fh2dGreater10_20GeVCorrelationN1N2",v[0],v[1]);
2083  if(jetpt>20 && jetpt <30) FillHist("fh2dGreater20_30GeVCorrelationN1N2",v[0],v[1]);
2084  if(jetpt>30 && jetpt <100) FillHist("fh2dGreater30_100GeVCorrelationN1N2",v[0],v[1]);
2085 
2086  }
2087  if (bn[0] && bn[2]) {
2088  FillHist("fh2dInclusiveCorrelationN1N3",v[1],v[2]);
2089  if(jetpt>10 && jetpt <20) FillHist("fh2dGreater10_20GeVCorrelationN1N3",v[0],v[2]);
2090  if(jetpt>20 && jetpt <30) FillHist("fh2dGreater20_30GeVCorrelationN1N3",v[0],v[2]);
2091  if(jetpt>30 && jetpt <100) FillHist("fh2dGreater30_100GeVCorrelationN1N3",v[0],v[2]);
2092  FillHist("fh2dInclusiveCorrelationN2N3",v[0],v[2]);
2093  if(jetpt>10 && jetpt <20) FillHist("fh2dGreater10_20GeVCorrelationN2N3",v[1],v[2]);
2094  if(jetpt>20 && jetpt <30) FillHist("fh2dGreater20_30GeVCorrelationN2N3",v[1],v[2]);
2095  if(jetpt>30 && jetpt <100) FillHist("fh2dGreater30_100GeVCorrelationN2N3",v[1],v[2]);
2096  }
2097  //Fill if possible mix distributions
2098  bool n3wasReady = false;
2099  double storedn3=-999;
2100  if(bn[0]){
2101 
2102  if(fIsMixSignalReady_n2) {
2103  double n2 =0;
2104  if(GetMixDCA(2,n2)) {
2105  FillHist("fh2dInclusiveCorrelationN1N2mix",v[0],n2);
2106  if(jetpt>10 && jetpt <20) FillHist("fh2dGreater10_20GeVCorrelationN1N2mix",v[0],n2);
2107  if(jetpt>20 && jetpt <30) FillHist("fh2dGreater20_30GeVCorrelationN1N2mix",v[0],n2);
2108  if(jetpt>30 && jetpt <100) FillHist("fh2dGreater30_100GeVCorrelationN1N2mix",v[0],n2);
2109  }
2110  }
2111  if(fIsMixSignalReady_n3) {
2112  double n3 =0;
2113  if(GetMixDCA(3,n3)) {
2114  n3wasReady=true;
2115  storedn3=n3;
2116  FillHist("fh2dInclusiveCorrelationN1N3mix",v[0],n3);
2117  if(jetpt>10 && jetpt <20) FillHist("fh2dGreater10_20GeVCorrelationN1N3mix",v[0],n3);
2118  if(jetpt>20 && jetpt <30) FillHist("fh2dGreater20_30GeVCorrelationN1N3mix",v[0],n3);
2119  if(jetpt>30 && jetpt <100) FillHist("fh2dGreater30_100GeVCorrelationN1N3mix",v[0],n3);
2120  }
2121  }
2122  }
2123 
2124  if(bn[1]) {
2125  if(n3wasReady) {
2126  double n3 =0;
2127  n3=storedn3;
2128  FillHist("fh2dInclusiveCorrelationN2N3mix",v[1],n3);
2129  if(jetpt>10 && jetpt <20) FillHist("fh2dGreater10_20GeVCorrelationN2N3mix",v[1],n3);
2130  if(jetpt>20 && jetpt <30) FillHist("fh2dGreater20_30GeVCorrelationN2N3mix",v[1],n3);
2131  if(jetpt>30 && jetpt <100) FillHist("fh2dGreater30_100GeVCorrelationN2N3mix",v[1],n3);
2132  }
2133  }
2134  }
2135  return;
2136  }
2137 
2138 
2139 
2140 void AliAnalysisTaskHFJetIPQA::UserCreateOutputObjects(){
2141  Printf("Analysing Jets with Radius: R=%f\n",fJetRadius);
2142 
2143  TString BJetCuts[25] = {
2144  "#sigma_{Dia}", //0
2145  "#sigma_{z}", //1
2146  "#sigma_{y}", //2
2147  "RelError_{z}", //3
2148  "RelError_{y}", //4
2149  "N_{Cont}", //5
2150  "MaxVtxZ", //6
2151  "d_{JetTrack}", //7
2152  "d_{z}", //8
2153  "d_{xy}", //9
2154  "DecayLength", //10
2155  "d_{chi2/ndf}", //11
2156  "p_{T,Track}^{min}", //12
2157  "p_{T,TrackMC}^{min}", //13
2158  "MinTPCClus", //14
2159  "n_{ITS Hits}", //15
2160  "chi2_{track}", //16
2161  "p_{T,Jet}^{min}",//17
2162  "p_{T,Jet}^{max}",//18
2163  "#eta_{Jet}^{min}", //19
2164  "#eta_{Jet}^{max}", //20
2165  "SPD Hits", //21
2166  "Kink",//22
2167  "TPC Refit",//23
2168  "ITS Refit" //24
2169  };
2170 
2171  TString sV0Cuts[18] = {
2172  "all"/*0*/,
2173  "reconstr."/*1*/,
2174  "tracks TPC"/*2*/,
2175  "Daugh pt"/*3*/,
2176  "DCA(Daugh<->PV)"/*4*/,
2177  "DCA(Daugh1<->Daugh2)"/*5*/,
2178  "cos(PA)"/*6*/,
2179  "R Decay"/*7*/,
2180  "Daugh Eta"/*8*/,
2181  "V0 rap, "/*9*/,
2182  "lifetime"/*10*/,
2183  "sigma dEdx"/*11*/,
2184  "pt Arm"/*12*/,
2185  "mass wind"/*13*/,
2186  "in jet"/*14*/,
2187  "a"/*15*/,
2188  "b"/*16*/,
2189  "c"/*17*/,
2190  };
2191 
2192  TString sTemp[7]={"Unidentified","udsg","c","b","udsgV0","cV0",""};
2193  for(int iS=0;iS<7;iS++){
2194  sTemplateFlavour.push_back(sTemp[iS]);
2195  }
2196 
2197  fIsMixSignalReady_n1=kFALSE;
2198  fIsMixSignalReady_n2=kFALSE;
2199  fIsMixSignalReady_n3=kFALSE;
2200  fn1_mix =-1;
2201  fn2_mix =-1;
2202  fn3_mix =-1;
2203  fIsSameEvent_n1=kFALSE;
2204  fIsSameEvent_n2=kFALSE;
2205  fIsSameEvent_n3=kFALSE;
2206 
2207  AliAnalysisTaskEmcal::UserCreateOutputObjects();
2208  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
2209  if(!mgr) AliError("Analysis manager not found!");
2210  AliVEventHandler *evhand = mgr->GetInputEventHandler();
2211  if (!evhand) AliError("Event handler not found!");
2212  if (evhand->InheritsFrom("AliESDInputHandler")) fIsEsd = kTRUE;
2213  else fIsEsd = kFALSE;
2214  OpenFile(1);
2215 
2216  Double_t lowIPxy =-1.; //ranges of xy axis of fh2dTracksImpParXY and fh2dTracksImpParZ
2217  Double_t highIPxy =1.;
2218 
2219  OpenFile(1);
2220  if(fIsPythia){
2221 
2222  if(!fPidResponse){
2223  AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
2224  AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
2225  fPidResponse = inputHandler->GetPIDResponse();
2226  }
2227  if (!fPidResponse) {
2228  AliFatal("NULL PID response");
2229  }
2230  if(!fCombined) fCombined = new AliPIDCombined();
2231  }/*
2232  //Graphs currently not in use
2233  const Int_t gfProtonN = 9;
2234  const Int_t gfAntiProtonN = 18;
2235  const Int_t gfAntiLambdaN = 34;
2236  const Int_t gfLambdaN = 2;
2237  const Int_t gfKMinusN =13 ;
2238  Double_t gfProtonX [gfProtonN] ={0,0.534483,1.29741,2.21552,3.0819,3.92241,4.5819,5.39655,1000};
2239  Double_t gfProtonY [gfProtonN] ={0.990964,0.990964,0.990964,0.990964,0.990964,0.990964,0.990964,0.990964,0.990964};
2240  Double_t gfAntiProtonX [gfAntiProtonN] = {0,0.806034,0.922414,1.09052,1.28448,1.5431,1.73707,1.89224,2.17672,2.43534,2.74569,3.06897,
2241  3.52155,3.88362,4.38793,5.03448,5.38362, 1000};
2242  Double_t gfAntiProtonY [gfAntiProtonN] = {0.922892,0.922892, 0.930723, 0.939157,0.94397,0.95241,0.956627,0.959639,0.964458,
2243  0.966867,0.971084,0.974096,0.978313,0.98012,0.983735,0.986747,0.989157,0.989157};
2244  Double_t gfAntiLambdaX [gfAntiLambdaN] = {0.,0.55555,0.64646,0.75757, 0.84848,0.94949,1.06061,1.15152,1.24242,1.35354,1.44444,
2245  1.54545,1.66667,1.75758,1.84848,1.9596,2.09091,2.30303,2.50505,2.68687,2.90909,3.11111,
2246  3.31313,3.51515,3.69697,3.89899,4.20202,4.66667,5.21212,5.74747,6.50505,7.51515,9.0101,1000};
2247  Double_t gfAntiLambdaY [gfAntiLambdaN] = {0.864925,0.864925,0.895896,0.908209,0.915672,0.921269,0.926866,0.931343,0.935821,0.938806,0.942164,
2248  0.945149,0.947761,0.95,0.952612,0.954478,0.957836,0.960821,0.96306,0.965672,0.968657,0.970149,
2249  0.972015,0.973507,0.975,0.976493,0.978358,0.981343,0.983955,0.986194,0.988433,0.991045,0.991045,0.991045};
2250  Double_t gfLambdaX [gfLambdaN] = {0.,1000};
2251  Double_t gfLambdaY [gfLambdaN] = {0.991045,0.991045};
2252  Double_t gfKMinusX [gfKMinusN] = {0,0.54741,0.74137,1.03879,1.36207,1.96983,2.52586,3.0819,3.67672,4.19397,5.03448,5.44828,1000};
2253  Double_t gfKMinusY [gfKMinusN] = {0,0.979518,0.983133,0.987349,0.989759,0.992169,0.993976,0.996386,0.995783,0.998193,0.99759,1,1000};
2254  fGeant3FlukaProton = new TGraph(gfProtonN,gfProtonX,gfProtonY);
2255  fGeant3FlukaAntiProton = new TGraph(gfAntiProtonN,gfAntiProtonX,gfAntiProtonY);
2256  fGeant3FlukaLambda = new TGraph(gfLambdaN,gfLambdaX,gfLambdaY);
2257  fGeant3FlukaAntiLambda = new TGraph(gfAntiLambdaN,gfAntiLambdaX,gfAntiLambdaY);
2258  fGeant3FlukaKMinus = new TGraph(gfKMinusN,gfKMinusX,gfKMinusY);
2259 
2260 
2261 
2262  //General Information
2263  /*fh1dEventRejectionRDHFCuts = (TH1D*)AddHistogramm("fh1dEventRejectionRDHFCuts","fh1dEventRejectionRDHFCuts;reason;count",12,0,12);
2264  fh1dEventRejectionRDHFCuts->GetXaxis()->SetBinLabel(1,"Event accepted");
2265  fh1dEventRejectionRDHFCuts->GetXaxis()->SetBinLabel(2,"Event rejected");
2266  fh1dEventRejectionRDHFCuts->GetXaxis()->SetBinLabel(3,"Wrong physics selection");
2267  fh1dEventRejectionRDHFCuts->GetXaxis()->SetBinLabel(4,"No vertex");
2268  fh1dEventRejectionRDHFCuts->GetXaxis()->SetBinLabel(5,"No contributors");
2269  fh1dEventRejectionRDHFCuts->GetXaxis()->SetBinLabel(6,"Less than 10 contributors");
2270  fh1dEventRejectionRDHFCuts->GetXaxis()->SetBinLabel(7,">10cm vertex Z distance");
2271  fh1dEventRejectionRDHFCuts->GetXaxis()->SetBinLabel(8,"Bad diamond X distance");
2272  fh1dEventRejectionRDHFCuts->GetXaxis()->SetBinLabel(9,"Bad diamond Y distance");
2273  fh1dEventRejectionRDHFCuts->GetXaxis()->SetBinLabel(10,"Bad diamond Z distance");
2274  fh1dEventRejectionRDHFCuts->GetXaxis()->SetBinLabel(11,"Chi2 vtx >1.5 ");*/
2275  //AddHistogramm("fh1dTracksAccepeted","# tracks before/after cuts;;",3,0,3);
2276 
2277  //****************************************
2278  //QA Plots
2279  fh1dCutsPrinted =(TH1D*)AddHistogramm("fh1dCutsPrinted","",3,0,3);
2280 
2281  fh1dTracksAccepeted =(TH1D*)AddHistogramm("fh1dTracksAccepeted","# tracks before/after cuts;;",3,0,3);
2282  fh1dTracksAccepeted->GetXaxis()->SetBinLabel(1,"total");
2283  fh1dTracksAccepeted->GetXaxis()->SetBinLabel(2,"accepted");
2284  fh1dTracksAccepeted->GetXaxis()->SetBinLabel(3,"rejected");
2285  //Event Properties
2286  fHistManager.CreateTH1("fh1dNoParticlesPerEvent","fh1dNoParticlesvsEvent;#;No Particles/Event",5000, 0, 5000,"s");
2287  fHistManager.CreateTH1("fh1dNoJetsPerEvent","fh1dNoJetsPerEvent;#;No Jets/Event",400, 0, 100,"s");
2288  fHistManager.CreateTH1("fh1dEventsAcceptedInRun","fh1dEventsAcceptedInRun;Events Accepted;count",1,0,1,"s");
2289  fHistManager.CreateTH1("fh1dPtHardMonitor","fh1dPtHardMonitor;ptHard;",500,0,250,"s");
2290  //Jet Properties
2291  fHistManager.CreateTH2("fh1dJetRecEtaPhiAccepted","detector level jet;#eta;phi",1,-0.5,0.5,1,0.,TMath::TwoPi(),"s");
2292  fHistManager.CreateTH2("fh2dAcceptedTracksEtaPhi","accepted tracks;#eta;phi",200,-0.9,0.9,200,0.,TMath::TwoPi(),"s");
2293  fHistManager.CreateTH1("fh1dJetRecPt","detector level jets;pt (GeV/c); count",500,0,250,"s");
2294  fHistManager.CreateTH1("fh1dJetRecPtAccepted","accepted detector level jets;pt (GeV/c); count",500,0,250,"s");
2295  fHistManager.CreateTH1("fh1dJetArea","fh1dJetArea;# Jet Area",100,0,1,"s");
2296  fHistManager.CreateTH1("fh1dParticlesPerJet","fh1dParticlesPerJet;#, Particles/Jet",100,0,100,"s");
2297  fHistManager.CreateTH2("fh2dNoAcceptedTracksvsJetArea","fh2dNoAcceptedTracksvsJetArea;No Accepted Tracks;JetArea",20,0,20,100,0,1);
2298  //MC properties
2299  if(fIsPythia){
2300  /*fHistManager.CreateTH1("fh1dJetGenPt","generator level jets;pt (GeV/c); count",250,0,250,"s");
2301  fHistManager.CreateTH1("fh1dJetGenPtUnidentified","generator level jets (no flavour assigned);pt (GeV/c); count",250,0,250,"s");
2302  fHistManager.CreateTH1("fh1dJetGenPtudsg","generator level udsg jets;pt (GeV/c); count",250,0,250,"s");
2303  fHistManager.CreateTH1("fh1dJetGenPtc","generator level c jets;pt (GeV/c); count",250,0,250,"s");
2304  fHistManager.CreateTH1("fh1dJetGenPtb","generator level b jets;pt (GeV/c); count",250,0,250,"s");
2305  fHistManager.CreateTH1("fh1dJetGenPts","generator level s jets;pt (GeV/c); count",250,0,250,"s");
2306  fHistManager.CreateTH2("fh2dJetGenPtVsJetRecPt","detector momentum response;gen pt;rec pt",500,0,250,500,0,250,"s");*/
2307  fHistManager.CreateTH1("fh1dJetRecPtudsg","detector level jets;pt (GeV/c); count",500,0,250,"s");
2308  fHistManager.CreateTH1("fh1dJetRecPtUnidentified","detector level jets;pt (GeV/c); count",500,0,250,"s");
2309  fHistManager.CreateTH1("fh1dJetRecPtc","detector level jets;pt (GeV/c); count",500,0,250,"s");
2310  fHistManager.CreateTH1("fh1dJetRecPtb","detector level jets;pt (GeV/c); count",500,0,250,"s");
2311  fHistManager.CreateTH1("fh1dJetRecPts","detector level jets;pt (GeV/c); count",500,0,250,"s");
2312  }
2313 
2314  fh1DCutInclusive=(TH1D*)AddHistogramm("fh1DCutInclusive","fh1DCutInclusive",30,0,30);
2315  fh1dCutudg=(TH1D*)AddHistogramm("fh1dCutudg","fh1dCutudg",30,0,30);
2316  fh1dCutc=(TH1D*)AddHistogramm("fh1dCutc","fh1dCutc",30,0,30);
2317  fh1dCutb=(TH1D*)AddHistogramm("fh1dCutb","fh1dCutb",30,0,30);
2318  fh1dCuts=(TH1D*)AddHistogramm("fh1dCuts","fh1dCuts",30,0,30);
2319  fh1dCuts->GetXaxis()->LabelsOption("v");
2320 
2321  for(Int_t iBin = 0; iBin < 25; iBin++){
2322  fh1DCutInclusive->GetXaxis()->SetBinLabel(iBin + 1, BJetCuts[iBin].Data());
2323  if(fIsPythia){
2324  fh1dCutudg->GetXaxis()->SetBinLabel(iBin + 1, BJetCuts[iBin].Data());
2325  fh1dCutb->GetXaxis()->SetBinLabel(iBin + 1, BJetCuts[iBin].Data());
2326  fh1dCutc->GetXaxis()->SetBinLabel(iBin + 1, BJetCuts[iBin].Data());
2327  fh1dCuts->GetXaxis()->SetBinLabel(iBin + 1, BJetCuts[iBin].Data());
2328  }
2329  }
2330 
2331  //****************************************
2332  //Lund Plane
2333  const Int_t dimSpec = 6;
2334  const Int_t nBinsSpec[6] = {50,100,100,20,100,2};
2335  const Double_t lowBinSpec[6] = {0.,-10,0,0,0,0};
2336  const Double_t hiBinSpec[6] = {5.,10.,100,20,100,2};
2337  fHLundIterative = new THnSparseF("fHLundIterative",
2338  "LundIterativePlot [log(1/theta),log(z*theta),pTjet,algo]",
2339  dimSpec,nBinsSpec,lowBinSpec,hiBinSpec);
2340  if(fDoLundPlane)fOutput->Add(fHLundIterative);
2341 
2342 
2343  //****************************************
2344  //Histograms for Probability Tagging
2345  const char * tagtype[6] = {"Full","Single1st","Single2nd","Single3rd","Double","Triple"};
2346 
2347  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
2348  for(int iType=0;iType<6;iType++){
2349  if(fDoJetProb){
2350  fHistManager.CreateTH2(Form("h2DProbDistsTag_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";; #",200, 0, 1,500, 0, 250);
2351  fHistManager.CreateTH2(Form("h2DLNProbDistsTag_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";; #",400, 0, 20,500, 0, 250);
2352  }
2353  if(fDoTCTagging!=TCNo)fHistManager.CreateTH1(Form("h1DTagged_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";; #",500, 0, 250);
2354  }
2355  }
2356  if(fDoJetProb){
2357  fHistManager.CreateTH2("h2DProbDists","h2DProbDistsAll",200, 0, 1,500, 0, 250);
2358  fHistManager.CreateTH2("h2DLNProbDists","h2DProbDistsAll",400, 0, 20,500, 0, 250);
2359 
2360  fHistManager.CreateTH2("h2DProb1Above0","h2DProb1Above0", 400,0,20,500,0,250);
2361  fHistManager.CreateTH2("h2DProb2Above0","h2DProb2Above0", 400,0,20,500,0,250);
2362  fHistManager.CreateTH2("h2DProb3Above0","h2DProb3Above0", 400,0,20,500,0,250);
2363 
2364  if(fIsPythia){
2365  fHistManager.CreateTH2("h2DProb1Above0_UDSG","h2DProb1Above0_UDSG", 400,0,20,500,0,250);
2366  fHistManager.CreateTH2("h2DProb2Above0_UDSG","h2DProb2Above0_UDSG", 400,0,20,500,0,250);
2367  fHistManager.CreateTH2("h2DProb3Above0_UDSG","h2DProb3Above0_UDSG", 400,0,20,500,0,250);
2368  fHistManager.CreateTH2("h2DProb1Above0_C","h2DProb1Above0_C", 400,0,20,500,0,250);
2369  fHistManager.CreateTH2("h2DProb2Above0_C","h2DProb2Above0_C", 400,0,20,500,0,250);
2370  fHistManager.CreateTH2("h2DProb3Above0_C","h2DProb3Above0_C", 400,0,20,500,0,250);
2371  fHistManager.CreateTH2("h2DProb1Above0_B","h2DProb1Above0_B", 400,0,20,500,0,250);
2372  fHistManager.CreateTH2("h2DProb2Above0_B","h2DProb2Above0_B", 400,0,20,500,0,250);
2373  fHistManager.CreateTH2("h2DProb3Above0_B","h2DProb3Above0_B", 400,0,20,500,0,250);
2374  fHistManager.CreateTH2("h2DProb1Above0_V0","h2DProb1Above0_V0", 400,0,20,500,0,250);
2375  fHistManager.CreateTH2("h2DProb2Above0_V0","h2DProb2Above0_V0", 400,0,20,500,0,250);
2376  fHistManager.CreateTH2("h2DProb3Above0_V0","h2DProb3Above0_V0", 400,0,20,500,0,250);
2377  }
2378 
2379  fHistManager.CreateTH2("h2DProb1AboveThresh","h2DProb1AboveThresh", 400,0,20,500,0,250);
2380  fHistManager.CreateTH2("h2DProb1AbThresh1Ab0","h2DProb1AbThresh1Ab0", 400,0,20,500,0,250);
2381  fHistManager.CreateTH2("h2DProb1AbThresh2Ab0","h2DProb1AbThresh2Ab0", 400,0,20,500,0,250);
2382 
2383  fHistManager.CreateTH2("h2DProb2AboveThresh","h2DProb2AboveThresh", 400,0,20,500,0,250);
2384  fHistManager.CreateTH2("h2DProb2AbThresh1Ab0","h2DProb2AbThresh1Ab0", 400,0,20,500,0,250);
2385  fHistManager.CreateTH2("h2DProb2AbThresh2Ab0","h2DProb2AbThresh2Ab0", 400,0,20,500,0,250);
2386  }
2387  if (fIsPythia){
2388  if(fDoMCEffs){
2389  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
2390  for(int iType=0;iType<6;iType++){
2391  fHistManager.CreateTH1(Form("h1DTrueBTagged_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";jet pt; #",500,0,250);
2392  fHistManager.CreateTH1(Form("h1DFalseCTagged_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";jet pt; #",500,0,250);
2393  fHistManager.CreateTH1(Form("h1DFalseUDSGTagged_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";jet pt; #",500,0,250);
2394  fHistManager.CreateTH1(Form("h1DFalseV0Tagged_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";jet pt; #",500,0,250);
2395  }
2396  }
2397  }
2398  if(fDoJetProb){
2399  fHistManager.CreateTH2("h2DLNProbDists_B",";; #",400, 0, 20,500, 0, 250);
2400  fHistManager.CreateTH2("h2DLNProbDists_C",";; #",400, 0, 20,500, 0, 250);
2401  fHistManager.CreateTH2("h2DLNProbDists_UDSG",";; #",400, 0, 20,500, 0, 250);
2402  fHistManager.CreateTH2("h2DLNProbDists_V0",";; #",400, 0, 20,500, 0, 250);
2403 
2404  fHistManager.CreateTH2("h2DProbDists_B",";; #",200, 0, 1,500, 0, 250);
2405  fHistManager.CreateTH2("h2DProbDists_C",";; #",200, 0, 1,500, 0, 250);
2406  fHistManager.CreateTH2("h2DProbDists_UDSG",";; #",200, 0, 1,500, 0, 250);
2407  fHistManager.CreateTH2("h2DProbDists_V0",";; #",200, 0, 1,500, 0, 250);
2408 
2409  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
2410  for(int iType=0;iType<6;iType++){
2411  fHistManager.CreateTH2(Form("h2DLNProbDistsTag_B_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";; #",400, 0, 20,500, 0, 250);
2412  fHistManager.CreateTH2(Form("h2DLNProbDistsTag_C_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";; #",400, 0, 20,500, 0, 250);
2413  fHistManager.CreateTH2(Form("h2DLNProbDistsTag_UDSG_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";; #",400, 0, 20,500, 0, 250);
2414  fHistManager.CreateTH2(Form("h2DLNProbDistsTag_V0_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";; #",400, 0, 20,500, 0, 250);
2415 
2416  fHistManager.CreateTH2(Form("h2DProbDistsTag_B_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";; #",200, 0, 1,500, 0, 250);
2417  fHistManager.CreateTH2(Form("h2DProbDistsTag_C_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";; #",200, 0, 1,500, 0, 250);
2418  fHistManager.CreateTH2(Form("h2DProbDistsTag_UDSG_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";; #",200, 0, 1,500, 0, 250);
2419  fHistManager.CreateTH2(Form("h2DProbDistsTag_V0_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";; #",200, 0, 1,500, 0, 250);
2420  }
2421  }
2422  }
2423  }
2424 
2425  //****************************************
2426  //Track Impact Parameter Distributions
2427  fHistManager.CreateTH2("fh2dTracksImpParXY","radial imp. parameter ;impact parameter xy (cm);a.u.",2000,lowIPxy,highIPxy,500,0,100.,"s");
2428  fHistManager.CreateTH2("fh2dTracksImpParXYZ","XYZ imp. parameter ;impact parameter xy (cm);a.u.",2000,-1,1,500,0,100.,"s");
2429  //fHistManager.CreateTH2("fh2dTracksImpParXYZSignificance","XYZ imp. parameter ;impact parameter xy (cm);a.u.",2000,-30,30,500,0,100.,"s");
2430  fHistManager.CreateTH2("fh2dTracksImpParZ","z imp. parameter ;impact parameter xy (cm);a.u.",2000,lowIPxy,highIPxy,500,0,10.,"s");
2431  fHistManager.CreateTH2("fh2dTracksImpParXYSignificance","radial imp. parameter sig;impact parameter xy (cm);a.u.",2000,-30,30,500,0,100.,"s");
2432  fHistManager.CreateTH2("fh2dTracksImpParZSignificance","z imp. parameter ;impact parameter xy (cm);a.u.",2000,-30,30,500,0,100.,"s");
2433  fHistManager.CreateTH1("fh1dTracksImpParXY","2d imp. parameter ;impact parameter 2d (cm);a.u.",400,-0.2,0.2,"s");
2434  fHistManager.CreateTH1("fh1dTracksImpParXYZ","3d imp. parameter ;impact parameter 3d (cm);a.u.",2000,0,1.,"s");
2435  fHistManager.CreateTH1("fh1dTracksImpParXYSignificance","radial imp. parameter ;impact parameter xy significance;a.u.",200,-30,30,"s");
2436  //fHistManager.CreateTH1 ("fh1dTracksImpParXYZSignificance","3d imp. parameter ;impact parameter 3d significance;a.u.",2000,0.,100.,"s");
2437 
2438  fHistManager.CreateTH2("fh1dTracksIPvsPt_B","Track IP vs Track Pt; p_{T,Track} (GeV/c); d_{0}",500,0,250,300,0,30);
2439  fHistManager.CreateTH2("fh1dTracksIPvsPt_V0JetTracks","Track IP vs Track Pt; p_{T,Track} (GeV/c); d_{0}",500,0,250,300,0,30);
2440  fHistManager.CreateTH2("fh1dTracksIPvsPt_V0inV0Jet","Track IP vs Track Pt; p_{T,Track} (GeV/c); d_{0}",500,0,250,300,0,30);
2441  fHistManager.CreateTH2("fh1dTracksIPvsPt_V0inBJet","Track IP vs Track Pt; p_{T,Track} (GeV/c); d_{0}",500,0,250,300,0,30);
2442 
2443 
2444 
2445  //****************************************
2446  //Pt Distributions for N1,N2,N3 Tracks
2447  if(fIsPythia){
2448  for(int iN=1;iN<=3;iN++){
2449  fHistManager.CreateTH1(Form("fh1dJetRecPt_n_%i_all_Accepted",iN),"detector level jets;pt (GeV/c); count",500,0,250,"s");
2450  fHistManager.CreateTH1(Form("fh1dTrackPt_n_%i_all_Accepted",iN),"detector level jets;pt (GeV/c); count",500,0,200,"s");
2451  for(long unsigned iFlav=0;iFlav<sTemplateFlavour.size();iFlav++){
2452  fHistManager.CreateTH1(Form("fh1dJetRecPt_n_%i_%s_Accepted",iN,sTemplateFlavour[iFlav].Data()),"detector level jets;pt (GeV/c); count",500,0,250,"s");
2453  }
2454  }
2455  }
2456 
2457  //V0Cuts
2458  if(fApplyV0Rej){
2459  //V0s from reconstruction
2460  const Int_t iNDimInJC = 5;
2461  Int_t binsKInJC[iNDimInJC] = {200, 200, 200, 200,5};
2462  Double_t xminKInJC[iNDimInJC] = {0.35, 0., -1., 0.,0};
2463  Double_t xmaxKInJC[iNDimInJC] = {0.65, 50., 1., 200.,5};
2464  Int_t binsLInJC[iNDimInJC] = {200, 200, 200, 200,5};
2465  Double_t xminLInJC[iNDimInJC] = {1.05, 0., -1., 0.,0};
2466  Double_t xmaxLInJC[iNDimInJC] = {1.25, 50., 1., 200.,5};
2467 
2468  fh2dKshortMassVsPt=(TH2D*)AddHistogramm("fh2dKshortMassVsPt","KShort Mass Vs Pt;p_{T} (GeV/c);Mass (GeV)",200,0,50,200,0.35, 0.65);
2469  fh2dLamdaMassVsPt =(TH2D*)AddHistogramm("fh2dLamdaMassVsPt","Lamda Mass Vs Pt;p_{T} (GeV/c);Mass (GeV)",200,0,50,200,1.05,1.25);
2470  fh2dAnLamdaMassVsPt =(TH2D*)AddHistogramm("fh2dAnLamdaMassVsPt","Anti Lamda Mass Vs Pt;p_{T} (GeV/c);Mass (GeV)",200,0,50,200,1.05,1.25);
2471 
2472  fhnV0InJetK0s = new THnSparseD("fhnV0InJetK0s", "K0s: Mass vs Pt in jets;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c}); IsMCTrueK0", iNDimInJC, binsKInJC, xminKInJC, xmaxKInJC);
2473  fhnV0InJetLambda = new THnSparseD("fhnV0InJetLambda", "Lambda: Mass vs Pt in jets;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c}); IsMCTrueLamda", iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
2474  fhnV0InJetALambda = new THnSparseD("fhnV0InJetALambda", "ALambda: Mass vs Pt in jets;#it{m}_{inv} (GeV/#it{c}^{2});#it{p}_{T}^{V0} (GeV/#it{c});#it{#eta}_{V0};#it{p}_{T}^{jet} (GeV/#it{c}); IsMCTrueALamda", iNDimInJC, binsLInJC, xminLInJC, xmaxLInJC);
2475 
2476  fh1V0CounterCentK0s = (TH1D*)AddHistogramm(Form("fh1V0CounterCentK0s_%d", 0), Form("Number of K0s candidates after cuts, cent %s;cut;counts","0-100%"), 18, 0, 18);
2477  fh1V0CounterCentLambda = (TH1D*)AddHistogramm(Form("fh1V0CounterCentLambda_%d", 0), Form("Number of Lambda candidates after cuts, cent %s;cut;counts", "0-100%"), 18, 0, 18);
2478  fh1V0CounterCentALambda = (TH1D*)AddHistogramm(Form("fh1V0CounterCentALambda_%d", 0), Form("Number of ALambda candidates after cuts, cent %s;cut;counts", "0-100%"), 18, 0, 18);
2479 
2480  fOutput->Add(fhnV0InJetK0s);
2481  fOutput->Add(fhnV0InJetLambda);
2482  fOutput->Add(fhnV0InJetALambda);
2483 
2484  for(Int_t j = 0; j < 18; j++){
2485  fh1V0CounterCentK0s->GetXaxis()->SetBinLabel(j + 1, sV0Cuts[j].Data());
2486  fh1V0CounterCentLambda->GetXaxis()->SetBinLabel(j + 1, sV0Cuts[j].Data());
2487  fh1V0CounterCentALambda->GetXaxis()->SetBinLabel(j + 1, sV0Cuts[j].Data());
2488  }
2489  if(fV0CandidateArray == NULL){
2490  fV0CandidateArray = new TClonesArray("AliAODv0",100);
2491  }
2492  fV0CandidateArray->Delete();//Reset the TClonesArray
2493 
2494  //V0s from MC
2495  fh1dKshortPtMC = new TH1D("fh1dKshortPtMC","KShort Pt MC;p_{T} (GeV/c)",200,0,50);
2496  fh1dLamdaPtMC = new TH1D("fh1dLamdaPtMC","Lamda Pt MC;p_{T} (GeV/c)",200,0,50);
2497  fh1dAnLamdaPtMC = new TH1D("fh1dAnLamdaPtMC","Anti Lamda Pt MC;p_{T} (GeV/c)",200,0,50);
2498  fh2dKshortPtVsJetPtMC = new TH2D("fh2dKshortPtVsJetPtMC","KShort Pt Vs Jet Pt MC;p_{T,V0} (GeV/c);#it{p}_{T,jet} (GeV/#it{c})",200,0,50,200,0, 200);
2499  fh2dLamdaPtVsJetPtMC = new TH2D("fh2dLamdaPtVsJetPtMC","Lamda Pt Vs Jet Pt MC;p_{T,V0} (GeV/c);#it{p}_{T,jet} (GeV/#it{c})",200,0,50,200,0,200);
2500  fh2dAnLamdaPtVsJetPtMC = new TH2D("fh2dAnLamdaPtVsJetPtMC","Anti Lamda Pt Vs Jet Pt MC;p_{T,V0} (GeV/c);#it{p}_{T,jet} (GeV/#it{c})",200,0,50,200,0,200);
2501 
2502  //V0 Tagging efficiency
2503  h1DV0FalseRec=(TH1D*)AddHistogramm("h1DV0FalseRec",";jetpt (GeV/c); N_{V0}",500,0,250);
2504  h1DV0TrueDataDef =(TH1D*)AddHistogramm("h1DV0TrueDataDef",";jetpt;N_{V0}",500,0,250);
2505  h1DV0TrueMCDef =(TH1D*)AddHistogramm("h1DV0TrueMCDef",";jetpt;N_{V0}",500,0,250);
2506  h1DV0TrueRec =(TH1D*)AddHistogramm("h1DV0TrueRec",";jetpt;N_{V0}",500,0,250);
2507 
2508  fOutput->Add(fh1dKshortPtMC);
2509  fOutput->Add(fh1dLamdaPtMC);
2510  fOutput->Add(fh1dAnLamdaPtMC);
2511  fOutput->Add(fh2dKshortPtVsJetPtMC);
2512  fOutput->Add(fh2dLamdaPtVsJetPtMC);
2513  fOutput->Add(fh2dAnLamdaPtVsJetPtMC);
2514  }
2515 
2516  //Final Tagged JetPt Spectra
2517  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
2518  for(int iType=0;iType<6;iType++){
2519  fHistManager.CreateTH1(Form("h1DTaggedJetPt_%s_%0.2f",tagtype[iType],fFracs[iThresh]),";jet pt; #",500, 0, 250);
2520  }
2521  }
2522 
2523  //Impact Parameter Template Generation
2524  const char * base = "fh2dJetSignedImpPar";
2525  const char * dim[2] = {"XY","XYZ"};
2526  const char * typ[2] = {"","Significance"};
2527  const char * ordpar [4] = {"","First","Second","Third"};
2528  const char * special [1] = {"",/*"McCorr"*/};
2529 
2530  Int_t ptbins = 250;
2531  Double_t ptlow = 0;
2532  Double_t pthigh = 250;
2533  Int_t ipbins = 1000;
2534  Double_t iplow = -.5;
2535  Double_t iphigh = .5;
2536  for (Int_t id = 0;id<1;++id) // XY or XY/
2537  for (unsigned int ifl = 0;ifl<sTemplateFlavour.size();++ifl) //flavour
2538  for (Int_t io = 0;io<4;++io) //order parameter
2539  for (Int_t is = 0;is<1;++is) //special comment
2540  for (Int_t it = 1;it<2;++it){ //significance or not
2541  if(it==1) {
2542  iplow=-30;
2543  iphigh=30; //from 30
2544  if(io==0 && ifl==4) ipbins = 1000;//2000;
2545  else ipbins =1000;//2000;\
2546  }else {
2547  iplow=-0.5;
2548  iphigh=0.5;
2549  ipbins =1000;//;2000;
2550  }
2551  if(id==0) ipbins =1000;//2000;
2552  if((fIsPythia||(!fIsPythia && ifl==6))){
2553  fHistManager.CreateTH2(Form("%s%s%s%s%s%s",base,dim[id],typ[it],sTemplateFlavour[ifl].Data(),ordpar[io],special[is]),
2554  Form("%s%s%s%s%s%s;;",base,dim[id],typ[it],sTemplateFlavour[ifl].Data(),ordpar[io],special[is]),
2555  ptbins,ptlow,pthigh,ipbins,iplow,iphigh,"s");
2556  //printf("Generating%s%s%s%s%s%s",base,dim[id],typ[it],flavour[ifl],ordpar[io],special[is]);
2557 
2558  }
2559  }
2560 
2561  TIter next(fHistManager.GetListOfHistograms());
2562  TObject* obj = 0;
2563  while ((obj = next())) {
2564  printf("Adding Object %s\n",obj->GetName());
2565  fOutput->Add(obj);
2566  }
2567 
2568  PostData(1, fOutput);
2569 }
2570 
2571 void AliAnalysisTaskHFJetIPQA::UserExecOnce(){
2572  AliJetContainer * jetconrec = static_cast<AliJetContainer*>(fJetCollArray.At(0));
2573  fAnalysisCuts[bAnalysisCut_MinJetPt]=jetconrec->GetJetPtCut();
2574  fAnalysisCuts[bAnalysisCut_MaxJetPt]=jetconrec->GetJetPtCutMax();
2575  fAnalysisCuts[bAnalysisCut_MinJetEta]=jetconrec->GetMinEta();
2576  fAnalysisCuts[bAnalysisCut_MaxJetEta]=jetconrec->GetMaxEta();
2577 
2578  PrintSettings();
2579  PrintV0Settings();
2580 }
2581 
2582 void AliAnalysisTaskHFJetIPQA::PrintV0Settings(){
2583  TString v0cuts="";
2584  Int_t version=1;
2585  int sForm[24]={1,2,2,0,0,0,0,1,1,0,
2586  0,1,0,1,0,3,3,0,0,0,3,3,2,2};
2587 
2588  v0cuts=version;
2589  for(int iV0Cut=0;iV0Cut<24;iV0Cut++){
2590  v0cuts+="+";
2591  switch (sForm[iV0Cut]){
2592  case 0:
2593  v0cuts+=Form("%0.f",fV0Cuts[iV0Cut]);
2594  break;
2595 
2596  case 1:
2597  v0cuts+=Form("%0.1f",fV0Cuts[iV0Cut]);
2598  break;
2599 
2600  case 2:
2601  v0cuts+=Form("%0.2f",fV0Cuts[iV0Cut]);
2602  break;
2603 
2604  case 3:
2605  v0cuts+=Form("%0.3f",fV0Cuts[iV0Cut]);
2606  break;
2607  }
2608  }
2609 
2610  fh1V0CounterCentK0s->SetTitle(v0cuts.Data());
2611 }
2612 
2613 void AliAnalysisTaskHFJetIPQA::PrintSettings(){
2614  TString jetcuts="";
2615  TString trackcuts="";
2616  TString vertexcuts="";
2617  Int_t version=3;
2618 
2619 
2620  jetcuts+=version;
2621  jetcuts+="+";
2622  jetcuts+=fAnalysisCuts[bAnalysisCut_MinJetPt];
2623  jetcuts+="+";
2624  jetcuts+=fAnalysisCuts[bAnalysisCut_MaxJetPt];
2625  jetcuts+="+";
2626  jetcuts+=fAnalysisCuts[bAnalysisCut_MinJetEta];
2627  jetcuts+="+";
2628  jetcuts+=fAnalysisCuts[bAnalysisCut_MaxJetEta];
2629  jetcuts+="+";
2630  jetcuts+=fNoJetConstituents;
2631  jetcuts+="+";
2632  jetcuts+=fDaughtersRadius;
2633  jetcuts+="+";
2634  jetcuts+=Form("%0.1f",fJetRadius);
2635 
2636  printf("Cut Track Settings: %s\n",jetcuts.Data());
2637 
2638  trackcuts+=version;
2639  trackcuts+="+";
2640  trackcuts+=Form("%0.2f",fAnalysisCuts[bAnalysisCut_DCAJetTrack]);
2641  trackcuts+="+";
2642  trackcuts+=fAnalysisCuts[bAnalysisCut_MaxDecayLength];
2643  trackcuts+="+";
2644  trackcuts+=fAnalysisCuts[bAnalysisCut_MinTrackPt];
2645  trackcuts+="+";
2646  trackcuts+=fAnalysisCuts[bAnalysisCut_MinTrackPtMC];
2647  trackcuts+="+";
2648  trackcuts+=fAnalysisCuts[bAnalysisCut_MaxDCA_Z];
2649  trackcuts+="+";
2650  trackcuts+=fAnalysisCuts[bAnalysisCut_MaxDCA_XY];
2651  trackcuts+="+";
2652  trackcuts+=fAnalysisCuts[bAnalysisCut_MinTPCClus];
2653  trackcuts+="+";
2654  trackcuts+=fAnalysisCuts[bAnalysisCut_MinITSLayersHit];
2655  trackcuts+="+";
2656  trackcuts+=fAnalysisCuts[bAnalysisCut_MinTrackChi2];
2657  trackcuts+="+";
2658  trackcuts+=fAnalysisCuts[bAnalysisCut_HasSDD];
2659  trackcuts+="+";
2660  trackcuts+=fAnalysisCuts[bAnalysisCut_KinkCand];
2661  trackcuts+="+";
2662  trackcuts+=fAnalysisCuts[bAnalysisCut_HasTPCrefit];
2663  trackcuts+="+";
2664  trackcuts+=fAnalysisCuts[bAnalysisCut_HasITSrefit];
2665 
2666 
2667  printf("Cut Vertex Settings %s\n", trackcuts.Data());
2668 
2669  vertexcuts+=version;
2670  vertexcuts+="+";
2671  vertexcuts+=Form("%0.f",fAnalysisCuts[bAnalysisCut_NContibutors]);
2672  vertexcuts+="+";
2673  vertexcuts+=Form("%0.f",fAnalysisCuts[bAnalysisCut_Sigma_Z]);
2674  vertexcuts+="+";
2675  vertexcuts+=Form("%0.f",fAnalysisCuts[bAnalysisCut_Sigma_Y]);
2676  vertexcuts+="+";
2677  vertexcuts+=Form("%0.f",fAnalysisCuts[bAnalysisCut_RelError_Z]);
2678  vertexcuts+="+";
2679  vertexcuts+=Form("%0.f",fAnalysisCuts[bAnalysisCut_RelError_Y]);
2680  vertexcuts+="+";
2681  vertexcuts+=fAnalysisCuts[bAnalysisCut_SigmaDiamond];
2682  vertexcuts+="+";
2683  vertexcuts+=fAnalysisCuts[bAnalysisCut_PtHardAndJetPtFactor];
2684  vertexcuts+="+";
2685  vertexcuts+=fAnalysisCuts[bAnalysisCut_MinNewVertexContrib];
2686  vertexcuts+="+";
2687  vertexcuts+=fAnalysisCuts[bAnalysisCut_MaxVtxZ];
2688  vertexcuts+="+";
2689  vertexcuts+=fAnalysisCuts[bAnalysisCut_Z_Chi2perNDF];
2690  vertexcuts+="+";
2691  vertexcuts+=fVertexRecalcMinPt;
2692  vertexcuts+="+";
2693  vertexcuts+=fDoMCCorrection;
2694  vertexcuts+="+";
2695  vertexcuts+=fRunSmearing;
2696  vertexcuts+="+";
2697  vertexcuts+=fDoUnderlyingEventSub;
2698  vertexcuts+="+";
2699  vertexcuts+=fDoFlavourMatching;
2700  vertexcuts+="+";
2701  vertexcuts+=fDoTCTagging;
2702  vertexcuts+="+";
2703  vertexcuts+=fDoProbTagging;
2704  vertexcuts+="+";
2705  vertexcuts+=fUseSignificance;
2706  vertexcuts+="+";
2707  vertexcuts+=fDoJetProb;
2708  vertexcuts+="+";
2709  vertexcuts+=Form("%0.3f",fTCThresholdPtFixed);
2710 
2711  fh1dCutsPrinted->SetTitle(jetcuts.Data());
2712  fh1dCutsPrinted->GetXaxis()->SetTitle(trackcuts.Data());
2713  fh1dCutsPrinted->GetYaxis()->SetTitle(vertexcuts.Data());
2714 
2715  printf("Vertex Cuts: %s\n",vertexcuts.Data());
2716 }
2717 
2718 //NotInUse
2719 /*void AliAnalysisTaskHFJetIPQA::GetMaxImpactParameterCutR(const AliVTrack * const track, Double_t &maximpactRcut){
2720  //
2721  // Get max impact parameter cut r (pt dependent)
2722  //
2723  Double_t pt = track->Pt();
2724  if(pt > 0.15) {
2725  maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01); // abs R cut
2726  }
2727  else maximpactRcut = 9999999999.0;
2728  }*/
2729 
2730 
2731  bool AliAnalysisTaskHFJetIPQA::GetPIDCombined(AliAODTrack * track, double * prob, int &nDetectors,UInt_t &usedDet ,AliPID::EParticleType &MostProbablePID, bool setTrackPID ){
2732  AliPIDResponse::EDetPidStatus status[AliPIDResponse::kNdetectors] = {AliPIDResponse::kDetNoSignal,AliPIDResponse::kDetNoSignal,AliPIDResponse::kDetNoSignal,AliPIDResponse::kDetNoSignal,AliPIDResponse::kDetNoSignal,AliPIDResponse::kDetNoSignal,AliPIDResponse::kDetNoSignal};
2733  unsigned int nGoodDet = 0;
2734  for (int j =0; j<AliPIDResponse::kNdetectors;j++)
2735  {
2736  for (int i =AliPID::kElectron; i<AliPID::kSPECIES;i++)
2737  {
2738  double val = 0;
2739  status[j] = fPidResponse->NumberOfSigmas(static_cast <AliPIDResponse::EDetector>(j), track, static_cast <AliPID::EParticleType>(i), val);
2740  if (status[j] == AliPIDResponse::kDetPidOk ){
2741  nGoodDet++;}
2742  }
2743  }
2744  if( nGoodDet/7 <2 ) return false;
2745  nDetectors = nGoodDet/7;
2746  Double_t probTPCTOF[AliPID::kSPECIES]={-1.,-1.,-1.,-1.,-1.};
2747  fCombined->SetDefaultTPCPriors();
2748  fCombined->SetDetectorMask(AliPIDResponse::kDetITS|AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF|AliPIDResponse::kDetTRD|AliPIDResponse::kDetEMCAL);
2749  usedDet = fCombined->ComputeProbabilities((AliVTrack*)track, fPidResponse , probTPCTOF);
2750  int maxpid=14;
2751  double maxpidv=0;
2752  for (int j =0 ;j< AliPID::kSPECIES;++j ){
2753  prob[j] =probTPCTOF[j];
2754  if(prob[j] >maxpidv ){
2755  maxpid=j;
2756  maxpidv=prob[j] ;
2757  }
2758  }
2759  MostProbablePID=static_cast <AliPID::EParticleType>(maxpid);
2760  return true;
2761  }
2762 
2763 
2764  bool AliAnalysisTaskHFJetIPQA::IsFromElectron(AliAODTrack*track){
2765 
2766  Double_t p[5] ={0};
2767  Int_t nDet=0;
2768  UInt_t nDetUCom=0;
2769  AliPID::EParticleType mpPID=AliPID::kUnknown;
2770  if(GetPIDCombined( track, p, nDet,nDetUCom ,mpPID, kFALSE )){
2771  if(mpPID==AliPID::kElectron && p[0] >0.90 && nDet >1) return kTRUE;
2772  }
2773  return false;
2774  }
2775 
2776  bool AliAnalysisTaskHFJetIPQA::IsFromPion(AliAODTrack*track){
2777  Double_t p[5] ={0};
2778  Int_t nDet=0;
2779  UInt_t nDetUCom=0;
2780  AliPID::EParticleType mpPID=AliPID::kUnknown;
2781 
2782  if(GetPIDCombined( track, p, nDet,nDetUCom ,mpPID, kFALSE )){
2783  if(mpPID==AliPID::kPion && p[2] >0.90 && nDet >1) return kTRUE;
2784  }
2785  return false;
2786  }
2787  bool AliAnalysisTaskHFJetIPQA::IsFromKaon(AliAODTrack*track){
2788  Double_t p[5] ={0};
2789  Int_t nDet=0;
2790  UInt_t nDetUCom=0;
2791  AliPID::EParticleType mpPID=AliPID::kUnknown;
2792  if(GetPIDCombined( track, p, nDet,nDetUCom ,mpPID, kFALSE )){
2793  if(mpPID==AliPID::kKaon && p[3] >0.90 && nDet >1) return kTRUE;
2794  }
2795  return false;
2796  }
2797  bool AliAnalysisTaskHFJetIPQA::IsFromProton(AliAODTrack*track){
2798  Double_t p[5] ={0};
2799  Int_t nDet=0;
2800  UInt_t nDetUCom=0;
2801  AliPID::EParticleType mpPID=AliPID::kUnknown;
2802  if(GetPIDCombined( track, p, nDet,nDetUCom ,mpPID, kFALSE )){
2803  if(mpPID==AliPID::kProton && p[3] >0.90 && nDet >1) return kTRUE;
2804  }
2805  return false;
2806  }
2807 
2808  int AliAnalysisTaskHFJetIPQA::DetermineUnsuitableVtxTracks(int *skipped, AliAODEvent * const aod, AliVTrack * const track){
2809  Int_t nTracks=aod->GetNumberOfTracks();
2810  AliAODTrack * t = nullptr;
2811  AliExternalTrackParam etp_at_r39_old; etp_at_r39_old.CopyFromVTrack(track);
2812  etp_at_r39_old.PropagateTo(3.9,InputEvent()->GetMagneticField());
2813  double angle0 = TMath::ATan2(etp_at_r39_old.Yv(),etp_at_r39_old.Xv());
2814  double zz0 = etp_at_r39_old.GetZ();
2815  int nTrksToSkip=1;
2816 
2817  for(Int_t i=0; i<nTracks; i++){
2818  t = (AliAODTrack *)(aod->GetTrack(i));
2819  if(!((((AliAODTrack*)t)->TestFilterBit(4))))continue;
2820  int id = (Int_t)t->GetID();
2821  AliExternalTrackParam etp_at_r39; etp_at_r39.CopyFromVTrack(t);
2822  etp_at_r39.PropagateTo(3.9,InputEvent()->GetMagneticField());
2823  double angle = TMath::ATan2(etp_at_r39.Yv(),etp_at_r39.Xv());
2824  double zz = etp_at_r39.GetZ();
2825  bool doskip=false;
2826  if(t->Pt()<fVertexRecalcMinPt) doskip=true;
2827  if(fabs(TVector2::Phi_mpi_pi(angle-angle0))>TMath::Pi()/6.) {
2828  doskip=true;
2829  }
2830  if(fabs(zz-zz0)>0.5) {
2831  doskip=true;
2832  }
2833  if(doskip && !fGlobalVertex){
2834  skipped[nTrksToSkip++] = id;
2835  }
2836  }
2837  return nTrksToSkip;
2838  }
2839 
2840 AliAODVertex *AliAnalysisTaskHFJetIPQA::RemoveDaughtersFromPrimaryVtx( const AliVTrack * const track) {
2841  //Initialisation of vertexer
2842  const AliAODEvent * aod = ((AliAODEvent*)InputEvent());
2843  AliAODVertex *vtxAOD =aod ->GetPrimaryVertex();
2844  if(!vtxAOD) return 0;
2845  //printf("Before remove:\n");
2846  //UvtxAOD->Print();
2847 
2848  TString title=vtxAOD->GetTitle();
2849  if(!title.Contains("VertexerTracks")) return 0;
2850 
2851  AliVertexerTracks vertexer(aod->GetMagneticField());
2852  vertexer.SetITSMode();
2853  vertexer.SetMinClusters(3);
2854  if(title.Contains("WithConstraint")) {
2855  Float_t diamondcovxy[3];
2856  aod->GetDiamondCovXY(diamondcovxy);
2857  Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
2858  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
2859  AliESDVertex diamond(pos,cov,1.,1);
2860  vertexer.SetVtxStart(&diamond);
2861  }
2862 
2863  //_____________________________
2864  //Determination of unsuited tracks
2865  Int_t skipped[5000]; for(Int_t i=0;i<5000;i++) skipped[i]=-1;
2866  Int_t id = (Int_t)track->GetID();
2867  if(!(id<0)) skipped[0] = id; //remove track under investigation from vertex
2868  int nTrksToSkip=1;
2869  //nTrksToSkip=DetermineUnsuitableVtxTracks(skipped, aod, track);
2870  vertexer.SetSkipTracks(nTrksToSkip,skipped);
2871 
2872  //________________________________
2873  //Determination of new ESD vertex
2874  AliESDVertex *vtxESDNew = vertexer.FindPrimaryVertex(aod);
2875  if(!vtxESDNew) return 0;
2876  Int_t nContrib =vtxESDNew->GetNContributors();
2877  if(vtxESDNew->GetNContributors()<fAnalysisCuts[bAnalysisCut_MinNewVertexContrib]) {
2878  //printf("Thrown away with %i contributors\n",vtxESDNew->GetNContributors());
2879  delete vtxESDNew; vtxESDNew=nullptr;
2880  return 0;
2881  }
2882  if(vtxESDNew->GetChi2toNDF()>fAnalysisCuts[bAnalysisCut_Z_Chi2perNDF]) {
2883  //printf("Thrown away with chi2 = %f\n",fAnalysisCuts[bAnalysisCut_Z_Chi2perNDF]);
2884  delete vtxESDNew; vtxESDNew=nullptr;
2885  return 0;
2886  }
2887 
2888  //________________________________
2889  //Conversion to AOD vertex
2890  Double_t pos[3];
2891  Double_t cov[6];
2892  Double_t chi2perNDF;
2893  vtxESDNew->GetXYZ(pos); // position
2894  vtxESDNew->GetCovMatrix(cov); //covariance matrix
2895  chi2perNDF = vtxESDNew->GetChi2toNDF(); //chisquare
2896  if(vtxESDNew) delete vtxESDNew;
2897  vtxESDNew=NULL;
2898  AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);
2899  vtxAODNew->SetNContributors(nContrib); //contributors
2900  return vtxAODNew;
2901 }
2902 void AliAnalysisTaskHFJetIPQA::FillParticleCompositionSpectra(AliEmcalJet * jet,const char * histname ){
2903  if(!jet) return;
2904  AliVTrack* tr=0x0;
2905  for(Int_t j = 0; j < jet->GetNumberOfTracks(); ++j) {
2906  tr = (AliVTrack*)GetParticleContainer(0)->GetParticle((jet->TrackAt(j)));
2907  if(!tr) continue;
2908  double pT=0x0;
2909  Int_t pCorr_indx=-1;
2910  GetMonteCarloCorrectionFactor(tr,pCorr_indx,pT);
2911  if(pCorr_indx<0 ) continue;
2912  FillHist(histname,pCorr_indx+0.5,pT, 1); //*this->fXsectionWeightingFactor );
2913  }
2914  return;
2915 }
2916 
2922 void AliAnalysisTaskHFJetIPQA::FillParticleCompositionEvent( ){
2923  AliVTrack* tr=0x0;
2924  for(Int_t j = 0; j < InputEvent()->GetNumberOfTracks(); ++j) {
2925  tr = (AliVTrack*)(InputEvent()->GetTrack(j));
2926  if(!tr) continue;
2927  double pT=0x0;
2928  Int_t pCorr_indx=-1;
2929  GetMonteCarloCorrectionFactor(tr,pCorr_indx,pT);
2930  if(pCorr_indx<0) continue;
2931  FillHist("fh2dParticleSpectra_Event",pCorr_indx+0.5,pT, 1); //*this->fXsectionWeightingFactor );
2932  }
2933  return;
2934 }
2935 
2936 
2937 void AliAnalysisTaskHFJetIPQA::GetOutOfJetParticleComposition(AliEmcalJet * jet, int flavour){
2938  if(!jet) return;
2939  AliEmcalJet * perp_jet1 = GetPerpendicularPseudoJet(jet,false);
2940  AliEmcalJet * perp_jet2 = GetPerpendicularPseudoJet(jet,true);
2941  FillParticleCompositionSpectra(jet,"fh2dParticleSpectra_InCone");
2942  if(flavour==3) FillParticleCompositionSpectra(jet,"fh2dParticleSpectra_InCone_bjet");
2943  else if(flavour==2) FillParticleCompositionSpectra(jet,"fh2dParticleSpectra_InCone_cjet");
2944  else if(flavour==1) FillParticleCompositionSpectra(jet,"fh2dParticleSpectra_InCone_lfjet");
2945  FillParticleCompositionSpectra(perp_jet1,"fh2dParticleSpectra_OutOfCone");
2946  FillParticleCompositionSpectra(perp_jet2,"fh2dParticleSpectra_OutOfCone");
2947  if(perp_jet1) delete perp_jet1;
2948  if(perp_jet2) delete perp_jet2;
2949 }
2950 
2951 AliEmcalJet * AliAnalysisTaskHFJetIPQA::GetPerpendicularPseudoJet (AliEmcalJet *jet_in , bool rev ){
2952  TVector3 j(jet_in->Px(), jet_in->Py(), jet_in->Pz());
2953  TVector3 p1(j);
2954  std::vector <int > track_inc;
2955  p1.RotateZ(rev ? -1*TMath::Pi()/2. :TMath::Pi()/2. );
2956  Double_t sumAllPt1 = 0;
2957  int nconst_1 =0;
2958  std::vector <int> const_idx1;
2959  for(long itrack= 0; itrack<GetParticleContainer(0)->GetNParticles();++itrack){
2960  AliVTrack * tr = static_cast<AliVTrack*>(GetParticleContainer(0)->GetParticle(itrack));
2961  if(!tr) continue;
2962  TVector3 v(tr->Px(), tr->Py(), tr->Pz());
2963  Double_t dR1 = v.DrEtaPhi(p1);
2964  if(v.Pt()>0.150){
2965  if(dR1 < fJetRadius) {
2966  sumAllPt1+=v.Pt();
2967  nconst_1++;
2968  const_idx1.push_back(itrack);
2969  }
2970  }
2971 }
2972 AliEmcalJet* jet1 =0;
2973 
2974 
2975 if (sumAllPt1>0) {
2976  jet1 = new AliEmcalJet(sumAllPt1, p1.Eta(), TVector2::Phi_0_2pi (p1.Phi()), 0);
2977  jet1->SetArea(fJetRadius*fJetRadius*TMath::Pi());
2978  jet1->SetNumberOfTracks(nconst_1);
2979  jet1->SetNumberOfClusters(0);
2980  for (int i = 0 ; i < (int) const_idx1.size();++i) {
2981  jet1->AddTrackAt(const_idx1.at(i), i);
2982  }}
2983 
2984 
2985  return jet1;
2986 }
2987 
2988 
2989 
2990 
2991 
2992 
2993 
2994 
2995 
2996 
3003 Double_t AliAnalysisTaskHFJetIPQA::GetValImpactParameter(TTypeImpPar type,Double_t *impar, Double_t * cov)
3004 {
3005  Float_t result = -999999;
3006  Float_t dFdx = 0;
3007  Float_t dFdy = 0;
3008  switch(type){
3009  case kXY:
3010  result = impar[0];
3011  break;
3012  case kXYSig:
3013  result = impar[0]/sqrt(cov[0]);
3014  break;
3015  case kXYZ:
3016  result = TMath::Sqrt(impar[0]*impar[0]+impar[1]*impar[1]);
3017  break;
3018  case kXYZSig:
3019  result = TMath::Sqrt(impar[0]*impar[0]+impar[1]*impar[1]);
3020  dFdx = impar[0]/result;
3021  dFdy = impar[1]/result;
3022  result /=TMath::Sqrt(cov[0]*dFdx*dFdx + cov[2]*dFdy*dFdy + 2* cov[1] *dFdx*dFdy);
3023  break;
3024  case kZSig:
3025  result = impar[1];
3026  result /=TMath::Sqrt(cov[2]);
3027  break;
3028  case kXYZSigmaOnly:
3029  result = TMath::Sqrt(impar[0]*impar[0]+impar[1]*impar[1]);
3030  dFdx = impar[0]/result;
3031  dFdy = impar[1]/result;
3032  result =TMath::Sqrt(cov[0]*dFdx*dFdx + cov[2]*dFdy*dFdy + 2* cov[1] *dFdx*dFdy);
3033  break;
3034  default:
3035  break;
3036  }
3037  return result;
3038 }
3039 
3040 //____________________________________________________
3041 void AliAnalysisTaskHFJetIPQA::FillCandidateJet(Int_t CutIndex, Int_t JetFlavor){
3042  if(JetFlavor>0) fh1DCutInclusive->Fill(CutIndex);
3043  if(fIsPythia){
3044  if(JetFlavor==1)fh1dCutudg->Fill(CutIndex);
3045  if(JetFlavor==2)fh1dCutc->Fill(CutIndex);
3046  if(JetFlavor==3)fh1dCutb->Fill(CutIndex);
3047  if(JetFlavor==4)fh1dCuts->Fill(CutIndex);
3048  }
3049 
3050 }
3051 
3052 Bool_t AliAnalysisTaskHFJetIPQA::IsTrackAccepted(AliVTrack* track , int jetflavour){
3053  if(!track) return kFALSE;
3054  if(fIsEsd){
3055  fESDTrackCut->SetMinNClustersITS(fAnalysisCuts[bAnalysisCut_MinITSLayersHit]);
3056  fESDTrackCut->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
3057  if(!(fESDTrackCut->AcceptTrack((AliESDtrack*)track))) return kFALSE;
3058  return kTRUE;
3059  }
3060  else {
3061  //HasMatchedGoodTracklet((AliAODTrack*)track);
3062  if(!(((AliAODTrack*)track)->TestFilterBit(9) || ((AliAODTrack*)track)->TestFilterBit(4)))return kFALSE;
3063 
3064  if(TMath::Abs(track->Eta())>0.9) return kFALSE;
3065 
3066  //SPD hits
3067  if(fAnalysisCuts[bAnalysisCut_HasSDD]){
3068  if(!(((AliAODTrack*)track->HasPointOnITSLayer(0))&&(AliAODTrack*)track->HasPointOnITSLayer(1))){
3069  FillCandidateJet(bAnalysisCut_HasSDD,jetflavour);
3070  return kFALSE;
3071  }
3072  }
3073 
3074  //n=hits in ITS layers
3075  Int_t SPDSSDHits = (int) track->HasPointOnITSLayer(0) + (int) track->HasPointOnITSLayer(1) + (int) track->HasPointOnITSLayer(4) + (int) track->HasPointOnITSLayer(5);
3076  if(SPDSSDHits<abs(fAnalysisCuts[bAnalysisCut_MinITSLayersHit])){
3077  //printf("Throw away due flav=%i ssd points 0 = %i, 1=%i, 2=%i, 3=%i, SPDSSDHits=%i cutvalue=%f\n",jetflavour,track->HasPointOnITSLayer(0),track->HasPointOnITSLayer(1),track->HasPointOnITSLayer(4),track->HasPointOnITSLayer(5),SPDSSDHits, abs(fAnalysisCuts[bAnalysisCut_MinITSLayersHit]));
3078  FillCandidateJet(bAnalysisCut_MinITSLayersHit,jetflavour);
3079  return kFALSE;
3080  }
3081 
3082  //TPC clusters
3083  if(((AliAODTrack*)track)->GetNcls(1)<fAnalysisCuts[bAnalysisCut_MinTPCClus]){
3084  //printf("Throw away due flav=%i TPCClus %i, cutvalue=%f\n",jetflavour,((AliAODTrack*)track)->GetNcls(1),fAnalysisCuts[bAnalysisCut_MinTPCClus]);
3085  FillCandidateJet(bAnalysisCut_MinTPCClus,jetflavour);
3086  return kFALSE;
3087  }
3088 
3089  if(track->Pt()<fAnalysisCuts[bAnalysisCut_MinTrackPt]){
3090  //printf("Throw away due flav=%i pt %f, cutvalue=%f\n",jetflavour,track->Pt(),fAnalysisCuts[bAnalysisCut_MinTrackPt]);
3091  FillCandidateJet(bAnalysisCut_MinTrackPt,jetflavour);
3092  return kFALSE;
3093  }
3094 
3095  AliAODVertex *aodvertex = (( AliAODTrack *)track)->GetProdVertex();
3096  if(!aodvertex) return kFALSE;
3097  if(fAnalysisCuts[bAnalysisCut_KinkCand]){
3098  if(aodvertex->GetType()==AliAODVertex::kKink){
3099  FillCandidateJet(bAnalysisCut_KinkCand,jetflavour);
3100  return kFALSE;
3101  }
3102  }
3103 
3104  ULong_t status=track->GetStatus();
3105  if(fAnalysisCuts[bAnalysisCut_HasTPCrefit]){
3106  if(!(status & AliAODTrack::kTPCrefit)){
3107  FillCandidateJet(bAnalysisCut_HasTPCrefit,jetflavour);
3108  return kFALSE;
3109  }
3110  }
3111 
3112  if(fAnalysisCuts[bAnalysisCut_HasITSrefit]){
3113  if(!(status & AliAODTrack::kITSrefit)){
3114  FillCandidateJet(bAnalysisCut_HasITSrefit,jetflavour);
3115  return kFALSE;
3116  }
3117  }
3118  if(((AliAODTrack*)track)->Chi2perNDF()>=fAnalysisCuts[bAnalysisCut_MinTrackChi2]){
3119  FillCandidateJet(bAnalysisCut_MinTrackChi2,jetflavour);
3120  //printf("Throw away due flav=%i chi2 %f, cutvalue=%f\n",jetflavour,((AliAODTrack*)track)->Chi2perNDF(),fAnalysisCuts[bAnalysisCut_MinTrackChi2]);
3121  return kFALSE;
3122  }
3123 
3124  return kTRUE;
3125  }
3126  return kTRUE;
3127 }
3128 
3129 
3130 bool AliAnalysisTaskHFJetIPQA::IsDCAAccepted(double decaylength, double ipwrtjet, Double_t * dca, int jetflavour){
3131  if(dca[0]>fAnalysisCuts[bAnalysisCut_MaxDCA_XY]){
3132  FillCandidateJet(bAnalysisCut_MaxDCA_XY,jetflavour);
3133  //printf("Throw away due to xy cut = %f, cutvalue=%f",dca[0],fAnalysisCuts[bAnalysisCut_MaxDCA_XY]);
3134  return kFALSE;
3135  }
3136  if(dca[1]>fAnalysisCuts[bAnalysisCut_MaxDCA_Z]){
3137  FillCandidateJet(bAnalysisCut_MaxDCA_Z,jetflavour);
3138  //printf("Throw away due to z cut = %f, cutvalue=%f",dca[1],fAnalysisCuts[bAnalysisCut_MaxDCA_Z]);
3139  return kFALSE;
3140  }
3141  if(decaylength>fAnalysisCuts[bAnalysisCut_MaxDecayLength]){
3142  FillCandidateJet(bAnalysisCut_MaxDecayLength,jetflavour);
3143  //printf("Throw away due to decaylength cut = %f, cutvalue=%f",decaylength,fAnalysisCuts[bAnalysisCut_MaxDecayLength]);
3144  return kFALSE;
3145  }
3146  if(ipwrtjet>fAnalysisCuts[bAnalysisCut_DCAJetTrack]){
3147  FillCandidateJet(bAnalysisCut_DCAJetTrack,jetflavour);
3148  //printf("Throw away due to dcajettrack = %f, cutvalue=%f",ipwrtjet,fAnalysisCuts[bAnalysisCut_DCAJetTrack]);
3149  return kFALSE;
3150  }
3151 
3152  return kTRUE;
3153 }
3154 
3166 Bool_t AliAnalysisTaskHFJetIPQA::MatchJetsGeometricDefault()
3167 {
3168  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
3169  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
3170  Double_t matchingpar1 =0.3;
3171  Double_t matchingpar2 =0.3;
3172  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return kFALSE;
3173  DoJetLoop();
3174  AliEmcalJet* jet1 = 0;
3175  jets1->ResetCurrentID();
3176  while ((jet1 = jets1->GetNextJet())) {
3177  AliEmcalJet *jet2 = jet1->ClosestJet();
3178  if (!jet2) continue;
3179  if (jet2->ClosestJet() != jet1) continue;
3180  if (jet1->ClosestJetDistance() > matchingpar1 || jet2->ClosestJetDistance() > matchingpar2) continue;
3181  // Matched jet found
3182  jet1->SetMatchedToClosest(1);
3183  jet2->SetMatchedToClosest(1);
3184  }
3185  return kTRUE;
3186 }
3187 
3195 void AliAnalysisTaskHFJetIPQA::DoJetLoop()
3196 {
3197  // Do the jet loop.
3198  Double_t minjetpt =1.;
3199  AliJetContainer *jets1 = static_cast<AliJetContainer*>(fJetCollArray.At(0));
3200  AliJetContainer *jets2 = static_cast<AliJetContainer*>(fJetCollArray.At(1));
3201  if (!jets1 || !jets1->GetArray() || !jets2 || !jets2->GetArray()) return;
3202  AliEmcalJet* jet1 = 0;
3203  AliEmcalJet* jet2 = 0;
3204  jets2->ResetCurrentID();
3205  while ((jet2 = jets2->GetNextJet())) jet2->ResetMatching();
3206  jets1->ResetCurrentID();
3207  while ((jet1 = jets1->GetNextJet())) {
3208  jet1->ResetMatching();
3209  if (jet1->MCPt() < minjetpt) continue;
3210  jets2->ResetCurrentID();
3211  while ((jet2 = jets2->GetNextJet())) {
3212  SetMatchingLevel(jet1, jet2, 1);
3213  } // jet2 loop
3214  } // jet1 loop
3215  }
3221  Bool_t AliAnalysisTaskHFJetIPQA::IsTruePrimary(AliVParticle * mcpart){
3222  if(!mcpart) return kFALSE;
3223  AliVParticle * mcmother = GetVParticleMother(mcpart);
3224  if(!mcmother) return kTRUE;
3225  Int_t istatus =-1;
3226  if(!fIsEsd) {
3227  istatus = ( (AliMCParticle*)mcmother)->Particle()->GetStatusCode();
3228  }
3229  else {
3230  istatus = ((AliAODMCParticle*)mcmother)->GetStatus();
3231  }
3232  if(istatus >11)return kTRUE;
3233  return kFALSE;
3234  }
3241  Double_t AliAnalysisTaskHFJetIPQA::GetWeightFactor( AliVTrack * track,Int_t &pCorr_indx, double &ppt){
3242  AliMCParticle *pMCESD = nullptr;
3243  AliAODMCParticle *pMCAOD = nullptr;
3244  if(track->GetLabel()< 0) return 1;
3245  if(!fIsPythia) return 1;
3246  if(fIsEsd){
3247  pMCESD = ((AliMCParticle*)MCEvent()->GetTrack(abs(track->GetLabel())));
3248  if(!(pMCESD)) return 1;
3249  }
3250  else {
3251  pMCAOD = static_cast<AliAODMCParticle*>(fMCArray->At(abs(track->GetLabel())));
3252  if(!(pMCAOD)) return 1;
3253  }
3254 
3255  AliVParticle * mcpart = fIsEsd ? ( AliVParticle * )pMCESD:( AliVParticle * )pMCAOD;
3256  Bool_t _particlesourcefound(kFALSE);
3257  Int_t _particlesourcepdg(mcpart->PdgCode());
3258  Int_t _particlesourceidx(25);
3259  Double_t _particlesourcept(0);
3260 
3261  AliVParticle * mcpartclone = mcpart;
3262  while(mcpart){//Strangenss
3263  if((abs(mcpart->PdgCode()) >0 && abs(mcpart->PdgCode()) <7)|| (abs(mcpart->PdgCode()) == 21)) break;
3264  _particlesourcept = mcpart->Pt();
3265  _particlesourcepdg = abs(mcpart->PdgCode());
3266  if (IsSelectionParticleStrange(mcpart,_particlesourcepdg,_particlesourcept,_particlesourceidx)){
3267  _particlesourcefound = kTRUE;
3268  break;
3269  }
3270  mcpart = GetVParticleMother(mcpart);
3271  }
3272  if (!_particlesourcefound) { //heavy mesons to improve templates
3273  mcpart = mcpartclone;
3274 
3275  while(mcpart){
3276  if((abs(mcpart->PdgCode()) >0 && abs(mcpart->PdgCode()) <7)|| (abs(mcpart->PdgCode()) == 21)) break;
3277  _particlesourcept = mcpart->Pt();
3278  _particlesourcepdg = abs(mcpart->PdgCode());
3279  if (IsSelectionParticleMeson(mcpart,_particlesourcepdg,_particlesourcept,_particlesourceidx)){
3280  _particlesourcefound = kTRUE;
3281  break;
3282  }
3283  mcpart = GetVParticleMother(mcpart);
3284  }
3285  }
3286  if (!_particlesourcefound) { //charged hadrons
3287  mcpart = mcpartclone;
3288  while(mcpart){
3289  if((abs(mcpart->PdgCode()) >0 && abs(mcpart->PdgCode()) <7)|| (abs(mcpart->PdgCode()) == 21)) break;
3290  _particlesourcept = mcpart->Pt();
3291  _particlesourcepdg = abs(mcpart->PdgCode());
3292  if (IsSelectionParticleOmegaXiSigmaP(mcpart,_particlesourcepdg,_particlesourcept,_particlesourceidx)){
3293  _particlesourcefound = kTRUE;
3294  break;
3295  }
3296  mcpart = GetVParticleMother(mcpart);
3297  }
3298  }
3299  if (!_particlesourcefound) {
3300  mcpart = mcpartclone;
3301  while(mcpart){
3302  if((abs(mcpart->PdgCode()) >0 && abs(mcpart->PdgCode()) <7)|| (abs(mcpart->PdgCode()) == 21)) break;
3303  _particlesourcept = mcpart->Pt();
3304  _particlesourcepdg = abs(mcpart->PdgCode());
3305  if (IsSelectionParticle(mcpart,_particlesourcepdg,_particlesourcept,_particlesourceidx)){
3306  _particlesourcefound = kTRUE;
3307  break;
3308  }
3309  mcpart = GetVParticleMother(mcpart);
3310  }
3311  }
3312  if (!_particlesourcefound) {
3313  mcpart = mcpartclone;
3314  while(mcpart){
3315  if((abs(mcpart->PdgCode()) >0 && abs(mcpart->PdgCode()) <7)|| (abs(mcpart->PdgCode()) == 21)) break;
3316  _particlesourcept = mcpart->Pt();
3317  if (IsSelectionParticleALICE(mcpart,_particlesourcepdg,_particlesourcept,_particlesourceidx)){
3318  _particlesourcefound = kTRUE;
3319  break;
3320  }
3321  mcpart = GetVParticleMother(mcpart);
3322  }
3323  }
3324 
3325  if (!_particlesourcefound) return 1.;
3326  // Do the weighting
3327  Double_t factor = 1;
3328  if (_particlesourceidx <0) return 1;
3329  if (_particlesourceidx >19) return 1;
3330  Double_t wpos = ((_particlesourcept - 0.15)/ 0.05);
3331  Double_t fractpart, intpart;
3332  fractpart = modf (wpos , &intpart);
3333  if (fractpart > 0) intpart = intpart + 1;
3334  Int_t bin = floor(intpart);
3335  if (bin > 497) bin = 497;// above weight definition
3336  if (_particlesourcept < 0.1+ 1E-5) bin = 0; //below weight definition
3337  if((_particlesourceidx == bIdxSigmaMinus) || (_particlesourceidx == bIdxSigmaPlus)) factor = fBackgroundFactorLinus[bIdxLambda][bin];
3338  else factor = fBackgroundFactorLinus[_particlesourceidx][bin];
3339  pCorr_indx = _particlesourceidx;
3340  Double_t flucafactor = 1;
3341 
3342  switch(mcpart->PdgCode())
3343  {
3344  case -bPhi:
3345  factor=1;
3346  flucafactor =1./fPhi->Eval(_particlesourcept);
3347  pCorr_indx=bIdxPhi;
3348  break;
3349  case -bK0S892:
3350  factor=1;
3351  flucafactor =1./fK0Star->Eval(_particlesourcept);
3352  pCorr_indx=bIdxK0S892;
3353  break;
3354  case -bK0S892plus:
3355  factor=1;
3356  flucafactor =1./fK0Star->Eval(_particlesourcept);
3357  pCorr_indx=bK0S892plus;
3358  break;
3359  case -bOmegaBaryon:
3360  pCorr_indx=bIdxOmegaBaryon;
3361  factor=1;
3362  flucafactor =1./fGraphOmega->Eval(_particlesourcept);
3363  break;
3364  case -bXiBaryon:
3365  pCorr_indx=bIdxXiBaryon;
3366  factor=1;
3367  flucafactor =1./fGraphXi->Eval(_particlesourcept);
3368  break;
3369  case bPhi:
3370  factor=1;
3371  flucafactor =1./fPhi->Eval(_particlesourcept);
3372  pCorr_indx=bIdxPhi;
3373  break;
3374  case bK0S892:
3375  factor=1;
3376  pCorr_indx=bIdxK0S892;
3377  flucafactor =1./fK0Star->Eval(_particlesourcept);
3378  break;
3379  case bK0S892plus:
3380  pCorr_indx=bK0S892plus;
3381  factor=1;
3382  flucafactor =1./fK0Star->Eval(_particlesourcept);
3383  break;
3384  case bOmegaBaryon:
3385  pCorr_indx=bIdxOmegaBaryon;
3386  factor=1;
3387  flucafactor =fGraphOmega->Eval(_particlesourcept);
3388  break;
3389  case bXiBaryon:
3390  pCorr_indx=bIdxXiBaryon;
3391  factor=1;
3392  flucafactor =fGraphXi->Eval(_particlesourcept);
3393  break;
3394 
3395  default:
3396  break;
3397  }
3398  factor*=flucafactor;
3399  if (factor <= 0 || factor > 100.) {
3400  return 1;
3401  }
3402  ppt = _particlesourcept;
3403  return factor ;
3404 }
3405 
3410 Bool_t AliAnalysisTaskHFJetIPQA::GetBMesonWeight( AliVParticle * mcpart ,Int_t &pdg,Double_t &pT,Int_t &idx )
3411 {
3412  pT = mcpart->Pt();
3413  switch(pdg){
3414  case bBPlus:
3415  idx = bIdxBPlus;
3416  return kTRUE;
3417  case bB0:
3418  idx = bIdxB0;
3419  return kTRUE;
3420  case bLambdaB:
3421  idx = bIdxLambdaB;
3422  return kTRUE;
3423  break;
3424  case bBStarPlus:
3425  idx = bIdxBStarPlus;
3426  return kTRUE;
3427  break;
3428  default:
3429  break;
3430  }
3431  return kFALSE;
3432 }
3437 Bool_t AliAnalysisTaskHFJetIPQA::IsSelectionParticle( AliVParticle * mcpart ,Int_t &pdg,Double_t &pT,Int_t &idx ){
3438  pT = mcpart->Pt();
3439  Int_t pdg2 = abs(mcpart->PdgCode());
3440  idx = -1;
3441 
3442  switch(pdg2){
3443  case bPi0:
3444  idx = bIdxPi0;
3445  if(!IsSecondaryFromWeakDecay(mcpart))return kTRUE;
3446  break;
3447  case bEta:
3448  idx = bIdxEta;
3449  if(!IsSecondaryFromWeakDecay(mcpart))return kTRUE;
3450  break;
3451  case bEtaPrime:
3452  idx = bIdxEtaPrime;
3453  if(!IsSecondaryFromWeakDecay(mcpart))return kTRUE;
3454  break;
3455  case bOmega:
3456  idx = bIdxOmega;
3457  if(!IsSecondaryFromWeakDecay(mcpart))return kTRUE;
3458  break;
3459  case bPhi:
3460  idx = bIdxPhi;
3461  if(!IsSecondaryFromWeakDecay(mcpart))return kTRUE;
3462  break;
3463  case bRho:
3464  idx = bIdxRho;
3465  if(!IsSecondaryFromWeakDecay(mcpart))return kTRUE;
3466  break;
3467  case bRhoPlus:
3468  idx = bIdxRho;
3469  if(!IsSecondaryFromWeakDecay(mcpart))return kTRUE;
3470  break;
3471  default:
3472  break;
3473  }
3474  return kFALSE;
3475 }
3480 Bool_t AliAnalysisTaskHFJetIPQA::IsSelectionParticleALICE( AliVParticle * mcpart ,Int_t &pdg,Double_t &pT,Int_t &idx ){
3481  pT = mcpart->Pt();
3482  AliVParticle * mother = nullptr;
3483  idx = -1;
3484  pdg = abs(mcpart->PdgCode());
3485  Bool_t pIsSecStrANGE= IsSecondaryFromWeakDecay(mcpart);
3486  if(pIsSecStrANGE ) return kFALSE;
3487  if(!IsPhysicalPrimary(mcpart))return kFALSE;
3488 
3489  switch(pdg){
3490  case bProton:
3491  mother = GetVParticleMother(mcpart);
3492  if(mother){
3493  if((abs(mother->PdgCode()) == 3222) )
3494  return kFALSE;
3495  }
3496  idx=bIdxProton;
3497  return kTRUE;
3498  break;
3499  case bPi:
3500  idx=bIdxPi;
3501  return kTRUE;
3502  break;
3503  case bKaon:
3504  idx=bIdxKaon;
3505  return kTRUE;
3506  break;
3507  default:
3508  return kFALSE;
3509  break;
3510  }
3511  return kFALSE;
3512 }
3517 Bool_t AliAnalysisTaskHFJetIPQA::IsSelectionParticleMeson( AliVParticle * mcpart ,Int_t &pdg,Double_t &pT,Int_t &idx ){
3518  pT = mcpart->Pt();
3519  idx = -1;
3520  pdg = abs(mcpart->PdgCode());
3521  if(TMath::Abs(mcpart->Y()) >=.5) return kFALSE;
3522  if(IsSecondaryFromWeakDecay(mcpart))return kFALSE;
3523  switch(pdg){
3524  case bD0:
3525  idx = bIdxD0;
3526  if(IsPromptDMeson(mcpart))return kTRUE;
3527  break;
3528  case bDPlus:
3529  idx = bIdxDPlus;
3530  if(IsPromptDMeson(mcpart))return kTRUE;
3531  break;
3532  case bDSPlus:
3533  idx = bIdxDSPlus;
3534  if(IsPromptDMeson(mcpart))return kTRUE;
3535  break;
3536  case bDStarPlus:
3537  idx = bIdxDStarPlus;
3538  if(IsPromptDMeson(mcpart))return kTRUE;
3539  break;
3540  case bLambdaC:
3541  idx = bIdxLambdaC;
3542  if(IsPromptDMeson(mcpart))return kTRUE;
3543  break;
3544  case bBPlus:
3545  idx = bIdxBPlus;
3546  if(IsPromptBMeson(mcpart))return kTRUE;
3547  break;
3548  case bB0:
3549  idx = bIdxB0;
3550  if(IsPromptBMeson(mcpart))return kTRUE;
3551  break;
3552  case bLambdaB:
3553  idx = bIdxLambdaB;
3554  if(IsPromptBMeson(mcpart))return kTRUE;
3555  break;
3556  case bBStarPlus:
3557  idx = bIdxBStarPlus;
3558  if(IsPromptBMeson(mcpart))return kTRUE;
3559  break;
3560  default:
3561  return kFALSE;
3562  break;
3563  }
3564  return kTRUE;
3565 }
3570 Bool_t AliAnalysisTaskHFJetIPQA::IsSelectionParticleOmegaXiSigmaP( AliVParticle * mcpart ,Int_t &pdg,Double_t &pT,Int_t &idx ){
3571  pT = mcpart->Pt();
3572  idx = -1;
3573  pdg = abs(mcpart->PdgCode());
3574 
3575  if (!IsPhysicalPrimary(mcpart)) return kFALSE;
3576  switch(pdg){
3577  case bSigmaMinus:
3578  idx = bIdxSigmaMinus; // use lambda as proxy
3579  return kTRUE;
3580  break;
3581  case bSigmaPlus:
3582  idx = bIdxSigmaPlus; //use lambda as proxy
3583  return kTRUE;
3584  break;
3585  case bXiBaryon:
3586  idx = bIdxBStarPlus;//dummy! Externally sotlved
3587  return kTRUE;
3588  break;
3589  case bOmegaBaryon:
3590  idx = bIdxBStarPlus;//dummy! Externally solved
3591  return kTRUE;
3592  break;
3593  default:
3594  return kFALSE;
3595  break;
3596  }
3597  return kFALSE;
3598  }
3603  AliVParticle * AliAnalysisTaskHFJetIPQA::GetVParticleMother(AliVParticle * part){
3604  AliVParticle * mother = nullptr;
3605  if (part->GetMother()<0) return nullptr;
3606  if(fIsEsd){
3607  mother = static_cast <AliMCParticle*>(MCEvent()->GetTrack(part->GetMother()));
3608  }
3609  else{
3610  mother =static_cast<AliAODMCParticle*>(fMCArray->At(part->GetMother()));
3611  }
3612  if(!mother) return nullptr;
3613  return mother;
3614  }
3619  Bool_t AliAnalysisTaskHFJetIPQA::IsPhysicalPrimary( AliVParticle * part){
3620  return fIsEsd ? MCEvent()->IsPhysicalPrimary(part->GetLabel()) : static_cast<AliAODMCParticle*>(part)->IsPhysicalPrimary() ;
3621  }
3626  Bool_t AliAnalysisTaskHFJetIPQA::IsSecondaryFromWeakDecay( AliVParticle * part){
3627  return fIsEsd ? MCEvent()->IsSecondaryFromWeakDecay(part->GetLabel()) : static_cast<AliAODMCParticle*>(part)->IsSecondaryFromWeakDecay() ;
3628 
3629  }
3634  Bool_t AliAnalysisTaskHFJetIPQA::IsSelectionParticleStrange( AliVParticle * mcpart ,Int_t &pdg,Double_t &pT,Int_t &idx ){
3635  pT = mcpart->Pt();
3636  AliVParticle * mother = nullptr;
3637  idx = -1;
3638  pdg = abs(mcpart->PdgCode());
3639  if(!IsPhysicalPrimary(mcpart))return kFALSE;
3640  switch(pdg){
3641  case bPhi:
3642  idx = bIdxPhi;//dummy will be overwritten in posrpos
3643  return kTRUE;
3644  break;
3645  case bK0S892:
3646  idx = bIdxK0s;//dummy will be overwritten in posrpos
3647  return kTRUE;
3648  break;
3649  case bK0S892plus:
3650  idx = bIdxK0s; //dummy will be overwritten in posrpos
3651  return kTRUE;
3652  break;
3653  case bK0s:
3654  idx = bIdxK0s;
3655  mother = GetVParticleMother(mcpart);
3656  if(mother){
3657  if((abs(mother->PdgCode()) == bPhi))
3658  return kFALSE;
3659  }
3660  return kTRUE;
3661  break;
3662  case bK0l:
3663  idx = bIdxK0s;
3664  mother = GetVParticleMother(mcpart);
3665  if(mother){
3666  if((abs(mother->PdgCode()) == bPhi))
3667  return kFALSE;
3668  }
3669  return kTRUE;
3670  break;
3671 
3672  case bLambda:
3673  idx = bIdxLambda;
3674  mother = GetVParticleMother(mcpart);
3675  if(mother){
3676  if((abs(mother->PdgCode()) == 3312) || (abs(mother->PdgCode()) == 3322) || (abs(mother->PdgCode()) == 3334))
3677  return kFALSE;
3678  }
3679  return kTRUE;
3680  break;
3681  default:
3682  return kFALSE;
3683  break;
3684  }
3685  return kFALSE;
3686  }
3691  Bool_t AliAnalysisTaskHFJetIPQA::IsPromptBMeson(AliVParticle * part )
3692  {
3693  if(!part) return kFALSE;
3694  Int_t pdg = TMath::Abs(part->PdgCode());
3695  if ((pdg >= 500 && pdg < 600 )||(pdg >= 5000 && pdg < 6000 ))
3696  {
3697  AliVParticle* pm = GetVParticleMother(part);
3698  Int_t mpdg = TMath::Abs(pm->PdgCode());
3699  if (!(mpdg >5000 && mpdg <6000) && !(mpdg >500 && mpdg <600))
3700  return kTRUE;
3701  }
3702  return kFALSE;
3703  }
3708  Bool_t AliAnalysisTaskHFJetIPQA::IsPromptDMeson(AliVParticle * part )
3709  {
3710  if(!part) return kFALSE;
3711  Int_t pdg = TMath::Abs(part->PdgCode());
3712  if ((pdg >= 400 && pdg < 500 )||(pdg >= 4000 && pdg < 5000 ))
3713  {
3714  AliVParticle* pm = GetVParticleMother(part);
3715  if(!pm) return kTRUE;
3716  Int_t mpdg = TMath::Abs(pm->PdgCode());
3717  if (!(mpdg >4000 && mpdg <6000) && !(mpdg >400 && mpdg <600))
3718  return kTRUE;
3719  }
3720 
3721  return kFALSE;
3722  }
3727  Bool_t AliAnalysisTaskHFJetIPQA::ParticleIsPossibleSource(Int_t pdg){
3728  Int_t pos[22] = {bPi0,bEta,bEtaPrime,bPhi,bRho,bOmega,bK0s,bLambda,bOmegaBaryon,bXiBaryon,bD0,bPi,bKaon,bProton,bDPlus,bDStarPlus,bDSPlus,bLambdaB,bLambdaC,bBPlus,bB0,bBStarPlus};
3729  for (Int_t i =0 ;i<22 ;++i){
3730  if (abs(pdg)==pos[i] ) return kTRUE;
3731  }
3732  return kFALSE;
3733  }
3739  void AliAnalysisTaskHFJetIPQA::SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Int_t matching)
3740  {
3741  Double_t d1 = - 1;
3742  Double_t d2 = -1;
3743 
3744  switch (matching) {
3745  case 1:
3746  GetGeometricalMatchingLevel(jet1,jet2,d1);
3747  d2 = d1;
3748  break;
3749  default:
3750  break;
3751  }
3752  if (d1 >= 0) {
3753  if (d1 < jet1->ClosestJetDistance()) {
3754  jet1->SetSecondClosestJet(jet1->ClosestJet(), jet1->ClosestJetDistance());
3755  jet1->SetClosestJet(jet2, d1);
3756  }
3757  else if (d1 < jet1->SecondClosestJetDistance()) {
3758  jet1->SetSecondClosestJet(jet2, d1);
3759  }
3760  }
3761  if (d2 >= 0) {
3762  if (d2 < jet2->ClosestJetDistance()) {
3763  jet2->SetSecondClosestJet(jet2->ClosestJet(), jet2->ClosestJetDistance());
3764  jet2->SetClosestJet(jet1, d2);
3765  }
3766  else if (d2 < jet2->SecondClosestJetDistance()) {
3767  jet2->SetSecondClosestJet(jet1, d2);
3768  }
3769  }
3770  }
3775  void AliAnalysisTaskHFJetIPQA::GetGeometricalMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Double_t &d) const
3776  {
3777  Double_t deta = jet2->Eta() - jet1->Eta();
3778  Double_t dphi = jet2->Phi() - jet1->Phi();
3779  dphi = TVector2::Phi_mpi_pi(dphi);
3780  d = sqrt(deta * deta + dphi * dphi);
3781  }
3786  Double_t AliAnalysisTaskHFJetIPQA::GetMonteCarloCorrectionFactor(AliVTrack* track,Int_t &pCorr_indx, double &ppt){
3787  printf("Doing MC Correction.\n");
3788  double val= GetWeightFactor(track,pCorr_indx,ppt);
3789  if(val > 0 ) return val;
3790  return 1.;
3791  }
3796  Bool_t AliAnalysisTaskHFJetIPQA::mysort(const SJetIpPati& i, const SJetIpPati& j)
3797  {
3798  if(i.first <= j.first)
3799  return kFALSE;
3800  else
3801  return kTRUE;
3802  }
3807  Double_t AliAnalysisTaskHFJetIPQA::GetPtCorrected(const AliEmcalJet *jet)
3808  {
3809  AliJetContainer * jetconrec = nullptr;
3810  jetconrec = static_cast<AliJetContainer*>(fJetCollArray.At(0));
3811  if(jet && jetconrec&&fDoUnderlyingEventSub){
3812  printf("Correct for Underlying Event.\n");
3813  return jet->Pt() - jetconrec->GetRhoVal() * jet->Area();
3814  }
3815  return -1.;
3816  }
3821  Double_t AliAnalysisTaskHFJetIPQA::GetPtCorrectedMC(const AliEmcalJet *jet)
3822  {
3823  AliJetContainer * jetconrec = nullptr;
3824  jetconrec = static_cast<AliJetContainer*>(fJetCollArray.At(1));
3825  if(jet && jetconrec&&fDoUnderlyingEventSub){
3826  printf("Correct for Underlying Event.\n");
3827  return jet->Pt() - jetconrec->GetRhoVal() * jet->Area();
3828  }
3829  return -1.;
3830  }
3831 
3832 
3837  Bool_t AliAnalysisTaskHFJetIPQA::IsJetTaggedTC(Int_t n, Double_t thres)
3838  {
3839  return kTRUE;
3840  }
3841 
3842 
3847  Bool_t AliAnalysisTaskHFJetIPQA::IsParton(int pdg){
3848  return ((pdg==1)||(pdg==2)||(pdg==3)||(pdg==4)||(pdg==5)||(pdg==21));
3849  }
3850 
3855  Bool_t AliAnalysisTaskHFJetIPQA::IsJetTaggedJetProb(Double_t thresProb)
3856  {
3857  return kTRUE;
3858  }
3863  Int_t AliAnalysisTaskHFJetIPQA::IsMCJetPartonFast(const AliEmcalJet *jet, Double_t radius,Bool_t &is_udg)
3864  {
3865  fJetCont.clear();
3866  fPUdsgJet.clear();
3867  fPSJet.clear();
3868  fPCJet.clear();
3869  fPBJet.clear();
3870  daughtermother.clear();
3871 
3872  double p_udsg_max=-999;
3873  double p_s_max=-999;
3874  AliVParticle *vp=0x0;
3875  Int_t pdg=0;
3876  Double_t p=0;
3877  int kJetOrigin=-999;
3878 
3879  //printf("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n");
3880  //printf("Starting loop\n");
3881  if(!jet) return 0;
3882  if(!(jet->GetNumberOfTracks()>fNoJetConstituents)){
3883  //printf("Throwing away jets with only too few constituent!\n");
3884  return 0;
3885  }
3886 
3887  if(fDoFlavourMatching){
3888  for(UInt_t i = 0; i < jet->GetNumberOfTracks(); i++) {//start trackloop jet
3889  vp = static_cast<AliVParticle*>(jet->Track(i));
3890  if (!vp){
3891  AliError("AliVParticle associated to constituent not found\n");
3892  continue;
3893  }
3894 
3895  AliAODMCParticle * part = static_cast<AliAODMCParticle*>(vp);
3896 
3897  if(!part){
3898  AliError("Finding no Part!\n");
3899  return 0;
3900  } // if(!part->IsPrimary()) continue;
3901  pdg = (abs(part->PdgCode()));
3902 
3903  fJetCont.push_back(part->Label());
3904  //printf("Daugther pdg=%i, Label=%i, Mother =%i, p=%f, MCStatusCode=%i\n",pdg, part->GetLabel(), part->GetMother(), p, part->MCStatusCode());
3905  }//end trackloop jet
3906  }
3907 
3908  for(Int_t iPrim = 0 ; iPrim<fMCArray->GetEntriesFast();iPrim++){//start trackloop MC
3909 
3910  AliAODMCParticle * part = static_cast<AliAODMCParticle*>(fMCArray->At(iPrim));
3911  if(!part) return 0;
3912  if(!part->IsPrimary()) continue;
3913  Double_t eta = part->Eta();
3914  Double_t phi= part->Phi();
3915  Double_t etajet = jet->Eta();
3916  Double_t phijet = jet->Phi();
3917  p=part->P();
3918 
3919  Int_t pdg = (abs(part->PdgCode()));
3920  Double_t deta = etajet - eta;
3921  Double_t dphi = phijet-phi;
3922  dphi = TVector2::Phi_mpi_pi(dphi);
3923  Double_t d = sqrt(deta * deta + dphi * dphi);
3924  // printf("LINE %i, deta%f, dphi%f, d=%f, fDoFlavourMatching=%i\n",__LINE__, deta,dphi,d,fDoFlavourMatching);
3925  if(!fDoFlavourMatching) {
3926  //if(!((part->GetStatus()==11) ||(part->GetStatus()==12))) continue;
3927  if(!IsParton(pdg)) continue;
3928  if(d > radius) continue;
3929  kJetOrigin=pdg;
3930  }
3931  else{
3932 
3933  if(IsParton(pdg)){
3934  if(d > fDaughtersRadius) continue;
3935  }
3936 
3937  //printf("i=%i, Mother Label %i, Mother PdG %i,Daughter:%i, Last Daughter: %i, MCStatusCode=%i, d=%f, p=%f\n", iPrim, part->Label(), part->PdgCode(), part->GetDaughterLabel(0),part->GetDaughterLabel(1),part->MCStatusCode(), d,p);
3938 
3939  int kFirstDaugh=part->GetDaughterLabel(0);
3940  int NDaugh=part->GetNDaughters();
3941  for(int idaugh=0;idaugh<NDaugh;idaugh++){
3942  int kDaughLabel=kFirstDaugh+idaugh;
3943  //printf("Dauglabel=%i, kFirstDaugh=%i, kLastDaugh=%i\n",kDaughLabel, kFirstDaugh,kLastDaugh);
3944 
3945  bool IsDaughter=std::find(fJetCont.begin(), fJetCont.end(),kDaughLabel) != fJetCont.end();
3946  if(IsDaughter){
3947  if(IsParton(pdg)){
3948  //printf("Directly matched %i with daughter =%i\n",part->GetLabel(), kDaughLabel);
3949  kJetOrigin=part->PdgCode();
3950  }
3951  else{
3952  bool Is2ndDaughter=daughtermother.find(part->Label()) != daughtermother.end();
3953  if(Is2ndDaughter){
3954  kJetOrigin=daughtermother.find(part->Label())->second;
3955  //printf("Matched with %i with 2nd daughter =%i\n",part->GetLabel(), kDaughLabel);
3956  }
3957  }
3958  }//end is jet daughter
3959  else{
3960  if(IsParton(pdg)){
3961  //printf("Writing Quarks in map: partlabel=%i, daughlabel=%i, status=%i\n", part->GetLabel(),kDaughLabel, part->MCStatusCode());
3962  daughtermother.emplace(kDaughLabel, part->PdgCode());
3963  }
3964  else{
3965  bool Is2ndDaughter=daughtermother.find(part->Label()) != daughtermother.end();
3966  if(Is2ndDaughter){
3967  //printf("Writing Daughters in map: partlabel=%i, daughlabel=%i, status=%i\n", part->GetLabel(),kDaughLabel, part->MCStatusCode());
3968  int kOriginalQuark=daughtermother.find(part->Label())->second;
3969  daughtermother.emplace(kDaughLabel, kOriginalQuark);
3970  }
3971  //printf("Asking for daughlabel%i\n",part->Label());
3972  }
3973  }//end related to jet?
3974  }//end daughterloop
3975  }//end else DoMatchFlavours
3976 
3977  //printf("i=%i, Mother Label %i, Mother PdG %i,Daughter:%i, Last Daughter: %i, MCStatusCode=%i, d=%f, p=%f\n", iPrim, part->Label(), part->PdgCode(), part->GetDaughterLabel(0),part->GetLastDaughter(),part->MCStatusCode(), d,p);
3978 
3979  if(abs(kJetOrigin) == 5) {
3980  fPBJet.push_back(p);
3981  }
3982  else if(abs(kJetOrigin)== 4) {
3983  fPCJet.push_back(p);
3984  }
3985  else if(abs(kJetOrigin) == 3 ) {
3986  fPSJet.push_back(p);
3987  //printf("Strange pushed with p=%f",p);
3988  }
3989  else if(abs(kJetOrigin)== 1 ||abs(kJetOrigin)== 2 || abs(kJetOrigin) == 21) {
3990  fPUdsgJet.push_back(p);
3991  //printf("Light pushed with p=%f",p);
3992  }
3993 
3994  }//end trackloop MC
3995 
3996  /*printf("Inside JetCont:\n");
3997  for(int i=0;i<fJetCont.size();i++){
3998  printf("%f\n", fJetCont[i]);
3999  }
4000  printf("Inside map:\n");
4001  for (auto& x: daughtermother) {
4002  std::cout << x.first << ": " << x.second << '\n';
4003  }*/
4004  if(fPCJet.size() ==0&& fPBJet.size()==0&& fPSJet.size()==0&&fPUdsgJet.size()==0) return 0; //udsg
4005  //check for b jet
4006  for (Int_t icj = 0 ; icj <(Int_t)fPBJet.size();++icj ){
4007  //printf("Bottom Flavour Jet!\n");
4008  return B;
4009  }
4010  //check for c jet
4011  for (Int_t icj = 0 ; icj <(Int_t)fPCJet.size();++icj ){
4012  //printf("Charm Flavour Jet!\n");
4013  return C;
4014  }
4015  //check for s and light jet
4016  if(fPUdsgJet.size()!=0){
4017  std::sort(fPUdsgJet.begin(), fPUdsgJet.end());
4018  p_udsg_max=fPUdsgJet[fPUdsgJet.size()-1];
4019  }
4020  if(fPSJet.size()!=0){
4021  std::sort(fPSJet.begin(), fPSJet.end());
4022  p_s_max=fPSJet[fPSJet.size()-1];
4023  }
4024 
4025  if(p_s_max>p_udsg_max){
4026  //printf("S prefered with psmax=%f, pudsgmax=%f\n", p_s_max,p_udsg_max);
4027  return UDSG;
4028  }
4029  else{
4030  if(fPUdsgJet.size()!=0){
4031  //printf("Light prefered with psmax=%f, pudsgmax=%f\n", p_s_max,p_udsg_max);
4032  is_udg =kTRUE;
4033  return UDSG;
4034  }
4035  }
4036  return 0;
4037  }
4038 
4039 
4044 //_________________________________________________________________________
4045 void AliAnalysisTaskHFJetIPQA::RecursiveParents(AliEmcalJet *fJet,AliJetContainer *fJetCont){
4046 
4047  std::vector<fastjet::PseudoJet> fInputVectors;
4048  fInputVectors.clear();
4049  fastjet::PseudoJet PseudoTracks;
4050 
4051  AliParticleContainer *fTrackCont = fJetCont->GetParticleContainer();
4052 
4053  if (fTrackCont) for (Int_t i=0; i<fJet->GetNumberOfTracks(); i++) {
4054  AliVParticle *fTrk = fJet->TrackAt(i, fTrackCont->GetArray());
4055  if (!fTrk) continue;
4056  //if(fDoTwoTrack==kTRUE && CheckClosePartner(i,fJet,fTrk,fTrackCont)) continue;
4057  PseudoTracks.reset(fTrk->Px(), fTrk->Py(), fTrk->Pz(),fTrk->E());
4058  PseudoTracks.set_user_index(fJet->TrackAt(i)+100);
4059  fInputVectors.push_back(PseudoTracks);
4060 
4061  }
4062  fastjet::JetAlgorithm jetalgo(fastjet::cambridge_algorithm);
4063 
4064 
4065 
4066  fastjet::JetDefinition fJetDef(jetalgo, 1., static_cast<fastjet::RecombinationScheme>(0), fastjet::BestFJ30 );
4067 
4068  try {
4069  fastjet::ClusterSequence fClustSeqSA(fInputVectors, fJetDef);
4070  std::vector<fastjet::PseudoJet> fOutputJets;
4071  fOutputJets.clear();
4072  fOutputJets=fClustSeqSA.inclusive_jets(0);
4073 
4074  fastjet::PseudoJet jj;
4075  fastjet::PseudoJet j1;
4076  fastjet::PseudoJet j2;
4077  jj=fOutputJets[0];
4078  double ktaverage=0;
4079  double thetaverage=0;
4080  double nall=0;
4081  double flagSubjet=0;
4082  while(jj.has_parents(j1,j2)){
4083  nall=nall+1;
4084  if(j1.perp() < j2.perp()) swap(j1,j2);
4085  flagSubjet=0;
4086  double delta_R=j1.delta_R(j2);
4087  double z=j2.perp()/(j1.perp()+j2.perp());
4088  double y =log(1.0/delta_R);
4089  double lnpt_rel=log(j2.perp()*delta_R);
4090  double yh=j1.e()+j2.e();
4091  vector < fastjet::PseudoJet > constitj1 = sorted_by_pt(j1.constituents());
4092  if(constitj1[0].perp()>fAnalysisCuts[bAnalysisCut_MinTrackPt]) flagSubjet=1;
4093  if(z>fHardCutOff){
4094  ktaverage=ktaverage+lnpt_rel;
4095  thetaverage=thetaverage+delta_R;
4096  Double_t LundEntries[6] = {y,lnpt_rel,fOutputJets[0].perp(),nall,yh,flagSubjet};
4097  fHLundIterative->Fill(LundEntries);}
4098  jj=j1;}
4099  } catch (fastjet::Error) {
4100  AliError(" [w] FJ Exception caught.");
4101  //return -1;
4102  }
4103  return;
4104 }
4109 void AliAnalysisTaskHFJetIPQA::FillHist(const char *name, Double_t x, Double_t w){
4110  TH1D * h1 =GetHist1D(name);
4111  if(h1) h1->Fill(x,w);
4112 }
4117 void AliAnalysisTaskHFJetIPQA::FillHist(const char *name, Double_t x, Double_t y, Double_t w){
4118  TH2D * h2 =GetHist2D(name);
4119  if(h2) h2->Fill(x,y,w);
4120 }
4125 void AliAnalysisTaskHFJetIPQA::IncHist(const char *name, Int_t bin){
4126  TH1D * h1 =GetHist1D(name);
4127  h1->SetBinContent(bin,h1->GetBinContent(bin)+1);
4128 }
4133 TH1 *AliAnalysisTaskHFJetIPQA::AddHistogramm(const char *name, const char *title, Int_t x, Double_t xlow, Double_t xhigh, Int_t y, Double_t ylow, Double_t yhigh){
4134  TObject * res = nullptr;
4135  res = fOutput->FindObject(name);
4136  if((res)) return nullptr;
4137 
4138  TH1 * phist=nullptr;
4139  if(y==0){ //TH1D*
4140  phist = new TH1D (name,title,x,xlow,xhigh);
4141  }
4142  else {
4143  phist = new TH2D(name,title,x,xlow,xhigh,y,ylow,yhigh);
4144  }
4145  phist->Sumw2();
4146 
4147  TString sName(name);
4148  if(sName.EqualTo("fh1dCutsPrinted")){
4149  fOutput->AddFirst(phist);
4150  }
4151  else{
4152  fOutput->Add(phist);
4153  }
4154  return (TH1*)phist;
4155 }
4156 
4157 void AliAnalysisTaskHFJetIPQA::setFFillCorrelations(const Bool_t &value)
4158 {
4159  fFillCorrelations = value;
4160 }
4161 
4162 
4163 
4164 void AliAnalysisTaskHFJetIPQA::setFMCglobalDCASmear(const Double_t value)
4165 {
4166  fMCglobalDCASmear = value;
4167 }
4168 
4169 Double_t AliAnalysisTaskHFJetIPQA::getFVertexRecalcMinPt() const
4170 {
4171  return fVertexRecalcMinPt;
4172 }
4173 
4174 void AliAnalysisTaskHFJetIPQA::setFVertexRecalcMinPt(const Double_t &value)
4175 {
4176  fVertexRecalcMinPt = value;
4177 }
4178 
4179 Double_t AliAnalysisTaskHFJetIPQA::getFMCglobalDCAxyShift() const
4180 {
4181  return fMCglobalDCAxyShift;
4182 }
4183 
4184 void AliAnalysisTaskHFJetIPQA::setFMCglobalDCAxyShift(const Double_t &value)
4185 {
4186  fMCglobalDCAxyShift = value;
4187 }
4188 
4189 
4190 
4191 
4192 
4199 //Currently not in use
4200 /*void AliAnalysisTaskHFJetIPQA::SubtractMean(Double_t val[], AliVTrack *track){
4201  Double_t deltamean=fGraphMean->Eval(track->Pt() < 3. ? track->Pt() : 3 );//(fCurrentMeanFactors[0]-fCurrentMeanFactors[1]*TMath::Exp(-1*fCurrentMeanFactors[2]*(track->Pt()-fCurrentMeanFactors[3]))) *1e-4;
4202 
4203  val[0] -=deltamean*1e-4;;
4204 }
4205 //Helpers*/
4206 
4207 
4208 Bool_t AliAnalysisTaskHFJetIPQA::GetImpactParameter(const AliAODTrack *track, const AliAODEvent *event, Double_t *dca, Double_t *cov, Double_t *XYZatDCA)
4209 {
4210  if(!track || !event) return kFALSE;
4211  if(dca==0 || cov ==0 ||XYZatDCA ==0 ) return kFALSE;
4212  AliExternalTrackParam etp;
4213  etp.CopyFromVTrack((AliVTrack*)track);
4214 
4215  const Double_t kBeampiperadius=3; //maximal dca used for track propagation
4216  fEventVertex = RemoveDaughtersFromPrimaryVtx(track);
4217  if(!fEventVertex)return kFALSE;
4218  //printf("After remove:\n");
4219  //fEventVertex->Print();
4220  const AliVVertex *vtxESDSkip =fEventVertex;//RemoveDaughtersFromPrimaryVtx(track);
4221  if(!vtxESDSkip) return kFALSE;
4222  if(!etp.PropagateTo(1,event->GetMagneticField())) return kFALSE;
4223  Bool_t success = kFALSE;
4224  Double_t x_at_dca[3];
4225  Double_t p_at_dca[3];
4226 
4227  //Classic dca calculation
4228  if(etp.PropagateToDCA(vtxESDSkip, event->GetMagneticField(), kBeampiperadius, dca, cov)){
4229  success = kTRUE;
4230  etp.GetXYZ(XYZatDCA);
4231  etp.GetXYZ(x_at_dca);
4232  etp.GetPxPyPz(p_at_dca);
4233  // if(fIsPythia) dca[0] *= fMCglobalDCASmear;
4234  // if(fIsPythia) dca[0] += fMCglobalDCAxyShift; // generic mean offset in LHC10e default is 0.007 == 7 µm
4235 
4236  } else return kFALSE;
4237  return success;
4238 }
4239 
4240 /*AliExternalTrackParam AliAnalysisTaskHFJetIPQA::GetExternalParamFromJet(const AliEmcalJet *jet, const AliAODEvent *event)
4241 {
4242  double vtx[3]= {0.};
4243  double cov [21] = {0.};
4244  double pxpypz[3] = {0.};
4245  jet->PxPyPz(pxpypz);
4246  (event->GetPrimaryVertex())->GetXYZ(vtx);
4247  AliExternalTrackParam etp (vtx, pxpypz, cov, (Short_t)0);
4248  return etp;
4249 }*/
4250 
4251 
4252 
4253 Bool_t AliAnalysisTaskHFJetIPQA::getJetVtxMass(AliEmcalJet *jet,double &value ){
4254 
4255  return true;
4256 
4257 
4258 
4259  const AliVVertex * vertex = InputEvent()->GetPrimaryVertex();
4260 
4261  Printf("Primary vertex %e %e %e",vertex->GetX(),vertex->GetY(),vertex->GetY() );
4262  /*
4263  if(!amvf) amvf = new AliHFAdaptiveMVF(InputEvent()->GetMagneticField());
4264  amvf->_fitter->_seeder->_vertex_penalty((AliAODVertex*)InputEvent()->GetPrimaryVertex());
4265  amvf->_fitter->_r_fitter_jet(jet,(AliAODEvent*)InputEvent());
4266  return kTRUE;*/
4267 }
4268 
4269 
4270 Bool_t AliAnalysisTaskHFJetIPQA::GetImpactParameterWrtToJet(const AliAODTrack *track, const AliAODEvent *event, const AliEmcalJet *jet, Double_t *dca, Double_t *cov, Double_t *XYZatDCA, Double_t &jetsign, int jetflavour)
4271 {
4272  if(!track || !event || !jet)return kFALSE;
4273  if(dca==0 || cov ==0 ||XYZatDCA ==0 ) return kFALSE;
4274 
4275  if(!GetImpactParameter(track,event,dca,cov,XYZatDCA)) return kFALSE;
4276  //vertex properties
4277  Double_t VxVyVz[3]= {0.,0.,0.};
4278  const AliVVertex *vtxESDSkip =fEventVertex; //GetImpactParameter does set fEventVertex to the recalculated one
4279  if(!fEventVertex) return kFALSE;
4280  vtxESDSkip->GetXYZ(VxVyVz);
4281  //printf("Vertex in wrt jet:\n");
4282  //vtxESDSkip->Print();
4283 
4284  //jet properties
4285  double jetp[3];
4286  jet->PxPyPz(jetp);
4287  TVector3 jetP3(jetp);
4288  double covjet [21] = {0.};
4289  double pxpypz[3] = {0.};
4290  jet->PxPyPz(pxpypz);
4291  AliExternalTrackParam etp_jet(VxVyVz, pxpypz, covjet, (Short_t)0);
4292 
4293  //Calculation of sign
4294  TVector3 JetDir =jetP3.Unit();
4295  TVector3 D0(XYZatDCA);
4296  TVector3 vertex(VxVyVz);
4297  TVector3 DD0(D0.x()-vertex.x(),D0.y()-vertex.y(),0.); //track impact parameter
4298  double ps =DD0.Dot(JetDir);
4299  double value = DD0.Mag()*(ps/fabs(ps)); //absolut track impact parameter
4300  jetsign = TMath::Sign(1.,value); //sign
4301  TVector3 dd0 =DD0.Unit();
4302 
4303  //track properties
4304  AliExternalTrackParam etp_track; etp_track.CopyFromVTrack(track);
4305  Double_t xa,xb,xyz_jet_global[3],xyz_track_global[3];
4306 
4307  etp_jet.GetDCA(&etp_track, event->GetMagneticField(), xa, xb);
4308  etp_jet.GetXYZAt(xa, event->GetMagneticField(),xyz_jet_global);
4309  etp_track.GetXYZAt(xb, event->GetMagneticField(),xyz_track_global);
4310  etp_track.PropagateTo(xb,event->GetMagneticField());
4311 
4312  if(fEventVertex) {
4313  delete fEventVertex;
4314  fEventVertex =nullptr;
4315 
4316  }
4317 
4318  double val = ((VxVyVz[0] - xyz_track_global[0]) * (VxVyVz[0] - xyz_track_global[0]) +
4319  (VxVyVz[1] - xyz_track_global[1]) * (VxVyVz[1] - xyz_track_global[1])+
4320  (VxVyVz[2] - xyz_track_global[2]) * (VxVyVz[2] - xyz_track_global[2]));
4321 
4322 
4323  double bdecaylength = val >0 ? sqrt(val) : 1000;
4324  //printf("decaylength:\n");
4325  //for(int i=0;i<3;i++){
4326  // printf("VxVyVZ=%f, xyztrackglobal=%f\n",VxVyVz[i], xyz_track_global[i]);
4327  //}
4328 
4329  double dcajetrack = sqrt((xyz_jet_global[0] - xyz_track_global[0]) * (xyz_jet_global[0] - xyz_track_global[0]) +
4330  (xyz_jet_global[1] - xyz_track_global[1]) * (xyz_jet_global[1] - xyz_track_global[1])+
4331  (xyz_jet_global[2] - xyz_track_global[2]) * (xyz_jet_global[2]- xyz_track_global[2]));
4332 
4333  //printf("decaylength:\n");
4334  //for(int i=0;i<3;i++){
4335  // printf("xyzjetglobal=%f, xyztrackglobal=%f\n",xyz_jet_global[i], xyz_track_global[i]);
4336  //}
4337 
4338  if(!(IsDCAAccepted(bdecaylength, dcajetrack, dca, jetflavour))) return kFALSE;
4339  //printf("decaylength=%f, pi=%f\n", bdecaylength, dcajetrack);
4340 
4341  return kTRUE;
4342 }
4343 
4344 
4345 void AliAnalysisTaskHFJetIPQA::SetTCThresholds(TObjArray** &threshs){
4346  if(fDoTCTagging==TCIPFixedPt){
4347  fNThresholds=1;
4348  }
4349  for(int iProbSet=0;iProbSet<fNThresholds;iProbSet++){
4350  TObjArray* oa=(TObjArray*)threshs[iProbSet];
4351 
4352  //printf("Pointer oa=%p\n",oa);
4353 
4354  h1DThresholdsFirst.push_back((TH1D*)oa->At(0));
4355  h1DThresholdsSecond.push_back((TH1D*)oa->At(1));
4356  h1DThresholdsThird.push_back((TH1D*)oa->At(2));
4357 
4358  TString sFrac=h1DThresholdsFirst[iProbSet]->GetTitle();
4359  fFracs.push_back(sFrac.Atof());
4360  }
4361  //checking fFracs
4362  //for(int iFrac=0;iFrac<fFracs.size();iFrac++){
4363  // printf("iFrac=%i, %f\n",iFrac,fFracs[iFrac]);
4364  //}
4365  /* int nPoints=h1DThresholdsFirst[0]->GetNbinsX();
4366 
4367  printf("LargestIP: 2nd Histogram bins:\n");
4368  for(int iPoint=0;iPoint<nPoints;iPoint++){
4369  printf("iPoint=%i, xval=%f, yval=%f\n",iPoint, h1DThresholdsFirst[1]->GetBinCenter(iPoint),h1DThresholdsFirst[1]->GetBinContent(iPoint));
4370  }
4371  printf("2ndLargestIP: 3rd Histogram bins:\n");
4372  for(int iPoint=0;iPoint<nPoints;iPoint++)
4373  printf("iPoint=%i, xval=%f, yval=%f\n",iPoint, h1DThresholdsSecond[2]->GetBinCenter(iPoint),h1DThresholdsSecond[2]->GetBinContent(iPoint));
4374  }
4375  printf("3rdLargestIP: 1st Histogram bins:\n");
4376  for(int iPoint=0;iPoint<nPoints;iPoint++){
4377  printf("iPoint=%i, xval=%f, yval=%f\n",iPoint, h1DThresholdsThird[0]->GetBinCenter(iPoint),h1DThresholdsThird[0]->GetBinContent(iPoint));
4378  }*/
4379 }
4380 
4381 void AliAnalysisTaskHFJetIPQA::SetProbThresholds(TObjArray** &threshs){
4382  for(int iProbSet=0;iProbSet<fNThresholds;iProbSet++){
4383  TObjArray* oa=(TObjArray*)threshs[iProbSet];
4384  if(!oa) AliError(Form(" No %i'th Probability Threshold object array!\n",iProbSet));
4385  printf("Pointer oa=%p\n",oa);
4386  h1DProbThresholds.push_back((TH1D*)oa->At(0));
4387  if(!h1DProbThresholds.back()) AliError(Form(" No %i'th Probability Threshold hist!\n",iProbSet));
4388  }
4389 
4390  /*int nPoints=h1DProbThresholds[0]->GetNbinsX();
4391  for(int iPoint=0;iPoint<nPoints;iPoint++){
4392  printf("iPoint=%i, xval=%f, yval=%f\n",iPoint, h1DProbThresholds[0]->GetXaxis()->GetBinLowEdge(iPoint),h1DProbThresholds[0]->GetBinContent(iPoint));
4393  }*/
4394 }
4395 
4396 // Read Threshold Histograms
4397 //==============================================================================
4398 void AliAnalysisTaskHFJetIPQA::ReadThresholdHists(TString PathToThresholds, TString taskname, int nTCThresh){
4399  TFile* fileThresholds=TFile::Open(PathToThresholds.Data());
4400  if(!fileThresholds ||(fileThresholds&& !fileThresholds->IsOpen())){AliError(Form("%s :: File with threshold values not found",taskname.Data()));}
4401 
4402  Printf("%s :: File %s successfully loaded, setting up threshold functions.",taskname.Data(),PathToThresholds.Data());
4403 
4404  if(fileThresholds){
4405  printf("Reading threshold histograms for track counting...\n");
4406 
4407  //TC Thresholds
4408  TObjArray** oaTCThresh=new TObjArray*[nTCThresh];
4409  for(int iThresh=0;iThresh<nTCThresh;iThresh++){
4410  fileThresholds->GetObject(Form("TCThres_%i",iThresh),oaTCThresh[iThresh]);
4411  }
4412 
4413  //ProbLookup hists
4414  TObjArray* oLookup;
4415  fileThresholds->GetObject("ProbLookup",oLookup);
4416 
4417  /* TObjArray** oaProbThresh=new TObjArray*[nTCThresh];
4418  for(int iThresh=0;iThresh<nTCThresh;iThresh++){
4419  fileThresholds->GetObject(Form("ProbThres_%i",iThresh),oaProbThresh[iThresh]);
4420  }*/
4421 
4422  this->setfNThresholds(nTCThresh);
4423  this->SetTCThresholds(oaTCThresh);
4424  //this->SetProbThresholds(oaProbThresh);
4425  this->ReadProbvsIPLookup(oLookup);
4426  }
4427 }
4428 
4429 void AliAnalysisTaskHFJetIPQA::ReadProbvsIPLookup(TObjArray*& oLookup){
4430 
4431  for(int iN=0;iN<3;iN++){
4432  h2DProbLookup.push_back((TH2D*)oLookup->At(iN));
4433  }
4434 }
4435 
4436 void AliAnalysisTaskHFJetIPQA::DoTCTagging(double jetpt, bool* hasIPs, double* ipval, bool **kTagDec){
4437  //printf("Start TCTagging!\n");
4438  //threshold values for tracks with largest, second and third largest IP
4439  int iJetPtBin=h1DThresholdsFirst[0]->FindBin(jetpt);
4440  double IPthresN1[fNThresholds]; //IP threshold values for individual separation power
4441  double IPthresN2[fNThresholds];
4442  double IPthresN3[fNThresholds];
4443 
4444  if(fDoTCTagging==TCIPSig){
4445  for(int iN=0;iN<fNThresholds;iN++){
4446  IPthresN1[iN]=h1DThresholdsFirst[iN]->GetBinContent(iJetPtBin);
4447  IPthresN2[iN]=h1DThresholdsSecond[iN]->GetBinContent(iJetPtBin);
4448  IPthresN3[iN]=h1DThresholdsThird[iN]->GetBinContent(iJetPtBin);
4449  }
4450  }
4451  if(fDoTCTagging==TCIPFixedPt){
4452  for(int iN=0;iN<fNThresholds;iN++){
4453  IPthresN1[iN]=fTCThresholdPtFixed;
4454  IPthresN2[iN]=fTCThresholdPtFixed;
4455  IPthresN3[iN]=fTCThresholdPtFixed;
4456  }
4457  }
4458 
4459  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
4460  if(!hasIPs[0]) continue;
4461  //printf("DoTCTagging:\n");
4462  //printf(" iJetPtBin=%i, IPthresN1=%f, IPthresN2=%f, IPthresN3=%f\n", iJetPtBin, IPthresN1[iThresh],IPthresN2[iThresh], IPthresN3[iThresh]);
4463 
4464 
4465  if(hasIPs[2]){
4466  //tripple tag
4467  //printf("ipval[0]=%f, ipval[1]=%f, ipval[2]=%f\n", ipval[0],ipval[1],ipval[2]);
4468  if(ipval[0]>IPthresN1[iThresh]&&ipval[1]>IPthresN2[iThresh]&&ipval[2]>IPthresN3[iThresh]){
4469  //printf("Triple %f!\n",fFracs[iThresh]);
4470  kTagDec[iThresh][Full]=kTRUE; kTagDec[iThresh][Triple]=kTRUE;
4471  }
4472 
4473  //single tag
4474  if(kTagLevel<2){
4475  //printf("Single catch\n");
4476  if(ipval[2]>IPthresN3[iThresh]) {
4477  // printf("Single3rd %f!\n",fFracs[iThresh]);
4478  kTagDec[iThresh][Full]=kTRUE; kTagDec[iThresh][Single3rd]=kTRUE;
4479  }
4480  }
4481  }
4482 
4483  if(hasIPs[1]){
4484  //printf("ipval[0]=%f, ipval[1]=%f\n", ipval[0],ipval[1]);
4485  //double tag
4486  if(kTagLevel<3){
4487  //printf("Double catch\n");
4488  if(ipval[0]>IPthresN1[iThresh]&&ipval[1]>IPthresN2[iThresh]) {
4489  //printf("Double %f!\n",fFracs[iThresh]);
4490  kTagDec[iThresh][Full]=kTRUE; kTagDec[iThresh][Double]=kTRUE;
4491  }
4492  }
4493  //single tag
4494  if(kTagLevel<2){
4495  //printf("Single catch\n");
4496  if(ipval[1]>IPthresN2[iThresh]) {
4497  //printf("Single2nd %f!\n",fFracs[iThresh]);
4498  kTagDec[iThresh][Full]=kTRUE; kTagDec[iThresh][Single2nd]=kTRUE;
4499  }
4500  }
4501  }
4502 
4503  //single tag
4504  if(hasIPs[0]){
4505  //printf("ipval[0]=%f", ipval[0]);
4506  if(kTagLevel<2){
4507  //printf("Single catch\n");
4508  if(ipval[0]>IPthresN1[iThresh]) {
4509  //printf("Single1st %f!\n",fFracs[iThresh]);
4510  kTagDec[iThresh][Full]=kTRUE; kTagDec[iThresh][Single1st]=kTRUE;
4511  }
4512  }
4513  }
4514  /*printf("Testing kTagLevel:: Line %i\n ",__LINE__);
4515  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
4516  for(int iType=0;iType<6;iType++){
4517  printf("iThresh=%f, %i, kTagDec=%i\n",fFracs[iThresh],iType,kTagDec[iThresh][iType]);
4518  }
4519  }*/
4520  }
4521 }
4522 
4523 void AliAnalysisTaskHFJetIPQA::DoProbTagging(double probval, double jetpt,bool **kTagDec){
4524  int iJetPtBin=h1DProbThresholds[0]->FindBin(jetpt);
4525  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
4526  double threshval=h1DProbThresholds[iThresh]->GetBinContent(iJetPtBin);
4527  //printf(" iThres=%i, iJetPtBin=%i, jetpt=%f, iThreshold=%f, probval=%f", iThresh, iJetPtBin, jetpt, threshval, probval);
4528 
4529  if(probval>threshval){
4530  kTagDec[iThresh][Full]=kTRUE;
4531  //printf("Tagging condition fullfilled %i!\n",jetflavour);
4532  }
4533  }
4534 }
4535 
4536 void AliAnalysisTaskHFJetIPQA::FillEfficiencyHists(bool** kTagDec, int jetflavour, double jetpt, bool hasIPs){
4537  //printf("Receiving BTagged decision: %i\n", kTagDec[Full]);
4538  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
4539  //printf("kTagDec=%i, jetflavour=%i, hasIPs=%i\n",kTagDec[iThresh][Full],jetflavour, hasIPs);
4540  if(!fIsPythia&&kTagDec[iThresh][Full]){
4541  FillHist(Form("h1DTagged_Full_%0.2f",fFracs[iThresh]),jetpt,1);
4542  if(kTagDec[iThresh][Single1st]) FillHist(Form("h1DTagged_Single1st_%0.2f",fFracs[iThresh]), jetpt, 1);
4543  if(kTagDec[iThresh][Single2nd]) FillHist(Form("h1DTagged_Single2nd_%0.2f",fFracs[iThresh]), jetpt, 1);
4544  if(kTagDec[iThresh][Single3rd])FillHist(Form("h1DTagged_Single3rd_%0.2f",fFracs[iThresh]), jetpt, 1);
4545  if(kTagDec[iThresh][Double])FillHist(Form("h1DTagged_Double_%0.2f",fFracs[iThresh]), jetpt, 1);
4546  if(kTagDec[iThresh][Triple])FillHist(Form("h1DTagged_Triple_%0.2f",fFracs[iThresh]), jetpt, 1);
4547  }
4548 
4549  if(fIsPythia){
4550  if(kTagDec[iThresh][Full]&&(jetflavour!=Unid)){
4551  FillHist(Form("h1DTagged_Full_%0.2f",fFracs[iThresh]),jetpt,1);
4552  if(kTagDec[iThresh][Single1st]) FillHist(Form("h1DTagged_Single1st_%0.2f",fFracs[iThresh]), jetpt, 1);
4553  if(kTagDec[iThresh][Single2nd]) FillHist(Form("h1DTagged_Single2nd_%0.2f",fFracs[iThresh]), jetpt, 1);
4554  if(kTagDec[iThresh][Single3rd])FillHist(Form("h1DTagged_Single3rd_%0.2f",fFracs[iThresh]), jetpt, 1);
4555  if(kTagDec[iThresh][Double])FillHist(Form("h1DTagged_Double_%0.2f",fFracs[iThresh]), jetpt, 1);
4556  if(kTagDec[iThresh][Triple])FillHist(Form("h1DTagged_Triple_%0.2f",fFracs[iThresh]), jetpt, 1);
4557  }
4558  if(kTagDec[iThresh][Full]&&(jetflavour==B)&&hasIPs){
4559  //printf("################################ Before: FoundJet with tagindex=%i!\n",kTagDec[iThresh][Full]);
4560 
4561  FillHist(Form("h1DTrueBTagged_Full_%0.2f",fFracs[iThresh]), jetpt, 1);
4562  if(kTagDec[iThresh][Single1st]) FillHist(Form("h1DTrueBTagged_Single1st_%0.2f",fFracs[iThresh]), jetpt, 1);
4563  if(kTagDec[iThresh][Single2nd]) FillHist(Form("h1DTrueBTagged_Single2nd_%0.2f",fFracs[iThresh]), jetpt, 1);
4564  if(kTagDec[iThresh][Single3rd])FillHist(Form("h1DTrueBTagged_Single3rd_%0.2f",fFracs[iThresh]), jetpt, 1);
4565  if(kTagDec[iThresh][Double])FillHist(Form("h1DTrueBTagged_Double_%0.2f",fFracs[iThresh]), jetpt, 1);
4566  if(kTagDec[iThresh][Triple])FillHist(Form("h1DTrueBTagged_Triple_%0.2f",fFracs[iThresh]), jetpt, 1);
4567 
4568  //printf("################################ FoundJet with tagindex=%i!\n",kTagDec[iThresh][Full]);
4569  }
4570  if(kTagDec[iThresh][Full]&&(jetflavour==C)&&hasIPs){
4571  FillHist(Form("h1DFalseCTagged_Full_%0.2f",fFracs[iThresh]), jetpt, 1);
4572  if(kTagDec[iThresh][Single1st]) FillHist(Form("h1DFalseCTagged_Single1st_%0.2f",fFracs[iThresh]), jetpt, 1);
4573  if(kTagDec[iThresh][Single2nd]) FillHist(Form("h1DFalseCTagged_Single2nd_%0.2f",fFracs[iThresh]), jetpt, 1);
4574  if(kTagDec[iThresh][Single3rd]) FillHist(Form("h1DFalseCTagged_Single3rd_%0.2f",fFracs[iThresh]), jetpt, 1);
4575  if(kTagDec[iThresh][Double])FillHist(Form("h1DFalseCTagged_Double_%0.2f",fFracs[iThresh]), jetpt, 1);
4576  if(kTagDec[iThresh][Triple])FillHist(Form("h1DFalseCTagged_Triple_%0.2f",fFracs[iThresh]), jetpt, 1);
4577  //printf("################################ CMistagged: flavour is=%i with tagindex=%i!\n",jetflavour,kTagDec[iThresh][Full]);
4578  }
4579  if(kTagDec[iThresh][Full]&&(jetflavour==UDSG)&&hasIPs){
4580  FillHist(Form("h1DFalseUDSGTagged_Full_%0.2f",fFracs[iThresh]), jetpt, 1);
4581  if(kTagDec[iThresh][Single1st]) FillHist(Form("h1DFalseUDSGTagged_Single1st_%0.2f",fFracs[iThresh]), jetpt, 1);
4582  if(kTagDec[iThresh][Single2nd]) FillHist(Form("h1DFalseUDSGTagged_Single2nd_%0.2f",fFracs[iThresh]), jetpt, 1);
4583  if(kTagDec[iThresh][Single3rd]) FillHist(Form("h1DFalseUDSGTagged_Single3rd_%0.2f",fFracs[iThresh]), jetpt, 1);
4584  if(kTagDec[iThresh][Double])FillHist(Form("h1DFalseUDSGTagged_Double_%0.2f",fFracs[iThresh]), jetpt, 1);
4585  if(kTagDec[iThresh][Triple])FillHist(Form("h1DFalseUDSGTagged_Triple_%0.2f",fFracs[iThresh]), jetpt, 1);
4586  //printf("################################ LightMistagged: flavour is=%i with tagindex=%i!\n",jetflavour,kTagDec[iThresh][Full]);
4587  }
4588  if((kTagDec[iThresh][Full])&&((jetflavour==UDSGV0)||(jetflavour==CV0))&&hasIPs){
4589  FillHist(Form("h1DFalseV0Tagged_Full_%0.2f",fFracs[iThresh]), jetpt, 1);
4590  if(kTagDec[iThresh][Single1st]) FillHist(Form("h1DFalseV0Tagged_Single1st_%0.2f",fFracs[iThresh]), jetpt, 1);
4591  if(kTagDec[iThresh][Single2nd]) FillHist(Form("h1DFalseV0Tagged_Single2nd_%0.2f",fFracs[iThresh]), jetpt, 1);
4592  if(kTagDec[iThresh][Single3rd]) FillHist(Form("h1DFalseV0Tagged_Single3rd_%0.2f",fFracs[iThresh]), jetpt, 1);
4593  if(kTagDec[iThresh][Double])FillHist(Form("h1DFalseV0Tagged_Double_%0.2f",fFracs[iThresh]), jetpt, 1);
4594  if(kTagDec[iThresh][Triple])FillHist(Form("h1DFalseV0Tagged_Triple_%0.2f",fFracs[iThresh]), jetpt, 1);
4595  //printf("################################ v0Mistagged: flavour is=%i with tagindex=%i!\n",jetflavour,kTagDec[iThresh][Full]);
4596  }
4597  if(!kTagDec[iThresh][Full]&&jetflavour==B&&hasIPs){
4598  //printf("################################ Missed one: flavour is=%i\n", jetflavour);
4599  }
4600  }
4601  }
4602 }
4603 
4604 void AliAnalysisTaskHFJetIPQA::FillTaggedJetPtDistribution(bool** kTagDec, double jetpt){
4605  const char * tagtype[6] = {"Full","Single1st","Single2nd","Single3rd","Double","Triple"};
4606 
4607  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
4608  for(int iType=0;iType<6;iType++){
4609  if(kTagDec[iThresh][iType]){
4610  //printf("Filling fracs=%f, iType=%s\n",fFracs[iThresh],tagtype[iType]);
4611  FillHist(Form("h1DTaggedJetPt_%s_%0.2f",tagtype[iType],fFracs[iThresh]),jetpt,1);
4612  }
4613  }
4614  }
4615 }
4616 
4617 void AliAnalysisTaskHFJetIPQA::FillV0EfficiencyHists(int isV0, int &jetflavour, double jetpt, bool &isV0Jet){
4618  //printf("isV0=%i, jetflavour=%i, jetpt=%f\n",isV0, jetflavour, jetpt);
4619  switch (isV0){
4620  case V0Rec:
4621  FillHist(Form("h1DV0FalseRec"), jetpt, 1);
4622  //printf("Found false Tag\n");
4623  break;
4624 
4625  case V0MC:
4626  FillHist(Form("h1DV0TrueDataDef"), jetpt, 1);
4627  if((jetflavour!=B)){
4628  isV0Jet=kTRUE;
4629  FillHist(Form("h1DV0TrueMCDef"), jetpt, 1);
4630  //printf("Found MC def true V0\n");
4631  }
4632  //printf("Found MC true V0 Jet: jetflavour=%i\n",jetflavour);
4633  if(jetflavour==UDSG) jetflavour=UDSGV0;
4634  if(jetflavour==C) jetflavour=CV0;
4635  break;
4636 
4637  case V0TrueRec:
4638  FillHist(Form("h1DV0TrueDataDef"), jetpt, 1);
4639  FillHist(Form("h1DV0TrueRec"), jetpt, 1);
4640  if(jetflavour==UDSG) jetflavour=UDSGV0;
4641  if(jetflavour==C) jetflavour=CV0;
4642  if(jetflavour!=B){
4643  FillHist(Form("h1DV0TrueMCDef"), jetpt, 1);
4644  isV0Jet=kTRUE;
4645  //printf("Found MC def true V0\n");
4646  }
4647  //printf("Found Rec true V0 Jet: jetflavour=%i\n",jetflavour);
4648  break;
4649  }
4650 }
4651 
4652 double AliAnalysisTaskHFJetIPQA::IntegrateIP(int iJetPtBin, int iIPBin, int iN){
4653  int iZeroIPBin=h2DProbLookup[iN]->GetXaxis()->FindBin(0.);
4654  int iStartIPBin=h2DProbLookup[iN]->GetXaxis()->FindBin(-25);
4655 
4656  double prob=h2DProbLookup[iN]->Integral(iStartIPBin,iIPBin,iJetPtBin,iJetPtBin);
4657  prob=prob/(h2DProbLookup[iN]->Integral(iStartIPBin,iZeroIPBin,iJetPtBin,iJetPtBin));
4658 
4659  //printf("Integrate: 0x=%f, lowx=%f, upx=%f, lowy=%f, upy=%f, prob=%f\n", h2DProbLookup[iN]->GetXaxis()->GetBinLowEdge(iZeroIPBin), h2DProbLookup[iN]->GetXaxis()->GetBinLowEdge(iStartIPBin),h2DProbLookup[iN]->GetXaxis()->GetBinLowEdge(iIPBin),h2DProbLookup[iN]->GetYaxis()->GetBinLowEdge(iJetPtBin),h2DProbLookup[iN]->GetYaxis()->GetBinLowEdge(iJetPtBin+1),prob);
4660 
4661  return prob;
4662 }
4663 
4664 double AliAnalysisTaskHFJetIPQA::GetTrackProbability(double jetpt, bool* hasIPs, double* ipval){
4665  //printf("Printing h2DProbLook\n");
4666  double prob=1;
4667  double probval[3]={0};
4668  int iJetPtBin=h2DProbLookup[0]->GetYaxis()->FindBin(jetpt);;
4669  int iIPBin[3]={0};
4670  //printf("ipval1=%f, ipval2=%f, ipval3=%f\n",ipval[0],ipval[1],ipval[2]);
4671 
4672  for(int iN=0;iN<3;iN++){
4673  if(!hasIPs[iN])continue;
4674  if(ipval[iN]<0) continue;
4675  iIPBin[iN]=h2DProbLookup[iN]->GetXaxis()->FindBin(-ipval[iN]);
4676  probval[iN]=IntegrateIP(iJetPtBin,iIPBin[iN], iN);
4677  //probval[iN]=h2DProbLookup[iN]->GetBinContent(iIPBin[iN],iJetPtBin);
4678  //printf("iN=%i, iIPBin=%i, ipval=%f, lowerIP=%f, higherIP=%f, || iJetPtBin=%i, jetpt=%f, lowerjetpt=%f, higherjetpt=%f, prob=%f\n", iN, iIPBin[iN],ipval[iN],h2DProbLookup[iN]->GetXaxis()->GetBinLowEdge(iIPBin[iN]),
4679  // h2DProbLookup[iN]->GetXaxis()->GetBinLowEdge(iIPBin[iN]+1), iJetPtBin,jetpt, h2DProbLookup[iN]->GetYaxis()->GetBinLowEdge(iJetPtBin), h2DProbLookup[iN]->GetYaxis()->GetBinLowEdge(iJetPtBin+1),probval[iN]);
4680  prob=prob*probval[iN];
4681  }
4682 
4683  double prob3=prob-prob*TMath::Log(prob)+prob*(TMath::Log(prob)*TMath::Log(prob))*0.5;
4684  double prob2=prob-prob*TMath::Log(prob);
4685  double prob1=prob;
4686 
4687 
4688  if(hasIPs[2]&&ipval[2]>0){
4689  //printf("3 tracks with prob1=%f, prob2=%f, prob3=%f, prob=%f, prob3=%f\n",probval[0],probval[1],probval[2],prob,prob3);
4690  return prob3;
4691  }
4692  if(hasIPs[1]&&ipval[1]>0){
4693  //printf("2 tracks with prob1=%f, prob2=%f prob=%f, prob2=%f\n",probval[0],probval[1],prob,prob2);
4694  return prob2;
4695  }
4696  if(hasIPs[0]&&ipval[0]>0){
4697  //printf("1 track with prob1=%f, prob=%f, prob1=%f\n",probval[0],prob,prob1);
4698  return prob1;
4699  }
4700  return 0;
4701 }
4702 
4703 void AliAnalysisTaskHFJetIPQA::FillProbabilityHists(double jetpt,double probval,int jetflavour,bool **kTagDec){
4704  // printf("Filling iflavou=%i, jetpt=%f, probval=%f into histogram\n", jetflavour,jetpt,probval);
4705  double lnprobval=-TMath::Log(probval);
4706 
4707  //printf("logval=%f\n",lnprobval);
4708  if(fIsPythia){
4709  FillHist(Form("h2DProbDists%s",sTemplateFlavour[jetflavour].Data()),probval,jetpt,1); //*this->fXsectionWeightingFactor );
4710  FillHist(Form("h2DLNProbDists%s",sTemplateFlavour[jetflavour].Data()),lnprobval,jetpt,1); //*this->fXsectionWeightingFactor );
4711  }
4712  //untagged Probability hist
4713  FillHist(Form("h2DProbDists"),probval,jetpt,1); //*this->fXsectionWeightingFactor );
4714  FillHist(Form("h2DLNProbDists"),lnprobval,jetpt,1); //*this->fXsectionWeightingFactor );
4715  const char * tagtype[6] = {"Full","Single1st","Single2nd","Single3rd","Double","Triple"};
4716 
4717  //tagged Probability hists
4718  if(fDoTCTagging!=TCNo){
4719  for(int iThresh=0;iThresh<fNThresholds;iThresh++){
4720  for(int iType=0;iType<6;iType++){
4721  if(kTagDec[iThresh][iType]){
4722  FillHist(Form("h2DProbDistsTag_%s_%0.2f",tagtype[iType],fFracs[iThresh]),probval,jetpt,1);
4723  FillHist(Form("h2DLNProbDistsTag_%s_%0.2f",tagtype[iType],fFracs[iThresh]),lnprobval,jetpt,1);
4724  }
4725  switch(jetflavour){
4726  case UDSG:
4727  if(iThresh==0&&iType==0){
4728  FillHist(Form("h2DProbDists_UDSG"),probval,jetpt,1);
4729  FillHist(Form("h2DLNProbDists_UDSG"),lnprobval,jetpt,1);
4730  //printf("Filling total UDSG %i\n",jetflavour);
4731  }
4732  if(kTagDec[iThresh][iType]){
4733  FillHist(Form("h2DProbDistsTag_UDSG_%s_%0.2f",tagtype[iType],fFracs[iThresh]),probval,jetpt,1);
4734  FillHist(Form("h2DLNProbDistsTag_UDSG_%s_%0.2f",tagtype[iType],fFracs[iThresh]),lnprobval,jetpt,1);
4735  //printf("Filling UDSG %i %s %0.2f\n", jetflavour, tagtype[iType],fFracs[iThresh]);
4736  }
4737  break;
4738 
4739  case B:
4740  if(iThresh==0&&iType==0){
4741  FillHist(Form("h2DProbDists_B"),probval,jetpt,1);
4742  FillHist(Form("h2DLNProbDists_B"),lnprobval,jetpt,1);
4743  //printf("Filling total B %i\n",jetflavour);
4744  }
4745  if(kTagDec[iThresh][iType]){
4746  FillHist(Form("h2DProbDistsTag_B_%s_%0.2f",tagtype[iType],fFracs[iThresh]),probval,jetpt,1);
4747  FillHist(Form("h2DLNProbDistsTag_B_%s_%0.2f",tagtype[iType],fFracs[iThresh]),lnprobval,jetpt,1);
4748  //printf("Filling B %i %s %0.2f\n,", jetflavour,tagtype[iType],fFracs[iThresh]);
4749  }
4750  break;
4751 
4752  case C:
4753  if(iThresh==0&&iType==0){
4754  FillHist(Form("h2DProbDists_C"),probval,jetpt,1);
4755  FillHist(Form("h2DLNProbDists_C"),lnprobval,jetpt,1);
4756  //printf("Filling total C %i\n",jetflavour);
4757  }
4758  if(kTagDec[iThresh][iType]){
4759  FillHist(Form("h2DProbDistsTag_C_%s_%0.2f",tagtype[iType],fFracs[iThresh]),probval,jetpt,1);
4760  FillHist(Form("h2DLNProbDistsTag_C_%s_%0.2f",tagtype[iType],fFracs[iThresh]),lnprobval,jetpt,1);
4761  //printf("Filling C %i %s %0.2f\n",jetflavour, tagtype[iType],fFracs[iThresh]);
4762  }
4763  break;
4764 
4765  case UDSGV0:
4766  if(iThresh==0&&iType==0){
4767  FillHist(Form("h2DProbDists_V0"),probval,jetpt,1);
4768  FillHist(Form("h2DLNProbDists_V0"),lnprobval,jetpt,1);
4769  //printf("Filling total V0 %i\n",jetflavour);
4770  }
4771  if(kTagDec[iThresh][iType]){
4772  FillHist(Form("h2DProbDistsTag_V0_%s_%0.2f",tagtype[iType],fFracs[iThresh]),probval,jetpt,1);
4773  FillHist(Form("h2DLNProbDistsTag_V0_%s_%0.2f",tagtype[iType],fFracs[iThresh]),lnprobval,jetpt,1);
4774  //printf("Filling V0 %i %s %0.2f\n",jetflavour, tagtype[iType],fFracs[iThresh]);
4775  }
4776  break;
4777 
4778  case CV0:
4779  if(iThresh==0&&iType==0){
4780  FillHist(Form("h2DProbDists_V0"),probval,jetpt,1);
4781  FillHist(Form("h2DLNProbDists_V0"),lnprobval,jetpt,1);
4782  //printf("Filling total V0 %i\n",jetflavour);
4783  }
4784  if(kTagDec[iThresh][iType]){
4785  FillHist(Form("h2DProbDistsTag_V0_%s_%0.2f",tagtype[iType],fFracs[iThresh]),probval,jetpt,1);
4786  FillHist(Form("h2DLNProbDistsTag_V0_%s_%0.2f",tagtype[iType],fFracs[iThresh]),lnprobval,jetpt,1);
4787  //printf("Filling V0 %i %s %0.2f\n", jetflavour,tagtype[iType],fFracs[iThresh]);
4788  }
4789  break;
4790  }
4791  }
4792  }
4793  }
4794 }
4795 
4796 void AliAnalysisTaskHFJetIPQA::FillProbThreshHists(double probval, double* ipval, double jetpt, int jetflavour, bool* hasIPs, bool** kTagDec){
4797  double lnprobval=-TMath::Log(probval);
4798  TString sFlavour[6]={"Unid","UDSG","C","B","V0","V0"};
4799 
4800  //printf("ipval0=%f, ipval1=%f, ipval2=%f, single1st=%i, double=%i\n", ipval[0],ipval[1],ipval[2],kTagDec[0][Single1st],kTagDec[0][Double]);
4801 
4802  //Untagged
4803  if(ipval[0]>0){FillHist("h2DProb1Above0",lnprobval,jetpt,1);}
4804  if((ipval[0]>0)&&(ipval[1]>0)){FillHist("h2DProb2Above0",lnprobval,jetpt,1);}
4805  if((ipval[0]>0)&&(ipval[1]>0)&&(ipval[2]>0)){FillHist("h2DProb3Above0",lnprobval,jetpt,1);}
4806  if(fIsPythia){
4807  if(ipval[0]>0) {FillHist(Form("h2DProb1Above0_%s",sFlavour[jetflavour].Data()),lnprobval,jetpt,1);}
4808  if((ipval[0]>0)&&(ipval[1]>0)) {FillHist(Form("h2DProb2Above0_%s",sFlavour[jetflavour].Data()),lnprobval,jetpt,1);}
4809  if((ipval[0]>0)&&(ipval[1]>0)&&(ipval[2]>0)) {FillHist(Form("h2DProb3Above0_%s",sFlavour[jetflavour].Data()),lnprobval,jetpt,1);}
4810  }
4811  //Single1st
4812  if(kTagDec[0][Single1st]){FillHist("h2DProb1AboveThresh",lnprobval,jetpt,1);}
4813  if(kTagDec[0][Single1st]&&(ipval[1]>0)) {FillHist("h2DProb1AbThresh1Ab0",lnprobval,jetpt,1);}
4814  if(kTagDec[0][Single1st]&&(ipval[1]>0)&&(ipval[2]>0)) {FillHist("h2DProb1AbThresh2Ab0",lnprobval,jetpt,1);}
4815 
4816  //Single2nd
4817  if(kTagDec[0][Double]) {FillHist("h2DProb2AboveThresh",lnprobval,jetpt,1);}
4818  if(kTagDec[0][Double]&&(ipval[2]>0)) { FillHist("h2DProb2AbThresh1Ab0",lnprobval,jetpt,1);}
4819  if(kTagDec[0][Double]&&(ipval[2]>0)&&(ipval[3]>0)) { FillHist("h2DProb2AbThresh2Ab0",lnprobval,jetpt,1);}
4820 }
4821 
4822 void AliAnalysisTaskHFJetIPQA::Terminate(Option_t *){
4823 
4824  printf("\n*********************************\n");
4825  printf("Corrections:\n");
4826  printf(" MC Corrections (Data/MC+Fluka):%i\n",fDoMCCorrection);
4827  printf(" Track Smearing:%i\n",fRunSmearing );
4828  printf(" Underlying Event Subtraction:%i\n", fDoUnderlyingEventSub);
4829  printf("*********************************\n");
4830 }
AliAnalysisTaskHFJetIPQA::GetTrackCurvature
Double_t GetTrackCurvature(AliAODTrack *track)
Definition: AliAnalysisTaskHFJetIPQA.cxx:1960
AliEmcalJet::SetSecondClosestJet
void SetSecondClosestJet(AliEmcalJet *j, Double_t d)
Definition: AliEmcalJet.h:338
pdg
Int_t pdg
Definition: charmCutsOptimization.C:34
filename
const char * filename
Definition: TestFCM.C:1
AliAnalysisTaskHFJetIPQA::ReadThresholdHists
void ReadThresholdHists(TString PathToThresholds, TString taskname, int nTCThresh)
Definition: AliAnalysisTaskHFJetIPQA.cxx:4398
AliEmcalJet::Area
Double_t Area() const
Definition: AliEmcalJet.h:130
AliEmcalJet::AddTrackAt
void AddTrackAt(Int_t track, Int_t idx)
Definition: AliEmcalJet.h:275
AliAnalysisTaskHFJetIPQA::fPSJet
std::vector< Double_t > fPSJet
Definition: AliAnalysisTaskHFJetIPQA.h:602
AliJetContainer::GetRhoVal
Double_t GetRhoVal() const
Definition: AliJetContainer.h:165
AliAnalysisTaskHFJetIPQA.h
kElectron
Definition: EstimateSingleTrackPIDsyst.C:43
AliAnalysisTaskEmcal::fIsPythia
Bool_t fIsPythia
trigger, if it is a PYTHIA production
Definition: AliAnalysisTaskEmcal.h:1194
AliAnalysisTaskHFJetIPQA::GetLocalAlphaAOD
Double_t GetLocalAlphaAOD(AliAODTrack *track)
Definition: AliAnalysisTaskHFJetIPQA.cxx:1939
AliAnalysisTaskHFJetIPQA::DetermineUnsuitableVtxTracks
int DetermineUnsuitableVtxTracks(int *skipped, AliAODEvent *const aod, AliVTrack *const track)
Definition: AliAnalysisTaskHFJetIPQA.cxx:2808
AliAnalysisTaskHFJetIPQA::GetVParticleMother
AliVParticle * GetVParticleMother(AliVParticle *part)
GetVParticleMother.
Definition: AliAnalysisTaskHFJetIPQA.cxx:3603
Double_t
double Double_t
Definition: External.C:58
AliAnalysisTaskHFJetIPQA::Terminate
virtual void Terminate(Option_t *option="")
Definition: AliAnalysisTaskHFJetIPQA.cxx:4822
AliEmcalJet::Xv
Double_t Xv() const
Definition: AliEmcalJet.h:112
AliAnalysisTaskHFJetIPQA::setFFillCorrelations
void setFFillCorrelations(const Bool_t &value)
Definition: AliAnalysisTaskHFJetIPQA.cxx:4157
AliAnalysisTaskHFJetIPQA::GetMCTrack
AliAODMCParticle * GetMCTrack(int iLabel)
Definition: AliAnalysisTaskHFJetIPQA.cxx:1135
AliAnalysisTaskHFJetIPQA::SetResFunctionPID
Bool_t SetResFunctionPID(const char *filename)
SetResFunction.
Definition: AliAnalysisTaskHFJetIPQA.cxx:2030
AliAnalysisTaskHFJetIPQA::IsFromElectron
bool IsFromElectron(AliAODTrack *track)
Definition: AliAnalysisTaskHFJetIPQA.cxx:2764
AliAnalysisTaskHFJetIPQA::GetOutOfJetParticleComposition
void GetOutOfJetParticleComposition(AliEmcalJet *jet, int flavour)
Definition: AliAnalysisTaskHFJetIPQA.cxx:2937
AliAnalysisTaskEmcalJet.h
title
const char * title
Definition: MakeQAPdf.C:27
AliAnalysisTaskHFJetIPQA::GetBMesonWeight
Bool_t GetBMesonWeight(AliVParticle *mcpart, Int_t &pdg, Double_t &pT, Int_t &idx)
GetBMesonWeight.
Definition: AliAnalysisTaskHFJetIPQA.cxx:3410
AliEmcalJet::MCPt
Double_t MCPt() const
Definition: AliEmcalJet.h:153
AliAnalysisTaskHFJetIPQA::fBackgroundFactorLinus
Double_t fBackgroundFactorLinus[21][498]
Definition: AliAnalysisTaskHFJetIPQA.h:600
AliEmcalJet::ClosestJet
AliEmcalJet * ClosestJet() const
Definition: AliEmcalJet.h:341
AliAnalysisTaskHFJetIPQA::FillProbThreshHists
void FillProbThreshHists(double proval, double *ipval, double jetpt, int jetflavour, bool *hasIPs, bool **kTagDec)
Definition: AliAnalysisTaskHFJetIPQA.cxx:4796
AliAnalysisTaskHFJetIPQA::GetPtCorrected
Double_t GetPtCorrected(const AliEmcalJet *jet)
GetPtCorrected.
Definition: AliAnalysisTaskHFJetIPQA.cxx:3807
AliAnalysisTaskHFJetIPQA::fPCJet
std::vector< Double_t > fPCJet
Definition: AliAnalysisTaskHFJetIPQA.h:603
AliEmcalJet::ClosestJetDistance
Double_t ClosestJetDistance() const
Definition: AliEmcalJet.h:342
AliAnalysisTaskHFJetIPQA::IsJetTaggedJetProb
Bool_t IsJetTaggedJetProb(double thresProb=0.90)
IsJetTaggedJetProb.
Definition: AliAnalysisTaskHFJetIPQA.cxx:3855
AliAnalysisTaskHFJetIPQA::GetHist2D
TH2D * GetHist2D(const char *name)
Definition: AliAnalysisTaskHFJetIPQA.h:467
AliAnalysisTaskHFJetIPQA::SetMatchingLevel
void SetMatchingLevel(AliEmcalJet *jet1, AliEmcalJet *jet2, Int_t matching=0)
SetMatchingLevel.
Definition: AliAnalysisTaskHFJetIPQA.cxx:3739
AliEmcalJet::Eta
Double_t Eta() const
Definition: AliEmcalJet.h:121
AliAnalysisTaskHFJetIPQA::FillCorrelations
void FillCorrelations(bool bn[3], double v[3], double jetpt)
Definition: AliAnalysisTaskHFJetIPQA.cxx:2062
AliAnalysisTaskHFJetIPQA::fh1V0CounterCentLambda
TH1D * fh1V0CounterCentLambda
number of K0s candidates after various cuts
Definition: AliAnalysisTaskHFJetIPQA.h:574
AliEmcalContainer::GetMaxEta
Double_t GetMaxEta() const
Definition: AliEmcalContainer.h:202
AliAnalysisTaskHFJetIPQA::GetMCTruth
int GetMCTruth(AliAODTrack *track, int &motherpdg)
Definition: AliAnalysisTaskHFJetIPQA.cxx:491