AliPhysics  31210d0 (31210d0)
AliAnalysisTaskSED0MassNonPromptFraction.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 
17 /* $Id$ */
18 
20 //
21 // AliAnalysisTaskSE for D0 candidates invariant mass histogram
22 // and comparison with the MC truth and cut variables distributions.
23 //
24 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
25 // Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
26 // Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign)
27 // Jeremy Wilkinson, jwilkinson@physi.uni-heidelberg.de (weighted Bayesian
29 
30 #include <Riostream.h>
31 #include <TClonesArray.h>
32 #include <TCanvas.h>
33 #include <TNtuple.h>
34 #include <TTree.h>
35 #include <TList.h>
36 #include <TH1F.h>
37 #include <TH2F.h>
38 #include <TDatabasePDG.h>
39 #include <THnSparse.h>
40 #include "AliVertexingHFUtils.h"
41 #include <TVector3.h>
42 #include <TVector2.h>
43 #include <AliAnalysisDataSlot.h>
44 #include <AliAnalysisDataContainer.h>
45 #include "AliAnalysisManager.h"
46 #include "AliESDtrack.h"
47 #include "AliVertexerTracks.h"
48 #include "AliAODHandler.h"
49 #include "AliAODEvent.h"
50 #include "AliAODVertex.h"
51 #include "AliAODTrack.h"
52 #include "AliAODMCHeader.h"
53 #include "AliAODMCParticle.h"
54 #include "AliAODRecoDecayHF2Prong.h"
55 #include "AliAODRecoCascadeHF.h"
56 #include "AliAnalysisVertexingHF.h"
57 #include "AliAnalysisTaskSE.h"
58 #include "AliAnalysisTaskSED0MassNonPromptFraction.h"
59 #include "AliNormalizationCounter.h"
60 
61 using std::cout;
62 using std::endl;
63 
65 ClassImp(AliAnalysisTaskSED0MassNonPromptFraction);
67 
68 //________________________________________________________________________
69 AliAnalysisTaskSED0MassNonPromptFraction::AliAnalysisTaskSED0MassNonPromptFraction():
70 AliAnalysisTaskSE(),
71  fOutputMass(0),
72  fOutputMassPt(0),
73  fOutputMassY(0),
74  fDistr(0),
75  fNentries(0),
76  fCuts(0),
77  fArray(0),
78  fReadMC(0),
79  fCutOnDistr(0),
80  fUsePid4Distr(0),
81  fCounter(0),
82  fNPtBins(1),
83  fLsNormalization(1.),
84  fFillOnlyD0D0bar(0),
85  fDaughterTracks(),
86  fIsSelectedCandidate(0),
87  fFillVarHists(kTRUE),
88  fFill(0),
89  fSys(0),
90  fIsRejectSDDClusters(0),
91  fFillPtHist(kTRUE),
92  fFillYHist(kFALSE),
93  fFillImpParHist(kFALSE),
94  fUseSelectionBit(kTRUE),
95  fAODProtection(1),
96  fFillAfterCuts(0),
97  fWriteVariableTree(kFALSE),
98  fVariablesTree(0),
99  fCandidateVariables(),
100  fPIDCheck(kFALSE),
101  fDrawDetSignal(kFALSE),
102  fUseQuarkTagInKine(kTRUE),
103  fDetSignal(0),
104  fhMultVZEROTPCoutTrackCorrNoCut(0x0),
105  fhMultVZEROTPCoutTrackCorr(0x0),
106  fEnablePileupRejVZEROTPCout(kFALSE)
107 {
109  for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
110 
111 }
112 
113 //________________________________________________________________________
114 AliAnalysisTaskSED0MassNonPromptFraction::AliAnalysisTaskSED0MassNonPromptFraction(const char *name,AliRDHFCutsD0toKpi* cuts):
115  AliAnalysisTaskSE(name),
116  fOutputMass(0),
117  fOutputMassPt(0),
118  fOutputMassY(0),
119  fDistr(0),
120  fNentries(0),
121  fCuts(0),
122  fArray(0),
123  fReadMC(0),
124  fCutOnDistr(0),
125  fUsePid4Distr(0),
126  fCounter(0),
127  fNPtBins(1),
128  fLsNormalization(1.),
129  fFillOnlyD0D0bar(0),
130  fDaughterTracks(),
131  fIsSelectedCandidate(0),
132  fFillVarHists(kTRUE),
133  fFill(0),
134  fSys(0),
135  fIsRejectSDDClusters(0),
136  fFillPtHist(kTRUE),
137  fFillYHist(kFALSE),
138  fFillImpParHist(kFALSE),
139  fUseSelectionBit(kTRUE),
140  fAODProtection(1),
141  fFillAfterCuts(0),
142  fWriteVariableTree(kFALSE),
143  fVariablesTree(0),
144  fCandidateVariables(),
145  fPIDCheck(kFALSE),
146  fDrawDetSignal(kFALSE),
147  fUseQuarkTagInKine(kTRUE),
148  fDetSignal(0),
149  fhMultVZEROTPCoutTrackCorrNoCut(0x0),
150  fhMultVZEROTPCoutTrackCorr(0x0),
151  fEnablePileupRejVZEROTPCout(kFALSE)
152 {
154 
155  fNPtBins=cuts->GetNPtBins();
156 
157  fCuts=cuts;
158  for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
159 
160  // Output slot #1 writes into a TList container (mass with cuts)
161  DefineOutput(1,TList::Class()); //My private output
162  // Output slot #2 writes into a TList container (distributions)
163  DefineOutput(2,TList::Class()); //My private output
164  // Output slot #3 writes into a TH1F container (number of events)
165  DefineOutput(3,TH1F::Class()); //My private output
166  // Output slot #4 writes into a TList container (cuts)
167  DefineOutput(4,AliRDHFCutsD0toKpi::Class()); //My private output
168  // Output slot #5 writes Normalization Counter
169  DefineOutput(5,AliNormalizationCounter::Class());
170  // Output slot #6 stores the mass vs pt and impact parameter distributions
171  DefineOutput(6,TList::Class()); //My private output
172  // Output slot #7 keeps a tree of the candidate variables after track selection
173  DefineOutput(7,TTree::Class()); //My private outpu
174  // Output slot #8 writes into a TList container (Detector signals)
175  DefineOutput(8, TList::Class()); //My private output
176  // Output slot #9 stores the mass vs rapidity (y) distributions
177  DefineOutput(9, TList::Class()); //My private output
178 }
179 
180 //________________________________________________________________________
181 AliAnalysisTaskSED0MassNonPromptFraction::~AliAnalysisTaskSED0MassNonPromptFraction()
182 {
183  if (fOutputMass) {
184  delete fOutputMass;
185  fOutputMass = 0;
186  }
187  if (fOutputMassPt) {
188  delete fOutputMassPt;
189  fOutputMassPt = 0;
190  }
191  if (fOutputMassY) {
192  delete fOutputMassY;
193  fOutputMassY = 0;
194  }
195  if (fDistr) {
196  delete fDistr;
197  fDistr = 0;
198  }
199  if (fCuts) {
200  delete fCuts;
201  fCuts = 0;
202  }
203  for(Int_t i=0; i<5; i++){
204  if(fHistMassPtImpParTC[i]) delete fHistMassPtImpParTC[i];
205  }
206  if (fNentries){
207  delete fNentries;
208  fNentries = 0;
209  }
210  if(fCounter){
211  delete fCounter;
212  fCounter=0;
213  }
214  if(fVariablesTree){
215  delete fVariablesTree;
216  fVariablesTree = 0;
217  }
218  if (fDetSignal) {
219  delete fDetSignal;
220  fDetSignal = 0;
221  }
222 
223  if(fhMultVZEROTPCoutTrackCorrNoCut){
224  delete fhMultVZEROTPCoutTrackCorrNoCut;
225  fhMultVZEROTPCoutTrackCorrNoCut = 0;
226  }
227  if (fhMultVZEROTPCoutTrackCorr) {
228  delete fhMultVZEROTPCoutTrackCorr;
229  fhMultVZEROTPCoutTrackCorr = 0;
230  }
231 
232 }
233 
234 //________________________________________________________________________
235 void AliAnalysisTaskSED0MassNonPromptFraction::Init()
236 {
238 
239  if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
240 
241 
242  AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
243  const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
244  copyfCuts->SetName(nameoutput);
245  // Post the data
246  PostData(4,copyfCuts);
247 
248 
249  return;
250 }
251 
252 //________________________________________________________________________
253 void AliAnalysisTaskSED0MassNonPromptFraction::UserCreateOutputObjects()
254 {
255 
257  //
258  if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
259 
260  // Several histograms are more conveniently managed in a TList
261  fOutputMass = new TList();
262  fOutputMass->SetOwner();
263  fOutputMass->SetName("listMass");
264 
265  fOutputMassPt = new TList();
266  fOutputMassPt->SetOwner();
267  fOutputMassPt->SetName("listMassPt");
268 
269  fOutputMassY = new TList();
270  fOutputMassY->SetOwner();
271  fOutputMassY->SetName("listMassY");
272 
273  fDistr = new TList();
274  fDistr->SetOwner();
275  fDistr->SetName("distributionslist");
276 
277  fDetSignal = new TList();
278  fDetSignal->SetOwner();
279  fDetSignal->SetName("detectorsignals");
280 
281  TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
282  TString nameMassPt="", nameSgnPt="", nameBkgPt="", nameRflPt="";
283  TString nameMassY="", nameSgnY="", nameBkgY="", nameRflY="";
284  Int_t nbins2dPt=60; Float_t binInPt=0., binFinPt=30.;
285  Int_t nbins2dY=60; Float_t binInY=-1.5, binFinY=1.5;
286 
287  for(Int_t i=0;i<fCuts->GetNPtBins();i++){
288 
289  nameMass="histMass_";
290  nameMass+=i;
291  nameSgn27="histSgn27_";
292  nameSgn27+=i;
293  nameSgn="histSgn_";
294  nameSgn+=i;
295 
296  nameBkg="histBkg_";
297  nameBkg+=i;
298 
299  nameRfl="histRfl_";
300  nameRfl+=i;
301  nameMassNocutsS="hMassS_";
302  nameMassNocutsS+=i;
303  nameMassNocutsB="hMassB_";
304  nameMassNocutsB+=i;
305 
306  //histograms of cut variable distributions
307  if(fFillVarHists){
308  if(fReadMC){
309 
310  namedistr="hNclsD0vsptS_";
311  namedistr+=i;
312  TH2F *hNclsD0vsptS = new TH2F(namedistr.Data(),"N cls distrubution [S];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
313  namedistr="hNclsD0barvsptS_";
314  namedistr+=i;
315  TH2F *hNclsD0barvsptS = new TH2F(namedistr.Data(),"N cls distrubution [S];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
316 
317  namedistr="hNITSpointsD0vsptS_";
318  namedistr+=i;
319  TH2F *hNITSpointsD0vsptS = new TH2F(namedistr.Data(),"N ITS points distrubution [S];p_{T} [GeV/c];N points",200,0.,20.,7,0.,7.);
320 
321  namedistr="hNSPDpointsD0S_";
322  namedistr+=i;
323  TH1I *hNSPDpointsD0S = new TH1I(namedistr.Data(),"N SPD points distrubution [S]; ;N tracks",4,0,4);
324  hNSPDpointsD0S->GetXaxis()->SetBinLabel(1, "no SPD");
325  hNSPDpointsD0S->GetXaxis()->SetBinLabel(2, "kOnlyFirst");
326  hNSPDpointsD0S->GetXaxis()->SetBinLabel(3, "kOnlySecond");
327  hNSPDpointsD0S->GetXaxis()->SetBinLabel(4, "kBoth");
328 
329  namedistr="hptD0S_";
330  namedistr+=i;
331  TH1F *hptD0S = new TH1F(namedistr.Data(), "p_{T} distribution [S];p_{T} [GeV/c]",200,0.,20.);
332  namedistr="hptD0barS_";
333  namedistr+=i;
334  TH1F *hptD0barS = new TH1F(namedistr.Data(), "p_{T} distribution [S];p_{T} [GeV/c]",200,0.,20.);
335 
336  namedistr="hphiD0S_";
337  namedistr+=i;
338  TH1F *hphiD0S = new TH1F(namedistr.Data(), "#phi distribution [S];#phi [rad]",100,0.,2*TMath::Pi());
339  namedistr="hphiD0barS_";
340  namedistr+=i;
341  TH1F *hphiD0barS = new TH1F(namedistr.Data(), "#phi distribution [S];#phi [rad]",100,0.,2*TMath::Pi());
342 
343 
344  namedistr="hetaphiD0candidateS_";
345  namedistr+=i;
346  TH2F *hetaphiD0candidateS = new TH2F(namedistr.Data(), "D^{0} candidates #eta #phi distribution [S];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
347  namedistr="hetaphiD0barcandidateS_";
348  namedistr+=i;
349  TH2F *hetaphiD0barcandidateS = new TH2F(namedistr.Data(), "anti-D^{0} candidates #eta #phi distribution [S];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
350 
351  namedistr="hetaphiD0candidatesignalregionS_";
352  namedistr+=i;
353  TH2F *hetaphiD0candidatesignalregionS = new TH2F(namedistr.Data(), "D^{0} candidates #eta #phi distribution [S] [mass cut];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
354  namedistr="hetaphiD0barcandidatesignalregionS_";
355  namedistr+=i;
356  TH2F *hetaphiD0barcandidatesignalregionS = new TH2F(namedistr.Data(), "anti-D^{0} candidates #eta #phi distribution [S] [mass cut];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
357 
358  // dca
359  namedistr="hdcaS_";
360  namedistr+=i;
361  TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
362  // impact parameter
363  namedistr="hd0piS_";
364  namedistr+=i;
365  TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
366 
367  namedistr="hd0KS_";
368  namedistr+=i;
369  TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
370 
371  namedistr="hd0d0S_";
372  namedistr+=i;
373  TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
374 
375  //decay lenght
376  namedistr="hdeclS_";
377  namedistr+=i;
378  TH1F *hdeclengthS = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm]",200,0,0.015);
379 
380  namedistr="hnormdeclS_";
381  namedistr+=i;
382  TH1F *hnormdeclengthS = new TH1F(namedistr.Data(), "Normalized Decay Length^{2} distribution;(Decay Length/Err)^{2} ",200,0,12.);
383 
384  namedistr="hdeclxyS_";
385  namedistr+=i;
386  TH1F* hdeclxyS=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
387 
388  namedistr="hnormdeclxyS_";
389  namedistr+=i;
390  TH1F* hnormdeclxyS=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
391 
392 
393  namedistr="hCorrCosX_";
394  namedistr+=i;
395  TH2F* hCorrCosPpdX =new TH2F(namedistr.Data(), "Correlation CosPointAngle vs Pseudo-proper decay-length",200, -1, 1, 1000,-1.5,1.5);
396 
397 
398 
399  namedistr="hCorrCosXtrue_";
400  namedistr+=i;
401  TH2F* hCorrCosPpdXtrue =new TH2F(namedistr.Data(),"Correlation CosPointAngle vs Pseudo-proper decay-length",200, -1, 1, 1000,-1.5,1.5);
402 
403  namedistr="hCorrCosX_prim_";
404  namedistr+=i;
405  TH2F* hCorrCosPpdX_prim =new TH2F(namedistr.Data(),"Correlation CosPointAngle vs Pseudo-proper decay-length",200, -1, 1, 1000,-1.5,1.5);
406 
407  namedistr="hCorrCosX_sec_";
408  namedistr+=i;
409  TH2F* hCorrCosPpdX_sec =new TH2F(namedistr.Data(),"Correlation CosPointAngle vs Pseudo-proper decay-length",200, -1, 1, 1000,-1.5,1.5);
410 
411 
412 
413 
414  namedistr="hCorrDistX_";
415  namedistr+=i;
416  TH2F* hCorrDistPpdX =new TH2F(namedistr.Data(),"Distance vs pseudo-proper decay-length correlation",1000, 0, 1.5, 1000,-1.5,1.5);
417 
418 
419  namedistr="hCorrDistXtrue_";
420  namedistr+=i;
421  TH2F* hCorrDistPpdXtrue =new TH2F(namedistr.Data(),"Distance vs pseudo-proper decay-length correlation",1000, 0, 1.5, 1000,-1.5,1.5);
422 
423 
424  namedistr="hCorrDistX_prim_";
425  namedistr+=i;
426  TH2F* hCorrDistPpdX_prim =new TH2F(namedistr.Data(),"Distance vs pseudo-proper decay-length correlation",1000, 0, 1.5, 1000,-1.5,1.5);
427 
428 
429  namedistr="hCorrDistX_sec_";
430  namedistr+=i;
431  TH2F* hCorrDistPpdX_sec =new TH2F(namedistr.Data(),"Distance vs pseudo-proper decay-length correlation",1000, 0, 1.5, 1000,-1.5,1.5);
432 
433  namedistr="hCorrNormLxyXtrue_";
434  namedistr+=i;
435  TH2F* hcorrnormLxyXtrue =new TH2F(namedistr.Data(),"Normalized Lxy vs pseudo-proper decay-length correlation",1000, 0, 1.5, 1000,-1.5,1.5);
436 
437 
438  namedistr="hCorrNormLxyX_prim_";
439  namedistr+=i;
440  TH2F* hcorrnormLxyX_prim =new TH2F(namedistr.Data(),"Normalized Lxy vs pseudo-proper decay-length correlatio",1000, 0, 1.5, 1000,-1.5,1.5);
441 
442 
443  namedistr="hCorrNormLxyX_sec_";
444  namedistr+=i;
445  TH2F* hcorrnormLxyX_sec =new TH2F(namedistr.Data(),"Normalized Lxy vs pseudo-proper decay-length correlatio",1000, 0, 1.5, 1000,-1.5,1.5);
446 
447 
448 
449 
450  namedistr="hpseudoProperDLRec_";
451  namedistr+=i;
452  TH1F* hpseudoProperDecLen =new TH1F(namedistr.Data(),"Pseudo-proper decay-length (reconstructed)",1000,-0.1,0.5);
453 
454 
455  namedistr="hpseudoProperGENERATE_";
456  namedistr+=i;
457  TH1F* hpseudoProperGENERATED =new TH1F(namedistr.Data(),"Pseudo proper decay length (MC-generated)",1000,-1.5,1.5);
458 
459 
460  namedistr="hpseudoProperPrimGenerate_";
461  namedistr+=i;
462  TH1F* hXPrimGenerate =new TH1F(namedistr.Data(),"Pseudo proper decay length (MC-generated)",1500,-1.5,1.5);
463 
464  namedistr="hpseudoProperSecGenerate_";
465  namedistr+=i;
466  TH1F* hXSecGenerate =new TH1F(namedistr.Data(),"Pseudo proper decay length (MC-generated)",1000,-1.5,1.5);
467 
468 
469 
470  namedistr="hGenMatchRec_";
471  namedistr+=i;
472  TH1F* hGenReconstructed =new TH1F(namedistr.Data(),"Pseudo proper decay length (MC-generated matched reconstructed)",1000,-1.5,1.5);
473 
474 
475  namedistr="hGenMatchRecPrim_";
476  namedistr+=i;
477  TH1F* hGenRecPrimaries =new TH1F(namedistr.Data(),"Pseudo proper decay length (MC-generated matched reconstructed)",2000,-1.5,1.5);
478 
479  namedistr="hGenMatchRecSec_";
480  namedistr+=i;
481  TH1F* hGenRecSecondaries =new TH1F(namedistr.Data(),"Pseudo proper decay length (MC-generated matched reconstructed)",1000,-1.5,1.5);
482 
483  namedistr="hresidualX_";
484  namedistr+=i;
485  TH1F* hresidualX =new TH1F(namedistr.Data(),"Residual x - x_{gen}",1000,-1.5,1.5);
486 
487 
488 
489  namedistr="hpseudoProperDLRecSec_";
490  namedistr+=i;
491  TH1F* hpseudoProperDecLen2 =new TH1F(namedistr.Data(),"Pseudo proper decay length (reconstructed)",1000,-0.1,0.5);
492 
493  namedistr="hresidualXSec_";
494  namedistr+=i;
495  TH1F* hresidualpseudoXsec =new TH1F(namedistr.Data(),"Residual x - x_{gen}",1000,-1.5,1.5);
496 
497 
498  namedistr="hresidual_Primary_";
499  namedistr+=i;
500  TH1F* hresidualXprim =new TH1F(namedistr.Data(),"Residual x - x_{gen}",1000,-1.5,1.5);
501 
502  namedistr="hpseudoProperDLRecPrim_";
503  namedistr+=i;
504  TH1F* hpseudoProperDecLenprim =new TH1F(namedistr.Data(),"Pseudo proper decay length (reconstructed)",1000,-0.1,0.5);
505 
506 
507 
508  namedistr="hctauS_";
509  namedistr+=i;
510  TH1F* hdctauS = new TH1F(namedistr.Data(),"Proper Decay time distribution (zeta component); ctau [cm]",200,0,0.5);
511 
512 
513 
514 
515  namedistr="hctauD0_";
516  namedistr+=i;
517  TH1F* hctauD0meson1=new TH1F(namedistr.Data(),"CtauD0 [cm]",200,0,0.5);
518 
519  namedistr="hctauD0S_";
520  namedistr+=i;
521  TH1F* hctauD0meson1S=new TH1F(namedistr.Data(),"CtauD0 [cm]",200,0,0.5);
522 
523 
524 
525 
526  namedistr="hctauXYctauZ_";
527  namedistr+=i;
528  TH2F* hctauXYvsZ=new TH2F(namedistr.Data(),"ctauD0-XY plane vs Ctau-Z [cm]",200,0,0.5,200,0,0.05);
529 
530 
531  namedistr="hctauXYctauZS_";
532  namedistr+=i;
533  TH2F* hctauXYvsZS=new TH2F(namedistr.Data(),"ctauD0-XY plane vs Ctau-Z [cm]",200,0,0.5,200,0,0.05);
534 
535 
536 
537 
538  namedistr="hctauD0prim_";
539  namedistr+=i;
540  TH1F* hctauD0primary=new TH1F(namedistr.Data(),"ctauD0 [cm]",200,0,0.5);
541 
542 
543 
544 
545  namedistr="hctauD0sec_";
546  namedistr+=i;
547  TH1F* hctauD0secondary =new TH1F(namedistr.Data(),"ctauD0 [cm]",200,0,0.5);
548 
549 
550 
551 
552  namedistr="hdeclxyd0d0S_";
553  namedistr+=i;
554  TH2F* hdeclxyd0d0S=new TH2F(namedistr.Data(),"Correlation decay Length XY [cm] vs d_{0}#times d_{0} [cm^{2}]",200,0,0.15,200,-0.001,0.001);
555 
556  namedistr="hnormdeclxyd0d0S_";
557  namedistr+=i;
558  TH2F* hnormdeclxyd0d0S=new TH2F(namedistr.Data(),"Correlation normalized decay Length XY - d_{0}#times d_{0};Decay Length XY/Err;d_{0}#times d_{0}[cm^{2}]",200,0,6,200,-0.001,0.001);
559 
560 
561 
562  //correlation ctau XY vs d0 D0
563  namedistr="hCorrXd0D0_";
564  namedistr+=i;
565  TH2F* hcorrXd0D0 = new TH2F(namedistr.Data(),"Correlation Proper Decay time XY - D0 d_{0} ; Proper decay time XY [cm]; d_{0}[cm]",200,-1.5, 1.5, 1000,-1000,1000);
566 
567 
568  namedistr="hCorrXd0D0true_";
569  namedistr+=i;
570  TH2F* hcorrXd0D0true = new TH2F(namedistr.Data(),"Correlation Proper Decay time XY - D0 d_{0} ; Proper decay time XY [cm]; d_{0}[cm]",200,-1.5, 1.5, 1000,-1000,1000);
571 
572 
573  namedistr="hCorrXd0D0_sec_";
574  namedistr+=i;
575  TH2F* hcorrXd0D0sec = new TH2F(namedistr.Data(),"Correlation Proper Decay time XY - D0 d_{0} ; Proper decay time XY [cm]; d_{0}[cm]",200,-1.5, 1.5, 1000,-1000,1000);
576 
577 
578  namedistr="hCorrXd0D0_prim_";
579  namedistr+=i;
580  TH2F* hcorrXd0D0prim = new TH2F(namedistr.Data(),"Correlation Proper Decay time XY - D0 d_{0} ; Proper decay time XY [cm]; d_{0}[cm]",200,-1.5, 1.5, 1000,-1000,1000);
581 
582 
583 
584  namedistr="hImpParXYtrue_";
585  namedistr+=i;
586  TH1F* hImpactParXYtrue = new TH1F(namedistr.Data()," D0 d_{0} ; Proper decay time XY [cm]; d_{0}[cm]", 1000,-1000,1000);
587 
588 
589  namedistr="hImpParXY_sec_";
590  namedistr+=i;
591  TH1F* hImpParXYsec = new TH1F(namedistr.Data(),"D0 d_{0} ; Proper decay time XY [cm]; d_{0}[cm]",1000,-1000,1000);
592 
593 
594  namedistr="hImpParXY_prim_";
595  namedistr+=i;
596  TH1F* hImpParXYprim = new TH1F(namedistr.Data(),"D0 d_{0} ; Proper decay time XY [cm]; d_{0}[cm]",1000,-1000,1000);
597 
598 
599 
600 
601  // costhetapoint
602  namedistr="hcosthetapointS_";
603  namedistr+=i;
604  TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
605 
606  namedistr="hcosthetapointxyS_";
607  namedistr+=i;
608  TH1F *hcosthetapointxyS = new TH1F(namedistr.Data(), "cos#theta_{Point} XYdistribution;cos#theta_{Point}",300,0.,1.);
609 
610  TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.5648,2.1648); //range (MD0-300MeV, mD0 + 300MeV)
611  tmpMS->Sumw2();
612 
613  fDistr->Add(hNclsD0vsptS);
614  fDistr->Add(hNclsD0barvsptS);
615  fDistr->Add(hNITSpointsD0vsptS);
616  fDistr->Add(hNSPDpointsD0S);
617  fDistr->Add(hptD0S);
618  fDistr->Add(hphiD0S);
619  fDistr->Add(hptD0barS);
620  fDistr->Add(hphiD0barS);
621  fDistr->Add(hetaphiD0candidateS);
622  fDistr->Add(hetaphiD0candidatesignalregionS);
623  fDistr->Add(hetaphiD0barcandidateS);
624  fDistr->Add(hetaphiD0barcandidatesignalregionS);
625  fDistr->Add(hdcaS);
626  fDistr->Add(hd0piS);
627  fDistr->Add(hd0KS);
628  fDistr->Add(hd0d0S);
629  fDistr->Add(hcosthetapointS);
630  fDistr->Add(hcosthetapointxyS);
631  fDistr->Add(hdeclengthS);
632  fDistr->Add(hnormdeclengthS);
633  fDistr->Add(hdeclxyS);
634  fDistr->Add(hdctauS);
635  fDistr->Add(hctauD0meson1);
636  fDistr->Add(hctauD0meson1S);
637  fDistr->Add(hctauXYvsZ);
638  fDistr->Add(hctauXYvsZS);
639  fDistr->Add(hctauD0primary);
640  fDistr->Add(hctauD0secondary);
641  fDistr->Add(hnormdeclxyS);
642  fDistr->Add(hCorrCosPpdX);
643  fDistr->Add(hCorrCosPpdXtrue);
644  fDistr->Add(hCorrCosPpdX_prim);
645  fDistr->Add(hCorrCosPpdX_sec);
646  fDistr->Add(hCorrDistPpdX);
647  fDistr->Add(hCorrDistPpdXtrue);
648  fDistr->Add(hCorrDistPpdX_prim);
649  fDistr->Add(hCorrDistPpdX_sec);
650  fDistr->Add(hcorrnormLxyXtrue);
651  fDistr->Add(hcorrnormLxyX_prim);
652  fDistr->Add(hcorrnormLxyX_sec);
653  fDistr->Add(hpseudoProperDecLen);
654  fDistr->Add(hpseudoProperGENERATED);
655  fDistr->Add(hXSecGenerate);
656  fDistr->Add(hXPrimGenerate);
657  fDistr->Add(hGenReconstructed);
658  fDistr->Add(hGenRecPrimaries);
659  fDistr->Add(hGenRecSecondaries);
660  fDistr->Add(hpseudoProperDecLen2);
661  fDistr->Add(hresidualX);
662  fDistr->Add(hresidualpseudoXsec);
663  fDistr->Add(hresidualXprim);
664  fDistr->Add(hpseudoProperDecLenprim);
665  fDistr->Add(hdeclxyd0d0S);
666  fDistr->Add(hcorrXd0D0);
667  fDistr->Add(hcorrXd0D0true);
668  fDistr->Add(hcorrXd0D0sec);
669  fDistr->Add(hcorrXd0D0prim);
670  fDistr->Add(hImpactParXYtrue);
671  fDistr->Add(hImpParXYsec);
672  fDistr->Add(hImpParXYprim);
673  fDistr->Add(hnormdeclxyd0d0S);
674  fDistr->Add(tmpMS);
675  }
676 
677 
678  //Ncls, phi, pt distrubutions
679 
680  namedistr="hNclsD0vsptB_";
681  namedistr+=i;
682  TH2F *hNclsD0vsptB = new TH2F(namedistr.Data(),"N cls distrubution [B];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
683  namedistr="hNclsD0barvsptB_";
684  namedistr+=i;
685  TH2F *hNclsD0barvsptB = new TH2F(namedistr.Data(),"N cls distrubution [B];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
686 
687  namedistr="hNITSpointsD0vsptB_";
688  namedistr+=i;
689  TH2F *hNITSpointsD0vsptB = new TH2F(namedistr.Data(),"N ITS points distrubution [B];p_{T} [GeV/c];N points",200,0.,20.,7,0.,7.);
690 
691  namedistr="hNSPDpointsD0B_";
692  namedistr+=i;
693  TH1I *hNSPDpointsD0B = new TH1I(namedistr.Data(),"N SPD points distrubution [B]; ;N tracks",4,0,4);
694  hNSPDpointsD0B->GetXaxis()->SetBinLabel(1, "no SPD");
695  hNSPDpointsD0B->GetXaxis()->SetBinLabel(2, "kOnlyFirst");
696  hNSPDpointsD0B->GetXaxis()->SetBinLabel(3, "kOnlySecond");
697  hNSPDpointsD0B->GetXaxis()->SetBinLabel(4, "kBoth");
698 
699  namedistr="hptD0B_";
700  namedistr+=i;
701  TH1F *hptD0B = new TH1F(namedistr.Data(), "p_{T} distribution [B];p_{T} [GeV/c]",200,0.,20.);
702  namedistr="hptD0barB_";
703  namedistr+=i;
704  TH1F *hptD0barB = new TH1F(namedistr.Data(), "p_{T} distribution [B];p_{T} [GeV/c]",200,0.,20.);
705 
706  namedistr="hphiD0B_";
707  namedistr+=i;
708  TH1F *hphiD0B = new TH1F(namedistr.Data(), "#phi distribution [B];#phi [rad]",100,0.,2*TMath::Pi());
709  namedistr="hphiD0barB_";
710  namedistr+=i;
711  TH1F *hphiD0barB = new TH1F(namedistr.Data(), "#phi distribution [B];#phi [rad]",100,0.,2*TMath::Pi());
712 
713  namedistr="hetaphiD0candidateB_";
714  namedistr+=i;
715  TH2F *hetaphiD0candidateB = new TH2F(namedistr.Data(), "D^{0} candidates #eta #phi distribution [B];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
716  namedistr="hetaphiD0barcandidateB_";
717  namedistr+=i;
718  TH2F *hetaphiD0barcandidateB = new TH2F(namedistr.Data(), "anti-D^{0} candidates #eta #phi distribution [B];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
719 
720  namedistr="hetaphiD0candidatesignalregionB_";
721  namedistr+=i;
722  TH2F *hetaphiD0candidatesignalregionB = new TH2F(namedistr.Data(), "D^{0} candidates #eta #phi distribution [B] [mass cut];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
723  namedistr="hetaphiD0barcandidatesignalregionB_";
724  namedistr+=i;
725  TH2F *hetaphiD0barcandidatesignalregionB = new TH2F(namedistr.Data(), "anti-D^{0} candidates #eta #phi distribution [B] [mass cut];#eta;#phi [rad]",100, -1.5, 1.5, 100, 0.,2*TMath::Pi());
726 
727  // dca
728  namedistr="hdcaB_";
729  namedistr+=i;
730  TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
731 
732  // impact parameter
733  namedistr="hd0B_";
734  namedistr+=i;
735  TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
736 
737  namedistr="hd0d0B_";
738  namedistr+=i;
739  TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
740 
741  //decay lenght
742  namedistr="hdeclB_";
743  namedistr+=i;
744  TH1F *hdeclengthB = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm^{2}]",200,0,0.015);
745 
746  namedistr="hnormdeclB_";
747  namedistr+=i;
748  TH1F *hnormdeclengthB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;(Decay Length/Err)^{2} ",200,0,12.);
749 
750 
751  namedistr="hctau_";
752  namedistr+=i;
753  TH1F *hctauD = new TH1F(namedistr.Data(), "Ctau D meson;(ctau);[cm] ",200,0,0.15);
754 
755 
756 
757 
758 
759  namedistr="hdeclxyB_";
760  namedistr+=i;
761  TH1F* hdeclxyB=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
762  namedistr="hnormdeclxyB_";
763  namedistr+=i;
764  TH1F* hnormdeclxyB=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
765 
766  namedistr="hdeclxyd0d0B_";
767  namedistr+=i;
768  TH2F* hdeclxyd0d0B=new TH2F(namedistr.Data(),"Correlation decay Length XY - d_{0}#times d_{0};Decay Length XY [cm];d_{0}#times d_{0}[cm^{2}]",200,0,0.15,200,-0.001,0.001);
769 
770  namedistr="hnormdeclxyd0d0B_";
771  namedistr+=i;
772  TH2F* hnormdeclxyd0d0B=new TH2F(namedistr.Data(),"Correlation normalized decay Length XY - d_{0}#times d_{0};Decay Length XY/Err;d_{0}#times d_{0}[cm^{2}]",200,0,6,200,-0.001,0.001);
773 
774  // costhetapoint
775  namedistr="hcosthetapointB_";
776  namedistr+=i;
777  TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
778 
779  namedistr="hcosthetapointxyB_";
780  namedistr+=i;
781  TH1F *hcosthetapointxyB = new TH1F(namedistr.Data(), "cos#theta_{Point} XY distribution;cos#theta_{Point} XY",300,0.,1.);
782 
783  TH1F* tmpMB = new TH1F(nameMassNocutsB.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,1.5648,2.1648); //range (MD0-300MeV, mD0 + 300MeV)
784  tmpMB->Sumw2();
785 
786 
787  fDistr->Add(hNclsD0vsptB);
788  fDistr->Add(hNclsD0barvsptB);
789  fDistr->Add(hNITSpointsD0vsptB);
790  fDistr->Add(hNSPDpointsD0B);
791  fDistr->Add(hptD0B);
792  fDistr->Add(hphiD0B);
793  fDistr->Add(hptD0barB);
794  fDistr->Add(hphiD0barB);
795  fDistr->Add(hetaphiD0candidateB);
796  fDistr->Add(hetaphiD0candidatesignalregionB);
797  fDistr->Add(hetaphiD0barcandidateB);
798  fDistr->Add(hetaphiD0barcandidatesignalregionB);
799  fDistr->Add(hdcaB);
800  fDistr->Add(hd0B);
801  fDistr->Add(hd0d0B);
802  fDistr->Add(hcosthetapointB);
803  fDistr->Add(hcosthetapointxyB);
804  fDistr->Add(hdeclengthB);
805  fDistr->Add(hnormdeclengthB);
806  fDistr->Add(hctauD);
807  fDistr->Add(hdeclxyB);
808  fDistr->Add(hnormdeclxyB);
809  fDistr->Add(hdeclxyd0d0B);
810  fDistr->Add(hnormdeclxyd0d0B);
811  fDistr->Add(tmpMB);
812 
813  //histograms filled only when the secondary vertex is recalculated w/o the daughter tracks (as requested in the cut object)
814 
815  if(fCuts->GetIsPrimaryWithoutDaughters()){
816  if(fReadMC){
817  namedistr="hd0vpiS_";
818  namedistr+=i;
819  TH1F *hd0vpiS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions)(vtx w/o these tracks);d0(#pi) [cm]",200,-0.1,0.1);
820 
821  namedistr="hd0vKS_";
822  namedistr+=i;
823  TH1F *hd0vKS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons) (vtx w/o these tracks);d0(K) [cm]",200,-0.1,0.1);
824 
825  namedistr="hd0d0vS_";
826  namedistr+=i;
827  TH1F *hd0d0vS = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution (vtx w/o these tracks);d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
828 
829  namedistr="hdeclvS_";
830  namedistr+=i;
831  TH1F *hdeclengthvS = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
832 
833  namedistr="hnormdeclvS_";
834  namedistr+=i;
835  TH1F *hnormdeclengthvS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
836 
837  fDistr->Add(hd0vpiS);
838  fDistr->Add(hd0vKS);
839  fDistr->Add(hd0d0vS);
840  fDistr->Add(hdeclengthvS);
841  fDistr->Add(hnormdeclengthvS);
842 
843  }
844 
845  namedistr="hd0vmoresB_";
846  namedistr+=i;
847  TH1F *hd0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
848 
849  namedistr="hd0d0vmoresB_";
850  namedistr+=i;
851  TH1F *hd0d0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
852 
853 
854  namedistr="hd0vB_";
855  namedistr+=i;
856  TH1F *hd0vB = new TH1F(namedistr.Data(), "Impact parameter distribution (vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
857 
858  namedistr="hd0vp0B_";
859  namedistr+=i;
860  TH1F *hd0vp0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong + ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
861 
862  namedistr="hd0vp1B_";
863  namedistr+=i;
864  TH1F *hd0vp1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong - ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
865 
866  namedistr="hd0d0vB_";
867  namedistr+=i;
868  TH1F *hd0d0vB = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution (vtx w/o these tracks);d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
869 
870  namedistr="hdeclvB_";
871  namedistr+=i;
872  TH1F *hdeclengthvB = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
873 
874  namedistr="hnormdeclvB_";
875  namedistr+=i;
876  TH1F *hnormdeclengthvB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
877 
878  fDistr->Add(hd0vB);
879  fDistr->Add(hd0vp0B);
880  fDistr->Add(hd0vp1B);
881  fDistr->Add(hd0vmoresB);
882 
883  fDistr->Add(hd0d0vB);
884  fDistr->Add(hd0d0vmoresB);
885 
886  fDistr->Add(hdeclengthvB);
887 
888  fDistr->Add(hnormdeclengthvB);
889  }
890 
891 
892  }
893  //histograms of invariant mass distributions
894 
895 
896  //MC signal
897  if(fReadMC){
898  TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",600,1.5648,2.1648);
899 
900  TH1F *tmpSl=(TH1F*)tmpSt->Clone();
901  tmpSt->Sumw2();
902  tmpSl->Sumw2();
903 
904  //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
905  TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",600,1.5648,2.1648);
906  //TH1F *tmpRl=(TH1F*)tmpRt->Clone();
907  TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",600,1.5648,2.1648);
908  //TH1F *tmpBl=(TH1F*)tmpBt->Clone();
909  tmpBt->Sumw2();
910  //tmpBl->Sumw2();
911  tmpRt->Sumw2();
912  //tmpRl->Sumw2();
913 
914  fOutputMass->Add(tmpSt);
915  fOutputMass->Add(tmpRt);
916  fOutputMass->Add(tmpBt);
917 
918  }
919 
920  //mass
921  TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",600,1.5648,2.1648);
922  //TH1F *tmpMl=(TH1F*)tmpMt->Clone();
923  tmpMt->Sumw2();
924  //tmpMl->Sumw2();
925  //distribution w/o cuts large range
926  // TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
927 
928  fOutputMass->Add(tmpMt);
929 
930  if(fSys==0){ //histograms filled only in pp to save time in PbPb
931  if(fFillVarHists){
932 
933  if(fReadMC){
934  // pT
935  namedistr="hptpiS_";
936  namedistr+=i;
937  TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
938 
939  namedistr="hptKS_";
940  namedistr+=i;
941  TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
942 
943  // costhetastar
944  namedistr="hcosthetastarS_";
945  namedistr+=i;
946  TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
947 
948  //pT no mass cut
949 
950  namedistr="hptpiSnoMcut_";
951  namedistr+=i;
952  TH1F *hptpiSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
953 
954  namedistr="hptKSnoMcut_";
955  namedistr+=i;
956  TH1F *hptKSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
957 
958  fDistr->Add(hptpiS);
959  fDistr->Add(hptKS);
960  fDistr->Add(hcosthetastarS);
961 
962  fDistr->Add(hptpiSnoMcut);
963  fDistr->Add(hptKSnoMcut);
964 
965  // costhetapoint vs d0 or d0d0
966  namedistr="hcosthpointd0S_";
967  namedistr+=i;
968  TH2F *hcosthpointd0S= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0};cos#theta_{Point};d_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
969  namedistr="hcosthpointd0d0S_";
970  namedistr+=i;
971  TH2F *hcosthpointd0d0S= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0}#timesd_{0};cos#theta_{Point};d_{0}#timesd_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
972 
973  fDistr->Add(hcosthpointd0S);
974  fDistr->Add(hcosthpointd0d0S);
975 
976  //to compare with AliAnalysisTaskCharmFraction
977  TH1F* tmpS27t = new TH1F(nameSgn27.Data(),"D^{0} invariant mass in M(D^{0}) +/- 27 MeV - MC; M [GeV]; Entries",600,1.5648,2.1648);
978  TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
979  tmpS27t->Sumw2();
980  tmpS27l->Sumw2();
981 
982  fOutputMass->Add(tmpS27t);
983  fOutputMass->Add(tmpS27l);
984 
985  }
986 
987  // pT
988  namedistr="hptB_";
989  namedistr+=i;
990  TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
991 
992  // costhetastar
993  namedistr="hcosthetastarB_";
994  namedistr+=i;
995  TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
996 
997  //pT no mass cut
998  namedistr="hptB1prongnoMcut_";
999  namedistr+=i;
1000  TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
1001 
1002  namedistr="hptB2prongsnoMcut_";
1003  namedistr+=i;
1004  TH1F *hptB2pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
1005 
1006  fDistr->Add(hptB);
1007  fDistr->Add(hcosthetastarB);
1008 
1009  fDistr->Add(hptB1pnoMcut);
1010  fDistr->Add(hptB2pnoMcut);
1011 
1012  //impact parameter of negative/positive track
1013  namedistr="hd0p0B_";
1014  namedistr+=i;
1015  TH1F *hd0p0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.1,0.1);
1016 
1017  namedistr="hd0p1B_";
1018  namedistr+=i;
1019  TH1F *hd0p1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong -);d0 [cm]",200,-0.1,0.1);
1020 
1021  //impact parameter corrected for strangeness
1022  namedistr="hd0moresB_";
1023  namedistr+=i;
1024  TH1F *hd0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
1025 
1026  namedistr="hd0d0moresB_";
1027  namedistr+=i;
1028  TH1F *hd0d0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
1029 
1030 
1031  namedistr="hcosthetapointmoresB_";
1032  namedistr+=i;
1033  TH1F *hcosthetapointmoresB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
1034 
1035  // costhetapoint vs d0 or d0d0
1036  namedistr="hcosthpointd0B_";
1037  namedistr+=i;
1038  TH2F *hcosthpointd0B= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0};cos#theta_{Point};d_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
1039 
1040  namedistr="hcosthpointd0d0B_";
1041  namedistr+=i;
1042  TH2F *hcosthpointd0d0B= new TH2F(namedistr.Data(),"Correlation cos#theta_{Point}-d_{0}#timesd_{0};cos#theta_{Point};d_{0}#timesd_{0} [cm^{2}]",200,0,1.,200,-0.001,0.001);
1043 
1044  fDistr->Add(hd0p0B);
1045  fDistr->Add(hd0p1B);
1046 
1047  fDistr->Add(hd0moresB);
1048  fDistr->Add(hd0d0moresB);
1049  fDistr->Add(hcosthetapointmoresB);
1050 
1051 
1052  fDistr->Add(hcosthpointd0B);
1053 
1054 
1055  fDistr->Add(hcosthpointd0d0B);
1056  }
1057 
1058  } //end pp histo
1059  }
1060 
1061 
1062  //for Like sign analysis
1063 
1064  if(fArray==1){
1065  namedistr="hpospair";
1066  TH1F* hpospair=new TH1F(namedistr.Data(),"Number of positive pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
1067  namedistr="hnegpair";
1068  TH1F* hnegpair=new TH1F(namedistr.Data(),"Number of negative pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
1069  fOutputMass->Add(hpospair);
1070  fOutputMass->Add(hnegpair);
1071  }
1072 
1073 
1074  // 2D Pt distributions and impact parameter histograms
1075  if(fFillPtHist) {
1076 
1077  nameMassPt="histMassPt";
1078  nameSgnPt="histSgnPt";
1079  nameBkgPt="histBkgPt";
1080  nameRflPt="histRflPt";
1081 
1082  //MC signal
1083  if(fReadMC){
1084  TH2F* tmpStPt = new TH2F(nameSgnPt.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.16484,nbins2dPt,binInPt,binFinPt);
1085  TH2F *tmpSlPt=(TH2F*)tmpStPt->Clone();
1086  tmpStPt->Sumw2();
1087  tmpSlPt->Sumw2();
1088 
1089  //Reflection: histo filled with D0MassV1 which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
1090  TH2F* tmpRtPt = new TH2F(nameRflPt.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
1091  TH2F* tmpBtPt = new TH2F(nameBkgPt.Data(), "Background invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
1092  tmpBtPt->Sumw2();
1093  tmpRtPt->Sumw2();
1094 
1095  fOutputMassPt->Add(tmpStPt);
1096  fOutputMassPt->Add(tmpRtPt);
1097  fOutputMassPt->Add(tmpBtPt);
1098 
1099  // cout<<endl<<endl<<"***************************************"<<endl;
1100  // cout << " created and added to the list "<< nameSgnPt.Data() <<" "<< tmpStPt <<
1101  // ", "<<nameRflPt.Data() <<" " << tmpRtPt<<", "<<nameBkgPt.Data()<< " " << tmpBtPt <<endl;
1102  // cout<<"***************************************"<<endl<<endl;
1103  }
1104 
1105  TH2F* tmpMtPt = new TH2F(nameMassPt.Data(),"D^{0} invariant mass; M [GeV]; Entries; Pt[GeV/c]",600,1.5648,2.1648,nbins2dPt,binInPt,binFinPt);
1106  tmpMtPt->Sumw2();
1107 
1108  fOutputMassPt->Add(tmpMtPt);
1109  }
1110 
1111  if(fFillImpParHist) CreateImpactParameterHistos();
1112 
1113  // 2D Y distributions
1114 
1115  if(fFillYHist) {
1116  for(Int_t i=0;i<fCuts->GetNPtBins();i++){
1117  nameMassY="histMassY_";
1118  nameMassY+=i;
1119  nameSgnY="histSgnY_";
1120  nameSgnY+=i;
1121  nameBkgY="histBkgY_";
1122  nameBkgY+=i;
1123  nameRflY="histRflY_";
1124  nameRflY+=i;
1125  //MC signal
1126  if(fReadMC){
1127  TH2F* tmpStY = new TH2F(nameSgnY.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries; y",600,1.5648,2.16484,nbins2dY,binInY,binFinY);
1128  tmpStY->Sumw2();
1129  //Reflection: histo filled with D0MassV1 which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
1130  TH2F* tmpRtY = new TH2F(nameRflY.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries; y",600,1.5648,2.1648,nbins2dY,binInY,binFinY);
1131  TH2F* tmpBtY = new TH2F(nameBkgY.Data(), "Background invariant mass - MC; M [GeV]; Entries; y",600,1.5648,2.1648,nbins2dY,binInY,binFinY);
1132  tmpBtY->Sumw2();
1133  tmpRtY->Sumw2();
1134 
1135  fOutputMassY->Add(tmpStY);
1136  fOutputMassY->Add(tmpRtY);
1137  fOutputMassY->Add(tmpBtY);
1138  }
1139  TH2F* tmpMtY = new TH2F(nameMassY.Data(),"D^{0} invariant mass; M [GeV]; Entries; y",600,1.5648,2.1648,nbins2dY,binInY,binFinY);
1140  tmpMtY->Sumw2();
1141  fOutputMassY->Add(tmpMtY);
1142  }
1143  }
1144 
1145 
1146  const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
1147 
1148  fNentries=new TH1F(nameoutput, "Integral(1,2) = number of AODs *** Integral(2,3) = number of candidates selected with cuts *** Integral(3,4) = number of D0 selected with cuts *** Integral(4,5) = events with good vertex *** Integral(5,6) = pt out of bounds", 19,-0.5,18.5);
1149 
1150  fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
1151  fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
1152  if(fReadMC) fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
1153  else fNentries->GetXaxis()->SetBinLabel(3,"Dstar<-D0");
1154  fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
1155  fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
1156  fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
1157  if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
1158  if(fFillVarHists || fPIDCheck){
1159  fNentries->GetXaxis()->SetBinLabel(8,"PID=0");
1160  fNentries->GetXaxis()->SetBinLabel(9,"PID=1");
1161  fNentries->GetXaxis()->SetBinLabel(10,"PID=2");
1162  fNentries->GetXaxis()->SetBinLabel(11,"PID=3");
1163  }
1164  if(fReadMC && fSys==0){
1165  fNentries->GetXaxis()->SetBinLabel(12,"K");
1166  fNentries->GetXaxis()->SetBinLabel(13,"Lambda");
1167  }
1168  fNentries->GetXaxis()->SetBinLabel(14,"Pile-up Rej");
1169  fNentries->GetXaxis()->SetBinLabel(15,"N. of 0SMH");
1170  if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
1171  if(fIsRejectSDDClusters) fNentries->GetXaxis()->SetBinLabel(17,"SDD-Cls Rej");
1172  fNentries->GetXaxis()->SetBinLabel(18,"Phys.Sel.Rej");
1173  fNentries->GetXaxis()->SetBinLabel(19,"D0 failed to be filled");
1174  fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
1175 
1176  fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
1177  fCounter->Init();
1178 
1179  //
1180  // Output slot 7 : tree of the candidate variables after track selection
1181  //
1182  nameoutput = GetOutputSlot(7)->GetContainer()->GetName();
1183  fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
1184  Int_t nVar = 15;
1185  fCandidateVariables = new Double_t [nVar];
1186  TString * fCandidateVariableNames = new TString[nVar];
1187  fCandidateVariableNames[0] = "massD0";
1188  fCandidateVariableNames[1] = "massD0bar";
1189  fCandidateVariableNames[2] = "pt";
1190  //fCandidateVariableNames[3] = "dca";
1191  //fCandidateVariableNames[4] = "costhsD0";
1192  //fCandidateVariableNames[5] = "costhsD0bar";
1193  //fCandidateVariableNames[6] = "ptk";
1194  //fCandidateVariableNames[7] = "ptpi";
1195  //fCandidateVariableNames[8] = "d0k";
1196  //fCandidateVariableNames[9] = "d0pi";
1197  //fCandidateVariableNames[10] = "d0xd0";
1198  //fCandidateVariableNames[11] = "costhp";
1199  //fCandidateVariableNames[12] = "costhpxy";
1200  //fCandidateVariableNames[13] = "lxy";
1201  //fCandidateVariableNames[14] = "specialcuts";
1202  fCandidateVariableNames[3] = "ctau";
1203  fCandidateVariableNames[4] = "PseudoProperDecayLength";
1204  fCandidateVariableNames[5] = "IsSelected";
1205 
1206  for(Int_t ivar=0; ivar<nVar; ivar++){
1207  fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/d",fCandidateVariableNames[ivar].Data()));
1208  }
1209 
1210  //
1211  // Output slot 8 : List for detector response histograms
1212  //
1213  if (fDrawDetSignal) {
1214  TH2F *TOFSigBefPID = new TH2F("TOFSigBefPID", "TOF signal of daughters before PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
1215  TH2F *TOFSigAftPID = new TH2F("TOFSigAftPID", "TOF signal after PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
1216 
1217  TH2F *TPCSigBefPID = new TH2F("TPCSigBefPID", "TPC dE/dx before PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
1218  TH2F *TPCSigAftPID = new TH2F("TPCSigAftPID", "TPC dE/dx after PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
1219 
1220  fDetSignal->Add(TOFSigBefPID);
1221  fDetSignal->Add(TOFSigAftPID);
1222  fDetSignal->Add(TPCSigBefPID);
1223  fDetSignal->Add(TPCSigAftPID);
1224  }
1225 
1226 
1227 
1228  fhMultVZEROTPCoutTrackCorrNoCut = new TH2F("hMultVZEROTPCoutTrackCorrNoCut", ";Tracks with kTPCout on;VZERO multiplicity", 1000, 0., 30000., 1000, 0., 40000.);
1229  fhMultVZEROTPCoutTrackCorr = new TH2F("hMultVZEROTPCoutTrackCorr", ";Tracks with kTPCout on;VZERO multiplicity", 1000, 0., 30000., 1000, 0., 40000.);
1230  fDistr->Add(fhMultVZEROTPCoutTrackCorrNoCut);
1231  fDistr->Add(fhMultVZEROTPCoutTrackCorr);
1232 
1233  // Post the data
1234  PostData(1,fOutputMass);
1235  PostData(2,fDistr);
1236  PostData(3,fNentries);
1237  PostData(5,fCounter);
1238  PostData(6,fOutputMassPt);
1239  PostData(7,fVariablesTree);
1240  PostData(8,fDetSignal);
1241  PostData(9,fOutputMassY);
1242 
1243  return;
1244 }
1245 
1246 //________________________________________________________________________
1247 
1248 
1249 
1250 void AliAnalysisTaskSED0MassNonPromptFraction::UserExec(Option_t */*option*/)
1251 {
1254 
1255  //cuts order
1256  // printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
1257  // printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
1258  // printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
1259  // printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
1260  // printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
1261  // printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
1262  // printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
1263  // printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
1264  // printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
1265 
1266  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
1267 
1268  if(fAODProtection>=0){
1269  // Protection against different number of events in the AOD and deltaAOD
1270  // In case of discrepancy the event is rejected.
1271  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
1272  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
1273  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
1274  fNentries->Fill(21);
1275  return;
1276  }
1277  fNentries->Fill(22);
1278  }
1279 
1280  TString bname;
1281  if(fArray==0){ //D0 candidates
1282  // load D0->Kpi candidates
1283  //cout<<"D0 candidates"<<endl;
1284  bname="D0toKpi";
1285  } else { //LikeSign candidates
1286  // load like sign candidates
1287  //cout<<"LS candidates"<<endl;
1288  bname="LikeSign2Prong";
1289  }
1290  TClonesArray *inputArray=0;
1291  if(!aod && AODEvent() && IsStandardAOD()) {
1292  // In case there is an AOD handler writing a standard AOD, use the AOD
1293  // event in memory rather than the input (ESD) event.
1294  aod = dynamic_cast<AliAODEvent*> (AODEvent());
1295  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
1296  // have to taken from the AOD event hold by the AliAODExtension
1297  AliAODHandler* aodHandler = (AliAODHandler*)
1298  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
1299 
1300 
1301 
1302 
1303  if(aodHandler->GetExtensions()) {
1304  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
1305  AliAODEvent* aodFromExt = ext->GetAOD();
1306  inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
1307  }
1308  } else if(aod) {
1309  inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
1310  }
1311 
1312  if(!inputArray || !aod) {
1313  printf("AliAnalysisTaskSED0MassNonPromptFraction::UserExec: input branch not found!\n");
1314  return;
1315  }
1316  // fix for temporary bug in ESDfilter
1317  // the AODs with null vertex pointer didn't pass the PhysSel
1318  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
1319 
1320 
1321  TClonesArray *mcArray = 0;
1322  AliAODMCHeader *mcHeader = 0;
1323 
1324  if(fReadMC) {
1325  // load MC particles
1326  mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1327  if(!mcArray) {
1328  printf("AliAnalysisTaskSED0MassNonPromptFraction::UserExec: MC particles branch not found!\n");
1329  return;
1330  }
1331 
1332  // load MC header
1333  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1334  if(!mcHeader) {
1335  printf("AliAnalysisTaskSED0MassNonPromptFraction::UserExec: MC header branch not found!\n");
1336  return;
1337  }
1338  }
1339 
1340  Int_t nTPCout=0;
1341  Float_t mTotV0=0;
1342  AliAODVZERO* v0data=(AliAODVZERO*)((AliAODEvent*)aod)->GetVZEROData();
1343  Float_t mTotV0A=v0data->GetMTotV0A();
1344  Float_t mTotV0C=v0data->GetMTotV0C();
1345  mTotV0=mTotV0A+mTotV0C;
1346  Int_t ntracksEv = aod->GetNumberOfTracks();
1347  for(Int_t itrack=0; itrack<ntracksEv; itrack++) { // loop on tacks
1348  // ... get the track
1349  AliAODTrack * track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrack));
1350  if(!track) {AliFatal("Not a standard AOD");}
1351  if(track->GetID()<0)continue;
1352  if((track->GetFlags())&(AliESDtrack::kTPCout)) nTPCout++;
1353  else continue;
1354  }
1355  if(fhMultVZEROTPCoutTrackCorrNoCut) fhMultVZEROTPCoutTrackCorrNoCut->Fill(nTPCout,mTotV0);
1356  Float_t mV0Cut=-2200.+(2.5*nTPCout)+(0.000012*nTPCout*nTPCout);
1357  if(fEnablePileupRejVZEROTPCout){
1358  if(mTotV0<mV0Cut) return;
1359  }
1360 
1361 
1362  if(fhMultVZEROTPCoutTrackCorr)fhMultVZEROTPCoutTrackCorr->Fill(nTPCout,mTotV0);
1363 
1364  //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
1365 
1366  //histogram filled with 1 for every AOD
1367  fNentries->Fill(0);
1368  fCounter->StoreEvent(aod,fCuts,fReadMC);
1369  //fCounter->StoreEvent(aod,fReadMC);
1370 
1371  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
1372  TString trigclass=aod->GetFiredTriggerClasses();
1373  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
1374 
1375  if(!fCuts->IsEventSelected(aod)) {
1376  if(fCuts->GetWhyRejection()==1) // rejected for pileup
1377  fNentries->Fill(13);
1378  if(fSys==1 && (fCuts->GetWhyRejection()==2 || fCuts->GetWhyRejection()==3)) fNentries->Fill(15);
1379  if(fCuts->GetWhyRejection()==7) fNentries->Fill(17);
1380  return;
1381  }
1382 
1383  // Check the Nb of SDD clusters
1384  if (fIsRejectSDDClusters) {
1385  Bool_t skipEvent = kFALSE;
1386  Int_t ntracks = 0;
1387  if (aod) ntracks = aod->GetNumberOfTracks();
1388  for(Int_t itrack=0; itrack<ntracks; itrack++) { // loop on tacks
1389  // ... get the track
1390  AliAODTrack * track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrack));
1391  if(!track) AliFatal("Not a standard AOD");
1392  if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
1393  skipEvent=kTRUE;
1394  fNentries->Fill(16);
1395  break;
1396  }
1397  }
1398  if (skipEvent) return;
1399  }
1400 
1401  // AOD primary vertex
1402  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
1403 
1404  Bool_t isGoodVtx=kFALSE;
1405 
1406  //vtx1->Print();
1407  TString primTitle = vtx1->GetTitle();
1408  if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
1409  isGoodVtx=kTRUE;
1410  fNentries->Fill(3);
1411  }
1412 
1413  // loop over candidates
1414  Int_t nInD0toKpi = inputArray->GetEntriesFast();
1415  if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
1416 
1417  // FILE *f=fopen("4display.txt","a");
1418  // printf("Number of D0->Kpi: %d\n",nInD0toKpi);
1419  Int_t nSelectedloose=0,nSelectedtight=0;
1420 
1421  // vHF object is needed to call the method that refills the missing info of the candidates
1422  // if they have been deleted in dAOD reconstruction phase
1423  // in order to reduce the size of the file
1424  AliAnalysisVertexingHF *vHF = new AliAnalysisVertexingHF();
1425 
1426  for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
1427  AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
1428 
1429  if(fUseSelectionBit && d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
1430  fNentries->Fill(2);
1431  continue; //skip the D0 from Dstar
1432  }
1433 
1434  if(!(vHF->FillRecoCand(aod,d))) {//Fill the data members of the candidate only if they are empty.
1435  fNentries->Fill(18); //monitor how often this fails
1436  continue;
1437  }
1438 
1439  if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
1440  nSelectedloose++;
1441  nSelectedtight++;
1442  if(fSys==0){
1443  if(fCuts->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
1444  }
1445  Int_t ptbin=fCuts->PtBin(d->Pt());
1446  if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
1447  fIsSelectedCandidate=fCuts->IsSelected(d,AliRDHFCuts::kAll,aod); //selected
1448  if(fFillVarHists) {
1449  //if(!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate)) {
1450  fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(0),0);
1451  fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(1),1);
1452  //check daughters
1453  if(!fDaughterTracks.UncheckedAt(0) || !fDaughterTracks.UncheckedAt(1)) {
1454  AliDebug(1,"at least one daughter not found!");
1455  fNentries->Fill(5);
1456  fDaughterTracks.Clear();
1457  continue;
1458  }
1459 
1460  // FillVarHists(aod,d,mcArray,fCuts,fDistr);
1461 
1462 
1463  if(!fFillAfterCuts) FillVarHists(aod,d,mcArray,fCuts,fDistr);
1464  if(fFillAfterCuts && fIsSelectedCandidate>=1) FillVarHists(aod,d,mcArray,fCuts,fDistr);
1465 
1466 
1467 
1468  }
1469 
1470  if (fDrawDetSignal) {
1471  DrawDetSignal(d, fDetSignal);
1472  }
1473 
1474  //filling invariant mass histograms
1475  FillMassHists(d,mcArray,mcHeader,fCuts,fOutputMass);
1476 
1477 
1478  if (fPIDCheck) {
1479  Int_t isSelectedPIDfill = 3;
1480  if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(d); //0 rejected,1 D0,2 D0bar, 3 both
1481 
1482  if (isSelectedPIDfill == 0)fNentries->Fill(7);
1483  if (isSelectedPIDfill == 1)fNentries->Fill(8);
1484  if (isSelectedPIDfill == 2)fNentries->Fill(9);
1485  if (isSelectedPIDfill == 3)fNentries->Fill(10);
1486  }
1487 
1488  }
1489 
1490  fDaughterTracks.Clear();
1491  //if(unsetvtx) d->UnsetOwnPrimaryVtx();
1492  } //end for prongs
1493 
1494  FillMCGenHists(aod, mcArray,fCuts,fDistr);
1495 
1496  fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1497  fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1498  delete vHF;
1499  // Post the data
1500  PostData(1,fOutputMass);
1501  PostData(2,fDistr);
1502  PostData(3,fNentries);
1503  PostData(5,fCounter);
1504  PostData(6,fOutputMassPt);
1505  PostData(7,fVariablesTree);
1506  PostData(8, fDetSignal);
1507 
1508  return;
1509 }
1510 
1511 //____________________________________________________________________________
1512 void AliAnalysisTaskSED0MassNonPromptFraction::FillMCGenHists(AliAODEvent* aod, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
1513 
1514 
1515  // cout << "numero di entries MC " << arrMC->GetEntries() << endl;
1516 
1517 
1518 
1519  for(Int_t i=0; i< arrMC->GetEntries(); i++){
1520 
1521  AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(i);
1522 
1523  Int_t ptbin=cuts->PtBin(partD0->Pt());
1524  Double_t pt = partD0->Pt();
1525 
1526  TString fillthisGEN ="";
1527  Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1528 
1529  if (TMath::Abs(((AliAODMCParticle*)partD0)->GetPdgCode() != 421)) continue;
1530 
1531  AliAODMCParticle *partDaughter = (AliAODMCParticle*)arrMC->At(partD0->GetDaughter(0));
1532  //if(!partDaughter) continue;
1533 
1534  if (!arrMC->At(partD0->GetDaughter(0))) continue;
1535 
1536  Double_t XvtxGen=0;
1537  Double_t YvtxGen=0;
1538  Double_t ZvtxGen=0;
1539 
1540 
1541  Double_t XpGen=0;
1542  Double_t YpGen=0;
1543  Double_t ZpGen=0;
1544 
1545  XvtxGen = partDaughter->Xv();
1546  YvtxGen = partDaughter->Yv();
1547  ZvtxGen = partDaughter->Zv();
1548 
1549 
1550 
1551  Double_t PvecGen[2]={-99,-99};
1552 
1553  PvecGen[0] = partD0->Px();
1554  PvecGen[1] = partD0->Py();
1555 
1556 
1557  Double_t ptravGen = partD0->Pt();
1558 
1559  AliAODMCHeader* mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1560  if(!mcHeader) {
1561  printf("AliAnalysisTaskSED0MassNonPromptFraction::UserExec: MC header branch not found!\n");
1562  return;
1563  }
1564 
1565  if (mcHeader){
1566  XpGen = mcHeader->GetVtxX();
1567  YpGen = mcHeader->GetVtxY();
1568  ZpGen = mcHeader->GetVtxZ();
1569  }
1570 
1571 
1572 
1573 
1574  TVector3 LGen ;
1575  LGen.SetXYZ (XvtxGen-XpGen, YvtxGen-YpGen, ZvtxGen-ZpGen);
1576 
1577 
1578  TVector3 PtvecGen;
1579 
1580  PtvecGen.SetXYZ (partD0->Px(), partD0->Py(), 0.);
1581 
1582 
1583 
1584  Double_t LxyGen;
1585  LxyGen = (LGen.Dot(PtvecGen))/ptravGen;
1586 
1587 
1588 
1589  Double_t xgenerate;
1590  xgenerate= LxyGen*mPDG/ptravGen;
1591 
1592 
1593 
1594  fillthisGEN= "hpseudoProperGENERATE_";
1595  fillthisGEN+=ptbin;
1596  ((TH1F*)listout->FindObject(fillthisGEN))->Fill(xgenerate);
1597 
1598 
1599  //old function if(CheckOrigin(arrMC,partD0)==5){
1600  if(AliVertexingHFUtils::CheckOrigin(arrMC,partD0,fUseQuarkTagInKine)==5){ //isPrimary=kFALSE;
1601  //isPrimary=kFALSE;
1602 
1603 
1604  fillthisGEN= "hpseudoProperSecGenerate_";
1605  fillthisGEN+=ptbin;
1606  ((TH1F*)listout->FindObject(fillthisGEN))->Fill(xgenerate);
1607 
1608 
1609 
1610  }
1611  else{
1612 
1613  fillthisGEN= "hpseudoProperPrimGenerate_";
1614  fillthisGEN+=ptbin;
1615  ((TH1F*)listout->FindObject(fillthisGEN))->Fill(xgenerate);
1616  }
1617 
1618 
1619 
1620 
1621  }//end ciclo for
1622 
1623 
1624 }//end void
1625 
1626 
1627 void AliAnalysisTaskSED0MassNonPromptFraction::DrawDetSignal(AliAODRecoDecayHF2Prong *part, TList *ListDetSignal)
1628 {
1629  //
1631  //
1632  fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(0), 0);
1633  fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(1), 1);
1634 
1635  AliESDtrack *esdtrack1 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1636  AliESDtrack *esdtrack2 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
1637 
1638 
1639  //For filling detector histograms
1640  Int_t isSelectedPIDfill = 3;
1641  if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
1642 
1643 
1644  //fill "before PID" histos with every daughter
1645  ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
1646  ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
1647  ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
1648  ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
1649 
1650  if (isSelectedPIDfill != 0) { //fill "After PID" with everything that's not rejected
1651  ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
1652  ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
1653  ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
1654  ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
1655 
1656  }
1657 
1658  delete esdtrack1;
1659  delete esdtrack2;
1660 
1661  return;
1662 }
1663 
1664 //____________________________________________________________________________
1665 void AliAnalysisTaskSED0MassNonPromptFraction::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
1666  //
1668  //
1669 
1670 
1671  Int_t pdgDgD0toKpi[2]={321,211};
1672  Int_t lab=-9999;
1673  if(fReadMC) lab=part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
1674  //Double_t pt = d->Pt(); //mother pt
1675  Int_t isSelectedPID=3;
1676  if(!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPID=cuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
1677  if (!fPIDCheck) { //if fPIDCheck, this is already filled elsewhere
1678  if (isSelectedPID==0)fNentries->Fill(7);
1679  if (isSelectedPID==1)fNentries->Fill(8);
1680  if (isSelectedPID==2)fNentries->Fill(9);
1681  if (isSelectedPID==3)fNentries->Fill(10);
1682  //fNentries->Fill(8+isSelectedPID);
1683  }
1684 
1685  if(fCutOnDistr && !fIsSelectedCandidate) return;
1686  //printf("\nif no cuts or cuts passed\n");
1687 
1688 
1689  //add distr here
1690  UInt_t pdgs[2];
1691 
1692  Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1693 
1694 
1695  Double_t Xvtx=0;
1696  Double_t Yvtx=0;
1697  Double_t Zvtx=0;
1698 
1699 
1700  Xvtx = part->GetSecondaryVtx()->GetX();
1701 
1702  Yvtx = part->GetSecondaryVtx()->GetY();
1703  Zvtx = part->GetSecondaryVtx()->GetZ();
1704 
1705  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
1706 
1707  Double_t Xp =0;
1708  Double_t Yp =0;
1709  Double_t Zp =0;
1710 
1711  Xp=aod->GetPrimaryVertex()->GetX();
1712  Yp=aod->GetPrimaryVertex()->GetY();
1713  Zp=aod->GetPrimaryVertex()->GetZ();
1714 
1715 
1716 
1717  Double_t CosinePointingAngle = part->CosPointingAngle();
1718 
1719  Double_t dSecToPrim=-99;
1720  dSecToPrim = TMath::Sqrt((Xvtx-Xp)*(Xvtx-Xp) +(Yvtx-Yp)*(Yvtx-Yp)+(Zvtx-Zp)*(Zvtx-Zp));
1721 
1722 
1723 
1724  Double_t Lvec[3]={-99,-99,-99};
1725  Double_t Pvec[2]={-99,-99};
1726 
1727  Pvec[0] = part->Px();
1728  Pvec[1] = part->Py();
1729 
1730  Double_t Pt = TMath::Sqrt(Pvec[0]*Pvec[0]+Pvec[1]*Pvec[1]);
1731 
1732 
1733  Double_t deltaX = Xp-Xvtx;
1734 
1735  Lvec[0]= Xp-Xvtx;
1736  Lvec[1]= Yp-Yvtx;
1737  Lvec[2]= Zp-Zvtx;
1738 
1739  TVector3 posMom;
1740 
1741  //posMom.SetXYZ(pTrack->Px(),pTrack->Py(),pTrack->Pz());
1742 
1743  TVector3 L ;
1744  L.SetXYZ (Xvtx-Xp, Yvtx-Yp, Zvtx-Zp);
1745 
1746  TVector3 L1 ;
1747  L1.SetXYZ (Xp-Xvtx, Yp-Yvtx, 0);
1748 
1749  TVector3 Ptvec;
1750 
1751  Ptvec.SetXYZ (part->Px(), part->Py(), 0);
1752 
1753 
1754 
1755  Double_t ptrav = part->Pt();
1756 
1757 
1758  Double_t Lxy;
1759  Lxy = (L.Dot(Ptvec))/ptrav;
1760 
1761  Double_t Lxy1;
1762  Lxy1 = (L1.Dot(Ptvec))/ptrav;
1763 
1764  Double_t x;
1765  x= Lxy*mPDG/ptrav;
1766 
1767 
1768 
1769 
1770 
1771 
1772  // Fill candidate variable Tree (track selection, no candidate selection)
1773  if( fWriteVariableTree && !part->HasBadDaughters() && fCuts->AreDaughtersSelected(part) && fCuts->IsSelectedPID(part) ){
1774 
1775 
1776 
1777  fCandidateVariables[0] = part->InvMassD0(); //ok
1778  fCandidateVariables[1] = part->InvMassD0bar(); //ok
1779  fCandidateVariables[2] = part->Pt(); //ok
1780 
1781  // fCandidateVariables[3] = part->Pt2Prong(0);
1782  // fCandidateVariables[4] = part->Pt2Prong(1);
1783  // fCandidateVariables[5] = part->Getd0Prong(0);
1784  // fCandidateVariables[6] = part->Getd0Prong(1);
1785  // fCandidateVariables[7] = part->Prodd0d0();
1786  // fCandidateVariables[8] = part->DecayLengthXY();
1787 
1788  fCandidateVariables[3] = part->CtD0();
1789  fCandidateVariables[4] = x;
1790  fCandidateVariables[5] = fIsSelectedCandidate;
1791  // fCandidateVariables[10]= part->ImpParXY()*10000.;
1792  // fCandidateVariables[14] = fCuts->IsSelectedSpecialCuts(part);
1793 
1794 
1795 
1796  fVariablesTree->Fill();
1797  }
1798 
1799 
1800 
1801  if (!fIsSelectedCandidate){
1802  //cout<<" cut " << cuts << " Rejected because "<<cuts->GetWhy()<<endl;
1803  return;
1804  }
1805 
1806  if(fDebug>2) cout << "Candidate selected" << endl;
1807 
1808 
1809  Double_t CtauD0 = -1;
1810 
1811  Double_t C= 2.99792458e10;
1812 
1813  //CtauD0 = (((part->DecayLength())*mPDG*C)/(TMath::Sqrt(part->P()))); //confronta con il valore CtD0
1814  CtauD0 = part->CtD0();
1815 
1816  pdgs[0]=211;
1817  pdgs[1]=321;
1818  Double_t minvD0 = part->InvMassD0();
1819  pdgs[1]=211;
1820  pdgs[0]=321;
1821  Double_t minvD0bar = part->InvMassD0bar();
1822  //cout<<"inside mass cut"<<endl;
1823  Double_t decayLengthxy = part->DecayLengthXY();
1824  Double_t normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1825 
1826 
1827  Double_t invmasscut=0.03;
1828 
1829  TString fillthispi="",fillthisK="",fillthis="", fillthispt="", fillthisetaphi="";
1830 
1831 
1832 
1833  Int_t ptbin=cuts->PtBin(part->Pt());
1834  Double_t pt = part->Pt();
1835 
1836  Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
1837  dz1[0]=-99; dz2[0]=-99;
1838  Double_t d0[2];
1839  Double_t decl[2]={-99,-99};
1840  Bool_t recalcvtx=kFALSE;
1841 
1842  fillthis= "hCorrCosX_";
1843  fillthis+=ptbin;
1844  ((TH2F*)listout->FindObject(fillthis))->Fill(CosinePointingAngle, x);
1845 
1846 
1847 
1848  fillthis= "hCorrDistX_";
1849  fillthis+=ptbin;
1850  ((TH2F*)listout->FindObject(fillthis))->Fill(dSecToPrim, x);
1851 
1852  Double_t impparXY=part->ImpParXY()*10000.;
1853 
1854  fillthis= "hCorrXd0D0_";
1855  fillthis+=ptbin;
1856  ((TH2F*)listout->FindObject(fillthis))->Fill(x,impparXY);
1857 
1858 
1859  fillthis="hctauD0_";
1860  fillthis+=ptbin;
1861  ((TH1F*)listout->FindObject(fillthis))->Fill(CtauD0);
1862 
1863  if(fCuts->GetIsPrimaryWithoutDaughters()){
1864  recalcvtx=kTRUE;
1865  //recalculate vertex
1866  AliAODVertex *vtxProp=0x0;
1867  vtxProp=GetPrimaryVtxSkipped(aod);
1868  if(vtxProp) {
1869  part->SetOwnPrimaryVtx(vtxProp);
1870  //Bool_t unsetvtx=kTRUE;
1871  //Calculate d0 for daughter tracks with Vtx Skipped
1872  AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1873  AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
1874  esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
1875  esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
1876  delete vtxProp; vtxProp=NULL;
1877  delete esdtrack1;
1878  delete esdtrack2;
1879  }
1880 
1881  d0[0]=dz1[0];
1882  d0[1]=dz2[0];
1883 
1884 
1885  decl[0]=part->DecayLength2();
1886  decl[1]=part->NormalizedDecayLength2();
1887  part->UnsetOwnPrimaryVtx();
1888 
1889  }
1890 
1891  Double_t cosThetaStarD0 = 99;
1892  Double_t cosThetaStarD0bar = 99;
1893  Double_t cosPointingAngle = 99;
1894  Double_t ctau = -1;
1895  // Double_t CtauD0 = -1;
1896  Double_t CtauD0z = -1;
1897  Double_t normalizedDecayLength2 = -1;
1898  Double_t decayLength2 = -1;
1899  Double_t ptProng[2]={-99,-99};
1900  Double_t pzprong0 = -99;
1901  Double_t pzprong1 = -99;
1902  Double_t d0Prong[2]={-99,-99};
1903  Double_t etaD = 99.;
1904  Double_t phiD = 99.;
1905  Double_t decayLengthz = -1;
1906 
1907  //disable the PID
1908  if(!fUsePid4Distr) isSelectedPID=0;
1909 
1910  Int_t labD0=-1;
1911  if (fReadMC) labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi);
1912  AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
1913  //Int_t pdgD0 = partD0->GetPdgCode();
1914  if(fReadMC){
1915 
1916 
1917  if(labD0>=0) {
1918 
1919  fillthis= "hpseudoProperDLRec_";
1920  fillthis+=ptbin;
1921  ((TH1F*)listout->FindObject(fillthis))->Fill(x);
1922 
1923 
1924  fillthis= "hCorrCosXtrue_";
1925  fillthis+=ptbin;
1926  ((TH2F*)listout->FindObject(fillthis))->Fill(CosinePointingAngle, x);
1927 
1928  fillthis= "hCorrDistXtrue_";
1929  fillthis+=ptbin;
1930  ((TH2F*)listout->FindObject(fillthis))->Fill(dSecToPrim, x);
1931 
1932 
1933  fillthis= "hCorrXd0D0true_";
1934  fillthis+=ptbin;
1935  ((TH2F*)listout->FindObject(fillthis))->Fill(x,impparXY);
1936 
1937 
1938  fillthis= "hImpParXYtrue_";
1939  fillthis+=ptbin;
1940  ((TH1F*)listout->FindObject(fillthis))->Fill(impparXY);
1941 
1942 
1943  fillthis= "hCorrNormLxyXtrue_";
1944  fillthis+=ptbin;
1945  ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy, x);
1946 
1947 
1948  fillthis="hctauD0S_";
1949  fillthis+=ptbin;
1950  ((TH1F*)listout->FindObject(fillthis))->Fill(CtauD0);
1951 
1952 
1953  Double_t XvtxGen=0;
1954  Double_t YvtxGen=0;
1955  Double_t ZvtxGen=0;
1956 
1957 
1958  Double_t XpGen=0;
1959  Double_t YpGen=0;
1960  Double_t ZpGen=0;
1961 
1962 
1963  AliAODMCParticle *partDaughter = (AliAODMCParticle*)arrMC->At(partD0->GetDaughter(0));
1964 
1965  if (!arrMC->At(partD0->GetDaughter(0))) return;
1966 
1967 
1968  XvtxGen = partDaughter->Xv();
1969  YvtxGen = partDaughter->Yv();
1970  ZvtxGen = partDaughter->Zv();
1971 
1972 
1973  Double_t PvecGen[2]={-99,-99};
1974 
1975  PvecGen[0] = partD0->Px();
1976  PvecGen[1] = partD0->Py();
1977 
1978 
1979  Double_t ptravGen = partD0->Pt();
1980 
1981  AliAODMCHeader* mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1982  if(!mcHeader) {
1983  printf("alysisTaskSED0MassNonPromptFraction::UserExec: MC header branch not found!\n");
1984  return;
1985  }
1986 
1987  if (mcHeader){
1988  XpGen = mcHeader->GetVtxX();
1989  YpGen = mcHeader->GetVtxY();
1990  ZpGen = mcHeader->GetVtxZ();}
1991 
1992 
1993 
1994  TVector3 LGen ;
1995  LGen.SetXYZ (XvtxGen-XpGen, YvtxGen-YpGen, ZvtxGen-ZpGen);
1996 
1997 
1998  TVector3 PtvecGen;
1999 
2000  PtvecGen.SetXYZ (partD0->Px(), partD0->Py(), 0);
2001 
2002 
2003  Double_t LxyGen;
2004  LxyGen = (LGen.Dot(PtvecGen))/ptravGen;
2005 
2006 
2007  Double_t xgen;
2008  xgen= LxyGen*mPDG/ptravGen;
2009 
2010 
2011  Double_t Residual;
2012  Residual= xgen -x;
2013 
2014  fillthis= "hGenMatchRec_";
2015  fillthis+=ptbin;
2016  ((TH1F*)listout->FindObject(fillthis))->Fill(xgen);
2017 
2018 
2019  fillthis= "hresidualX_";
2020  fillthis+=ptbin;
2021  ((TH1F*)listout->FindObject(fillthis))->Fill(Residual);
2022 
2023 
2024 
2025  if(AliVertexingHFUtils::CheckOrigin(arrMC,partD0,fUseQuarkTagInKine)==5){ //isPrimary=kFALSE;
2026 
2027  //old function if(CheckOrigin(arrMC,partD0)==5){
2028  //isPrimary=kFALSE; //la funzione checkorigin restituisce 5 se la D deriva da B, ovvero se è secondaria!
2029 
2030 
2031 
2032  fillthis= "hpseudoProperDLRecSec_";
2033  fillthis+=ptbin;
2034  ((TH1F*)listout->FindObject(fillthis))->Fill(x);
2035 
2036 
2037  fillthis= "hCorrDistX_sec_";
2038  fillthis+=ptbin;
2039  ((TH2F*)listout->FindObject(fillthis))->Fill(dSecToPrim, x);
2040 
2041  fillthis= "hCorrNormLxyX_sec_";
2042  fillthis+=ptbin;
2043  ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy, x);
2044 
2045  fillthis= "hCorrCosX_sec_";
2046  fillthis+=ptbin;
2047  ((TH2F*)listout->FindObject(fillthis))->Fill(CosinePointingAngle, x);
2048 
2049  fillthis= "hCorrXd0D0_sec_";
2050  fillthis+=ptbin;
2051  ((TH2F*)listout->FindObject(fillthis))->Fill(x,impparXY);
2052 
2053  fillthis= "hImpParXY_sec_";
2054  fillthis+=ptbin;
2055  ((TH1F*)listout->FindObject(fillthis))->Fill(impparXY);
2056 
2057 
2058 
2059  fillthis= "hresidualXSec_";
2060  fillthis+=ptbin;
2061  ((TH1F*)listout->FindObject(fillthis))->Fill(Residual);
2062 
2063  fillthis= "hGenMatchRecSec_";
2064  fillthis+=ptbin;
2065  ((TH1F*)listout->FindObject(fillthis))->Fill(xgen);
2066 
2067 
2068 
2069 
2070  fillthis="hctauD0sec_";
2071  fillthis+=ptbin;
2072  ((TH1F*)listout->FindObject(fillthis))->Fill(CtauD0);
2073 
2074 
2075 
2076 
2077 
2078  }
2079  else{
2080 
2081  fillthis= "hpseudoProperDLRecPrim_";
2082  fillthis+=ptbin;
2083  ((TH1F*)listout->FindObject(fillthis))->Fill(x);
2084 
2085 
2086  fillthis= "hresidual_Primary_";
2087  fillthis+=ptbin;
2088  ((TH1F*)listout->FindObject(fillthis))->Fill(Residual);
2089 
2090  fillthis= "hGenMatchRecPrim_";
2091  fillthis+=ptbin;
2092  ((TH1F*)listout->FindObject(fillthis))->Fill(xgen);
2093 
2094  fillthis= "hCorrNormLxyX_prim_";
2095  fillthis+=ptbin;
2096  ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy, x);
2097 
2098  fillthis= "hCorrDistX_prim_";
2099  fillthis+=ptbin;
2100  ((TH2F*)listout->FindObject(fillthis))->Fill(dSecToPrim, x);
2101 
2102  fillthis= "hCorrCosX_prim_";
2103  fillthis+=ptbin;
2104  ((TH2F*)listout->FindObject(fillthis))->Fill(CosinePointingAngle, x);
2105 
2106  fillthis= "hCorrXd0D0_prim_";
2107  fillthis+=ptbin;
2108  ((TH2F*)listout->FindObject(fillthis))->Fill(x,impparXY);
2109 
2110  fillthis= "hImpParXY_prim_";
2111  fillthis+=ptbin;
2112  ((TH1F*)listout->FindObject(fillthis))->Fill(impparXY);
2113 
2114 
2115  fillthis="hctauD0prim_";
2116  fillthis+=ptbin;
2117  ((TH1F*)listout->FindObject(fillthis))->Fill(CtauD0);
2118 
2119 
2120  }
2121 
2122  }
2123 
2124  }
2125 
2126 
2127 
2128  if((lab>=0 && fReadMC) || (!fReadMC && isSelectedPID)){ //signal (from MC or PID)
2129 
2130 
2131 
2132  //check pdg of the prongs
2133  AliAODTrack *prong0=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
2134  AliAODTrack *prong1=(AliAODTrack*)fDaughterTracks.UncheckedAt(1);
2135 
2136  if(!prong0 || !prong1) {
2137  return;
2138  }
2139  Int_t labprong[2];
2140  if(fReadMC){
2141  labprong[0]=prong0->GetLabel();
2142  labprong[1]=prong1->GetLabel();
2143  }
2144  AliAODMCParticle *mcprong=0;
2145  Int_t pdgProng[2]={0,0};
2146 
2147  for (Int_t iprong=0;iprong<2;iprong++){
2148  if(fReadMC && labprong[iprong]>=0) {
2149  mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
2150  pdgProng[iprong]=mcprong->GetPdgCode();
2151  }
2152  }
2153 
2154 
2155 
2156 
2157  if(fSys==0){
2158  //no mass cut ditributions: ptbis
2159 
2160  fillthispi="hptpiSnoMcut_";
2161  fillthispi+=ptbin;
2162 
2163  fillthisK="hptKSnoMcut_";
2164  fillthisK+=ptbin;
2165 
2166  if ((TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321)
2167  || (isSelectedPID==1 || isSelectedPID==3)){
2168  ((TH1F*)listout->FindObject(fillthispi))->Fill(prong0->Pt());
2169  ((TH1F*)listout->FindObject(fillthisK))->Fill(prong1->Pt());
2170  }
2171 
2172  if ((TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211)
2173  || (isSelectedPID==2 || isSelectedPID==3)){
2174  ((TH1F*)listout->FindObject(fillthisK))->Fill(prong0->Pt());
2175  ((TH1F*)listout->FindObject(fillthispi))->Fill(prong1->Pt());
2176  }
2177  }
2178 
2179  //no mass cut ditributions: mass
2180 
2181  etaD = part->Eta();
2182  phiD = part->Phi();
2183 
2184 
2185  fillthis="hMassS_";
2186  fillthis+=ptbin;
2187  fillthispt="histSgnPt";
2188 
2189  if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)
2190  || (!fReadMC && (isSelectedPID==1 || isSelectedPID==3))){//D0
2191  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
2192  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
2193 
2194  fillthisetaphi="hetaphiD0candidateS_";
2195  fillthisetaphi+=ptbin;
2196  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
2197 
2198  if(TMath::Abs(minvD0-mPDG)<0.05){
2199  fillthisetaphi="hetaphiD0candidatesignalregionS_";
2200  fillthisetaphi+=ptbin;
2201  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
2202  }
2203 
2204  }
2205  else { //D0bar
2206  if(fReadMC || (!fReadMC && isSelectedPID > 1)){
2207  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
2208  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
2209 
2210  fillthisetaphi="hetaphiD0barcandidateS_";
2211  fillthisetaphi+=ptbin;
2212  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
2213 
2214  if(TMath::Abs(minvD0bar-mPDG)<0.05){
2215  fillthisetaphi="hetaphiD0barcandidatesignalregionS_";
2216  fillthisetaphi+=ptbin;
2217  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
2218  }
2219 
2220  }
2221  }
2222 
2223  //apply cut on invariant mass on the pair
2224 
2225  //invmasscut
2226 
2227  if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
2228 
2229  cosThetaStarD0 = part->CosThetaStarD0();
2230  cosThetaStarD0bar = part->CosThetaStarD0bar();
2231  cosPointingAngle = part->CosPointingAngle();
2232  normalizedDecayLength2 = part->NormalizedDecayLength2();
2233  ctau=part->CtD0();
2234 
2235  decayLength2 = part->DecayLength2();
2236 
2237 
2238  decayLengthxy = part->DecayLengthXY();
2239 
2240  normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
2241 
2242  ptProng[0]=prong0->Pt(); ptProng[1]=prong1->Pt();
2243  d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
2244  decayLengthz = TMath::Sqrt(decayLength2-decayLengthxy*decayLengthxy);
2245  pzprong0 = prong0->Pz();
2246  pzprong1 = prong1->Pz();
2247 
2248 
2249 
2250 
2251  if(fArray==1) cout<<"LS signal: ERROR"<<endl;
2252  for (Int_t iprong=0; iprong<2; iprong++){
2253  AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
2254  if (fReadMC) labprong[iprong]=prong->GetLabel();
2255 
2256  //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
2257  Int_t pdgprong=0;
2258  if(fReadMC && labprong[iprong]>=0) {
2259  mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
2260  pdgprong=mcprong->GetPdgCode();
2261  }
2262 
2263  Bool_t isPionHere[2]={(isSelectedPID==1 || isSelectedPID==3) ? kTRUE : kFALSE,(isSelectedPID==2 || isSelectedPID==3) ? kTRUE : kFALSE};
2264 
2265  if(TMath::Abs(pdgprong)==211 || isPionHere[iprong]) {
2266  //cout<<"pi"<<endl;
2267 
2268  if(fSys==0){
2269  fillthispi="hptpiS_";
2270  fillthispi+=ptbin;
2271  ((TH1F*)listout->FindObject(fillthispi))->Fill(ptProng[iprong]);
2272 
2273 
2274 
2275  }
2276 
2277  fillthispi="hd0piS_";
2278  fillthispi+=ptbin;
2279  ((TH1F*)listout->FindObject(fillthispi))->Fill(d0Prong[iprong]);
2280  if(recalcvtx) {
2281 
2282  fillthispi="hd0vpiS_";
2283  fillthispi+=ptbin;
2284  ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
2285  }
2286 
2287  }
2288 
2289  if(TMath::Abs(pdgprong)==321 || !isPionHere[iprong]) {
2290  //cout<<"kappa"<<endl;
2291  if(fSys==0){
2292  fillthisK="hptKS_";
2293  fillthisK+=ptbin;
2294  ((TH1F*)listout->FindObject(fillthisK))->Fill(ptProng[iprong]);
2295 
2296 
2297  }
2298 
2299  fillthisK="hd0KS_";
2300  fillthisK+=ptbin;
2301  ((TH1F*)listout->FindObject(fillthisK))->Fill(d0Prong[iprong]);
2302  if (recalcvtx){
2303  fillthisK="hd0vKS_";
2304  fillthisK+=ptbin;
2305  ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
2306  }
2307  }
2308 
2309  if(fSys==0){
2310  fillthis="hcosthpointd0S_";
2311  fillthis+=ptbin;
2312  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[iprong]);
2313  }
2314  } //end loop on prongs
2315 
2316  fillthis="hdcaS_";
2317  fillthis+=ptbin;
2318  ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
2319 
2320  fillthis="hcosthetapointS_";
2321  fillthis+=ptbin;
2322  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
2323 
2324  fillthis="hcosthetapointxyS_";
2325  fillthis+=ptbin;
2326  ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
2327 
2328 
2329  fillthis="hd0d0S_";
2330  fillthis+=ptbin;
2331  ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
2332 
2333  fillthis="hdeclS_";
2334  fillthis+=ptbin;
2335  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
2336 
2337  fillthis="hnormdeclS_";
2338  fillthis+=ptbin;
2339  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
2340 
2341  //ctau DzeroZ
2342  fillthis="hctauS_";
2343  fillthis+=ptbin;
2344  ((TH1F*)listout->FindObject(fillthis))->Fill(CtauD0z);
2345 
2346 
2347 
2348 
2349 
2350 
2351  //correlation ctauXY ctauZ
2352  fillthis="hctauXYctauZS_" ;
2353  fillthis+=ptbin;
2354  ((TH2F*)listout->FindObject(fillthis))->Fill(CtauD0,CtauD0z);
2355 
2356 
2357 
2358 
2359 
2360 
2361 
2362  fillthis="hdeclxyS_";
2363  fillthis+=ptbin;
2364  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
2365 
2366  fillthis="hnormdeclxyS_";
2367  fillthis+=ptbin;
2368  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
2369 
2370  fillthis="hdeclxyd0d0S_";
2371  fillthis+=ptbin;
2372  ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
2373 
2374 
2375 
2376  fillthis="hnormdeclxyd0d0S_";
2377  fillthis+=ptbin;
2378  ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
2379 
2380  if(recalcvtx) {
2381  fillthis="hdeclvS_";
2382  fillthis+=ptbin;
2383  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
2384 
2385  fillthis="hnormdeclvS_";
2386  fillthis+=ptbin;
2387  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
2388 
2389  fillthis="hd0d0vS_";
2390  fillthis+=ptbin;
2391  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
2392  }
2393 
2394  if(fSys==0){
2395  fillthis="hcosthetastarS_";
2396  fillthis+=ptbin;
2397  if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)) ((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
2398  else {
2399  if (fReadMC || isSelectedPID>1)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
2400  if(isSelectedPID==1 || isSelectedPID==3)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
2401  }
2402 
2403  fillthis="hcosthpointd0d0S_";
2404  fillthis+=ptbin;
2405  ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
2406  }
2407 
2408  if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)){
2409  for(Int_t it=0; it<2; it++){
2410  fillthis="hptD0S_";
2411  fillthis+=ptbin;
2412  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
2413  fillthis="hphiD0S_";
2414  fillthis+=ptbin;
2415  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
2416  Int_t nPointsITS = 0;
2417  for (Int_t il=0; il<6; il++){
2418  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
2419  }
2420  fillthis="hNITSpointsD0vsptS_";
2421  fillthis+=ptbin;
2422  ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(),nPointsITS);
2423  fillthis="hNSPDpointsD0S_";
2424  fillthis+=ptbin;
2425  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
2426  ((TH1I*)listout->FindObject(fillthis))->Fill(0);
2427  }
2428  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
2429  ((TH1I*)listout->FindObject(fillthis))->Fill(1);
2430  }
2431  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
2432  ((TH1I*)listout->FindObject(fillthis))->Fill(2);
2433  }
2434  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
2435  ((TH1I*)listout->FindObject(fillthis))->Fill(3);
2436  }
2437  fillthis="hNclsD0vsptS_";
2438  fillthis+=ptbin;
2439  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
2440  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->GetTPCNcls();
2441  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
2442  }
2443  }
2444  else {
2445  if (fReadMC || isSelectedPID>1){
2446  for(Int_t it=0; it<2; it++){
2447  fillthis="hptD0barS_";
2448  fillthis+=ptbin;
2449  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
2450  fillthis="hphiD0barS_";
2451  fillthis+=ptbin;
2452  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
2453  fillthis="hNclsD0barvsptS_";
2454  fillthis+=ptbin;
2455  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
2456  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
2457  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
2458  }
2459  }
2460  if(isSelectedPID==1 || isSelectedPID==3){
2461  for(Int_t it=0; it<2; it++){
2462  fillthis="hptD0S_";
2463  fillthis+=ptbin;
2464  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
2465  fillthis="hphiD0S_";
2466  fillthis+=ptbin;
2467  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
2468  Int_t nPointsITS = 0;
2469  for (Int_t il=0; il<6; il++){
2470  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
2471  }
2472  fillthis="hNITSpointsD0vsptS_";
2473  fillthis+=ptbin;
2474  ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(), nPointsITS);
2475  fillthis="hNSPDpointsD0S_";
2476  fillthis+=ptbin;
2477  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
2478  ((TH1I*)listout->FindObject(fillthis))->Fill(0);
2479  }
2480  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
2481  ((TH1I*)listout->FindObject(fillthis))->Fill(1);
2482  }
2483  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
2484  ((TH1I*)listout->FindObject(fillthis))->Fill(2);
2485  }
2486  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
2487  ((TH1I*)listout->FindObject(fillthis))->Fill(3);
2488  }
2489  fillthis="hNclsD0vsptS_";
2490  fillthis+=ptbin;
2491  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
2492  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->GetTPCNcls();
2493  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
2494  }
2495  }
2496  }
2497 
2498 
2499  } //end mass cut
2500 
2501  } else{ //Background or LS
2502  //if(!fReadMC){
2503  //cout<<"is background"<<endl;
2504 
2505  etaD = part->Eta();
2506  phiD = part->Phi();
2507 
2508  //no mass cut distributions: mass, ptbis
2509  fillthis="hMassB_";
2510  fillthis+=ptbin;
2511  fillthispt="histBkgPt";
2512 
2513  if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3))) {
2514  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
2515  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
2516 
2517  fillthisetaphi="hetaphiD0candidateB_";
2518  fillthisetaphi+=ptbin;
2519  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
2520 
2521  if(TMath::Abs(minvD0-mPDG)<0.05){
2522  fillthisetaphi="hetaphiD0candidatesignalregionB_";
2523  fillthisetaphi+=ptbin;
2524  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
2525  }
2526  }
2527  if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) {
2528  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
2529  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
2530 
2531  fillthisetaphi="hetaphiD0barcandidateB_";
2532  fillthisetaphi+=ptbin;
2533  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
2534 
2535  if(TMath::Abs(minvD0bar-mPDG)<0.05){
2536  fillthisetaphi="hetaphiD0barcandidatesignalregionB_";
2537  fillthisetaphi+=ptbin;
2538  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
2539  }
2540 
2541  }
2542  if(fSys==0){
2543  fillthis="hptB1prongnoMcut_";
2544  fillthis+=ptbin;
2545 
2546  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
2547 
2548  fillthis="hptB2prongsnoMcut_";
2549  fillthis+=ptbin;
2550  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
2551  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
2552  }
2553 
2554 
2555  //apply cut on invariant mass on the pair
2556  if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
2557  if(fSys==0){
2558  ptProng[0]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt(); ptProng[1]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt();
2559  cosThetaStarD0 = part->CosThetaStarD0();
2560  cosThetaStarD0bar = part->CosThetaStarD0bar();
2561  }
2562 // AliAODTrack *prong0=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
2563 
2564 // AliAODTrack *prong1=(AliAODTrack*)fDaughterTracks.UncheckedAt(1);
2565 
2566  cosPointingAngle = part->CosPointingAngle();
2567  normalizedDecayLength2 = part->NormalizedDecayLength2();
2568  ctau = part->CtD0();
2569  Double_t pzprong0 = -99;
2570  Double_t pzprong1 = -99;
2571 
2572  decayLengthz = TMath::Sqrt(part->DecayLength2() - (part->DecayLengthXY())*(part->DecayLength2() - (part->DecayLengthXY())));
2573 
2574 
2575 
2576 
2577  Double_t C= 2.99792458e10;
2578 
2579  CtauD0z = ((decayLengthz*mPDG*C)/(TMath::Sqrt(part->Pz()* part->Pz())));
2580 
2581 
2582  decayLength2 = part->DecayLength2();
2583  decayLengthxy = part->DecayLengthXY();
2584  normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
2585  d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
2586 
2587 
2588  AliAODTrack *prongg=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
2589  if(!prongg) {
2590  if(fDebug>2) cout<<"No daughter found";
2591  return;
2592  }
2593  else{
2594  if(fArray==1){
2595  if(prongg->Charge()==1) {
2596  //fTotPosPairs[ptbin]++;
2597  ((TH1F*)fOutputMass->FindObject("hpospair"))->Fill(ptbin);
2598  } else {
2599  //fTotNegPairs[ptbin]++;
2600  ((TH1F*)fOutputMass->FindObject("hnegpair"))->Fill(ptbin);
2601  }
2602  }
2603  }
2604 
2605  //fill pt and phi distrib for prongs with M cut
2606 
2607  if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3))){
2608  for(Int_t it=0; it<2; it++){
2609  fillthis="hptD0B_";
2610  fillthis+=ptbin;
2611  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
2612  fillthis="hphiD0B_";
2613  fillthis+=ptbin;
2614  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
2615 
2616  Int_t nPointsITS = 0;
2617  for (Int_t il=0; il<6; il++){
2618  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
2619  }
2620  fillthis="hNITSpointsD0vsptB_";
2621  fillthis+=ptbin;
2622  ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(), nPointsITS);
2623  fillthis="hNSPDpointsD0B_";
2624  fillthis+=ptbin;
2625  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
2626  ((TH1I*)listout->FindObject(fillthis))->Fill(0);
2627  }
2628  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
2629  ((TH1I*)listout->FindObject(fillthis))->Fill(1);
2630  }
2631  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
2632  ((TH1I*)listout->FindObject(fillthis))->Fill(2);
2633  }
2634  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
2635  ((TH1I*)listout->FindObject(fillthis))->Fill(3);
2636  }
2637  fillthis="hNclsD0vsptB_";
2638  fillthis+=ptbin;
2639  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
2640  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
2641  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
2642  }
2643 
2644 
2645  }
2646 
2647  if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) {
2648  for(Int_t it=0; it<2; it++){
2649  fillthis="hptD0barB_";
2650  fillthis+=ptbin;
2651  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
2652  fillthis="hphiD0barB_";
2653  fillthis+=ptbin;
2654  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
2655  fillthis="hNclsD0barvsptB_";
2656  fillthis+=ptbin;
2657  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
2658  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
2659  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
2660  }
2661  }
2662 
2663  fillthis="hd0B_";
2664  fillthis+=ptbin;
2665  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
2666  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
2667 
2668  if(fReadMC){
2669  Int_t pdgMother[2]={0,0};
2670  Double_t factor[2]={1,1};
2671 
2672  for(Int_t iprong=0;iprong<2;iprong++){
2673  AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
2674  lab=prong->GetLabel();
2675  if(lab>=0){
2676  AliAODMCParticle* mcprong=(AliAODMCParticle*)arrMC->At(lab);
2677  if(mcprong){
2678  Int_t labmom=mcprong->GetMother();
2679  if(labmom>=0){
2680  AliAODMCParticle* mcmother=(AliAODMCParticle*)arrMC->At(labmom);
2681  if(mcmother) pdgMother[iprong]=mcmother->GetPdgCode();
2682  }
2683  }
2684  }
2685 
2686  if(fSys==0){
2687 
2688  fillthis="hd0moresB_";
2689  fillthis+=ptbin;
2690 
2691  if(TMath::Abs(pdgMother[iprong])==310 || TMath::Abs(pdgMother[iprong])==130 || TMath::Abs(pdgMother[iprong])==321){ //K^0_S, K^0_L, K^+-
2692  if(ptProng[iprong]<=1)factor[iprong]=1./.7;
2693  else factor[iprong]=1./.6;
2694  fNentries->Fill(11);
2695  }
2696 
2697  if(TMath::Abs(pdgMother[iprong])==3122) { //Lambda
2698  factor[iprong]=1./0.25;
2699  fNentries->Fill(12);
2700  }
2701  fillthis="hd0moresB_";
2702  fillthis+=ptbin;
2703 
2704  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[iprong],factor[iprong]);
2705 
2706  if(recalcvtx){
2707  fillthis="hd0vmoresB_";
2708  fillthis+=ptbin;
2709  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[iprong],factor[iprong]);
2710  }
2711  }
2712  } //loop on prongs
2713 
2714  if(fSys==0){
2715  fillthis="hd0d0moresB_";
2716  fillthis+=ptbin;
2717  ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0(),factor[0]*factor[1]);
2718 
2719  fillthis="hcosthetapointmoresB_";
2720  fillthis+=ptbin;
2721  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,factor[0]*factor[1]);
2722 
2723  if(recalcvtx){
2724  fillthis="hd0d0vmoresB_";
2725  fillthis+=ptbin;
2726  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1],factor[0]*factor[1]);
2727  }
2728  }
2729 
2730 
2731 
2732  fillthis="hctau_";
2733  fillthis+=ptbin;
2734  ((TH1F*)listout->FindObject(fillthis))->Fill(CtauD0z);
2735 
2736 
2737 
2738  //correlation ctauXY ctauZ
2739  fillthis="hctauXYctauZ_" ;
2740  fillthis+=ptbin;
2741  ((TH2F*)listout->FindObject(fillthis))->Fill(CtauD0,CtauD0z);
2742 
2743 
2744 
2745 
2746  Int_t pdgDgD0toKpi[2]={321,211};
2747  Int_t labD0=-1;
2748  Bool_t isPrimary=kTRUE;
2749  // Double_t impparXY=part->ImpParXY()*10000.;
2750  Double_t trueImpParXY=0.;
2751  Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
2752  Double_t arrayForSparse[3]={invmassD0,pt,impparXY};
2753  Double_t arrayForSparseTrue[3]={invmassD0,pt,trueImpParXY};
2754  if (fReadMC) labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
2755  // cout << "labD0 "<<labD0 <<endl;
2756  fNentries->Fill(1);
2757  //count true D0 selected by cuts
2758  if (fReadMC && labD0>=0) fNentries->Fill(2);
2759 
2760  if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
2761 
2762  // cout << "confronto ctau primarie e secondarie DOPO if" <<endl;
2763 
2764  //arrayForSparse[0]=invmassD0; arrayForSparse[2]=impparXY;
2765  if(fReadMC){
2766  if(labD0>=0) {
2767  if(fArray==1) cout<<"LS signal ERROR"<<endl;
2768 
2769  AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
2770  Int_t pdgD0 = partD0->GetPdgCode();
2771 
2772 
2773 
2774 
2775  }
2776  }
2777 
2778  }
2779 
2780 
2781  Double_t weigD0=1.;
2782  Double_t weigD0bar=1.;
2783  if (fCuts->GetCombPID() && (fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeight || fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeightNoFilter)) {
2784  weigD0=fCuts->GetWeightsNegative()[AliPID::kKaon] * fCuts->GetWeightsPositive()[AliPID::kPion];
2785  weigD0bar=fCuts->GetWeightsPositive()[AliPID::kKaon] * fCuts->GetWeightsNegative()[AliPID::kPion];
2786  if (weigD0 > 1.0 || weigD0 < 0.) {weigD0 = 0.;}
2787  if (weigD0bar > 1.0 || weigD0bar < 0.) {weigD0bar = 0.;} //Prevents filling with weight > 1, or < 0
2788  }
2789 
2790 
2791 
2792 
2793 
2794 
2795  } //readMC
2796 
2797  if(fSys==0){
2798  //normalise pt distr to half afterwards
2799  fillthis="hptB_";
2800  fillthis+=ptbin;
2801  ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[0]);
2802  ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[1]);
2803 
2804  fillthis="hcosthetastarB_";
2805  fillthis+=ptbin;
2806  if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3)))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
2807  if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
2808 
2809 
2810  fillthis="hd0p0B_";
2811  fillthis+=ptbin;
2812  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
2813  fillthis="hd0p1B_";
2814  fillthis+=ptbin;
2815  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
2816 
2817  fillthis="hcosthpointd0d0B_";
2818  fillthis+=ptbin;
2819  ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
2820 
2821  fillthis="hcosthpointd0B_";
2822  fillthis+=ptbin;
2823  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[0]);
2824  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[1]);
2825 
2826 
2827  if(recalcvtx){
2828 
2829  fillthis="hd0vp0B_";
2830  fillthis+=ptbin;
2831  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
2832  fillthis="hd0vp1B_";
2833  fillthis+=ptbin;
2834  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
2835 
2836  fillthis="hd0vB_";
2837  fillthis+=ptbin;
2838  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
2839  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
2840 
2841  }
2842 
2843  }
2844 
2845  fillthis="hdcaB_";
2846  fillthis+=ptbin;
2847  ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
2848 
2849  fillthis="hd0d0B_";
2850  fillthis+=ptbin;
2851  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]*d0Prong[1]);
2852 
2853  if(recalcvtx){
2854  fillthis="hd0d0vB_";
2855  fillthis+=ptbin;
2856  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
2857  }
2858 
2859  fillthis="hcosthetapointB_";
2860  fillthis+=ptbin;
2861  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
2862 
2863  fillthis="hcosthetapointxyB_";
2864  fillthis+=ptbin;
2865  ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
2866 
2867  fillthis="hdeclB_";
2868  fillthis+=ptbin;
2869  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
2870 
2871  fillthis="hnormdeclB_";
2872  fillthis+=ptbin;
2873  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
2874 
2875 
2876 
2877 
2878 
2879 
2880 
2882  fillthis="hdeclxyB_";
2883  fillthis+=ptbin;
2884  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
2885 
2886  fillthis="hnormdeclxyB_";
2887  fillthis+=ptbin;
2888  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
2889 
2890  fillthis="hdeclxyd0d0B_";
2891  fillthis+=ptbin;
2892  ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
2893 
2894  fillthis="hnormdeclxyd0d0B_";
2895  fillthis+=ptbin;
2896  ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
2897 
2898 
2899  if(recalcvtx) {
2900 
2901  fillthis="hdeclvB_";
2902  fillthis+=ptbin;
2903  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
2904 
2905  fillthis="hnormdeclvB_";
2906  fillthis+=ptbin;
2907  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
2908 
2909 
2910  }
2911  }//mass cut
2912  }//else (background)
2913 
2914  return;
2915 }
2916 
2917 //____________________________________________________________________________
2918 void AliAnalysisTaskSED0MassNonPromptFraction::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAODMCHeader *mcHeader, AliRDHFCutsD0toKpi* cuts, TList *listout){
2919  //
2921  //
2922 
2923 Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2924 
2925  /*
2926 
2927  //cout<<"is selected = "<<fIsSelectedCandidate<<endl;
2928 
2929  // Fill candidate variable Tree (track selection, no candidate selection)
2930  if( fWriteVariableTree && !part->HasBadDaughters()
2931  && fCuts->AreDaughtersSelected(part) && fCuts->IsSelectedPID(part) ){
2932  fCandidateVariables[0] = part->InvMassD0(); //ok
2933  fCandidateVariables[1] = part->InvMassD0bar(); //ok
2934  fCandidateVariables[2] = part->Pt(); //ok
2935 
2936  // fCandidateVariables[3] = part->Pt2Prong(0);
2937  // fCandidateVariables[4] = part->Pt2Prong(1);
2938  // fCandidateVariables[5] = part->Getd0Prong(0);
2939  // fCandidateVariables[6] = part->Getd0Prong(1);
2940  // fCandidateVariables[7] = part->Prodd0d0();
2941 
2942  // fCandidateVariables[8] = part->DecayLengthXY();
2943  // Double_t PseudoProperDecLen =
2944 
2945  fCandidateVariables[3] = part->CtD0();
2946  fCandidateVariables[4] = PseudoProperDecLen;
2947  // fCandidateVariables[10]= part->ImpParXY()*10000.;
2948  // fCandidateVariables[14] = fCuts->IsSelectedSpecialCuts(part);
2949  fVariablesTree->Fill();
2950  }
2951 
2952  */
2953 
2954 
2955  //cout<<"check cuts = "<<endl;
2956  //cuts->PrintAll();
2957  if (!fIsSelectedCandidate){
2958  //cout<<" cut " << cuts << " Rejected because "<<cuts->GetWhy()<<endl;
2959  return;
2960  }
2961 
2962 
2963  if(fDebug>2) cout<<"Candidate selected"<<endl;
2964 
2965  Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
2966  //printf("SELECTED\n");
2967  Int_t ptbin=cuts->PtBin(part->Pt());
2968  Double_t pt = part->Pt();
2969  Double_t y = part->YD0();
2970 
2971  Double_t impparXY=part->ImpParXY()*10000.;
2972  Double_t trueImpParXY=0.;
2973  Double_t arrayForSparse[3]={invmassD0,pt,impparXY};
2974  Double_t arrayForSparseTrue[3]={invmassD0,pt,trueImpParXY};
2975 
2976 
2977 
2978  TString fillthis="", fillthispt="", fillthismasspt="", fillthismassy="";
2979  Int_t pdgDgD0toKpi[2]={321,211};
2980  Int_t labD0=-1;
2981 
2982 
2983 
2984  Double_t ctau = -1;
2985  Double_t ptProng[2]={-99,-99};
2986  AliAODTrack *prong0=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
2987  AliAODTrack *prong1=(AliAODTrack*)fDaughterTracks.UncheckedAt(1);
2988  Double_t CtauD0 = -1;
2989 
2990 
2991 
2992 
2993  Bool_t isPrimary=kTRUE;
2994  if (fReadMC) labD0 = part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
2995 
2996  //Define weights for Bayesian (if appropriate)
2997 
2998  Double_t weigD0=1.;
2999  Double_t weigD0bar=1.;
3000  if (fCuts->GetCombPID() && (fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeight || fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeightNoFilter)) {
3001  weigD0=fCuts->GetWeightsNegative()[AliPID::kKaon] * fCuts->GetWeightsPositive()[AliPID::kPion];
3002  weigD0bar=fCuts->GetWeightsPositive()[AliPID::kKaon] * fCuts->GetWeightsNegative()[AliPID::kPion];
3003  if (weigD0 > 1.0 || weigD0 < 0.) {weigD0 = 0.;}
3004  if (weigD0bar > 1.0 || weigD0bar < 0.) {weigD0bar = 0.;} //Prevents filling with weight > 1, or < 0
3005  }
3006 
3007  //count candidates selected by cuts
3008  fNentries->Fill(1);
3009  //count true D0 selected by cuts
3010  if (fReadMC && labD0>=0) fNentries->Fill(2);
3011 
3012  if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
3013 
3014  arrayForSparse[0]=invmassD0; arrayForSparse[2]=impparXY;
3015 
3016  if(fReadMC){
3017  if(labD0>=0) {
3018  if(fArray==1) cout<<"LS signal ERROR"<<endl;
3019 
3020  AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
3021  Int_t pdgD0 = partD0->GetPdgCode();
3022  // cout<<"pdg = "<<pdgD0<<endl;
3023  // ctau=part->CtD0();
3024 
3025 
3026 
3027  // old function if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE; //la funzione checkorigin restituisce 5 se la D deriva da B, ovvero se è secondaria!
3028  if(AliVertexingHFUtils::CheckOrigin(arrMC,partD0,fUseQuarkTagInKine)==5) isPrimary=kFALSE;
3029 
3030 
3031  if(!isPrimary)
3032  trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
3033  arrayForSparseTrue[0]=invmassD0; arrayForSparseTrue[2]=trueImpParXY;
3034 
3035 
3036 
3037  if (pdgD0==421){ //D0
3038  // cout<<"Fill S with D0"<<endl;
3039  fillthis="histSgn_";
3040  fillthis+=ptbin;
3041  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
3042 
3043 
3044 
3045  if(fFillPtHist){
3046  fillthismasspt="histSgnPt";
3047  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
3048  }
3049  if(fFillImpParHist){
3050  if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0);
3051  else {
3052  fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0);
3053  fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0);
3054  }
3055  }
3056 
3057  if(fFillYHist){
3058  fillthismassy="histSgnY_";
3059  fillthismassy+=ptbin;
3060  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
3061  }
3062 
3063  if(fSys==0){
3064  if(TMath::Abs(invmassD0 - mPDG) < 0.027 && fFillVarHists){
3065  fillthis="histSgn27_";
3066  fillthis+=ptbin;
3067  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
3068  }
3069  }
3070  } else{ //it was a D0bar
3071  fillthis="histRfl_";
3072  fillthis+=ptbin;
3073  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
3074 
3075  if(fFillPtHist){
3076  fillthismasspt="histRflPt";
3077  // cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
3078  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
3079  }
3080 
3081  if(fFillYHist){
3082  fillthismassy="histRflY_";
3083  fillthismassy+=ptbin;
3084  // cout << " Filling "<<fillthismassy<<" D0bar"<<endl;
3085  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
3086  }
3087 
3088  }
3089  } else {//background
3090  fillthis="histBkg_";
3091  fillthis+=ptbin;
3092  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
3093 
3094  if(fFillPtHist){
3095  fillthismasspt="histBkgPt";
3096  // cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
3097  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
3098  }
3099  if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0);
3100 
3101  if(fFillYHist){
3102  fillthismassy="histBkgY_";
3103  fillthismassy+=ptbin;
3104  // cout << " Filling "<<fillthismassy<<" D0bar"<<endl;
3105  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
3106  }
3107 
3108  }
3109 
3110  }else{
3111  fillthis="histMass_";
3112  fillthis+=ptbin;
3113  // cout<<"Filling "<<fillthis<<endl;
3114 
3115  // printf("Fill mass with D0");
3116  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
3117 
3118 
3119  if(fFillPtHist){
3120  fillthismasspt="histMassPt";
3121  // cout<<"Filling "<<fillthismasspt<<endl;
3122  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
3123  }
3124  if(fFillImpParHist) {
3125  // cout << "Filling fHistMassPtImpParTC[0]"<<endl;
3126  fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0);
3127  }
3128 
3129  if(fFillYHist){
3130  fillthismassy="histMassY_";
3131  fillthismassy+=ptbin;
3132  // cout<<"Filling "<<fillthismassy<<endl;
3133  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
3134  }
3135 
3136  }
3137 
3138  }
3139  if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
3140 
3141  arrayForSparse[0]=invmassD0bar; arrayForSparse[2]=impparXY;
3142 
3143  if(fReadMC){
3144  if(labD0>=0) {
3145  if(fArray==1) cout<<"LS signal ERROR"<<endl;
3146  AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
3147  Int_t pdgD0 = partD0->GetPdgCode();
3148  // cout<<" pdg = "<<pdgD0<<endl;
3149 
3150 
3151  //old function if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE;
3152  if(AliVertexingHFUtils::CheckOrigin(arrMC,partD0,fUseQuarkTagInKine)==5) isPrimary=kFALSE;
3153  if(!isPrimary)
3154  trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
3155  arrayForSparseTrue[0]=invmassD0bar; arrayForSparseTrue[2]=trueImpParXY;
3156 
3157  if (pdgD0==-421){ //D0bar
3158  fillthis="histSgn_";
3159  fillthis+=ptbin;
3160  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
3161 
3162 
3163  if(fFillPtHist){
3164  fillthismasspt="histSgnPt";
3165  // cout<<" Filling "<< fillthismasspt << endl;
3166  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
3167  }
3168  if(fFillImpParHist){
3169  // cout << " Filling impact parameter thnsparse"<<endl;
3170  if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0bar);
3171  else {
3172  fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0bar);
3173  fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0bar);
3174  }
3175  }
3176 
3177  if(fFillYHist){
3178  fillthismassy="histSgnY_";
3179  fillthismassy+=ptbin;
3180  // cout<<" Filling "<< fillthismassy << endl;
3181  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
3182  }
3183 
3184  } else{
3185  fillthis="histRfl_";
3186  fillthis+=ptbin;
3187  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
3188  if(fFillPtHist){
3189  fillthismasspt="histRflPt";
3190  // cout << " Filling "<<fillthismasspt<<endl;
3191  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
3192  }
3193  if(fFillYHist){
3194  fillthismassy="histRflY_";
3195  fillthismassy+=ptbin;
3196  // cout << " Filling "<<fillthismassy<<endl;
3197  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
3198  }
3199  }
3200  } else {//background or LS
3201  fillthis="histBkg_";
3202  fillthis+=ptbin;
3203  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
3204 
3205  if(fFillPtHist){
3206  fillthismasspt="histBkgPt";
3207  // cout<<" Filling "<< fillthismasspt << endl;
3208  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
3209  }
3210  if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0bar);
3211  if(fFillYHist){
3212  fillthismassy="histBkgY_";
3213  fillthismassy+=ptbin;
3214  // cout<<" Filling "<< fillthismassy << endl;
3215  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
3216  }
3217  }
3218  }else{
3219  fillthis="histMass_";
3220  fillthis+=ptbin;
3221  // printf("Fill mass with D0bar");
3222 
3223  ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar,weigD0bar);
3224 
3225 
3226  if(fFillPtHist){
3227  fillthismasspt="histMassPt";
3228  // cout<<" Filling "<< fillthismasspt << endl;
3229  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
3230  }
3231  if(fFillImpParHist) fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0bar);
3232  if(fFillYHist){
3233  fillthismassy="histMassY_";
3234  fillthismassy+=ptbin;
3235  // cout<<" Filling "<< fillthismassy << endl;
3236  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
3237  }
3238  }
3239  }
3240 
3241  return;
3242 }
3243 
3244 //__________________________________________________________________________
3245 AliAODVertex* AliAnalysisTaskSED0MassNonPromptFraction::GetPrimaryVtxSkipped(AliAODEvent *aodev){
3247 
3248  Int_t skipped[2];
3249  Int_t nTrksToSkip=2;
3250  AliAODTrack *dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(0);
3251  if(!dgTrack){
3252  AliDebug(2,"no daughter found!");
3253  return 0x0;
3254  }
3255  skipped[0]=dgTrack->GetID();
3256  dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(1);
3257  if(!dgTrack){
3258  AliDebug(2,"no daughter found!");
3259  return 0x0;
3260  }
3261  skipped[1]=dgTrack->GetID();
3262 
3263  AliESDVertex *vertexESD=0x0;
3264  AliAODVertex *vertexAOD=0x0;
3265  AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
3266 
3267  //
3268  vertexer->SetSkipTracks(nTrksToSkip,skipped);
3269  vertexer->SetMinClusters(4);
3270  vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
3271  if(!vertexESD) return vertexAOD;
3272  if(vertexESD->GetNContributors()<=0) {
3273  AliDebug(2,"vertexing failed");
3274  delete vertexESD; vertexESD=NULL;
3275  return vertexAOD;
3276  }
3277 
3278  delete vertexer; vertexer=NULL;
3279 
3280 
3281  // convert to AliAODVertex
3282  Double_t pos[3],cov[6],chi2perNDF;
3283  vertexESD->GetXYZ(pos); // position
3284  vertexESD->GetCovMatrix(cov); //covariance matrix
3285  chi2perNDF = vertexESD->GetChi2toNDF();
3286  delete vertexESD; vertexESD=NULL;
3287 
3288  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
3289  return vertexAOD;
3290 
3291 }
3292 
3293 
3294 //________________________________________________________________________
3295 void AliAnalysisTaskSED0MassNonPromptFraction::Terminate(Option_t */*option*/)
3296 {
3298  //
3299  if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
3300 
3301 
3302  fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
3303  if (!fOutputMass) {
3304  printf("ERROR: fOutputMass not available\n");
3305  return;
3306  }
3307  fOutputMassPt = dynamic_cast<TList*> (GetOutputData(6));
3308  if ((fFillPtHist || fFillImpParHist) && !fOutputMassPt) {
3309  printf("ERROR: fOutputMass not available\n");
3310  return;
3311  }
3312 
3313  if(fFillVarHists){
3314  fDistr = dynamic_cast<TList*> (GetOutputData(2));
3315  if (!fDistr) {
3316  printf("ERROR: fDistr not available\n");
3317  return;
3318  }
3319  }
3320 
3321  fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
3322 
3323  if(!fNentries){
3324  printf("ERROR: fNEntries not available\n");
3325  return;
3326  }
3327  fCuts = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(4));
3328  if(!fCuts){
3329  printf("ERROR: fCuts not available\n");
3330  return;
3331  }
3332  fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(5));
3333  if (!fCounter) {
3334  printf("ERROR: fCounter not available\n");
3335  return;
3336  }
3337  if (fDrawDetSignal) {
3338  fDetSignal = dynamic_cast<TList*>(GetOutputData(8));
3339  if (!fDetSignal) {
3340  printf("ERROR: fDetSignal not available\n");
3341  return;
3342  }
3343  }
3344  if(fFillYHist){
3345  fOutputMassY = dynamic_cast<TList*> (GetOutputData(9));
3346  if (fFillYHist && !fOutputMassY) {
3347  printf("ERROR: fOutputMassY not available\n");
3348  return;
3349  }
3350  }
3351 
3352  Int_t nptbins=fCuts->GetNPtBins();
3353  for(Int_t ipt=0;ipt<nptbins;ipt++){
3354 
3355  if(fArray==1 && fFillVarHists){
3356  fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fOutputMass->FindObject("hpospair"))->Integral(nptbins+ipt+1,nptbins+ipt+2)*((TH1F*)fOutputMass->FindObject("hnegpair"))->Integral(nptbins+ipt+1,nptbins+ipt+2)); //after cuts
3357 
3358 
3359  if(fLsNormalization>1e-6) {
3360 
3361  TString massName="histMass_";
3362  massName+=ipt;
3363  ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
3364 
3365  }
3366 
3367 
3368  fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fOutputMass->FindObject("hpospair"))->Integral(ipt+1,ipt+2)*((TH1F*)fOutputMass->FindObject("hnegpair"))->Integral(ipt+1,ipt+2));
3369  //fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
3370 
3371  if(fLsNormalization>1e-6) {
3372 
3373  TString nameDistr="hdcaB_";
3374  nameDistr+=ipt;
3375  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
3376  nameDistr="hd0B_";
3377  nameDistr+=ipt;
3378  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
3379  nameDistr="hd0d0B_";
3380  nameDistr+=ipt;
3381  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
3382  nameDistr="hcosthetapointB_";
3383  nameDistr+=ipt;
3384  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
3385  if(fSys==0){
3386  nameDistr="hptB_";
3387  nameDistr+=ipt;
3388  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
3389  nameDistr="hcosthetastarB_";
3390  nameDistr+=ipt;
3391  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
3392  nameDistr="hcosthpointd0d0B_";
3393  nameDistr+=ipt;
3394  ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
3395  }
3396  }
3397  }
3398  }
3399  TString cvname,cstname;
3400 
3401  if (fArray==0){
3402  cvname="D0invmass";
3403  cstname="cstat0";
3404  } else {
3405  cvname="LSinvmass";
3406  cstname="cstat1";
3407  }
3408 
3409  TCanvas *cMass=new TCanvas(cvname,cvname);
3410  cMass->cd();
3411  ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
3412 
3413  TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
3414  cStat->cd();
3415  cStat->SetGridy();
3416  fNentries->Draw("htext0");
3417 
3418  // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
3419  // ccheck->cd();
3420 
3421  return;
3422 }
3423 
3424 
3425 //________________________________________________________________________
3426 void AliAnalysisTaskSED0MassNonPromptFraction::CreateImpactParameterHistos(){
3428 
3429  Int_t nmassbins=200;
3430  Double_t fLowmasslimit=1.5648, fUpmasslimit=2.1648;
3431  Int_t fNImpParBins=400;
3432  Double_t fLowerImpPar=-2000., fHigherImpPar=2000.;
3433  Int_t nbins[3]={nmassbins,200,fNImpParBins};
3434  Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar};
3435  Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar};
3436 
3437 
3438  fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
3439  "Mass vs. pt vs.imppar - All",
3440  3,nbins,xmin,xmax);
3441  fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
3442  "Mass vs. pt vs.imppar - promptD",
3443  3,nbins,xmin,xmax);
3444  fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
3445  "Mass vs. pt vs.imppar - DfromB",
3446  3,nbins,xmin,xmax);
3447  fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
3448  "Mass vs. pt vs.true imppar -DfromB",
3449  3,nbins,xmin,xmax);
3450  fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
3451  "Mass vs. pt vs.imppar - backgr.",
3452  3,nbins,xmin,xmax);
3453 
3454  for(Int_t i=0; i<5;i++){
3455  fOutputMassPt->Add(fHistMassPtImpParTC[i]);
3456  }
3457 }
3458 
3459 //_________________________________________________________________________________________________
3460 Float_t AliAnalysisTaskSED0MassNonPromptFraction::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD0) const {
3462 
3463  printf(" AliAnalysisTaskSED0MassNonPromptFractionV1::GetTrueImpactParameter() \n");
3464 
3465  Double_t vtxTrue[3];
3466  mcHeader->GetVertex(vtxTrue);
3467  Double_t origD[3];
3468  partD0->XvYvZv(origD);
3469  Short_t charge=partD0->Charge();
3470  Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
3471  for(Int_t iDau=0; iDau<2; iDau++){
3472  pXdauTrue[iDau]=0.;
3473  pYdauTrue[iDau]=0.;
3474  pZdauTrue[iDau]=0.;
3475  }
3476 
3477  // Int_t nDau=partD0->GetNDaughters();
3478  Int_t labelFirstDau = partD0->GetDaughter(0);
3479 
3480  for(Int_t iDau=0; iDau<2; iDau++){
3481  Int_t ind = labelFirstDau+iDau;
3482  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
3483  if(!part) continue;
3484  Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
3485  if(!part){
3486  AliError("Daughter particle not found in MC array");
3487  return 99999.;
3488  }
3489  if(pdgCode==211 || pdgCode==321){
3490  pXdauTrue[iDau]=part->Px();
3491  pYdauTrue[iDau]=part->Py();
3492  pZdauTrue[iDau]=part->Pz();
3493  }
3494  }
3495 
3496  Double_t d0dummy[2]={0.,0.};
3497  AliAODRecoDecayHF aodDzeroMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
3498  return aodDzeroMC.ImpParXY();
3499 
3500 }
3501 
3502 //_________________________________________________________________________________________________
3503 Int_t AliAnalysisTaskSED0MassNonPromptFraction::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
3504  //
3506  //
3507  //printf(" AliAnalysisTaskSED0MassNonPromptFraction V1::CheckOrigin() \n");
3508 
3509  Int_t pdgGranma = 0;
3510  Int_t mother = 0;
3511  mother = mcPartCandidate->GetMother();
3512  Int_t istep = 0;
3513  Int_t abspdgGranma =0;
3514  Bool_t isFromB=kFALSE;
3515  Bool_t isQuarkFound=kFALSE;
3516  while (mother >0 ){
3517  istep++;
3518  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
3519  if (mcGranma){
3520  pdgGranma = mcGranma->GetPdgCode();
3521  abspdgGranma = TMath::Abs(pdgGranma);
3522  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
3523  isFromB=kTRUE;
3524  }
3525  if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
3526  mother = mcGranma->GetMother();
3527  }else{
3528  AliError("Failed casting the mother particle!");
3529  break;
3530  }
3531  }
3532 
3533  if(isFromB) return 5;
3534  else return 4;
3535 }
charge
Int_t charge
Definition: MuonTrackingEfficiency.C:58
AliAnalysisTaskSED0MassNonPromptFraction::fDetSignal
TList * fDetSignal
flag to decide whether to draw the TPC dE/dx and TOF signal before/after PID
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:134
AliAnalysisTaskSED0MassNonPromptFraction::fFillYHist
Bool_t fFillYHist
flag to fill Pt and Impact Parameter Histograms
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:123
AliAnalysisTaskSED0MassNonPromptFraction::fCandidateVariables
Double_t * fCandidateVariables
! variables to be written to the tree
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:131
Double_t
double Double_t
Definition: External.C:58
AliAnalysisTaskSED0MassNonPromptFraction::fDrawDetSignal
Bool_t fDrawDetSignal
flag to decide whether to fill "PID = x" bins in fNentrie
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:133
Draw
void Draw(const char *filename, const char *title="", const char *others="ALL", const char *options="DEFAULT", const char *outFlg="ALL", UShort_t rebin=5, Float_t eff=0, const char *base="")
Definition: DrawdNdeta.C:3603
AliNormalizationCounter::StoreCandidates
void StoreCandidates(AliVEvent *, Int_t nCand=0, Bool_t flagFilter=kTRUE)
Definition: AliNormalizationCounter.cxx:289
TH2F
Definition: External.C:236
AliRDHFCutsD0toKpi::GetCombPID
Bool_t GetCombPID() const
Definition: AliRDHFCutsD0toKpi.h:107
AliAnalysisTaskSED0MassNonPromptFraction::fNentries
TH1F * fNentries
! histogram with number of events on output slot 3
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:104
AliAODRecoCascadeHF.h
AliAODRecoDecayHF::HasSelectionBit
Bool_t HasSelectionBit(Int_t i) const
Definition: AliAODRecoDecayHF.h:123
AliAnalysisVertexingHF.h
AliVertexingHFUtils.h
AliAODRecoDecayHF::NormalizedDecayLength2
Double_t NormalizedDecayLength2() const
Definition: AliAODRecoDecayHF.h:69
AliAODRecoDecayHF::ImpParXY
Double_t ImpParXY() const
Definition: AliAODRecoDecayHF.h:83
AliRDHFCuts::CheckMatchingAODdeltaAODevents
static Int_t CheckMatchingAODdeltaAODevents()
Definition: AliRDHFCuts.cxx:761
AliAnalysisTaskSED0MassNonPromptFraction::fFillVarHists
Bool_t fFillVarHists
flag for filling var hists only after cuts
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:118
AliAnalysisTaskSED0MassNonPromptFraction::fVariablesTree
TTree * fVariablesTree
flag to decide whether to write the candidate variables on a tree variables
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:130
AliAnalysisTaskSED0MassNonPromptFraction::fCounter
AliNormalizationCounter * fCounter
! AliNormalizationCounter on output slot 5
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:111
AliRDHFCuts::GetWhyRejection
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:315
AliAODRecoDecayHF::CosPointingAngleXY
Double_t CosPointingAngleXY() const
Definition: AliAODRecoDecayHF.h:81
AliAnalysisVertexingHF::FillRecoCand
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
Definition: AliAnalysisVertexingHF.cxx:1698
AliAODRecoDecayHF::GetSelectionMap
ULong_t GetSelectionMap() const
Definition: AliAODRecoDecayHF.h:124
AliRDHFCutsD0toKpi::IsInFiducialAcceptance
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
Definition: AliRDHFCutsD0toKpi.cxx:798
AliVertexingHFUtils::CheckOrigin
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Functions to check the decay tree.
Definition: AliVertexingHFUtils.cxx:637
AliAnalysisTaskSED0MassNonPromptFraction::GetTrueImpactParameter
Float_t GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray *arrayMC, AliAODMCParticle *partD0) const
Definition: AliAnalysisTaskSED0MassNonPromptFraction.cxx:3460
AliAnalysisTaskSED0MassNonPromptFraction::FillMCGenHists
void FillMCGenHists(AliAODEvent *aod, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout)
Definition: AliAnalysisTaskSED0MassNonPromptFraction.cxx:1512
AliAnalysisTaskSED0MassNonPromptFraction::fAODProtection
Int_t fAODProtection
flag to check or not the selection bit
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:126
AliRDHFCutsD0toKpi::GetBayesianStrategy
Int_t GetBayesianStrategy() const
Definition: AliRDHFCutsD0toKpi.h:88
AliAODRecoDecayHF::HasBadDaughters
Bool_t HasBadDaughters() const
Definition: AliAODRecoDecayHF.h:247
AliAnalysisTaskSED0MassNonPromptFraction::fDaughterTracks
TObjArray fDaughterTracks
flag to fill mass histogram with D0/D0bar only (0 = fill with both, 1 = fill with D0 only...
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:115
Int_t
int Int_t
Definition: External.C:63
TH1I
Definition: External.C:204
AliAnalysisTaskSED0MassNonPromptFraction.h
UInt_t
unsigned int UInt_t
Definition: External.C:33
AliAODRecoDecayHF2Prong::CosThetaStarD0bar
Double_t CosThetaStarD0bar() const
angle of K
Definition: AliAODRecoDecayHF2Prong.h:44
Float_t
float Float_t
Definition: External.C:68
AliAnalysisTaskSED0MassNonPromptFraction::fUsePid4Distr
Bool_t fUsePid4Distr
flag to decide if apply cut also on distributions: 0 no cuts, 1 looser cuts, 2 tighter/ cuts ...
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:110
AliAnalysisTaskSED0MassNonPromptFraction::fUseSelectionBit
Bool_t fUseSelectionBit
flag to fill Pt and Impact Parameter Histograms
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:125
AliAODRecoDecayHF2Prong.h
AliAnalysisTaskSED0MassNonPromptFraction::FillVarHists
void FillVarHists(AliAODEvent *aodev, AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout)
Definition: AliAnalysisTaskSED0MassNonPromptFraction.cxx:1665
AliAODRecoDecayHF::DecayLength2
Double_t DecayLength2() const
kinematics & topology
Definition: AliAODRecoDecayHF.h:61
Data
Bool_t Data(TH1F *h, Double_t *rangefit, Bool_t writefit, Double_t &sgn, Double_t &errsgn, Double_t &bkg, Double_t &errbkg, Double_t &sgnf, Double_t &errsgnf, Double_t &sigmafit, Int_t &status)
Definition: charmCutsOptimization.C:371
Short_t
short Short_t
Definition: External.C:23
AliAODRecoDecayHF::SetOwnPrimaryVtx
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Definition: AliAODRecoDecayHF.h:44
AliAnalysisTaskSED0MassNonPromptFraction::UserCreateOutputObjects
virtual void UserCreateOutputObjects()
Implementation of interface methods.
Definition: AliAnalysisTaskSED0MassNonPromptFraction.cxx:253
AliAODRecoDecayHF::DecayLengthXY
Double_t DecayLengthXY() const
Definition: AliAODRecoDecayHF.h:71
AliRDHFCuts::GetIsPrimaryWithoutDaughters
Bool_t GetIsPrimaryWithoutDaughters() const
Definition: AliRDHFCuts.h:278
AliRDHFCutsD0toKpi::GetWeightsPositive
Double_t * GetWeightsPositive() const
Definition: AliRDHFCutsD0toKpi.h:80
AliRDHFCuts::AreDaughtersSelected
Bool_t AreDaughtersSelected(AliAODRecoDecayHF *rd, const AliAODEvent *aod=0x0) const
Definition: AliRDHFCuts.cxx:731
AliRDHFCuts::IsEventSelected
Bool_t IsEventSelected(AliVEvent *event)
Definition: AliRDHFCuts.cxx:500
AliNormalizationCounter::StoreEvent
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
Definition: AliNormalizationCounter.cxx:179
AliAnalysisTaskSED0MassNonPromptFraction::fHistMassPtImpParTC
THnSparseF * fHistMassPtImpParTC[5]
! histograms for impact paramter studies
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:106
AliAnalysisTaskSED0MassNonPromptFraction::fFillPtHist
Bool_t fFillPtHist
flag to reject events with SDD clusters
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:122
nVar
const Int_t nVar
Definition: AliAnalysisTaskSubJetFraction.h:22
AliRDHFCutsD0toKpi::IsSelectedPID
virtual Int_t IsSelectedPID(AliAODRecoDecayHF *rd)
Definition: AliRDHFCutsD0toKpi.cxx:828
AliAnalysisTaskSED0MassNonPromptFraction::CheckOrigin
Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPartCandidate) const
Definition: AliAnalysisTaskSED0MassNonPromptFraction.cxx:3503
AliNormalizationCounter.h
AliAnalysisTaskSED0MassNonPromptFraction::fCutOnDistr
Bool_t fCutOnDistr
flag for MC array: kTRUE = read it, kFALSE = do not read it
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:109
AliAnalysisTaskSED0MassNonPromptFraction::FillMassHists
void FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAODMCHeader *mcHeader, AliRDHFCutsD0toKpi *cuts, TList *listout)
Definition: AliAnalysisTaskSED0MassNonPromptFraction.cxx:2918
Scale
TH2 * Scale(TH2 *h, TH1 *g)
Definition: SimpleCorrect.C:568
AliAnalysisTaskSED0MassNonPromptFraction::fFill
Bool_t fFill
flag to enable filling variable histos
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:119
AliRDHFCutsD0toKpi::IsSelected
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Definition: AliRDHFCutsD0toKpi.h:39
AliAnalysisTaskSED0MassNonPromptFraction::DrawDetSignal
void DrawDetSignal(AliAODRecoDecayHF2Prong *part, TList *ListDetSignal)
Definition: AliAnalysisTaskSED0MassNonPromptFraction.cxx:1627
Option_t
const char Option_t
Definition: External.C:48
AliRDHFCuts::GetNPtBins
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:249
nbins
const Int_t nbins
Definition: AverageDmesonRaa.C:63
Bool_t
bool Bool_t
Definition: External.C:53
AliAODRecoDecayHF::CosPointingAngle
Double_t CosPointingAngle() const
Definition: AliAODRecoDecayHF.h:79
AliAODRecoDecayHF::DecayLengthXYError
Double_t DecayLengthXYError() const
Definition: AliAODRecoDecayHF.h:73
AliAnalysisTaskSED0MassNonPromptFraction::fIsRejectSDDClusters
Bool_t fIsRejectSDDClusters
fSys=0 -> p-p; fSys=1 ->PbPb (in this case fFillVarHists=kFALSE by default: set it to kTRUE after if ...
Definition: AliAnalysisTaskSED0MassNonPromptFraction.h:121
AliRDHFCuts::PtBin
Int_t PtBin(Double_t pt) const
Definition: AliRDHFCuts.cxx:1306
nptbins
Int_t nptbins
Definition: CalculateAveragePt.C:56
AliRDHFCutsD0toKpi::GetWeightsNegative
Double_t * GetWeightsNegative() const
Definition: AliRDHFCutsD0toKpi.h:86