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