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