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