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