AliPhysics  251aa1e (251aa1e)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskSED0Mass.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 /* $Id$ */
17 
19 //
20 // AliAnalysisTaskSE for D0 candidates invariant mass histogram
21 // and comparison with the MC truth and cut variables distributions.
22 //
23 // Authors: A.Dainese, andrea.dainese@lnl.infn.it
24 // Chiara Bianchin, chiara.bianchin@pd.infn.it (invariant mass)
25 // Carmelo Di Giglio, carmelo.digiglio@ba.infn.it (like sign)
26 // Jeremy Wilkinson, jwilkinson@physi.uni-heidelberg.de (weighted Bayesian
28 
29 #include <Riostream.h>
30 #include <TClonesArray.h>
31 #include <TCanvas.h>
32 #include <TNtuple.h>
33 #include <TTree.h>
34 #include <TList.h>
35 #include <TH1F.h>
36 #include <TH2F.h>
37 #include <TDatabasePDG.h>
38 #include <THnSparse.h>
39 #include "AliVertexingHFUtils.h"
40 #include <AliAnalysisDataSlot.h>
41 #include <AliAnalysisDataContainer.h>
42 #include "AliAnalysisManager.h"
43 #include "AliESDtrack.h"
44 #include "AliVertexerTracks.h"
45 #include "AliAODHandler.h"
46 #include "AliAODEvent.h"
47 #include "AliAODVertex.h"
48 #include "AliAODTrack.h"
49 #include "AliAODMCHeader.h"
50 #include "AliAODMCParticle.h"
52 #include "AliAODRecoCascadeHF.h"
53 #include "AliAnalysisVertexingHF.h"
54 #include "AliAnalysisTaskSE.h"
57 
58 using std::cout;
59 using std::endl;
60 
64 
65 //________________________________________________________________________
68  fOutputMass(0),
69  fOutputMassPt(0),
70  fOutputMassY(0),
71  fDistr(0),
72  fNentries(0),
73  fMCAccPrompt(0),
74  fMCAccBFeed(0),
75  fStepMCAcc(kTRUE),
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  fSys(0),
89  fIsRejectSDDClusters(0),
90  fFillPtHist(kTRUE),
91  fFillYHist(kFALSE),
92  fFillImpParHist(kFALSE),
93  fUseSelectionBit(kTRUE),
94  fAODProtection(1),
95  fWriteVariableTree(kFALSE),
96  fVariablesTree(0),
97  fCandidateVariables(),
98  fPIDCheck(kFALSE),
99  fDrawDetSignal(kFALSE),
100  fUseQuarkTagInKine(kTRUE),
101  fhStudyImpParSingleTrackSign(0),
102  fhStudyImpParSingleTrackCand(0),
103  fhStudyImpParSingleTrackFd(0),
104  fDetSignal(0)
105 {
107  for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
108 
109 }
110 
111 //________________________________________________________________________
113  AliAnalysisTaskSE(name),
114  fOutputMass(0),
115  fOutputMassPt(0),
116  fOutputMassY(0),
117  fDistr(0),
118  fNentries(0),
119  fMCAccPrompt(0),
120  fMCAccBFeed(0),
121  fStepMCAcc(kTRUE),
122  fCuts(0),
123  fArray(0),
124  fReadMC(0),
125  fCutOnDistr(0),
126  fUsePid4Distr(0),
127  fCounter(0),
128  fNPtBins(1),
129  fLsNormalization(1.),
130  fFillOnlyD0D0bar(0),
131  fDaughterTracks(),
132  fIsSelectedCandidate(0),
133  fFillVarHists(kTRUE),
134  fSys(0),
135  fIsRejectSDDClusters(0),
136  fFillPtHist(kTRUE),
137  fFillYHist(kFALSE),
138  fFillImpParHist(kFALSE),
139  fUseSelectionBit(kTRUE),
140  fAODProtection(1),
141  fWriteVariableTree(kFALSE),
142  fVariablesTree(0),
143  fCandidateVariables(),
144  fPIDCheck(kFALSE),
145  fDrawDetSignal(kFALSE),
146  fUseQuarkTagInKine(kTRUE),
147  fhStudyImpParSingleTrackSign(0),
148  fhStudyImpParSingleTrackCand(0),
149  fhStudyImpParSingleTrackFd(0),
150  fDetSignal(0)
151 {
153 
154  fNPtBins=cuts->GetNPtBins();
155 
156  fCuts=cuts;
157  for(Int_t ih=0; ih<5; ih++) fHistMassPtImpParTC[ih]=0x0;
158 
159  // Output slot #1 writes into a TList container (mass with cuts)
160  DefineOutput(1,TList::Class()); //My private output
161  // Output slot #2 writes into a TList container (distributions)
162  DefineOutput(2,TList::Class()); //My private output
163  // Output slot #3 writes into a TH1F container (number of events)
164  DefineOutput(3,TH1F::Class()); //My private output
165  // Output slot #4 writes into a TList container (cuts)
166  DefineOutput(4,AliRDHFCutsD0toKpi::Class()); //My private output
167  // Output slot #5 writes Normalization Counter
168  DefineOutput(5,AliNormalizationCounter::Class());
169  // Output slot #6 stores the mass vs pt and impact parameter distributions
170  DefineOutput(6,TList::Class()); //My private output
171  // Output slot #7 keeps a tree of the candidate variables after track selection
172  DefineOutput(7,TTree::Class()); //My private outpu
173  // Output slot #8 writes into a TList container (Detector signals)
174  DefineOutput(8, TList::Class()); //My private output
175  // Output slot #9 stores the mass vs rapidity (y) distributions
176  DefineOutput(9, TList::Class()); //My private output
177 }
178 
179 //________________________________________________________________________
181 {
182  if (fOutputMass) {
183  delete fOutputMass;
184  fOutputMass = 0;
185  }
186  if (fOutputMassPt) {
187  delete fOutputMassPt;
188  fOutputMassPt = 0;
189  }
190  if (fOutputMassY) {
191  delete fOutputMassY;
192  fOutputMassY = 0;
193  }
194  if (fDistr) {
195  delete fDistr;
196  fDistr = 0;
197  }
198  if (fCuts) {
199  delete fCuts;
200  fCuts = 0;
201  }
202  for(Int_t i=0; i<5; i++){
203  if(fHistMassPtImpParTC[i]) delete fHistMassPtImpParTC[i];
204  }
205  if (fNentries){
206  delete fNentries;
207  fNentries = 0;
208  }
209  if(fCounter){
210  delete fCounter;
211  fCounter=0;
212  }
213  if(fVariablesTree){
214  delete fVariablesTree;
215  fVariablesTree = 0;
216  }
217  if (fDetSignal) {
218  delete fDetSignal;
219  fDetSignal = 0;
220  }
221  delete fMCAccPrompt;
222  delete fMCAccBFeed;
226 
227 }
228 
229 //________________________________________________________________________
231 {
233 
234  if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
235 
236 
238  const char* nameoutput=GetOutputSlot(4)->GetContainer()->GetName();
239  copyfCuts->SetName(nameoutput);
240  // Post the data
241  PostData(4,copyfCuts);
242 
243 
244  return;
245 }
246 
247 //________________________________________________________________________
249 {
250 
252  //
253  if(fDebug > 1) printf("AnalysisTaskSED0Mass::UserCreateOutputObjects() \n");
254 
255  // Several histograms are more conveniently managed in a TList
256  fOutputMass = new TList();
257  fOutputMass->SetOwner();
258  fOutputMass->SetName("listMass");
259 
260  fOutputMassPt = new TList();
261  fOutputMassPt->SetOwner();
262  fOutputMassPt->SetName("listMassPt");
263 
264  fOutputMassY = new TList();
265  fOutputMassY->SetOwner();
266  fOutputMassY->SetName("listMassY");
267 
268  fDistr = new TList();
269  fDistr->SetOwner();
270  fDistr->SetName("distributionslist");
271 
272  fDetSignal = new TList();
273  fDetSignal->SetOwner();
274  fDetSignal->SetName("detectorsignals");
275 
276  TString nameMass=" ",nameSgn27=" ",nameSgn=" ", nameBkg=" ", nameRfl=" ",nameMassNocutsS =" ",nameMassNocutsB =" ", namedistr=" ";
277  TString nameMassPt="", nameSgnPt="", nameBkgPt="", nameRflPt="";
278  TString nameMassY="", nameSgnY="", nameBkgY="", nameRflY="";
279  Int_t nbins2dPt=60; Float_t binInPt=0., binFinPt=30.;
280  Int_t nbins2dY=60; Float_t binInY=-1.5, binFinY=1.5;
281 
282  for(Int_t i=0;i<fCuts->GetNPtBins();i++){
283 
284  nameMass="histMass_";
285  nameMass+=i;
286  nameSgn27="histSgn27_";
287  nameSgn27+=i;
288  nameSgn="histSgn_";
289  nameSgn+=i;
290  nameBkg="histBkg_";
291  nameBkg+=i;
292  nameRfl="histRfl_";
293  nameRfl+=i;
294  nameMassNocutsS="hMassS_";
295  nameMassNocutsS+=i;
296  nameMassNocutsB="hMassB_";
297  nameMassNocutsB+=i;
298 
299  //histograms of cut variable distributions
300  if(fFillVarHists){
301  if(fReadMC){
302 
303  namedistr="hNclsD0vsptS_";
304  namedistr+=i;
305  TH2F *hNclsD0vsptS = new TH2F(namedistr.Data(),"N cls distrubution [S];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
306  namedistr="hNclsD0barvsptS_";
307  namedistr+=i;
308  TH2F *hNclsD0barvsptS = new TH2F(namedistr.Data(),"N cls distrubution [S];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
309 
310  namedistr="hNITSpointsD0vsptS_";
311  namedistr+=i;
312  TH2F *hNITSpointsD0vsptS = new TH2F(namedistr.Data(),"N ITS points distrubution [S];p_{T} [GeV/c];N points",200,0.,20.,7,0.,7.);
313 
314  namedistr="hNSPDpointsD0S_";
315  namedistr+=i;
316  TH1I *hNSPDpointsD0S = new TH1I(namedistr.Data(),"N SPD points distrubution [S]; ;N tracks",4,0,4);
317  hNSPDpointsD0S->GetXaxis()->SetBinLabel(1, "no SPD");
318  hNSPDpointsD0S->GetXaxis()->SetBinLabel(2, "kOnlyFirst");
319  hNSPDpointsD0S->GetXaxis()->SetBinLabel(3, "kOnlySecond");
320  hNSPDpointsD0S->GetXaxis()->SetBinLabel(4, "kBoth");
321 
322  namedistr="hptD0S_";
323  namedistr+=i;
324  TH1F *hptD0S = new TH1F(namedistr.Data(), "p_{T} distribution [S];p_{T} [GeV/c]",200,0.,20.);
325  namedistr="hptD0barS_";
326  namedistr+=i;
327  TH1F *hptD0barS = new TH1F(namedistr.Data(), "p_{T} distribution [S];p_{T} [GeV/c]",200,0.,20.);
328 
329  namedistr="hphiD0S_";
330  namedistr+=i;
331  TH1F *hphiD0S = new TH1F(namedistr.Data(), "#phi distribution [S];#phi [rad]",100,0.,2*TMath::Pi());
332  namedistr="hphiD0barS_";
333  namedistr+=i;
334  TH1F *hphiD0barS = new TH1F(namedistr.Data(), "#phi distribution [S];#phi [rad]",100,0.,2*TMath::Pi());
335 
336 
337  namedistr="hetaphiD0candidateS_";
338  namedistr+=i;
339  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());
340  namedistr="hetaphiD0barcandidateS_";
341  namedistr+=i;
342  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());
343 
344  namedistr="hetaphiD0candidatesignalregionS_";
345  namedistr+=i;
346  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());
347  namedistr="hetaphiD0barcandidatesignalregionS_";
348  namedistr+=i;
349  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());
350 
351  // dca
352  namedistr="hdcaS_";
353  namedistr+=i;
354  TH1F *hdcaS = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
355  // impact parameter
356  namedistr="hd0piS_";
357  namedistr+=i;
358  TH1F *hd0piS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions);d0(#pi) [cm]",200,-0.1,0.1);
359 
360  namedistr="hd0KS_";
361  namedistr+=i;
362  TH1F *hd0KS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons);d0(K) [cm]",200,-0.1,0.1);
363  namedistr="hd0d0S_";
364  namedistr+=i;
365  TH1F *hd0d0S = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
366 
367  //decay lenght
368  namedistr="hdeclS_";
369  namedistr+=i;
370  TH1F *hdeclengthS = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm]",200,0,0.015);
371 
372  namedistr="hnormdeclS_";
373  namedistr+=i;
374  TH1F *hnormdeclengthS = new TH1F(namedistr.Data(), "Normalized Decay Length^{2} distribution;(Decay Length/Err)^{2} ",200,0,12.);
375 
376  namedistr="hdeclxyS_";
377  namedistr+=i;
378  TH1F* hdeclxyS=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
379  namedistr="hnormdeclxyS_";
380  namedistr+=i;
381  TH1F* hnormdeclxyS=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
382 
383  namedistr="hdeclxyd0d0S_";
384  namedistr+=i;
385  TH2F* hdeclxyd0d0S=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);
386 
387  namedistr="hnormdeclxyd0d0S_";
388  namedistr+=i;
389  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);
390 
391  // costhetapoint
392  namedistr="hcosthetapointS_";
393  namedistr+=i;
394  TH1F *hcosthetapointS = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
395 
396  namedistr="hcosthetapointxyS_";
397  namedistr+=i;
398  TH1F *hcosthetapointxyS = new TH1F(namedistr.Data(), "cos#theta_{Point} XYdistribution;cos#theta_{Point}",300,0.,1.);
399 
400  TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",600,1.6248,2.2248); //range (MD0-300MeV, mD0 + 300MeV)
401  tmpMS->Sumw2();
402 
403  fDistr->Add(hNclsD0vsptS);
404  fDistr->Add(hNclsD0barvsptS);
405  fDistr->Add(hNITSpointsD0vsptS);
406  fDistr->Add(hNSPDpointsD0S);
407  fDistr->Add(hptD0S);
408  fDistr->Add(hphiD0S);
409  fDistr->Add(hptD0barS);
410  fDistr->Add(hphiD0barS);
411  fDistr->Add(hetaphiD0candidateS);
412  fDistr->Add(hetaphiD0candidatesignalregionS);
413  fDistr->Add(hetaphiD0barcandidateS);
414  fDistr->Add(hetaphiD0barcandidatesignalregionS);
415 
416  fDistr->Add(hdcaS);
417 
418  fDistr->Add(hd0piS);
419  fDistr->Add(hd0KS);
420 
421  fDistr->Add(hd0d0S);
422 
423  fDistr->Add(hcosthetapointS);
424 
425  fDistr->Add(hcosthetapointxyS);
426 
427  fDistr->Add(hdeclengthS);
428 
429  fDistr->Add(hnormdeclengthS);
430 
431  fDistr->Add(hdeclxyS);
432 
433  fDistr->Add(hnormdeclxyS);
434 
435  fDistr->Add(hdeclxyd0d0S);
436  fDistr->Add(hnormdeclxyd0d0S);
437 
438  fDistr->Add(tmpMS);
439  }
440 
441 
442  //Ncls, phi, pt distrubutions
443 
444  namedistr="hNclsD0vsptB_";
445  namedistr+=i;
446  TH2F *hNclsD0vsptB = new TH2F(namedistr.Data(),"N cls distrubution [B];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
447  namedistr="hNclsD0barvsptB_";
448  namedistr+=i;
449  TH2F *hNclsD0barvsptB = new TH2F(namedistr.Data(),"N cls distrubution [B];p_{T} [GeV/c];N cls",200,0.,20.,100,0.,200.);
450 
451  namedistr="hNITSpointsD0vsptB_";
452  namedistr+=i;
453  TH2F *hNITSpointsD0vsptB = new TH2F(namedistr.Data(),"N ITS points distrubution [B];p_{T} [GeV/c];N points",200,0.,20.,7,0.,7.);
454 
455  namedistr="hNSPDpointsD0B_";
456  namedistr+=i;
457  TH1I *hNSPDpointsD0B = new TH1I(namedistr.Data(),"N SPD points distrubution [B]; ;N tracks",4,0,4);
458  hNSPDpointsD0B->GetXaxis()->SetBinLabel(1, "no SPD");
459  hNSPDpointsD0B->GetXaxis()->SetBinLabel(2, "kOnlyFirst");
460  hNSPDpointsD0B->GetXaxis()->SetBinLabel(3, "kOnlySecond");
461  hNSPDpointsD0B->GetXaxis()->SetBinLabel(4, "kBoth");
462 
463  namedistr="hptD0B_";
464  namedistr+=i;
465  TH1F *hptD0B = new TH1F(namedistr.Data(), "p_{T} distribution [B];p_{T} [GeV/c]",200,0.,20.);
466  namedistr="hptD0barB_";
467  namedistr+=i;
468  TH1F *hptD0barB = new TH1F(namedistr.Data(), "p_{T} distribution [B];p_{T} [GeV/c]",200,0.,20.);
469 
470  namedistr="hphiD0B_";
471  namedistr+=i;
472  TH1F *hphiD0B = new TH1F(namedistr.Data(), "#phi distribution [B];#phi [rad]",100,0.,2*TMath::Pi());
473  namedistr="hphiD0barB_";
474  namedistr+=i;
475  TH1F *hphiD0barB = new TH1F(namedistr.Data(), "#phi distribution [B];#phi [rad]",100,0.,2*TMath::Pi());
476 
477  namedistr="hetaphiD0candidateB_";
478  namedistr+=i;
479  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());
480  namedistr="hetaphiD0barcandidateB_";
481  namedistr+=i;
482  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());
483 
484  namedistr="hetaphiD0candidatesignalregionB_";
485  namedistr+=i;
486  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());
487  namedistr="hetaphiD0barcandidatesignalregionB_";
488  namedistr+=i;
489  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());
490 
491  // dca
492  namedistr="hdcaB_";
493  namedistr+=i;
494  TH1F *hdcaB = new TH1F(namedistr.Data(), "DCA distribution;dca [cm]",200,0.,0.1);
495 
496  // impact parameter
497  namedistr="hd0B_";
498  namedistr+=i;
499  TH1F *hd0B = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
500 
501  namedistr="hd0d0B_";
502  namedistr+=i;
503  TH1F *hd0d0B = new TH1F(namedistr.Data(), "d_{0}#timesd_{0} distribution;d_{0}#timesd_{0} [cm^{2}]",200,-0.001,0.001);
504 
505  //decay lenght
506  namedistr="hdeclB_";
507  namedistr+=i;
508  TH1F *hdeclengthB = new TH1F(namedistr.Data(), "Decay Length^{2} distribution;Decay Length^{2} [cm^{2}]",200,0,0.015);
509 
510  namedistr="hnormdeclB_";
511  namedistr+=i;
512  TH1F *hnormdeclengthB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution;(Decay Length/Err)^{2} ",200,0,12.);
513 
514  namedistr="hdeclxyB_";
515  namedistr+=i;
516  TH1F* hdeclxyB=new TH1F(namedistr.Data(),"Decay Length XY distribution;Decay Length XY [cm]",200,0,0.15);
517  namedistr="hnormdeclxyB_";
518  namedistr+=i;
519  TH1F* hnormdeclxyB=new TH1F(namedistr.Data(),"Normalized decay Length XY distribution;Decay Length XY/Err",200,0,6.);
520 
521  namedistr="hdeclxyd0d0B_";
522  namedistr+=i;
523  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);
524 
525  namedistr="hnormdeclxyd0d0B_";
526  namedistr+=i;
527  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);
528 
529  // costhetapoint
530  namedistr="hcosthetapointB_";
531  namedistr+=i;
532  TH1F *hcosthetapointB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
533 
534  namedistr="hcosthetapointxyB_";
535  namedistr+=i;
536  TH1F *hcosthetapointxyB = new TH1F(namedistr.Data(), "cos#theta_{Point} XY distribution;cos#theta_{Point} XY",300,0.,1.);
537 
538  TH1F* tmpMB = new TH1F(nameMassNocutsB.Data(),"D^{0} invariant mass; M [GeV]; Entries",600,1.6248,2.2248); //range (MD0-300MeV, mD0 + 300MeV)
539  tmpMB->Sumw2();
540 
541 
542  fDistr->Add(hNclsD0vsptB);
543  fDistr->Add(hNclsD0barvsptB);
544  fDistr->Add(hNITSpointsD0vsptB);
545  fDistr->Add(hNSPDpointsD0B);
546  fDistr->Add(hptD0B);
547  fDistr->Add(hphiD0B);
548  fDistr->Add(hptD0barB);
549  fDistr->Add(hphiD0barB);
550  fDistr->Add(hetaphiD0candidateB);
551  fDistr->Add(hetaphiD0candidatesignalregionB);
552  fDistr->Add(hetaphiD0barcandidateB);
553  fDistr->Add(hetaphiD0barcandidatesignalregionB);
554 
555  fDistr->Add(hdcaB);
556 
557  fDistr->Add(hd0B);
558 
559  fDistr->Add(hd0d0B);
560 
561  fDistr->Add(hcosthetapointB);
562 
563  fDistr->Add(hcosthetapointxyB);
564 
565  fDistr->Add(hdeclengthB);
566 
567  fDistr->Add(hnormdeclengthB);
568 
569  fDistr->Add(hdeclxyB);
570 
571  fDistr->Add(hnormdeclxyB);
572 
573  fDistr->Add(hdeclxyd0d0B);
574  fDistr->Add(hnormdeclxyd0d0B);
575 
576  fDistr->Add(tmpMB);
577 
578  //histograms filled only when the secondary vertex is recalculated w/o the daughter tracks (as requested in the cut object)
579 
581  if(fReadMC){
582  namedistr="hd0vpiS_";
583  namedistr+=i;
584  TH1F *hd0vpiS = new TH1F(namedistr.Data(), "Impact parameter distribution (pions)(vtx w/o these tracks);d0(#pi) [cm]",200,-0.1,0.1);
585 
586  namedistr="hd0vKS_";
587  namedistr+=i;
588  TH1F *hd0vKS = new TH1F(namedistr.Data(), "Impact parameter distribution (kaons) (vtx w/o these tracks);d0(K) [cm]",200,-0.1,0.1);
589 
590  namedistr="hd0d0vS_";
591  namedistr+=i;
592  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);
593 
594  namedistr="hdeclvS_";
595  namedistr+=i;
596  TH1F *hdeclengthvS = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
597 
598  namedistr="hnormdeclvS_";
599  namedistr+=i;
600  TH1F *hnormdeclengthvS = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
601 
602  fDistr->Add(hd0vpiS);
603  fDistr->Add(hd0vKS);
604  fDistr->Add(hd0d0vS);
605  fDistr->Add(hdeclengthvS);
606  fDistr->Add(hnormdeclengthvS);
607 
608  }
609 
610  namedistr="hd0vmoresB_";
611  namedistr+=i;
612  TH1F *hd0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
613 
614  namedistr="hd0d0vmoresB_";
615  namedistr+=i;
616  TH1F *hd0d0vmoresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
617 
618 
619  namedistr="hd0vB_";
620  namedistr+=i;
621  TH1F *hd0vB = new TH1F(namedistr.Data(), "Impact parameter distribution (vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
622 
623  namedistr="hd0vp0B_";
624  namedistr+=i;
625  TH1F *hd0vp0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong + ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
626 
627  namedistr="hd0vp1B_";
628  namedistr+=i;
629  TH1F *hd0vp1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong - ** vtx w/o these tracks);d0 [cm]",200,-0.1,0.1);
630 
631  namedistr="hd0d0vB_";
632  namedistr+=i;
633  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);
634 
635  namedistr="hdeclvB_";
636  namedistr+=i;
637  TH1F *hdeclengthvB = new TH1F(namedistr.Data(), "Decay Length distribution (vtx w/o tracks);Decay Length [cm]",200,0,0.6);
638 
639  namedistr="hnormdeclvB_";
640  namedistr+=i;
641  TH1F *hnormdeclengthvB = new TH1F(namedistr.Data(), "Normalized Decay Length distribution (vtx w/o tracks);Decay Length/Err ",200,0,10.);
642 
643  fDistr->Add(hd0vB);
644  fDistr->Add(hd0vp0B);
645  fDistr->Add(hd0vp1B);
646  fDistr->Add(hd0vmoresB);
647 
648  fDistr->Add(hd0d0vB);
649  fDistr->Add(hd0d0vmoresB);
650 
651  fDistr->Add(hdeclengthvB);
652 
653  fDistr->Add(hnormdeclengthvB);
654  }
655 
656 
657  }
658  //histograms of invariant mass distributions
659 
660 
661  //MC signal
662  if(fReadMC){
663  TH1F* tmpSt = new TH1F(nameSgn.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries",600,1.6248,2.2248);
664 
665  TH1F *tmpSl=(TH1F*)tmpSt->Clone();
666  tmpSt->Sumw2();
667  tmpSl->Sumw2();
668 
669  //Reflection: histo filled with D0Mass which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
670  TH1F* tmpRt = new TH1F(nameRfl.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries",600,1.6248,2.2248);
671  //TH1F *tmpRl=(TH1F*)tmpRt->Clone();
672  TH1F* tmpBt = new TH1F(nameBkg.Data(), "Background invariant mass - MC; M [GeV]; Entries",600,1.6248,2.2248);
673  //TH1F *tmpBl=(TH1F*)tmpBt->Clone();
674  tmpBt->Sumw2();
675  //tmpBl->Sumw2();
676  tmpRt->Sumw2();
677  //tmpRl->Sumw2();
678 
679  fOutputMass->Add(tmpSt);
680  fOutputMass->Add(tmpRt);
681  fOutputMass->Add(tmpBt);
682 
683  }
684 
685  //mass
686  TH1F* tmpMt = new TH1F(nameMass.Data(),"D^{0} invariant mass; M [GeV]; Entries",600,1.6248,2.2248);
687  //TH1F *tmpMl=(TH1F*)tmpMt->Clone();
688  tmpMt->Sumw2();
689  //tmpMl->Sumw2();
690  //distribution w/o cuts large range
691  // TH1F* tmpMS = new TH1F(nameMassNocutsS.Data(),"D^{0} invariant mass; M [GeV]; Entries",300,0.7,3.);
692 
693  fOutputMass->Add(tmpMt);
694 
695  if(fSys==0){ //histograms filled only in pp to save time in PbPb
696  if(fFillVarHists){
697 
698  if(fReadMC){
699  // pT
700  namedistr="hptpiS_";
701  namedistr+=i;
702  TH1F *hptpiS = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
703 
704  namedistr="hptKS_";
705  namedistr+=i;
706  TH1F *hptKS = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
707 
708  // costhetastar
709  namedistr="hcosthetastarS_";
710  namedistr+=i;
711  TH1F *hcosthetastarS = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
712 
713  //pT no mass cut
714 
715  namedistr="hptpiSnoMcut_";
716  namedistr+=i;
717  TH1F *hptpiSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (pions);p_{T} [GeV/c]",200,0.,8.);
718 
719  namedistr="hptKSnoMcut_";
720  namedistr+=i;
721  TH1F *hptKSnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution (kaons);p_{T} [GeV/c]",200,0.,8.);
722 
723  fDistr->Add(hptpiS);
724  fDistr->Add(hptKS);
725  fDistr->Add(hcosthetastarS);
726 
727  fDistr->Add(hptpiSnoMcut);
728  fDistr->Add(hptKSnoMcut);
729 
730  // costhetapoint vs d0 or d0d0
731  namedistr="hcosthpointd0S_";
732  namedistr+=i;
733  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);
734  namedistr="hcosthpointd0d0S_";
735  namedistr+=i;
736  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);
737 
738  fDistr->Add(hcosthpointd0S);
739  fDistr->Add(hcosthpointd0d0S);
740 
741  //to compare with AliAnalysisTaskCharmFraction
742  TH1F* tmpS27t = new TH1F(nameSgn27.Data(),"D^{0} invariant mass in M(D^{0}) +/- 27 MeV - MC; M [GeV]; Entries",600,1.6248,2.2248);
743  TH1F *tmpS27l=(TH1F*)tmpS27t->Clone();
744  tmpS27t->Sumw2();
745  tmpS27l->Sumw2();
746 
747  fOutputMass->Add(tmpS27t);
748  fOutputMass->Add(tmpS27l);
749 
750  }
751 
752  // pT
753  namedistr="hptB_";
754  namedistr+=i;
755  TH1F *hptB = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
756 
757  // costhetastar
758  namedistr="hcosthetastarB_";
759  namedistr+=i;
760  TH1F *hcosthetastarB = new TH1F(namedistr.Data(), "cos#theta* distribution;cos#theta*",200,-1.,1.);
761 
762  //pT no mass cut
763  namedistr="hptB1prongnoMcut_";
764  namedistr+=i;
765  TH1F *hptB1pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
766 
767  namedistr="hptB2prongsnoMcut_";
768  namedistr+=i;
769  TH1F *hptB2pnoMcut = new TH1F(namedistr.Data(), "P_{T} distribution;p_{T} [GeV/c]",200,0.,8.);
770 
771  fDistr->Add(hptB);
772  fDistr->Add(hcosthetastarB);
773 
774  fDistr->Add(hptB1pnoMcut);
775  fDistr->Add(hptB2pnoMcut);
776 
777  //impact parameter of negative/positive track
778  namedistr="hd0p0B_";
779  namedistr+=i;
780  TH1F *hd0p0B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.1,0.1);
781 
782  namedistr="hd0p1B_";
783  namedistr+=i;
784  TH1F *hd0p1B = new TH1F(namedistr.Data(), "Impact parameter distribution (prong -);d0 [cm]",200,-0.1,0.1);
785 
786  //impact parameter corrected for strangeness
787  namedistr="hd0moresB_";
788  namedistr+=i;
789  TH1F *hd0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (both);d0 [cm]",200,-0.1,0.1);
790 
791  namedistr="hd0d0moresB_";
792  namedistr+=i;
793  TH1F *hd0d0moresB = new TH1F(namedistr.Data(), "Impact parameter distribution (prong +);d0 [cm]",200,-0.001,0.001);
794 
795 
796  namedistr="hcosthetapointmoresB_";
797  namedistr+=i;
798  TH1F *hcosthetapointmoresB = new TH1F(namedistr.Data(), "cos#theta_{Point} distribution;cos#theta_{Point}",200,0,1.);
799 
800  // costhetapoint vs d0 or d0d0
801  namedistr="hcosthpointd0B_";
802  namedistr+=i;
803  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);
804 
805  namedistr="hcosthpointd0d0B_";
806  namedistr+=i;
807  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);
808 
809  fDistr->Add(hd0p0B);
810  fDistr->Add(hd0p1B);
811 
812  fDistr->Add(hd0moresB);
813  fDistr->Add(hd0d0moresB);
814  fDistr->Add(hcosthetapointmoresB);
815 
816 
817  fDistr->Add(hcosthpointd0B);
818 
819 
820  fDistr->Add(hcosthpointd0d0B);
821  }
822 
823  } //end pp histo
824  }
825  //create THnSparse for impact parameter analysis in DATA sample.
826  //pt, normImpParTrk1, normImpParTrk2, decLXY, normDecLXY, massd0, cut, pid, D0D0bar
827  Int_t nbinsImpParStudy[9]= {50, 40, 40, 20, 15, 600, 3,4,2};
828  Double_t limitLowImpParStudy[9]={0, -5, -5, 0, 0,1.6248,1,0,0};
829  Double_t limitUpImpParStudy[9]= {50., 5, 5, 0.2,15,2.2248,4,4,2};
830  TString axTit[9]={"#it{p}_{T} (GeV/c)","normalized imp par residual, trk1","normalized imp par residual, trk2","#it{L}_{xy} (cm)","norm #it{L}_{xy}","MassD0_{k#pi} (GeV/#it{c}^{2})", "cutSel","PIDinfo","D0D0bar"};
831  fhStudyImpParSingleTrackCand=new THnSparseF("fhStudyImpParSingleTrackCand","fhStudyImpParSingleTrackCand",9,nbinsImpParStudy,limitLowImpParStudy,limitUpImpParStudy);
832  for(Int_t iax=0; iax<9; iax++) fhStudyImpParSingleTrackCand->GetAxis(iax)->SetTitle(axTit[iax].Data());
834  if(fReadMC){
835  //pt, ptB, normImpParTrk1, normImpParTrk2, decLXY, normDecLXY, iscut, ispid
836  Int_t nbinsImpParStudy[8]= {50,50,40, 40, 20, 15, 3, 4};
837  Double_t limitLowImpParStudy[8]={0, 0, -5,-5., 0., 0., 1.,0.};
838  Double_t limitUpImpParStudy[8]= {50.,50., 5, 5, 0.2, 15, 3.,4.};
839 
840  fhStudyImpParSingleTrackSign=new THnSparseF("fhStudyImpParSingleTrackSign","fhStudyImpParSingleTrackSign",8,nbinsImpParStudy,limitLowImpParStudy,limitUpImpParStudy);
841  TString axTitMC[8]={"#it{p}_{T} (GeV/c)","#it{p}_{T} (GeV/c)","normalized imp par residual, trk1","normalized imp par residual, trk2","#it{L}_{xy} (cm)","norm #it{L}_{xy}","cutSel","PIDinfo"};
842  for(Int_t iax=0; iax<8; iax++) fhStudyImpParSingleTrackSign->GetAxis(iax)->SetTitle(axTitMC[iax].Data());
844 
845 
846  fhStudyImpParSingleTrackFd=new THnSparseF("fhStudyImpParSingleTrackFd","fhStudyImpParSingleTrackFd",8,nbinsImpParStudy,limitLowImpParStudy,limitUpImpParStudy);
847  for(Int_t iax=0; iax<8; iax++) fhStudyImpParSingleTrackFd->GetAxis(iax)->SetTitle(axTitMC[iax].Data());
849 
851  }
852 
853 
854  //for Like sign analysis
855 
856  if(fArray==1){
857  namedistr="hpospair";
858  TH1F* hpospair=new TH1F(namedistr.Data(),"Number of positive pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
859  namedistr="hnegpair";
860  TH1F* hnegpair=new TH1F(namedistr.Data(),"Number of negative pairs",fCuts->GetNPtBins(),-0.5,fCuts->GetNPtBins()-0.5);
861  fOutputMass->Add(hpospair);
862  fOutputMass->Add(hnegpair);
863  }
864 
865 
866  // 2D Pt distributions and impact parameter histograms
867  if(fFillPtHist) {
868 
869  nameMassPt="histMassPt";
870  nameSgnPt="histSgnPt";
871  nameBkgPt="histBkgPt";
872  nameRflPt="histRflPt";
873 
874  //MC signal
875  if(fReadMC){
876  TH2F* tmpStPt = new TH2F(nameSgnPt.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.6248,2.2248,nbins2dPt,binInPt,binFinPt);
877  TH2F *tmpSlPt=(TH2F*)tmpStPt->Clone();
878  tmpStPt->Sumw2();
879  tmpSlPt->Sumw2();
880 
881  //Reflection: histo filled with D0MassV1 which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
882  TH2F* tmpRtPt = new TH2F(nameRflPt.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.6248,2.2248,nbins2dPt,binInPt,binFinPt);
883  TH2F* tmpBtPt = new TH2F(nameBkgPt.Data(), "Background invariant mass - MC; M [GeV]; Entries; Pt[GeV/c]",600,1.6248,2.2248,nbins2dPt,binInPt,binFinPt);
884  tmpBtPt->Sumw2();
885  tmpRtPt->Sumw2();
886 
887  fOutputMassPt->Add(tmpStPt);
888  fOutputMassPt->Add(tmpRtPt);
889  fOutputMassPt->Add(tmpBtPt);
890 
891  // cout<<endl<<endl<<"***************************************"<<endl;
892  // cout << " created and added to the list "<< nameSgnPt.Data() <<" "<< tmpStPt <<
893  // ", "<<nameRflPt.Data() <<" " << tmpRtPt<<", "<<nameBkgPt.Data()<< " " << tmpBtPt <<endl;
894  // cout<<"***************************************"<<endl<<endl;
895  }
896 
897  TH2F* tmpMtPt = new TH2F(nameMassPt.Data(),"D^{0} invariant mass; M [GeV]; Entries; Pt[GeV/c]",600,1.6248,2.2248,nbins2dPt,binInPt,binFinPt);
898  tmpMtPt->Sumw2();
899 
900  fOutputMassPt->Add(tmpMtPt);
901  }
902 
904 
905  // 2D Y distributions
906 
907  if(fFillYHist) {
908  for(Int_t i=0;i<fCuts->GetNPtBins();i++){
909  nameMassY="histMassY_";
910  nameMassY+=i;
911  nameSgnY="histSgnY_";
912  nameSgnY+=i;
913  nameBkgY="histBkgY_";
914  nameBkgY+=i;
915  nameRflY="histRflY_";
916  nameRflY+=i;
917  //MC signal
918  if(fReadMC){
919  TH2F* tmpStY = new TH2F(nameSgnY.Data(), "D^{0} invariant mass - MC; M [GeV]; Entries; y",600,1.6248,2.2248,nbins2dY,binInY,binFinY);
920  tmpStY->Sumw2();
921  //Reflection: histo filled with D0MassV1 which pass the cut (also) as D0bar and with D0bar which pass (also) the cut as D0
922  TH2F* tmpRtY = new TH2F(nameRflY.Data(), "Reflected signal invariant mass - MC; M [GeV]; Entries; y",600,1.6248,2.2248,nbins2dY,binInY,binFinY);
923  TH2F* tmpBtY = new TH2F(nameBkgY.Data(), "Background invariant mass - MC; M [GeV]; Entries; y",600,1.6248,2.2248,nbins2dY,binInY,binFinY);
924  tmpBtY->Sumw2();
925  tmpRtY->Sumw2();
926 
927  fOutputMassY->Add(tmpStY);
928  fOutputMassY->Add(tmpRtY);
929  fOutputMassY->Add(tmpBtY);
930  }
931  TH2F* tmpMtY = new TH2F(nameMassY.Data(),"D^{0} invariant mass; M [GeV]; Entries; y",600,1.6248,2.2248,nbins2dY,binInY,binFinY);
932  tmpMtY->Sumw2();
933  fOutputMassY->Add(tmpMtY);
934  }
935  }
936 
937 
938  const char* nameoutput=GetOutputSlot(3)->GetContainer()->GetName();
939 
940  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", 23,-0.5,22.5);
941 
942  fNentries->GetXaxis()->SetBinLabel(1,"nEventsAnal");
943  fNentries->GetXaxis()->SetBinLabel(2,"nCandSel(Cuts)");
944  if(fReadMC) fNentries->GetXaxis()->SetBinLabel(3,"nD0Selected");
945  else fNentries->GetXaxis()->SetBinLabel(3,"Dstar<-D0");
946  fNentries->GetXaxis()->SetBinLabel(4,"nEventsGoodVtxS");
947  fNentries->GetXaxis()->SetBinLabel(5,"ptbin = -1");
948  fNentries->GetXaxis()->SetBinLabel(6,"no daughter");
949  if(fSys==0) fNentries->GetXaxis()->SetBinLabel(7,"nCandSel(Tr)");
950  if(fFillVarHists || fPIDCheck){
951  fNentries->GetXaxis()->SetBinLabel(8,"PID=0");
952  fNentries->GetXaxis()->SetBinLabel(9,"PID=1");
953  fNentries->GetXaxis()->SetBinLabel(10,"PID=2");
954  fNentries->GetXaxis()->SetBinLabel(11,"PID=3");
955  }
956  if(fReadMC && fSys==0){
957  fNentries->GetXaxis()->SetBinLabel(12,"K");
958  fNentries->GetXaxis()->SetBinLabel(13,"Lambda");
959  }
960  fNentries->GetXaxis()->SetBinLabel(14,"Pile-up Rej");
961  fNentries->GetXaxis()->SetBinLabel(15,"N. of 0SMH");
962  if(fSys==1) fNentries->GetXaxis()->SetBinLabel(16,"Nev in centr");
963  if(fIsRejectSDDClusters) fNentries->GetXaxis()->SetBinLabel(17,"SDD-Cls Rej");
964  fNentries->GetXaxis()->SetBinLabel(18,"Phys.Sel.Rej");
965  fNentries->GetXaxis()->SetBinLabel(19,"D0 failed to be filled");
966  fNentries->GetXaxis()->SetBinLabel(20,"fisFilled is 0");
967  fNentries->GetXaxis()->SetBinLabel(21,"fisFilled is 1");
968  fNentries->GetXaxis()->SetBinLabel(22,"AOD/dAOD mismatch");
969  fNentries->GetXaxis()->SetBinLabel(23,"AOD/dAOD #events ok");
970  fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
971 
972  fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
973  fCounter->Init();
974 
975  //
976  // Output slot 7 : tree of the candidate variables after track selection
977  //
978  nameoutput = GetOutputSlot(7)->GetContainer()->GetName();
979  fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
980  Int_t nVar = 15;
982  TString * fCandidateVariableNames = new TString[nVar];
983  fCandidateVariableNames[0] = "massD0";
984  fCandidateVariableNames[1] = "massD0bar";
985  fCandidateVariableNames[2] = "pt";
986  fCandidateVariableNames[3] = "dca";
987  fCandidateVariableNames[4] = "costhsD0";
988  fCandidateVariableNames[5] = "costhsD0bar";
989  fCandidateVariableNames[6] = "ptk";
990  fCandidateVariableNames[7] = "ptpi";
991  fCandidateVariableNames[8] = "d0k";
992  fCandidateVariableNames[9] = "d0pi";
993  fCandidateVariableNames[10] = "d0xd0";
994  fCandidateVariableNames[11] = "costhp";
995  fCandidateVariableNames[12] = "costhpxy";
996  fCandidateVariableNames[13] = "lxy";
997  fCandidateVariableNames[14] = "specialcuts";
998  for(Int_t ivar=0; ivar<nVar; ivar++){
999  fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/d",fCandidateVariableNames[ivar].Data()));
1000  }
1001 
1002 
1003  //
1004  // Output slot 8 : List for detector response histograms
1005  //
1006  if (fDrawDetSignal) {
1007  TH2F *TOFSigBefPID = new TH2F("TOFSigBefPID", "TOF signal of daughters before PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
1008  TH2F *TOFSigAftPID = new TH2F("TOFSigAftPID", "TOF signal after PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
1009 
1010  TH2F *TPCSigBefPID = new TH2F("TPCSigBefPID", "TPC dE/dx before PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
1011  TH2F *TPCSigAftPID = new TH2F("TPCSigAftPID", "TPC dE/dx after PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
1012 
1013  fDetSignal->Add(TOFSigBefPID);
1014  fDetSignal->Add(TOFSigAftPID);
1015  fDetSignal->Add(TPCSigBefPID);
1016  fDetSignal->Add(TPCSigAftPID);
1017  }
1018 
1019  // Post the data
1020  PostData(1,fOutputMass);
1021  PostData(2,fDistr);
1022  PostData(3,fNentries);
1023  PostData(5,fCounter);
1024  PostData(6,fOutputMassPt);
1025  PostData(7,fVariablesTree);
1026  PostData(8, fDetSignal);
1027  PostData(9,fOutputMassY);
1028 
1029  return;
1030 }
1031 
1032 //________________________________________________________________________
1034 
1035 {
1038 
1039  //cuts order
1040  // printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
1041  // printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
1042  // printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
1043  // printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
1044  // printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
1045  // printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
1046  // printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
1047  // printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
1048  // printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
1049 
1050  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
1051 
1052  if(fAODProtection>=0){
1053  // Protection against different number of events in the AOD and deltaAOD
1054  // In case of discrepancy the event is rejected.
1055  Int_t matchingAODdeltaAODlevel = AliRDHFCuts::CheckMatchingAODdeltaAODevents();
1056  if (matchingAODdeltaAODlevel<0 || (matchingAODdeltaAODlevel==0 && fAODProtection==1)) {
1057  // AOD/deltaAOD trees have different number of entries || TProcessID do not match while it was required
1058  fNentries->Fill(21);
1059  return;
1060  }
1061  fNentries->Fill(22);
1062  }
1063 
1064  TString bname;
1065  if(fArray==0){ //D0 candidates
1066  // load D0->Kpi candidates
1067  //cout<<"D0 candidates"<<endl;
1068  bname="D0toKpi";
1069  } else { //LikeSign candidates
1070  // load like sign candidates
1071  //cout<<"LS candidates"<<endl;
1072  bname="LikeSign2Prong";
1073  }
1074  TClonesArray *inputArray=0;
1075  if(!aod && AODEvent() && IsStandardAOD()) {
1076  // In case there is an AOD handler writing a standard AOD, use the AOD
1077  // event in memory rather than the input (ESD) event.
1078  aod = dynamic_cast<AliAODEvent*> (AODEvent());
1079  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
1080  // have to taken from the AOD event hold by the AliAODExtension
1081  AliAODHandler* aodHandler = (AliAODHandler*)
1082  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
1083 
1084  if(aodHandler->GetExtensions()) {
1085  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
1086  AliAODEvent* aodFromExt = ext->GetAOD();
1087  inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
1088  }
1089  } else if(aod) {
1090  inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
1091  }
1092 
1093  if(!inputArray || !aod) {
1094  printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
1095  return;
1096  }
1097  // fix for temporary bug in ESDfilter
1098  // the AODs with null vertex pointer didn't pass the PhysSel
1099  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
1100 
1101  TClonesArray *mcArray = 0;
1102  AliAODMCHeader *mcHeader = 0;
1103 
1104  if(fReadMC) {
1105  // load MC particles
1106  mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1107  if(!mcArray) {
1108  printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
1109  return;
1110  }
1111 
1112  // load MC header
1113  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1114  if(!mcHeader) {
1115  printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
1116  return;
1117  }
1118  }
1119  //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
1120 
1121  //histogram filled with 1 for every AOD
1122  fNentries->Fill(0);
1123 
1124 
1126  //fCounter->StoreEvent(aod,fReadMC);
1127  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
1128  TString trigclass=aod->GetFiredTriggerClasses();
1129  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
1130  if(fReadMC && fStepMCAcc){
1131  FillMCAcceptanceHistos(mcArray, mcHeader);
1132  }
1133 
1134  if(!fCuts->IsEventSelected(aod)) {
1135  if(fCuts->GetWhyRejection()==1) // rejected for pileup
1136  fNentries->Fill(13);
1137  if(fSys==1 && (fCuts->GetWhyRejection()==2 || fCuts->GetWhyRejection()==3)) fNentries->Fill(15);
1138  if(fCuts->GetWhyRejection()==7) fNentries->Fill(17);
1139  return;
1140  }
1141  // Check the Nb of SDD clusters
1142  if (fIsRejectSDDClusters) {
1143  Bool_t skipEvent = kFALSE;
1144  Int_t ntracks = 0;
1145  if (aod) ntracks = aod->GetNumberOfTracks();
1146  for(Int_t itrack=0; itrack<ntracks; itrack++) { // loop on tacks
1147  // ... get the track
1148  AliAODTrack * track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrack));
1149  if(!track) AliFatal("Not a standard AOD");
1150  if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
1151  skipEvent=kTRUE;
1152  fNentries->Fill(16);
1153  break;
1154  }
1155  }
1156  if (skipEvent) return;
1157  }
1158 
1159  // AOD primary vertex
1160  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
1161 
1162  Bool_t isGoodVtx=kFALSE;
1163 
1164  //vtx1->Print();
1165  TString primTitle = vtx1->GetTitle();
1166  if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
1167  isGoodVtx=kTRUE;
1168  fNentries->Fill(3);
1169  }
1170 
1171  // loop over candidates
1172  Int_t nInD0toKpi = inputArray->GetEntriesFast();
1173  if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
1174 
1175  // FILE *f=fopen("4display.txt","a");
1176  // printf("Number of D0->Kpi: %d\n",nInD0toKpi);
1177  Int_t nSelectedloose=0,nSelectedtight=0;
1178 
1179  // vHF object is needed to call the method that refills the missing info of the candidates
1180  // if they have been deleted in dAOD reconstruction phase
1181  // in order to reduce the size of the file
1183 
1184  for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
1185  AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
1186 
1188  fNentries->Fill(2);
1189  continue; //skip the D0 from Dstar
1190  }
1191  if(d->GetIsFilled()==0)fNentries->Fill(19);//tmp check
1192  if(d->GetIsFilled()==1)fNentries->Fill(20);//tmp check
1193  if(!(vHF->FillRecoCand(aod,d))) {//Fill the data members of the candidate only if they are empty.
1194  fNentries->Fill(18); //monitor how often this fails
1195  continue;
1196  }
1197 
1198  if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
1199  nSelectedloose++;
1200  nSelectedtight++;
1201  if(fSys==0){
1202  if(fCuts->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
1203  }
1204  Int_t ptbin=fCuts->PtBin(d->Pt());
1205  if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
1207  if(fFillVarHists) {
1208  //if(!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate)) {
1209  fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(0),0);
1210  fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(1),1);
1211  //check daughters
1212  if(!fDaughterTracks.UncheckedAt(0) || !fDaughterTracks.UncheckedAt(1)) {
1213  AliDebug(1,"at least one daughter not found!");
1214  fNentries->Fill(5);
1215  fDaughterTracks.Clear();
1216  continue;
1217  }
1218  //}
1219  FillVarHists(aod,d,mcArray,fCuts,fDistr);
1220  }
1221 
1222  if (fDrawDetSignal) {
1224  }
1225  FillMassHists(d,mcArray,mcHeader,fCuts,fOutputMass);
1226  NormIPvar(aod, d,mcArray);
1227  if (fPIDCheck) {
1228  Int_t isSelectedPIDfill = 3;
1229  if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(d); //0 rejected,1 D0,2 Dobar, 3 both
1230 
1231  if (isSelectedPIDfill == 0)fNentries->Fill(7);
1232  if (isSelectedPIDfill == 1)fNentries->Fill(8);
1233  if (isSelectedPIDfill == 2)fNentries->Fill(9);
1234  if (isSelectedPIDfill == 3)fNentries->Fill(10);
1235  }
1236 
1237  }
1238 
1239  fDaughterTracks.Clear();
1240  //if(unsetvtx) d->UnsetOwnPrimaryVtx();
1241  } //end for prongs
1242  fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1243  fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1244  delete vHF;
1245  // Post the data
1246  PostData(1,fOutputMass);
1247  PostData(2,fDistr);
1248  PostData(3,fNentries);
1249  PostData(5,fCounter);
1250  PostData(6,fOutputMassPt);
1251  PostData(7,fVariablesTree);
1252  PostData(8, fDetSignal);
1253 
1254  return;
1255 }
1256 
1257 //____________________________________________________________________________
1259 {
1260  //
1262  //
1263  fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(0), 0);
1264  fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(1), 1);
1265 
1266  AliESDtrack *esdtrack1 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1267  AliESDtrack *esdtrack2 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
1268 
1269 
1270  //For filling detector histograms
1271  Int_t isSelectedPIDfill = 3;
1272  if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
1273 
1274 
1275  //fill "before PID" histos with every daughter
1276  ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
1277  ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
1278  ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
1279  ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
1280 
1281  if (isSelectedPIDfill != 0) { //fill "After PID" with everything that's not rejected
1282  ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
1283  ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
1284  ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
1285  ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
1286 
1287  }
1288 
1289  delete esdtrack1;
1290  delete esdtrack2;
1291 
1292  return;
1293 }
1294 
1295 //____________________________________________________________________________
1297  //
1299  //
1300 
1301 
1302  Int_t pdgDgD0toKpi[2]={321,211};
1303  Int_t lab=-9999;
1304  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)
1305  //Double_t pt = d->Pt(); //mother pt
1306  Int_t isSelectedPID=3;
1307  if(!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPID=cuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
1308  if (!fPIDCheck) { //if fPIDCheck, this is already filled elsewhere
1309  if (isSelectedPID==0)fNentries->Fill(7);
1310  if (isSelectedPID==1)fNentries->Fill(8);
1311  if (isSelectedPID==2)fNentries->Fill(9);
1312  if (isSelectedPID==3)fNentries->Fill(10);
1313  //fNentries->Fill(8+isSelectedPID);
1314  }
1315 
1316  if(fCutOnDistr && !fIsSelectedCandidate) return;
1317  //printf("\nif no cuts or cuts passed\n");
1318 
1319 
1320  //add distr here
1321  UInt_t pdgs[2];
1322 
1323  Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1324  pdgs[0]=211;
1325  pdgs[1]=321;
1326  Double_t minvD0 = part->InvMassD0();
1327  pdgs[1]=211;
1328  pdgs[0]=321;
1329  Double_t minvD0bar = part->InvMassD0bar();
1330  //cout<<"inside mass cut"<<endl;
1331 
1332  Double_t invmasscut=0.03;
1333 
1334  TString fillthispi="",fillthisK="",fillthis="", fillthispt="", fillthisetaphi="";
1335 
1336  Int_t ptbin=cuts->PtBin(part->Pt());
1337  Double_t pt = part->Pt();
1338 
1339  Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
1340  dz1[0]=-99; dz2[0]=-99;
1341  Double_t d0[2];
1342  Double_t decl[2]={-99,-99};
1343  Bool_t recalcvtx=kFALSE;
1344 
1345 
1346 
1348  recalcvtx=kTRUE;
1349  //recalculate vertex
1350  AliAODVertex *vtxProp=0x0;
1351  vtxProp=GetPrimaryVtxSkipped(aod);
1352  if(vtxProp) {
1353  part->SetOwnPrimaryVtx(vtxProp);
1354  //Bool_t unsetvtx=kTRUE;
1355  //Calculate d0 for daughter tracks with Vtx Skipped
1356  AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1357  AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
1358  esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
1359  esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
1360  delete vtxProp; vtxProp=NULL;
1361  delete esdtrack1;
1362  delete esdtrack2;
1363  }
1364 
1365  d0[0]=dz1[0];
1366  d0[1]=dz2[0];
1367 
1368  decl[0]=part->DecayLength2();
1369  decl[1]=part->NormalizedDecayLength2();
1370  part->UnsetOwnPrimaryVtx();
1371 
1372  }
1373 
1374  Double_t cosThetaStarD0 = 99;
1375  Double_t cosThetaStarD0bar = 99;
1376  Double_t cosPointingAngle = 99;
1377  Double_t normalizedDecayLength2 = -1, normalizedDecayLengthxy=-1;
1378  Double_t decayLength2 = -1, decayLengthxy=-1;
1379  Double_t ptProng[2]={-99,-99};
1380  Double_t d0Prong[2]={-99,-99};
1381  Double_t etaD = 99.;
1382  Double_t phiD = 99.;
1383 
1384 
1385  //disable the PID
1386  if(!fUsePid4Distr) isSelectedPID=0;
1387  if((lab>=0 && fReadMC) || (!fReadMC && isSelectedPID)){ //signal (from MC or PID)
1388 
1389  //check pdg of the prongs
1390  AliAODTrack *prong0=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1391  AliAODTrack *prong1=(AliAODTrack*)fDaughterTracks.UncheckedAt(1);
1392 
1393  if(!prong0 || !prong1) {
1394  return;
1395  }
1396  Int_t labprong[2];
1397  if(fReadMC){
1398  labprong[0]=prong0->GetLabel();
1399  labprong[1]=prong1->GetLabel();
1400  }
1401  AliAODMCParticle *mcprong=0;
1402  Int_t pdgProng[2]={0,0};
1403 
1404  for (Int_t iprong=0;iprong<2;iprong++){
1405  if(fReadMC && labprong[iprong]>=0) {
1406  mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1407  pdgProng[iprong]=mcprong->GetPdgCode();
1408  }
1409  }
1410 
1411  if(fSys==0){
1412  //no mass cut ditributions: ptbis
1413 
1414  fillthispi="hptpiSnoMcut_";
1415  fillthispi+=ptbin;
1416 
1417  fillthisK="hptKSnoMcut_";
1418  fillthisK+=ptbin;
1419 
1420  if ((TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321)
1421  || (isSelectedPID==1 || isSelectedPID==3)){
1422  ((TH1F*)listout->FindObject(fillthispi))->Fill(prong0->Pt());
1423  ((TH1F*)listout->FindObject(fillthisK))->Fill(prong1->Pt());
1424  }
1425 
1426  if ((TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211)
1427  || (isSelectedPID==2 || isSelectedPID==3)){
1428  ((TH1F*)listout->FindObject(fillthisK))->Fill(prong0->Pt());
1429  ((TH1F*)listout->FindObject(fillthispi))->Fill(prong1->Pt());
1430  }
1431  }
1432 
1433  //no mass cut ditributions: mass
1434 
1435  etaD = part->Eta();
1436  phiD = part->Phi();
1437 
1438 
1439  fillthis="hMassS_";
1440  fillthis+=ptbin;
1441  fillthispt="histSgnPt";
1442 
1443  if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)
1444  || (!fReadMC && (isSelectedPID==1 || isSelectedPID==3))){//D0
1445  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1446  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
1447 
1448  fillthisetaphi="hetaphiD0candidateS_";
1449  fillthisetaphi+=ptbin;
1450  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1451 
1452  if(TMath::Abs(minvD0-mPDG)<0.05){
1453  fillthisetaphi="hetaphiD0candidatesignalregionS_";
1454  fillthisetaphi+=ptbin;
1455  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1456  }
1457 
1458  }
1459  else { //D0bar
1460  if(fReadMC || (!fReadMC && isSelectedPID > 1)){
1461  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1462  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
1463 
1464  fillthisetaphi="hetaphiD0barcandidateS_";
1465  fillthisetaphi+=ptbin;
1466  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1467 
1468  if(TMath::Abs(minvD0bar-mPDG)<0.05){
1469  fillthisetaphi="hetaphiD0barcandidatesignalregionS_";
1470  fillthisetaphi+=ptbin;
1471  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1472  }
1473 
1474  }
1475  }
1476 
1477  //apply cut on invariant mass on the pair
1478  if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1479 
1480  cosThetaStarD0 = part->CosThetaStarD0();
1481  cosThetaStarD0bar = part->CosThetaStarD0bar();
1482  cosPointingAngle = part->CosPointingAngle();
1483  normalizedDecayLength2 = part->NormalizedDecayLength2();
1484  decayLength2 = part->DecayLength2();
1485  decayLengthxy = part->DecayLengthXY();
1486  normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1487 
1488  ptProng[0]=prong0->Pt(); ptProng[1]=prong1->Pt();
1489  d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
1490 
1491  if(fArray==1) cout<<"LS signal: ERROR"<<endl;
1492  for (Int_t iprong=0; iprong<2; iprong++){
1493  AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
1494  if (fReadMC) labprong[iprong]=prong->GetLabel();
1495 
1496  //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
1497  Int_t pdgprong=0;
1498  if(fReadMC && labprong[iprong]>=0) {
1499  mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1500  pdgprong=mcprong->GetPdgCode();
1501  }
1502 
1503  Bool_t isPionHere[2]={(isSelectedPID==1 || isSelectedPID==3) ? kTRUE : kFALSE,(isSelectedPID==2 || isSelectedPID==3) ? kTRUE : kFALSE};
1504 
1505  if(TMath::Abs(pdgprong)==211 || isPionHere[iprong]) {
1506  //cout<<"pi"<<endl;
1507 
1508  if(fSys==0){
1509  fillthispi="hptpiS_";
1510  fillthispi+=ptbin;
1511  ((TH1F*)listout->FindObject(fillthispi))->Fill(ptProng[iprong]);
1512  }
1513 
1514  fillthispi="hd0piS_";
1515  fillthispi+=ptbin;
1516  ((TH1F*)listout->FindObject(fillthispi))->Fill(d0Prong[iprong]);
1517  if(recalcvtx) {
1518 
1519  fillthispi="hd0vpiS_";
1520  fillthispi+=ptbin;
1521  ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
1522  }
1523 
1524  }
1525 
1526  if(TMath::Abs(pdgprong)==321 || !isPionHere[iprong]) {
1527  //cout<<"kappa"<<endl;
1528  if(fSys==0){
1529  fillthisK="hptKS_";
1530  fillthisK+=ptbin;
1531  ((TH1F*)listout->FindObject(fillthisK))->Fill(ptProng[iprong]);
1532  }
1533 
1534 
1535  fillthisK="hd0KS_";
1536  fillthisK+=ptbin;
1537  ((TH1F*)listout->FindObject(fillthisK))->Fill(d0Prong[iprong]);
1538  if (recalcvtx){
1539  fillthisK="hd0vKS_";
1540  fillthisK+=ptbin;
1541  ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
1542  }
1543  }
1544 
1545  if(fSys==0){
1546  fillthis="hcosthpointd0S_";
1547  fillthis+=ptbin;
1548  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[iprong]);
1549  }
1550  } //end loop on prongs
1551 
1552  fillthis="hdcaS_";
1553  fillthis+=ptbin;
1554  ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1555 
1556  fillthis="hcosthetapointS_";
1557  fillthis+=ptbin;
1558  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
1559 
1560  fillthis="hcosthetapointxyS_";
1561  fillthis+=ptbin;
1562  ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
1563 
1564 
1565  fillthis="hd0d0S_";
1566  fillthis+=ptbin;
1567  ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
1568 
1569  fillthis="hdeclS_";
1570  fillthis+=ptbin;
1571  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
1572 
1573  fillthis="hnormdeclS_";
1574  fillthis+=ptbin;
1575  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
1576 
1577  fillthis="hdeclxyS_";
1578  fillthis+=ptbin;
1579  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
1580 
1581  fillthis="hnormdeclxyS_";
1582  fillthis+=ptbin;
1583  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
1584 
1585  fillthis="hdeclxyd0d0S_";
1586  fillthis+=ptbin;
1587  ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
1588 
1589  fillthis="hnormdeclxyd0d0S_";
1590  fillthis+=ptbin;
1591  ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
1592 
1593  if(recalcvtx) {
1594  fillthis="hdeclvS_";
1595  fillthis+=ptbin;
1596  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1597 
1598  fillthis="hnormdeclvS_";
1599  fillthis+=ptbin;
1600  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1601 
1602  fillthis="hd0d0vS_";
1603  fillthis+=ptbin;
1604  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1605  }
1606 
1607  if(fSys==0){
1608  fillthis="hcosthetastarS_";
1609  fillthis+=ptbin;
1610  if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)) ((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1611  else {
1612  if (fReadMC || isSelectedPID>1)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
1613  if(isSelectedPID==1 || isSelectedPID==3)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1614  }
1615 
1616  fillthis="hcosthpointd0d0S_";
1617  fillthis+=ptbin;
1618  ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1619  }
1620 
1621  if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)){
1622  for(Int_t it=0; it<2; it++){
1623  fillthis="hptD0S_";
1624  fillthis+=ptbin;
1625  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1626  fillthis="hphiD0S_";
1627  fillthis+=ptbin;
1628  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1629  Int_t nPointsITS = 0;
1630  for (Int_t il=0; il<6; il++){
1631  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
1632  }
1633  fillthis="hNITSpointsD0vsptS_";
1634  fillthis+=ptbin;
1635  ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(),nPointsITS);
1636  fillthis="hNSPDpointsD0S_";
1637  fillthis+=ptbin;
1638  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
1639  ((TH1I*)listout->FindObject(fillthis))->Fill(0);
1640  }
1641  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
1642  ((TH1I*)listout->FindObject(fillthis))->Fill(1);
1643  }
1644  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
1645  ((TH1I*)listout->FindObject(fillthis))->Fill(2);
1646  }
1647  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
1648  ((TH1I*)listout->FindObject(fillthis))->Fill(3);
1649  }
1650  fillthis="hNclsD0vsptS_";
1651  fillthis+=ptbin;
1652  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1653  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->GetTPCNcls();
1654  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1655  }
1656  }
1657  else {
1658  if (fReadMC || isSelectedPID>1){
1659  for(Int_t it=0; it<2; it++){
1660  fillthis="hptD0barS_";
1661  fillthis+=ptbin;
1662  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1663  fillthis="hphiD0barS_";
1664  fillthis+=ptbin;
1665  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1666  fillthis="hNclsD0barvsptS_";
1667  fillthis+=ptbin;
1668  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1669  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
1670  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1671  }
1672  }
1673  if(isSelectedPID==1 || isSelectedPID==3){
1674  for(Int_t it=0; it<2; it++){
1675  fillthis="hptD0S_";
1676  fillthis+=ptbin;
1677  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1678  fillthis="hphiD0S_";
1679  fillthis+=ptbin;
1680  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1681  Int_t nPointsITS = 0;
1682  for (Int_t il=0; il<6; il++){
1683  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
1684  }
1685  fillthis="hNITSpointsD0vsptS_";
1686  fillthis+=ptbin;
1687  ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(), nPointsITS);
1688  fillthis="hNSPDpointsD0S_";
1689  fillthis+=ptbin;
1690  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
1691  ((TH1I*)listout->FindObject(fillthis))->Fill(0);
1692  }
1693  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
1694  ((TH1I*)listout->FindObject(fillthis))->Fill(1);
1695  }
1696  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
1697  ((TH1I*)listout->FindObject(fillthis))->Fill(2);
1698  }
1699  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
1700  ((TH1I*)listout->FindObject(fillthis))->Fill(3);
1701  }
1702  fillthis="hNclsD0vsptS_";
1703  fillthis+=ptbin;
1704  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1705  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->GetTPCNcls();
1706  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1707  }
1708  }
1709  }
1710 
1711 
1712  } //end mass cut
1713 
1714  } else{ //Background or LS
1715  //if(!fReadMC){
1716  //cout<<"is background"<<endl;
1717 
1718  etaD = part->Eta();
1719  phiD = part->Phi();
1720 
1721  //no mass cut distributions: mass, ptbis
1722  fillthis="hMassB_";
1723  fillthis+=ptbin;
1724  fillthispt="histBkgPt";
1725 
1727  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1728  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
1729 
1730  fillthisetaphi="hetaphiD0candidateB_";
1731  fillthisetaphi+=ptbin;
1732  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1733 
1734  if(TMath::Abs(minvD0-mPDG)<0.05){
1735  fillthisetaphi="hetaphiD0candidatesignalregionB_";
1736  fillthisetaphi+=ptbin;
1737  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1738  }
1739  }
1740  if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) {
1741  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1742  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
1743 
1744  fillthisetaphi="hetaphiD0barcandidateB_";
1745  fillthisetaphi+=ptbin;
1746  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1747 
1748  if(TMath::Abs(minvD0bar-mPDG)<0.05){
1749  fillthisetaphi="hetaphiD0barcandidatesignalregionB_";
1750  fillthisetaphi+=ptbin;
1751  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1752  }
1753 
1754  }
1755  if(fSys==0){
1756  fillthis="hptB1prongnoMcut_";
1757  fillthis+=ptbin;
1758 
1759  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1760 
1761  fillthis="hptB2prongsnoMcut_";
1762  fillthis+=ptbin;
1763  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1764  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1765  }
1766 
1767 
1768  //apply cut on invariant mass on the pair
1769  if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1770  if(fSys==0){
1771  ptProng[0]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt(); ptProng[1]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt();
1772  cosThetaStarD0 = part->CosThetaStarD0();
1773  cosThetaStarD0bar = part->CosThetaStarD0bar();
1774  }
1775 
1776  cosPointingAngle = part->CosPointingAngle();
1777  normalizedDecayLength2 = part->NormalizedDecayLength2();
1778  decayLength2 = part->DecayLength2();
1779  decayLengthxy = part->DecayLengthXY();
1780  normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1781  d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
1782 
1783 
1784  AliAODTrack *prongg=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1785  if(!prongg) {
1786  if(fDebug>2) cout<<"No daughter found";
1787  return;
1788  }
1789  else{
1790  if(fArray==1){
1791  if(prongg->Charge()==1) {
1792  //fTotPosPairs[ptbin]++;
1793  ((TH1F*)fOutputMass->FindObject("hpospair"))->Fill(ptbin);
1794  } else {
1795  //fTotNegPairs[ptbin]++;
1796  ((TH1F*)fOutputMass->FindObject("hnegpair"))->Fill(ptbin);
1797  }
1798  }
1799  }
1800 
1801  //fill pt and phi distrib for prongs with M cut
1802 
1804  for(Int_t it=0; it<2; it++){
1805  fillthis="hptD0B_";
1806  fillthis+=ptbin;
1807  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1808  fillthis="hphiD0B_";
1809  fillthis+=ptbin;
1810  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1811 
1812  Int_t nPointsITS = 0;
1813  for (Int_t il=0; il<6; il++){
1814  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
1815  }
1816  fillthis="hNITSpointsD0vsptB_";
1817  fillthis+=ptbin;
1818  ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(), nPointsITS);
1819  fillthis="hNSPDpointsD0B_";
1820  fillthis+=ptbin;
1821  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
1822  ((TH1I*)listout->FindObject(fillthis))->Fill(0);
1823  }
1824  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
1825  ((TH1I*)listout->FindObject(fillthis))->Fill(1);
1826  }
1827  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
1828  ((TH1I*)listout->FindObject(fillthis))->Fill(2);
1829  }
1830  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
1831  ((TH1I*)listout->FindObject(fillthis))->Fill(3);
1832  }
1833  fillthis="hNclsD0vsptB_";
1834  fillthis+=ptbin;
1835  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1836  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
1837  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1838  }
1839 
1840 
1841  }
1842 
1843  if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) {
1844  for(Int_t it=0; it<2; it++){
1845  fillthis="hptD0barB_";
1846  fillthis+=ptbin;
1847  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1848  fillthis="hphiD0barB_";
1849  fillthis+=ptbin;
1850  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1851  fillthis="hNclsD0barvsptB_";
1852  fillthis+=ptbin;
1853  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1854  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
1855  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1856  }
1857  }
1858 
1859  fillthis="hd0B_";
1860  fillthis+=ptbin;
1861  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1862  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
1863 
1864  if(fReadMC){
1865  Int_t pdgMother[2]={0,0};
1866  Double_t factor[2]={1,1};
1867 
1868  for(Int_t iprong=0;iprong<2;iprong++){
1869  AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
1870  lab=prong->GetLabel();
1871  if(lab>=0){
1872  AliAODMCParticle* mcprong=(AliAODMCParticle*)arrMC->At(lab);
1873  if(mcprong){
1874  Int_t labmom=mcprong->GetMother();
1875  if(labmom>=0){
1876  AliAODMCParticle* mcmother=(AliAODMCParticle*)arrMC->At(labmom);
1877  if(mcmother) pdgMother[iprong]=mcmother->GetPdgCode();
1878  }
1879  }
1880  }
1881 
1882  if(fSys==0){
1883 
1884  fillthis="hd0moresB_";
1885  fillthis+=ptbin;
1886 
1887  if(TMath::Abs(pdgMother[iprong])==310 || TMath::Abs(pdgMother[iprong])==130 || TMath::Abs(pdgMother[iprong])==321){ //K^0_S, K^0_L, K^+-
1888  if(ptProng[iprong]<=1)factor[iprong]=1./.7;
1889  else factor[iprong]=1./.6;
1890  fNentries->Fill(11);
1891  }
1892 
1893  if(TMath::Abs(pdgMother[iprong])==3122) { //Lambda
1894  factor[iprong]=1./0.25;
1895  fNentries->Fill(12);
1896  }
1897  fillthis="hd0moresB_";
1898  fillthis+=ptbin;
1899 
1900  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[iprong],factor[iprong]);
1901 
1902  if(recalcvtx){
1903  fillthis="hd0vmoresB_";
1904  fillthis+=ptbin;
1905  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[iprong],factor[iprong]);
1906  }
1907  }
1908  } //loop on prongs
1909 
1910  if(fSys==0){
1911  fillthis="hd0d0moresB_";
1912  fillthis+=ptbin;
1913  ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0(),factor[0]*factor[1]);
1914 
1915  fillthis="hcosthetapointmoresB_";
1916  fillthis+=ptbin;
1917  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,factor[0]*factor[1]);
1918 
1919  if(recalcvtx){
1920  fillthis="hd0d0vmoresB_";
1921  fillthis+=ptbin;
1922  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1],factor[0]*factor[1]);
1923  }
1924  }
1925  } //readMC
1926 
1927  if(fSys==0){
1928  //normalise pt distr to half afterwards
1929  fillthis="hptB_";
1930  fillthis+=ptbin;
1931  ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[0]);
1932  ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[1]);
1933 
1934  fillthis="hcosthetastarB_";
1935  fillthis+=ptbin;
1936  if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3)))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1937  if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
1938 
1939 
1940  fillthis="hd0p0B_";
1941  fillthis+=ptbin;
1942  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1943  fillthis="hd0p1B_";
1944  fillthis+=ptbin;
1945  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
1946 
1947  fillthis="hcosthpointd0d0B_";
1948  fillthis+=ptbin;
1949  ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1950 
1951  fillthis="hcosthpointd0B_";
1952  fillthis+=ptbin;
1953  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[0]);
1954  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[1]);
1955 
1956 
1957  if(recalcvtx){
1958 
1959  fillthis="hd0vp0B_";
1960  fillthis+=ptbin;
1961  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1962  fillthis="hd0vp1B_";
1963  fillthis+=ptbin;
1964  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1965 
1966  fillthis="hd0vB_";
1967  fillthis+=ptbin;
1968  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1969  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1970 
1971  }
1972 
1973  }
1974 
1975  fillthis="hdcaB_";
1976  fillthis+=ptbin;
1977  ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1978 
1979  fillthis="hd0d0B_";
1980  fillthis+=ptbin;
1981  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]*d0Prong[1]);
1982 
1983  if(recalcvtx){
1984  fillthis="hd0d0vB_";
1985  fillthis+=ptbin;
1986  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1987  }
1988 
1989  fillthis="hcosthetapointB_";
1990  fillthis+=ptbin;
1991  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
1992 
1993  fillthis="hcosthetapointxyB_";
1994  fillthis+=ptbin;
1995  ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
1996 
1997  fillthis="hdeclB_";
1998  fillthis+=ptbin;
1999  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
2000 
2001  fillthis="hnormdeclB_";
2002  fillthis+=ptbin;
2003  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
2004 
2005  fillthis="hdeclxyB_";
2006  fillthis+=ptbin;
2007  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
2008 
2009  fillthis="hnormdeclxyB_";
2010  fillthis+=ptbin;
2011  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
2012 
2013  fillthis="hdeclxyd0d0B_";
2014  fillthis+=ptbin;
2015  ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
2016 
2017  fillthis="hnormdeclxyd0d0B_";
2018  fillthis+=ptbin;
2019  ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
2020 
2021 
2022  if(recalcvtx) {
2023 
2024  fillthis="hdeclvB_";
2025  fillthis+=ptbin;
2026  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
2027 
2028  fillthis="hnormdeclvB_";
2029  fillthis+=ptbin;
2030  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
2031 
2032 
2033  }
2034  }//mass cut
2035  }//else (background)
2036 
2037  return;
2038 }
2039 
2040 //____________________________________________________________________________
2041 void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAODMCHeader *mcHeader, AliRDHFCutsD0toKpi* cuts, TList *listout){
2042  //
2044  //
2045 
2046 
2047  Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
2048 
2049  //cout<<"is selected = "<<fIsSelectedCandidate<<endl;
2050 
2051  // Fill candidate variable Tree (track selection, no candidate selection)
2052  if( fWriteVariableTree && !part->HasBadDaughters()
2053  && fCuts->AreDaughtersSelected(part) && fCuts->IsSelectedPID(part) ){
2054  fCandidateVariables[0] = part->InvMassD0();
2055  fCandidateVariables[1] = part->InvMassD0bar();
2056  fCandidateVariables[2] = part->Pt();
2057  fCandidateVariables[3] = part->GetDCA();
2058  Double_t ctsD0=0. ,ctsD0bar=0.; part->CosThetaStarD0(ctsD0,ctsD0bar);
2059  fCandidateVariables[4] = ctsD0;
2060  fCandidateVariables[5] = ctsD0bar;
2061  fCandidateVariables[6] = part->Pt2Prong(0);
2062  fCandidateVariables[7] = part->Pt2Prong(1);
2063  fCandidateVariables[8] = part->Getd0Prong(0);
2064  fCandidateVariables[9] = part->Getd0Prong(1);
2065  fCandidateVariables[10] = part->Prodd0d0();
2066  fCandidateVariables[11] = part->CosPointingAngle();
2070  fVariablesTree->Fill();
2071  }
2072 
2073  //cout<<"check cuts = "<<endl;
2074  //cuts->PrintAll();
2075  if (!fIsSelectedCandidate){
2076  //cout<<" cut " << cuts << " Rejected because "<<cuts->GetWhy()<<endl;
2077  return;
2078  }
2079 
2080 
2081  if(fDebug>2) cout<<"Candidate selected"<<endl;
2082 
2083  Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
2084  //printf("SELECTED\n");
2085  Int_t ptbin=cuts->PtBin(part->Pt());
2086  Double_t pt = part->Pt();
2087  Double_t y = part->YD0();
2088 
2089  Double_t impparXY=part->ImpParXY()*10000.;
2090  Double_t trueImpParXY=0.;
2091  Double_t arrayForSparse[3]={invmassD0,pt,impparXY};
2092  Double_t arrayForSparseTrue[3]={invmassD0,pt,trueImpParXY};
2093 
2094 
2095  // AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
2096  // if(!prong) {
2097  // AliDebug(2,"No daughter found");
2098  // return;
2099  // }
2100  // else{
2101  // if(prong->Charge()==1) {
2102  // ((TH1F*)fDistr->FindObject("hpospair"))->Fill(fCuts->GetNPtBins()+ptbin);
2103  // //fTotPosPairs[ptbin]++;
2104  // } else {
2105  // ((TH1F*)fDistr->FindObject("hnegpair"))->Fill(fCuts->GetNPtBins()+ptbin);
2106  // //fTotNegPairs[ptbin]++;
2107  // }
2108  // }
2109 
2110  // for(Int_t it=0;it<2;it++){
2111 
2112  // //request on spd points to be addes
2113  // if(/*nSPD==2 && */part->Pt() > 5. && (TMath::Abs(invmassD0-mPDG)<0.01 || TMath::Abs(invmassD0bar-mPDG)<0.01)){
2114  // FILE *f=fopen("4display.txt","a");
2115  // fprintf(f,"pt: %f \n Rapidity: %f \t Period Number: %x \t Run Number: %d \t BunchCrossNumb: %d \t OrbitNumber: %d\n",part->Pt(),part->Y(421),aod->GetPeriodNumber(),aod->GetRunNumber(),aod->GetBunchCrossNumber(),aod->GetOrbitNumber());
2116  // fclose(f);
2117  // //printf("PrimVtx NContributors: %d \n Prongs Rel Angle: %f \n \n",ncont,relangle);
2118  // }
2119  // }
2120 
2121  TString fillthis="", fillthispt="", fillthismasspt="", fillthismassy="";
2122  Int_t pdgDgD0toKpi[2]={321,211};
2123  Int_t labD0=-1;
2124  Bool_t isPrimary=kTRUE;
2125  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)
2126 
2127  //Define weights for Bayesian (if appropriate)
2128 
2129  Double_t weigD0=1.;
2130  Double_t weigD0bar=1.;
2132  weigD0=fCuts->GetWeightsNegative()[AliPID::kKaon] * fCuts->GetWeightsPositive()[AliPID::kPion];
2133  weigD0bar=fCuts->GetWeightsPositive()[AliPID::kKaon] * fCuts->GetWeightsNegative()[AliPID::kPion];
2134  if (weigD0 > 1.0 || weigD0 < 0.) {weigD0 = 0.;}
2135  if (weigD0bar > 1.0 || weigD0bar < 0.) {weigD0bar = 0.;} //Prevents filling with weight > 1, or < 0
2136  }
2137 
2138  //count candidates selected by cuts
2139  fNentries->Fill(1);
2140  //count true D0 selected by cuts
2141  if (fReadMC && labD0>=0) fNentries->Fill(2);
2142 
2143  if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
2144 
2145  arrayForSparse[0]=invmassD0; arrayForSparse[2]=impparXY;
2146 
2147  if(fReadMC){
2148  if(labD0>=0) {
2149  if(fArray==1) cout<<"LS signal ERROR"<<endl;
2150 
2151  AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
2152  Int_t pdgD0 = partD0->GetPdgCode();
2153  // cout<<"pdg = "<<pdgD0<<endl;
2154 
2155  // old function if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE;
2156  if(AliVertexingHFUtils::CheckOrigin(arrMC,partD0,fUseQuarkTagInKine)==5) isPrimary=kFALSE;
2157  if(!isPrimary)
2158  trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
2159  arrayForSparseTrue[0]=invmassD0; arrayForSparseTrue[2]=trueImpParXY;
2160 
2161  if (pdgD0==421){ //D0
2162  // cout<<"Fill S with D0"<<endl;
2163  fillthis="histSgn_";
2164  fillthis+=ptbin;
2165  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2166 
2167  if(fFillPtHist){
2168  fillthismasspt="histSgnPt";
2169  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2170  }
2171  if(fFillImpParHist){
2172  if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0);
2173  else {
2174  fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0);
2175  fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0);
2176  }
2177  }
2178 
2179  if(fFillYHist){
2180  fillthismassy="histSgnY_";
2181  fillthismassy+=ptbin;
2182  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2183  }
2184 
2185  if(fSys==0){
2186  if(TMath::Abs(invmassD0 - mPDG) < 0.027 && fFillVarHists){
2187  fillthis="histSgn27_";
2188  fillthis+=ptbin;
2189  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2190  }
2191  }
2192  } else{ //it was a D0bar
2193  fillthis="histRfl_";
2194  fillthis+=ptbin;
2195  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2196 
2197  if(fFillPtHist){
2198  fillthismasspt="histRflPt";
2199  // cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
2200  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2201  }
2202 
2203  if(fFillYHist){
2204  fillthismassy="histRflY_";
2205  fillthismassy+=ptbin;
2206  // cout << " Filling "<<fillthismassy<<" D0bar"<<endl;
2207  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2208  }
2209 
2210  }
2211  } else {//background
2212  fillthis="histBkg_";
2213  fillthis+=ptbin;
2214  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2215 
2216  if(fFillPtHist){
2217  fillthismasspt="histBkgPt";
2218  // cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
2219  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2220  }
2221  if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0);
2222 
2223  if(fFillYHist){
2224  fillthismassy="histBkgY_";
2225  fillthismassy+=ptbin;
2226  // cout << " Filling "<<fillthismassy<<" D0bar"<<endl;
2227  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2228  }
2229 
2230  }
2231 
2232  }else{
2233  fillthis="histMass_";
2234  fillthis+=ptbin;
2235  // cout<<"Filling "<<fillthis<<endl;
2236 
2237  // printf("Fill mass with D0");
2238  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2239 
2240 
2241  if(fFillPtHist){
2242  fillthismasspt="histMassPt";
2243  // cout<<"Filling "<<fillthismasspt<<endl;
2244  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2245  }
2246  if(fFillImpParHist) {
2247  // cout << "Filling fHistMassPtImpParTC[0]"<<endl;
2248  fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0);
2249  }
2250 
2251  if(fFillYHist){
2252  fillthismassy="histMassY_";
2253  fillthismassy+=ptbin;
2254  // cout<<"Filling "<<fillthismassy<<endl;
2255  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2256  }
2257 
2258  }
2259 
2260  }
2261  if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
2262 
2263  arrayForSparse[0]=invmassD0bar; arrayForSparse[2]=impparXY;
2264 
2265  if(fReadMC){
2266  if(labD0>=0) {
2267  if(fArray==1) cout<<"LS signal ERROR"<<endl;
2268  AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
2269  Int_t pdgD0 = partD0->GetPdgCode();
2270  // cout<<" pdg = "<<pdgD0<<endl;
2271 
2272  //old function if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE;
2273  if(AliVertexingHFUtils::CheckOrigin(arrMC,partD0,fUseQuarkTagInKine)==5) isPrimary=kFALSE;
2274  if(!isPrimary)
2275  trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
2276  arrayForSparseTrue[0]=invmassD0bar; arrayForSparseTrue[2]=trueImpParXY;
2277 
2278  if (pdgD0==-421){ //D0bar
2279  fillthis="histSgn_";
2280  fillthis+=ptbin;
2281  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
2282  // if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
2283  // fillthis="histSgn27_";
2284  // fillthis+=ptbin;
2285  // ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
2286  // }
2287 
2288  if(fFillPtHist){
2289  fillthismasspt="histSgnPt";
2290  // cout<<" Filling "<< fillthismasspt << endl;
2291  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2292  }
2293  if(fFillImpParHist){
2294  // cout << " Filling impact parameter thnsparse"<<endl;
2295  if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0bar);
2296  else {
2297  fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0bar);
2298  fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0bar);
2299  }
2300  }
2301 
2302  if(fFillYHist){
2303  fillthismassy="histSgnY_";
2304  fillthismassy+=ptbin;
2305  // cout<<" Filling "<< fillthismassy << endl;
2306  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2307  }
2308 
2309  } else{
2310  fillthis="histRfl_";
2311  fillthis+=ptbin;
2312  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
2313  if(fFillPtHist){
2314  fillthismasspt="histRflPt";
2315  // cout << " Filling "<<fillthismasspt<<endl;
2316  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2317  }
2318  if(fFillYHist){
2319  fillthismassy="histRflY_";
2320  fillthismassy+=ptbin;
2321  // cout << " Filling "<<fillthismassy<<endl;
2322  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2323  }
2324  }
2325  } else {//background or LS
2326  fillthis="histBkg_";
2327  fillthis+=ptbin;
2328  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
2329 
2330  if(fFillPtHist){
2331  fillthismasspt="histBkgPt";
2332  // cout<<" Filling "<< fillthismasspt << endl;
2333  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2334  }
2335  if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0bar);
2336  if(fFillYHist){
2337  fillthismassy="histBkgY_";
2338  fillthismassy+=ptbin;
2339  // cout<<" Filling "<< fillthismassy << endl;
2340  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2341  }
2342  }
2343  }else{
2344  fillthis="histMass_";
2345  fillthis+=ptbin;
2346  // printf("Fill mass with D0bar");
2347 
2348  ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar,weigD0bar);
2349 
2350 
2351  if(fFillPtHist){
2352  fillthismasspt="histMassPt";
2353  // cout<<" Filling "<< fillthismasspt << endl;
2354  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2355  }
2356  if(fFillImpParHist) fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0bar);
2357  if(fFillYHist){
2358  fillthismassy="histMassY_";
2359  fillthismassy+=ptbin;
2360  // cout<<" Filling "<< fillthismassy << endl;
2361  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2362  }
2363  }
2364  }
2365 
2366  return;
2367 }
2368 
2369 //__________________________________________________________________________
2372 
2373  Int_t skipped[2];
2374  Int_t nTrksToSkip=2;
2375  AliAODTrack *dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(0);
2376  if(!dgTrack){
2377  AliDebug(2,"no daughter found!");
2378  return 0x0;
2379  }
2380  skipped[0]=dgTrack->GetID();
2381  dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(1);
2382  if(!dgTrack){
2383  AliDebug(2,"no daughter found!");
2384  return 0x0;
2385  }
2386  skipped[1]=dgTrack->GetID();
2387 
2388  AliESDVertex *vertexESD=0x0;
2389  AliAODVertex *vertexAOD=0x0;
2390  AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
2391 
2392  //
2393  vertexer->SetSkipTracks(nTrksToSkip,skipped);
2394  vertexer->SetMinClusters(4);
2395  vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
2396  if(!vertexESD) return vertexAOD;
2397  if(vertexESD->GetNContributors()<=0) {
2398  AliDebug(2,"vertexing failed");
2399  delete vertexESD; vertexESD=NULL;
2400  return vertexAOD;
2401  }
2402 
2403  delete vertexer; vertexer=NULL;
2404 
2405 
2406  // convert to AliAODVertex
2407  Double_t pos[3],cov[6],chi2perNDF;
2408  vertexESD->GetXYZ(pos); // position
2409  vertexESD->GetCovMatrix(cov); //covariance matrix
2410  chi2perNDF = vertexESD->GetChi2toNDF();
2411  delete vertexESD; vertexESD=NULL;
2412 
2413  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
2414  return vertexAOD;
2415 
2416 }
2417 
2418 
2419 //________________________________________________________________________
2421 {
2423  //
2424  if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
2425 
2426 
2427  fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
2428  if (!fOutputMass) {
2429  printf("ERROR: fOutputMass not available\n");
2430  return;
2431  }
2432  fOutputMassPt = dynamic_cast<TList*> (GetOutputData(6));
2433  if ((fFillPtHist || fFillImpParHist) && !fOutputMassPt) {
2434  printf("ERROR: fOutputMass not available\n");
2435  return;
2436  }
2437 
2438  if(fFillVarHists){
2439  fDistr = dynamic_cast<TList*> (GetOutputData(2));
2440  if (!fDistr) {
2441  printf("ERROR: fDistr not available\n");
2442  return;
2443  }
2444  }
2445 
2446  fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
2447 
2448  if(!fNentries){
2449  printf("ERROR: fNEntries not available\n");
2450  return;
2451  }
2452  fCuts = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(4));
2453  if(!fCuts){
2454  printf("ERROR: fCuts not available\n");
2455  return;
2456  }
2457  fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(5));
2458  if (!fCounter) {
2459  printf("ERROR: fCounter not available\n");
2460  return;
2461  }
2462  if (fDrawDetSignal) {
2463  fDetSignal = dynamic_cast<TList*>(GetOutputData(8));
2464  if (!fDetSignal) {
2465  printf("ERROR: fDetSignal not available\n");
2466  return;
2467  }
2468  }
2469  if(fFillYHist){
2470  fOutputMassY = dynamic_cast<TList*> (GetOutputData(9));
2471  if (fFillYHist && !fOutputMassY) {
2472  printf("ERROR: fOutputMassY not available\n");
2473  return;
2474  }
2475  }
2476 
2478  for(Int_t ipt=0;ipt<nptbins;ipt++){
2479 
2480  if(fArray==1 && fFillVarHists){
2481  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
2482 
2483 
2484  if(fLsNormalization>1e-6) {
2485 
2486  TString massName="histMass_";
2487  massName+=ipt;
2488  ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
2489 
2490  }
2491 
2492 
2493  fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fOutputMass->FindObject("hpospair"))->Integral(ipt+1,ipt+2)*((TH1F*)fOutputMass->FindObject("hnegpair"))->Integral(ipt+1,ipt+2));
2494  //fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
2495 
2496  if(fLsNormalization>1e-6) {
2497 
2498  TString nameDistr="hdcaB_";
2499  nameDistr+=ipt;
2500  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2501  nameDistr="hd0B_";
2502  nameDistr+=ipt;
2503  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2504  nameDistr="hd0d0B_";
2505  nameDistr+=ipt;
2506  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2507  nameDistr="hcosthetapointB_";
2508  nameDistr+=ipt;
2509  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2510  if(fSys==0){
2511  nameDistr="hptB_";
2512  nameDistr+=ipt;
2513  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2514  nameDistr="hcosthetastarB_";
2515  nameDistr+=ipt;
2516  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2517  nameDistr="hcosthpointd0d0B_";
2518  nameDistr+=ipt;
2519  ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
2520  }
2521  }
2522  }
2523  }
2524  TString cvname,cstname;
2525 
2526  if (fArray==0){
2527  cvname="D0invmass";
2528  cstname="cstat0";
2529  } else {
2530  cvname="LSinvmass";
2531  cstname="cstat1";
2532  }
2533 
2534  TCanvas *cMass=new TCanvas(cvname,cvname);
2535  cMass->cd();
2536  ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
2537 
2538  TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
2539  cStat->cd();
2540  cStat->SetGridy();
2541  fNentries->Draw("htext0");
2542 
2543  // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
2544  // ccheck->cd();
2545 
2546  return;
2547 }
2548 
2549 
2550 //________________________________________________________________________
2553 
2554  Int_t nmassbins=200;
2555  Double_t fLowmasslimit=1.5648, fUpmasslimit=2.1648;
2556  Int_t fNImpParBins=400;
2557  Double_t fLowerImpPar=-2000., fHigherImpPar=2000.;
2558  Int_t nbins[3]={nmassbins,200,fNImpParBins};
2559  Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar};
2560  Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar};
2561 
2562 
2563  fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
2564  "Mass vs. pt vs.imppar - All",
2565  3,nbins,xmin,xmax);
2566  fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
2567  "Mass vs. pt vs.imppar - promptD",
2568  3,nbins,xmin,xmax);
2569  fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
2570  "Mass vs. pt vs.imppar - DfromB",
2571  3,nbins,xmin,xmax);
2572  fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
2573  "Mass vs. pt vs.true imppar -DfromB",
2574  3,nbins,xmin,xmax);
2575  fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
2576  "Mass vs. pt vs.imppar - backgr.",
2577  3,nbins,xmin,xmax);
2578 
2579  for(Int_t i=0; i<5;i++){
2581  }
2582 }
2583 
2584 //_________________________________________________________________________________________________
2585 Float_t AliAnalysisTaskSED0Mass::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD0) const {
2587 
2588  printf(" AliAnalysisTaskSED0MassV1::GetTrueImpactParameter() \n");
2589 
2590  Double_t vtxTrue[3];
2591  mcHeader->GetVertex(vtxTrue);
2592  Double_t origD[3];
2593  partD0->XvYvZv(origD);
2594  Short_t charge=partD0->Charge();
2595  Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
2596  for(Int_t iDau=0; iDau<2; iDau++){
2597  pXdauTrue[iDau]=0.;
2598  pYdauTrue[iDau]=0.;
2599  pZdauTrue[iDau]=0.;
2600  }
2601 
2602  // Int_t nDau=partD0->GetNDaughters();
2603  Int_t labelFirstDau = partD0->GetDaughter(0);
2604 
2605  for(Int_t iDau=0; iDau<2; iDau++){
2606  Int_t ind = labelFirstDau+iDau;
2607  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
2608  if(!part) continue;
2609  Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
2610  if(!part){
2611  AliError("Daughter particle not found in MC array");
2612  return 99999.;
2613  }
2614  if(pdgCode==211 || pdgCode==321){
2615  pXdauTrue[iDau]=part->Px();
2616  pYdauTrue[iDau]=part->Py();
2617  pZdauTrue[iDau]=part->Pz();
2618  }
2619  }
2620 
2621  Double_t d0dummy[2]={0.,0.};
2622  AliAODRecoDecayHF aodDzeroMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
2623  return aodDzeroMC.ImpParXY();
2624 
2625 }
2626 
2627 //_________________________________________________________________________________________________
2628 Int_t AliAnalysisTaskSED0Mass::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
2629  // obsolete method
2631  //
2632  printf(" AliAnalysisTaskSED0Mass V1::CheckOrigin() \n");
2633 
2634  Int_t pdgGranma = 0;
2635  Int_t mother = 0;
2636  mother = mcPartCandidate->GetMother();
2637  Int_t istep = 0;
2638  Int_t abspdgGranma =0;
2639  Bool_t isFromB=kFALSE;
2640  Bool_t isQuarkFound=kFALSE;
2641  while (mother >0 ){
2642  istep++;
2643  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
2644  if (mcGranma){
2645  pdgGranma = mcGranma->GetPdgCode();
2646  abspdgGranma = TMath::Abs(pdgGranma);
2647  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
2648  isFromB=kTRUE;
2649  }
2650  if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
2651  mother = mcGranma->GetMother();
2652  }else{
2653  AliError("Failed casting the mother particle!");
2654  break;
2655  }
2656  }
2657 
2658  if(isFromB) return 5;
2659  else return 4;
2660 }
2661 //_______________________________________
2664 
2665  const Int_t nVarPrompt = 2;
2666  const Int_t nVarFD = 3;
2667 
2668  Int_t nbinsPrompt[nVarPrompt]={200,100};
2669  Int_t nbinsFD[nVarFD]={200,100,200};
2670 
2671  Double_t xminPrompt[nVarPrompt] = {0.,-1.};
2672  Double_t xmaxPrompt[nVarPrompt] = {40.,1.};
2673 
2674  Double_t xminFD[nVarFD] = {0.,-1.,0.};
2675  Double_t xmaxFD[nVarFD] = {40.,1.,40.};
2676 
2677  //pt, y
2678  fMCAccPrompt = new THnSparseF("hMCAccPrompt","kStepMCAcceptance pt vs. y - promptD",nVarPrompt,nbinsPrompt,xminPrompt,xmaxPrompt);
2679  fMCAccPrompt->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
2680  fMCAccPrompt->GetAxis(1)->SetTitle("y");
2681 
2682  //pt,y,ptB
2683  fMCAccBFeed = new THnSparseF("hMCAccBFeed","kStepMCAcceptance pt vs. y vs. ptB - DfromB",nVarFD,nbinsFD,xminFD,xmaxFD);
2684  fMCAccBFeed->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
2685  fMCAccBFeed->GetAxis(1)->SetTitle("y");
2686  fMCAccBFeed->GetAxis(2)->SetTitle("p_{T}^{B} (GeV/c)");
2687 
2688  fOutputMass->Add(fMCAccPrompt);
2689  fOutputMass->Add(fMCAccBFeed);
2690 }
2691 //___________________________________________________________________________________________________
2692 void AliAnalysisTaskSED0Mass::FillMCAcceptanceHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader){
2694  const Int_t nProng = 2;
2695  Double_t zMCVertex = mcHeader->GetVtxZ(); //vertex MC
2696 
2697  for(Int_t iPart=0; iPart<arrayMC->GetEntriesFast(); iPart++){
2698  AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(iPart));
2699  if (TMath::Abs(mcPart->GetPdgCode()) == 421){
2700 
2701  Int_t orig=AliVertexingHFUtils::CheckOrigin(arrayMC,mcPart,fUseQuarkTagInKine);//Prompt = 4, FeedDown = 5
2702 
2703  Int_t deca = 0;
2704  Bool_t isGoodDecay=kFALSE;
2705  Int_t labDau[4]={-1,-1,-1,-1};
2706  Bool_t isInAcc = kFALSE;
2707  Bool_t isFidAcc = kFALSE;
2708 
2709  deca=AliVertexingHFUtils::CheckD0Decay(arrayMC,mcPart,labDau);
2710  if(deca > 0) isGoodDecay=kTRUE;
2711 
2712  if(labDau[0]==-1){
2713  continue; //protection against unfilled array of labels
2714  }
2715 
2716  isFidAcc=fCuts->IsInFiducialAcceptance(mcPart->Pt(),mcPart->Y());
2717  isInAcc=CheckAcc(arrayMC,nProng,labDau);
2718 
2719  if(isGoodDecay && TMath::Abs(zMCVertex) < fCuts->GetMaxVtxZ() && isFidAcc && isInAcc) {
2720  //for prompt
2721  if(orig == 4){
2722  //fill histo for prompt
2723  Double_t arrayMCprompt[2] = {mcPart->Pt(),mcPart->Y()};
2724  fMCAccPrompt->Fill(arrayMCprompt);
2725  }
2726  //for FD
2727  else if(orig == 5){
2728  Double_t ptB = AliVertexingHFUtils::GetBeautyMotherPt(arrayMC,mcPart);
2729  //fill histo for FD
2730  Double_t arrayMCFD[3] = {mcPart->Pt(),mcPart->Y(),ptB};
2731  fMCAccBFeed->Fill(arrayMCFD);
2732  }
2733  else
2734  continue;
2735  }
2736  }
2737  }
2738 }
2739 //______________________________________________________________________
2740 Bool_t AliAnalysisTaskSED0Mass::CheckAcc(TClonesArray* arrayMC,Int_t nProng, Int_t *labDau){
2742  for (Int_t iProng = 0; iProng<nProng; iProng++){
2743  AliAODMCParticle* mcPartDaughter=dynamic_cast<AliAODMCParticle*>(arrayMC->At(labDau[iProng]));
2744  if(!mcPartDaughter) return kFALSE;
2745  Double_t eta = mcPartDaughter->Eta();
2746  Double_t pt = mcPartDaughter->Pt();
2747  if (TMath::Abs(eta) > 0.9 || pt < 0.1) return kFALSE;
2748  }
2749  return kTRUE;
2750 }
2751 //________________________________________
2753  Double_t normIP[2];
2754  Double_t pointD0[9];
2755  Double_t pointD0bar[9];
2756  Double_t pointD0MC[9];
2757  Float_t ptB;
2758  Int_t signalType=0;//background
2759 
2760  //fill variables
2761  Double_t diffIP[2], errdiffIP[2];
2762  part->Getd0MeasMinusExpProng(0,aod->GetMagneticField(),diffIP[0],errdiffIP[0]);
2763  part->Getd0MeasMinusExpProng(1,aod->GetMagneticField(),diffIP[1],errdiffIP[1]);
2764  normIP[0]=diffIP[0]/errdiffIP[0];
2765  normIP[1]=diffIP[1]/errdiffIP[1];
2766 
2767  AliAODVertex* secvtx=part->GetSecondaryVtx();
2768  AliAODVertex *primvtx=part->GetPrimaryVtx();
2769  Double_t err2decaylength=secvtx->Error2DistanceXYToVertex(primvtx);
2770  Double_t lxy=part->AliAODRecoDecay::DecayLengthXY(primvtx);
2771  Bool_t ispid=fCuts->GetIsUsePID();
2772 
2773  if(fReadMC){
2774  // MC THnSparse
2775  Int_t pdgDgD0toKpi[2]={321,211};
2776  Int_t lab=-9999;
2777  lab=part->MatchToMC(421,arrMC,2,pdgDgD0toKpi); //return MC particle label if the array corresponds to a D0, -1 if not (cf. AliAODRecoDecay.cxx)
2778  if(lab>=0){
2779  signalType=1;//signal
2780  AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(lab);
2782  if(orig==4)signalType=1;//prompt
2783  else if(orig==5)signalType=2;//feed down
2784  //pt, ptB, normImpParTrk1, normImpParTrk2, decLXY, normDecLXY, iscut, ispid
2785  if(ispid)fCuts->SetUsePID(kFALSE);// if PID on, switch it off
2786  Int_t isCuts=fCuts->IsSelected(part,AliRDHFCuts::kAll,aod);
2787  if(ispid)fCuts->SetUsePID(kTRUE);//if PID was on, switch it on
2788  if(!isCuts) return;
2789  Int_t isPid= fCuts->IsSelectedPID(part);
2790  ptB=AliVertexingHFUtils::GetBeautyMotherPt(arrMC,partD0);
2791  Double_t arrayMC[8]={part->Pt(), ptB,normIP[0], normIP[1], lxy, lxy/TMath::Sqrt(err2decaylength),(Double_t)isCuts, (Double_t)isPid};
2792 
2793  //fill pointD0MC
2794  for(Int_t i=0; i<8; i++){
2795  pointD0MC[i]=arrayMC[i];
2796  }
2797  if(signalType==1)fhStudyImpParSingleTrackSign->Fill(pointD0MC);
2798  if(signalType==2){
2799  fhStudyImpParSingleTrackFd->Fill(pointD0MC);
2800  }
2801  }
2802  }else{
2803  // data
2804  if(ispid)fCuts->SetUsePID(kFALSE);// if PID on, switch it off
2805  Int_t IsSelectedPIDoff=fCuts->IsSelected(part,AliRDHFCuts::kAll,aod);
2806  if(ispid)fCuts->SetUsePID(kTRUE);//if PID was on, switch it on
2807  if(!IsSelectedPIDoff) return;
2808  Int_t pid=fCuts->IsSelectedPID(part);
2809  if((IsSelectedPIDoff==1 || IsSelectedPIDoff==3) && fFillOnlyD0D0bar<2){
2810  Double_t array[9]={part->Pt(),normIP[0],normIP[1],lxy,lxy/TMath::Sqrt(err2decaylength),part->InvMassD0(),(Double_t)IsSelectedPIDoff,(Double_t)pid,0.};
2811  for(Int_t i=0;i<9;i++){
2812  pointD0[i]=array[i];
2813  }
2814  fhStudyImpParSingleTrackCand->Fill(pointD0);
2815  }
2816  if (IsSelectedPIDoff>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)){
2817  Double_t arrayD0bar[9]={part->Pt(),normIP[0],normIP[1],lxy,lxy/TMath::Sqrt(err2decaylength),part->InvMassD0bar(),(Double_t)IsSelectedPIDoff,(Double_t)pid,1.};
2818  for(Int_t i=0;i<9;i++){
2819  pointD0bar[i]=arrayD0bar[i];
2820  }
2821  fhStudyImpParSingleTrackCand->Fill(pointD0bar);
2822  }
2823  }
2824 }
Int_t charge
THnSparseF * fMCAccBFeed
!histo for StepMCAcc for D0 FD (pt,y,ptB)
Double_t NormalizedDecayLengthXY() const
TList * fOutputMassPt
! list send on output slot 6
Bool_t fDrawDetSignal
flag to decide whether to fill "PID = x" bins in fNentrie
double Double_t
Definition: External.C:58
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
Bool_t fFillYHist
flag to fill Pt and Impact Parameter Histograms
void StoreCandidates(AliVEvent *, Int_t nCand=0, Bool_t flagFilter=kTRUE)
Definition: External.C:236
void Getd0MeasMinusExpProng(Int_t ip, Double_t magf, Double_t &d0diff, Double_t &errd0diff) const
void FillVarHists(AliAODEvent *aodev, AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout)
Bool_t GetCombPID() const
Double_t * fCandidateVariables
! variables to be written to the tree
Bool_t HasSelectionBit(Int_t i) const
Double_t NormalizedDecayLength2() const
Double_t ImpParXY() const
Int_t fAODProtection
flag to check or not the selection bit
static Int_t CheckMatchingAODdeltaAODevents()
TList * fDistr
! list send on output slot 2
void NormIPvar(AliAODEvent *aod, AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC)
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:299
TList * fDetSignal
!Detector signal histograms (on output slot 8)
Double_t CosPointingAngleXY() const
Bool_t fUseSelectionBit
flag to fill Pt and Impact Parameter Histograms
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
void FillMCAcceptanceHistos(TClonesArray *arrayMC, AliAODMCHeader *mcHeader)
Int_t fSys
flag to enable filling variable histos
ULong_t GetSelectionMap() const
TList * fOutputMass
! list send on output slot 1
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
Bool_t CheckAcc(TClonesArray *arrayMC, Int_t nProng, Int_t *labDau)
Bool_t fReadMC
can be D0 or Like Sign candidates
static Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPart, Bool_t searchUpToQuark=kTRUE)
Double_t GetMaxVtxZ() const
Definition: AliRDHFCuts.h:242
Int_t GetBayesianStrategy() const
Bool_t HasBadDaughters() const
THnSparseF * fHistMassPtImpParTC[5]
! histograms for impact paramter studies
TH1F * fNentries
! histogram with number of events on output slot 3
virtual void UserCreateOutputObjects()
Implementation of interface methods.
Bool_t fFillPtHist
flag to reject events with SDD clusters
Bool_t fFillImpParHist
flag to fill Y Histograms
int Int_t
Definition: External.C:63
Definition: External.C:204
unsigned int UInt_t
Definition: External.C:33
Int_t GetIsFilled() const
Double_t CosThetaStarD0bar() const
angle of K
float Float_t
Definition: External.C:68
Bool_t fUseQuarkTagInKine
flag to decide whether to draw the TPC dE/dx and TOF signal before/after PID
Bool_t fIsRejectSDDClusters
fSys=0 -> p-p; fSys=1 ->PbPb (in this case fFillVarHists=kFALSE by default: set it to kTRUE after if ...
static Int_t CheckD0Decay(AliStack *stack, Int_t label, Int_t *arrayDauLab)
Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPartCandidate) const
virtual void Terminate(Option_t *option)
Float_t GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray *arrayMC, AliAODMCParticle *partD0) const
void DrawDetSignal(AliAODRecoDecayHF2Prong *part, TList *ListDetSignal)
Double_t fLsNormalization
number of pt bins
Bool_t fUsePid4Distr
flag to decide if apply cut also on distributions: 0 no cuts, 1 looser cuts, 2 tighter/ cuts ...
Double_t DecayLength2() const
kinematics & topology
THnSparseF * fhStudyImpParSingleTrackFd
! sparse with imp par residual cuts for MC
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)
Bool_t fCutOnDistr
flag for MC array: kTRUE = read it, kFALSE = do not read it
AliAODVertex * GetPrimaryVtxSkipped(AliAODEvent *aodev)
short Short_t
Definition: External.C:23
Int_t IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
THnSparseF * fhStudyImpParSingleTrackSign
! sparse with imp par residual cuts for MC
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Double_t DecayLengthXY() const
static Double_t GetBeautyMotherPt(TClonesArray *arrayMC, AliAODMCParticle *mcPart)
void FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAODMCHeader *mcHeader, AliRDHFCutsD0toKpi *cuts, TList *listout)
Int_t fIsSelectedCandidate
keeps the daughter tracks
Bool_t GetIsPrimaryWithoutDaughters() const
Definition: AliRDHFCuts.h:263
Double_t * GetWeightsPositive() const
Bool_t AreDaughtersSelected(AliAODRecoDecayHF *rd, const AliAODEvent *aod=0x0) const
Bool_t IsEventSelected(AliVEvent *event)
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
void SetUsePID(Bool_t flag=kTRUE)
Definition: AliRDHFCuts.h:205
THnSparseF * fhStudyImpParSingleTrackCand
! sparse with imp par residual cuts for Data
const Int_t nVar
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
virtual Int_t IsSelectedPID(AliAODRecoDecayHF *rd)
TList * fOutputMassY
! list send on output slot 9
TTree * fVariablesTree
flag to decide whether to write the candidate variables on a tree variables
TH2 * Scale(TH2 *h, TH1 *g)
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Bool_t GetIsUsePID() const
Definition: AliRDHFCuts.h:255
const char Option_t
Definition: External.C:48
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:234
AliAODVertex * GetPrimaryVtx() const
virtual void UserExec(Option_t *option)
TObjArray fDaughterTracks
flag to fill mass histogram with D0/D0bar only (0 = fill with both, 1 = fill with D0 only...
const Int_t nbins
bool Bool_t
Definition: External.C:53
Double_t CosPointingAngle() const
Double_t DecayLengthXYError() const
Int_t fFillOnlyD0D0bar
normalization
THnSparseF * fMCAccPrompt
!histo for StepMCAcc for D0 prompt (pt,y,ptB)
Int_t PtBin(Double_t pt) const
Int_t nptbins
Bool_t fFillVarHists
selection outcome
AliNormalizationCounter * fCounter
! AliNormalizationCounter on output slot 5
Double_t * GetWeightsNegative() const