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