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