AliPhysics  7dff9e1 (7dff9e1)
 All Classes Namespaces Files Functions Variables Typedefs 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"
50 #include "AliAODRecoDecayHF2Prong.h"
51 #include "AliAODRecoCascadeHF.h"
52 #include "AliAnalysisVertexingHF.h"
53 #include "AliAnalysisTaskSE.h"
54 #include "AliAnalysisTaskSED0Mass.h"
55 #include "AliNormalizationCounter.h"
56 
57 using std::cout;
58 using std::endl;
59 
61 ClassImp(AliAnalysisTaskSED0Mass);
63 
64 //________________________________________________________________________
65 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass():
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 //________________________________________________________________________
103 AliAnalysisTaskSED0Mass::AliAnalysisTaskSED0Mass(const char *name,AliRDHFCutsD0toKpi* cuts):
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 //________________________________________________________________________
163 AliAnalysisTaskSED0Mass::~AliAnalysisTaskSED0Mass()
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 //________________________________________________________________________
208 void AliAnalysisTaskSED0Mass::Init()
209 {
211 
212  if(fDebug > 1) printf("AnalysisTaskSED0Mass::Init() \n");
213 
214 
215  AliRDHFCutsD0toKpi* copyfCuts=new AliRDHFCutsD0toKpi(*fCuts);
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 //________________________________________________________________________
226 void AliAnalysisTaskSED0Mass::UserCreateOutputObjects()
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 
558  if(fCuts->GetIsPrimaryWithoutDaughters()){
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 
854  if(fFillImpParHist) CreateImpactParameterHistos();
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", 19,-0.5,18.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()->SetBinLabel(19,"D0 failed to be filled");
917  fNentries->GetXaxis()->SetNdivisions(1,kFALSE);
918 
919  fCounter = new AliNormalizationCounter(Form("%s",GetOutputSlot(5)->GetContainer()->GetName()));
920  fCounter->Init();
921 
922  //
923  // Output slot 7 : tree of the candidate variables after track selection
924  //
925  nameoutput = GetOutputSlot(7)->GetContainer()->GetName();
926  fVariablesTree = new TTree(nameoutput,"Candidates variables tree");
927  Int_t nVar = 15;
928  fCandidateVariables = new Double_t [nVar];
929  TString * fCandidateVariableNames = new TString[nVar];
930  fCandidateVariableNames[0] = "massD0";
931  fCandidateVariableNames[1] = "massD0bar";
932  fCandidateVariableNames[2] = "pt";
933  fCandidateVariableNames[3] = "dca";
934  fCandidateVariableNames[4] = "costhsD0";
935  fCandidateVariableNames[5] = "costhsD0bar";
936  fCandidateVariableNames[6] = "ptk";
937  fCandidateVariableNames[7] = "ptpi";
938  fCandidateVariableNames[8] = "d0k";
939  fCandidateVariableNames[9] = "d0pi";
940  fCandidateVariableNames[10] = "d0xd0";
941  fCandidateVariableNames[11] = "costhp";
942  fCandidateVariableNames[12] = "costhpxy";
943  fCandidateVariableNames[13] = "lxy";
944  fCandidateVariableNames[14] = "specialcuts";
945  for(Int_t ivar=0; ivar<nVar; ivar++){
946  fVariablesTree->Branch(fCandidateVariableNames[ivar].Data(),&fCandidateVariables[ivar],Form("%s/d",fCandidateVariableNames[ivar].Data()));
947  }
948 
949  //
950  // Output slot 8 : List for detector response histograms
951  //
952  if (fDrawDetSignal) {
953  TH2F *TOFSigBefPID = new TH2F("TOFSigBefPID", "TOF signal of daughters before PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
954  TH2F *TOFSigAftPID = new TH2F("TOFSigAftPID", "TOF signal after PID;p(daught)(GeV/c);Signal", 500, 0, 10, 5000, 0, 50e3);
955 
956  TH2F *TPCSigBefPID = new TH2F("TPCSigBefPID", "TPC dE/dx before PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
957  TH2F *TPCSigAftPID = new TH2F("TPCSigAftPID", "TPC dE/dx after PID;p(daught)(GeV/c);dE/dx (arb. units)", 1000, 0, 10, 1000, 0, 500);
958 
959  fDetSignal->Add(TOFSigBefPID);
960  fDetSignal->Add(TOFSigAftPID);
961  fDetSignal->Add(TPCSigBefPID);
962  fDetSignal->Add(TPCSigAftPID);
963  }
964 
965  // Post the data
966  PostData(1,fOutputMass);
967  PostData(2,fDistr);
968  PostData(3,fNentries);
969  PostData(5,fCounter);
970  PostData(6,fOutputMassPt);
971  PostData(7,fVariablesTree);
972  PostData(8, fDetSignal);
973  PostData(9,fOutputMassY);
974 
975  return;
976 }
977 
978 //________________________________________________________________________
979 void AliAnalysisTaskSED0Mass::UserExec(Option_t */*option*/)
980 {
983 
984  //cuts order
985  // printf(" |M-MD0| [GeV] < %f\n",fD0toKpiCuts[0]);
986  // printf(" dca [cm] < %f\n",fD0toKpiCuts[1]);
987  // printf(" cosThetaStar < %f\n",fD0toKpiCuts[2]);
988  // printf(" pTK [GeV/c] > %f\n",fD0toKpiCuts[3]);
989  // printf(" pTpi [GeV/c] > %f\n",fD0toKpiCuts[4]);
990  // printf(" |d0K| [cm] < %f\n",fD0toKpiCuts[5]);
991  // printf(" |d0pi| [cm] < %f\n",fD0toKpiCuts[6]);
992  // printf(" d0d0 [cm^2] < %f\n",fD0toKpiCuts[7]);
993  // printf(" cosThetaPoint > %f\n",fD0toKpiCuts[8]);
994 
995  AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
996  TString bname;
997  if(fArray==0){ //D0 candidates
998  // load D0->Kpi candidates
999  //cout<<"D0 candidates"<<endl;
1000  bname="D0toKpi";
1001  } else { //LikeSign candidates
1002  // load like sign candidates
1003  //cout<<"LS candidates"<<endl;
1004  bname="LikeSign2Prong";
1005  }
1006  TClonesArray *inputArray=0;
1007  if(!aod && AODEvent() && IsStandardAOD()) {
1008  // In case there is an AOD handler writing a standard AOD, use the AOD
1009  // event in memory rather than the input (ESD) event.
1010  aod = dynamic_cast<AliAODEvent*> (AODEvent());
1011  // in this case the braches in the deltaAOD (AliAOD.VertexingHF.root)
1012  // have to taken from the AOD event hold by the AliAODExtension
1013  AliAODHandler* aodHandler = (AliAODHandler*)
1014  ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
1015 
1016  if(aodHandler->GetExtensions()) {
1017  AliAODExtension *ext = (AliAODExtension*)aodHandler->GetExtensions()->FindObject("AliAOD.VertexingHF.root");
1018  AliAODEvent* aodFromExt = ext->GetAOD();
1019  inputArray=(TClonesArray*)aodFromExt->GetList()->FindObject(bname.Data());
1020  }
1021  } else if(aod) {
1022  inputArray=(TClonesArray*)aod->GetList()->FindObject(bname.Data());
1023  }
1024 
1025  if(!inputArray || !aod) {
1026  printf("AliAnalysisTaskSED0Mass::UserExec: input branch not found!\n");
1027  return;
1028  }
1029  // fix for temporary bug in ESDfilter
1030  // the AODs with null vertex pointer didn't pass the PhysSel
1031  if(!aod->GetPrimaryVertex() || TMath::Abs(aod->GetMagneticField())<0.001) return;
1032 
1033 
1034  TClonesArray *mcArray = 0;
1035  AliAODMCHeader *mcHeader = 0;
1036 
1037  if(fReadMC) {
1038  // load MC particles
1039  mcArray = (TClonesArray*)aod->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1040  if(!mcArray) {
1041  printf("AliAnalysisTaskSED0Mass::UserExec: MC particles branch not found!\n");
1042  return;
1043  }
1044 
1045  // load MC header
1046  mcHeader = (AliAODMCHeader*)aod->GetList()->FindObject(AliAODMCHeader::StdBranchName());
1047  if(!mcHeader) {
1048  printf("AliAnalysisTaskSED0Mass::UserExec: MC header branch not found!\n");
1049  return;
1050  }
1051  }
1052 
1053  //printf("VERTEX Z %f %f\n",vtx1->GetZ(),mcHeader->GetVtxZ());
1054 
1055  //histogram filled with 1 for every AOD
1056  fNentries->Fill(0);
1057  fCounter->StoreEvent(aod,fCuts,fReadMC);
1058  //fCounter->StoreEvent(aod,fReadMC);
1059 
1060  // trigger class for PbPb C0SMH-B-NOPF-ALLNOTRD, C0SMH-B-NOPF-ALL
1061  TString trigclass=aod->GetFiredTriggerClasses();
1062  if(trigclass.Contains("C0SMH-B-NOPF-ALLNOTRD") || trigclass.Contains("C0SMH-B-NOPF-ALL")) fNentries->Fill(14);
1063 
1064  if(!fCuts->IsEventSelected(aod)) {
1065  if(fCuts->GetWhyRejection()==1) // rejected for pileup
1066  fNentries->Fill(13);
1067  if(fSys==1 && (fCuts->GetWhyRejection()==2 || fCuts->GetWhyRejection()==3)) fNentries->Fill(15);
1068  if(fCuts->GetWhyRejection()==7) fNentries->Fill(17);
1069  return;
1070  }
1071 
1072  // Check the Nb of SDD clusters
1073  if (fIsRejectSDDClusters) {
1074  Bool_t skipEvent = kFALSE;
1075  Int_t ntracks = 0;
1076  if (aod) ntracks = aod->GetNumberOfTracks();
1077  for(Int_t itrack=0; itrack<ntracks; itrack++) { // loop on tacks
1078  // ... get the track
1079  AliAODTrack * track = dynamic_cast<AliAODTrack*>(aod->GetTrack(itrack));
1080  if(!track) AliFatal("Not a standard AOD");
1081  if(TESTBIT(track->GetITSClusterMap(),2) || TESTBIT(track->GetITSClusterMap(),3) ){
1082  skipEvent=kTRUE;
1083  fNentries->Fill(16);
1084  break;
1085  }
1086  }
1087  if (skipEvent) return;
1088  }
1089 
1090  // AOD primary vertex
1091  AliAODVertex *vtx1 = (AliAODVertex*)aod->GetPrimaryVertex();
1092 
1093  Bool_t isGoodVtx=kFALSE;
1094 
1095  //vtx1->Print();
1096  TString primTitle = vtx1->GetTitle();
1097  if(primTitle.Contains("VertexerTracks") && vtx1->GetNContributors()>0) {
1098  isGoodVtx=kTRUE;
1099  fNentries->Fill(3);
1100  }
1101 
1102  // loop over candidates
1103  Int_t nInD0toKpi = inputArray->GetEntriesFast();
1104  if(fDebug>2) printf("Number of D0->Kpi: %d\n",nInD0toKpi);
1105 
1106  // FILE *f=fopen("4display.txt","a");
1107  // printf("Number of D0->Kpi: %d\n",nInD0toKpi);
1108  Int_t nSelectedloose=0,nSelectedtight=0;
1109 
1110  // vHF object is needed to call the method that refills the missing info of the candidates
1111  // if they have been deleted in dAOD reconstruction phase
1112  // in order to reduce the size of the file
1113  AliAnalysisVertexingHF *vHF = new AliAnalysisVertexingHF();
1114 
1115  for (Int_t iD0toKpi = 0; iD0toKpi < nInD0toKpi; iD0toKpi++) {
1116  AliAODRecoDecayHF2Prong *d = (AliAODRecoDecayHF2Prong*)inputArray->UncheckedAt(iD0toKpi);
1117 
1118  if(fUseSelectionBit && d->GetSelectionMap()) if(!d->HasSelectionBit(AliRDHFCuts::kD0toKpiCuts)){
1119  fNentries->Fill(2);
1120  continue; //skip the D0 from Dstar
1121  }
1122 
1123  if(!(vHF->FillRecoCand(aod,d))) {//Fill the data members of the candidate only if they are empty.
1124  fNentries->Fill(18); //monitor how often this fails
1125  continue;
1126  }
1127 
1128  if ( fCuts->IsInFiducialAcceptance(d->Pt(),d->Y(421)) ) {
1129  nSelectedloose++;
1130  nSelectedtight++;
1131  if(fSys==0){
1132  if(fCuts->IsSelected(d,AliRDHFCuts::kTracks,aod))fNentries->Fill(6);
1133  }
1134  Int_t ptbin=fCuts->PtBin(d->Pt());
1135  if(ptbin==-1) {fNentries->Fill(4); continue;} //out of bounds
1136  fIsSelectedCandidate=fCuts->IsSelected(d,AliRDHFCuts::kAll,aod); //selected
1137  if(fFillVarHists) {
1138  //if(!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate)) {
1139  fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(0),0);
1140  fDaughterTracks.AddAt((AliAODTrack*)d->GetDaughter(1),1);
1141  //check daughters
1142  if(!fDaughterTracks.UncheckedAt(0) || !fDaughterTracks.UncheckedAt(1)) {
1143  AliDebug(1,"at least one daughter not found!");
1144  fNentries->Fill(5);
1145  fDaughterTracks.Clear();
1146  continue;
1147  }
1148  //}
1149  FillVarHists(aod,d,mcArray,fCuts,fDistr);
1150  }
1151 
1152  if (fDrawDetSignal) {
1153  DrawDetSignal(d, fDetSignal);
1154  }
1155 
1156  FillMassHists(d,mcArray,mcHeader,fCuts,fOutputMass);
1157  if (fPIDCheck) {
1158  Int_t isSelectedPIDfill = 3;
1159  if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(d); //0 rejected,1 D0,2 Dobar, 3 both
1160 
1161  if (isSelectedPIDfill == 0)fNentries->Fill(7);
1162  if (isSelectedPIDfill == 1)fNentries->Fill(8);
1163  if (isSelectedPIDfill == 2)fNentries->Fill(9);
1164  if (isSelectedPIDfill == 3)fNentries->Fill(10);
1165  }
1166 
1167  }
1168 
1169  fDaughterTracks.Clear();
1170  //if(unsetvtx) d->UnsetOwnPrimaryVtx();
1171  } //end for prongs
1172  fCounter->StoreCandidates(aod,nSelectedloose,kTRUE);
1173  fCounter->StoreCandidates(aod,nSelectedtight,kFALSE);
1174  delete vHF;
1175  // Post the data
1176  PostData(1,fOutputMass);
1177  PostData(2,fDistr);
1178  PostData(3,fNentries);
1179  PostData(5,fCounter);
1180  PostData(6,fOutputMassPt);
1181  PostData(7,fVariablesTree);
1182  PostData(8, fDetSignal);
1183 
1184  return;
1185 }
1186 
1187 //____________________________________________________________________________
1188 void AliAnalysisTaskSED0Mass::DrawDetSignal(AliAODRecoDecayHF2Prong *part, TList *ListDetSignal)
1189 {
1190  //
1192  //
1193  fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(0), 0);
1194  fDaughterTracks.AddAt((AliAODTrack*)part->GetDaughter(1), 1);
1195 
1196  AliESDtrack *esdtrack1 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1197  AliESDtrack *esdtrack2 = new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
1198 
1199 
1200  //For filling detector histograms
1201  Int_t isSelectedPIDfill = 3;
1202  if (!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPIDfill = fCuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
1203 
1204 
1205  //fill "before PID" histos with every daughter
1206  ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
1207  ((TH2F*)ListDetSignal->FindObject("TOFSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
1208  ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
1209  ((TH2F*)ListDetSignal->FindObject("TPCSigBefPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
1210 
1211  if (isSelectedPIDfill != 0) { //fill "After PID" with everything that's not rejected
1212  ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTOFsignal());
1213  ((TH2F*)ListDetSignal->FindObject("TOFSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTOFsignal());
1214  ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack1->P(), esdtrack1->GetTPCsignal());
1215  ((TH2F*)ListDetSignal->FindObject("TPCSigAftPID"))->Fill(esdtrack2->P(), esdtrack2->GetTPCsignal());
1216 
1217  }
1218 
1219  delete esdtrack1;
1220  delete esdtrack2;
1221 
1222  return;
1223 }
1224 
1225 //____________________________________________________________________________
1226 void AliAnalysisTaskSED0Mass::FillVarHists(AliAODEvent* aod,AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout){
1227  //
1229  //
1230 
1231 
1232  Int_t pdgDgD0toKpi[2]={321,211};
1233  Int_t lab=-9999;
1234  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)
1235  //Double_t pt = d->Pt(); //mother pt
1236  Int_t isSelectedPID=3;
1237  if(!fReadMC || (fReadMC && fUsePid4Distr)) isSelectedPID=cuts->IsSelectedPID(part); //0 rejected,1 D0,2 Dobar, 3 both
1238  if (!fPIDCheck) { //if fPIDCheck, this is already filled elsewhere
1239  if (isSelectedPID==0)fNentries->Fill(7);
1240  if (isSelectedPID==1)fNentries->Fill(8);
1241  if (isSelectedPID==2)fNentries->Fill(9);
1242  if (isSelectedPID==3)fNentries->Fill(10);
1243  //fNentries->Fill(8+isSelectedPID);
1244  }
1245 
1246  if(fCutOnDistr && !fIsSelectedCandidate) return;
1247  //printf("\nif no cuts or cuts passed\n");
1248 
1249 
1250  //add distr here
1251  UInt_t pdgs[2];
1252 
1253  Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1254  pdgs[0]=211;
1255  pdgs[1]=321;
1256  Double_t minvD0 = part->InvMassD0();
1257  pdgs[1]=211;
1258  pdgs[0]=321;
1259  Double_t minvD0bar = part->InvMassD0bar();
1260  //cout<<"inside mass cut"<<endl;
1261 
1262  Double_t invmasscut=0.03;
1263 
1264  TString fillthispi="",fillthisK="",fillthis="", fillthispt="", fillthisetaphi="";
1265 
1266  Int_t ptbin=cuts->PtBin(part->Pt());
1267  Double_t pt = part->Pt();
1268 
1269  Double_t dz1[2],dz2[2],covar1[3],covar2[3];//,d0xd0proper,errd0xd0proper;
1270  dz1[0]=-99; dz2[0]=-99;
1271  Double_t d0[2];
1272  Double_t decl[2]={-99,-99};
1273  Bool_t recalcvtx=kFALSE;
1274 
1275 
1276 
1277  if(fCuts->GetIsPrimaryWithoutDaughters()){
1278  recalcvtx=kTRUE;
1279  //recalculate vertex
1280  AliAODVertex *vtxProp=0x0;
1281  vtxProp=GetPrimaryVtxSkipped(aod);
1282  if(vtxProp) {
1283  part->SetOwnPrimaryVtx(vtxProp);
1284  //Bool_t unsetvtx=kTRUE;
1285  //Calculate d0 for daughter tracks with Vtx Skipped
1286  AliESDtrack *esdtrack1=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(0));
1287  AliESDtrack *esdtrack2=new AliESDtrack((AliVTrack*)fDaughterTracks.UncheckedAt(1));
1288  esdtrack1->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz1,covar1);
1289  esdtrack2->PropagateToDCA(vtxProp,aod->GetMagneticField(),1.,dz2,covar2);
1290  delete vtxProp; vtxProp=NULL;
1291  delete esdtrack1;
1292  delete esdtrack2;
1293  }
1294 
1295  d0[0]=dz1[0];
1296  d0[1]=dz2[0];
1297 
1298  decl[0]=part->DecayLength2();
1299  decl[1]=part->NormalizedDecayLength2();
1300  part->UnsetOwnPrimaryVtx();
1301 
1302  }
1303 
1304  Double_t cosThetaStarD0 = 99;
1305  Double_t cosThetaStarD0bar = 99;
1306  Double_t cosPointingAngle = 99;
1307  Double_t normalizedDecayLength2 = -1, normalizedDecayLengthxy=-1;
1308  Double_t decayLength2 = -1, decayLengthxy=-1;
1309  Double_t ptProng[2]={-99,-99};
1310  Double_t d0Prong[2]={-99,-99};
1311  Double_t etaD = 99.;
1312  Double_t phiD = 99.;
1313 
1314 
1315  //disable the PID
1316  if(!fUsePid4Distr) isSelectedPID=0;
1317  if((lab>=0 && fReadMC) || (!fReadMC && isSelectedPID)){ //signal (from MC or PID)
1318 
1319  //check pdg of the prongs
1320  AliAODTrack *prong0=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1321  AliAODTrack *prong1=(AliAODTrack*)fDaughterTracks.UncheckedAt(1);
1322 
1323  if(!prong0 || !prong1) {
1324  return;
1325  }
1326  Int_t labprong[2];
1327  if(fReadMC){
1328  labprong[0]=prong0->GetLabel();
1329  labprong[1]=prong1->GetLabel();
1330  }
1331  AliAODMCParticle *mcprong=0;
1332  Int_t pdgProng[2]={0,0};
1333 
1334  for (Int_t iprong=0;iprong<2;iprong++){
1335  if(fReadMC && labprong[iprong]>=0) {
1336  mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1337  pdgProng[iprong]=mcprong->GetPdgCode();
1338  }
1339  }
1340 
1341  if(fSys==0){
1342  //no mass cut ditributions: ptbis
1343 
1344  fillthispi="hptpiSnoMcut_";
1345  fillthispi+=ptbin;
1346 
1347  fillthisK="hptKSnoMcut_";
1348  fillthisK+=ptbin;
1349 
1350  if ((TMath::Abs(pdgProng[0]) == 211 && TMath::Abs(pdgProng[1]) == 321)
1351  || (isSelectedPID==1 || isSelectedPID==3)){
1352  ((TH1F*)listout->FindObject(fillthispi))->Fill(prong0->Pt());
1353  ((TH1F*)listout->FindObject(fillthisK))->Fill(prong1->Pt());
1354  }
1355 
1356  if ((TMath::Abs(pdgProng[0]) == 321 && TMath::Abs(pdgProng[1]) == 211)
1357  || (isSelectedPID==2 || isSelectedPID==3)){
1358  ((TH1F*)listout->FindObject(fillthisK))->Fill(prong0->Pt());
1359  ((TH1F*)listout->FindObject(fillthispi))->Fill(prong1->Pt());
1360  }
1361  }
1362 
1363  //no mass cut ditributions: mass
1364 
1365  etaD = part->Eta();
1366  phiD = part->Phi();
1367 
1368 
1369  fillthis="hMassS_";
1370  fillthis+=ptbin;
1371  fillthispt="histSgnPt";
1372 
1373  if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)
1374  || (!fReadMC && (isSelectedPID==1 || isSelectedPID==3))){//D0
1375  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1376  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
1377 
1378  fillthisetaphi="hetaphiD0candidateS_";
1379  fillthisetaphi+=ptbin;
1380  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1381 
1382  if(TMath::Abs(minvD0-mPDG)<0.05){
1383  fillthisetaphi="hetaphiD0candidatesignalregionS_";
1384  fillthisetaphi+=ptbin;
1385  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1386  }
1387 
1388  }
1389  else { //D0bar
1390  if(fReadMC || (!fReadMC && isSelectedPID > 1)){
1391  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1392  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
1393 
1394  fillthisetaphi="hetaphiD0barcandidateS_";
1395  fillthisetaphi+=ptbin;
1396  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1397 
1398  if(TMath::Abs(minvD0bar-mPDG)<0.05){
1399  fillthisetaphi="hetaphiD0barcandidatesignalregionS_";
1400  fillthisetaphi+=ptbin;
1401  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1402  }
1403 
1404  }
1405  }
1406 
1407  //apply cut on invariant mass on the pair
1408  if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1409 
1410  cosThetaStarD0 = part->CosThetaStarD0();
1411  cosThetaStarD0bar = part->CosThetaStarD0bar();
1412  cosPointingAngle = part->CosPointingAngle();
1413  normalizedDecayLength2 = part->NormalizedDecayLength2();
1414  decayLength2 = part->DecayLength2();
1415  decayLengthxy = part->DecayLengthXY();
1416  normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1417 
1418  ptProng[0]=prong0->Pt(); ptProng[1]=prong1->Pt();
1419  d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
1420 
1421  if(fArray==1) cout<<"LS signal: ERROR"<<endl;
1422  for (Int_t iprong=0; iprong<2; iprong++){
1423  AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
1424  if (fReadMC) labprong[iprong]=prong->GetLabel();
1425 
1426  //cout<<"prong name = "<<prong->GetName()<<" label = "<<prong->GetLabel()<<endl;
1427  Int_t pdgprong=0;
1428  if(fReadMC && labprong[iprong]>=0) {
1429  mcprong= (AliAODMCParticle*)arrMC->At(labprong[iprong]);
1430  pdgprong=mcprong->GetPdgCode();
1431  }
1432 
1433  Bool_t isPionHere[2]={(isSelectedPID==1 || isSelectedPID==3) ? kTRUE : kFALSE,(isSelectedPID==2 || isSelectedPID==3) ? kTRUE : kFALSE};
1434 
1435  if(TMath::Abs(pdgprong)==211 || isPionHere[iprong]) {
1436  //cout<<"pi"<<endl;
1437 
1438  if(fSys==0){
1439  fillthispi="hptpiS_";
1440  fillthispi+=ptbin;
1441  ((TH1F*)listout->FindObject(fillthispi))->Fill(ptProng[iprong]);
1442  }
1443 
1444  fillthispi="hd0piS_";
1445  fillthispi+=ptbin;
1446  ((TH1F*)listout->FindObject(fillthispi))->Fill(d0Prong[iprong]);
1447  if(recalcvtx) {
1448 
1449  fillthispi="hd0vpiS_";
1450  fillthispi+=ptbin;
1451  ((TH1F*)listout->FindObject(fillthispi))->Fill(d0[iprong]);
1452  }
1453 
1454  }
1455 
1456  if(TMath::Abs(pdgprong)==321 || !isPionHere[iprong]) {
1457  //cout<<"kappa"<<endl;
1458  if(fSys==0){
1459  fillthisK="hptKS_";
1460  fillthisK+=ptbin;
1461  ((TH1F*)listout->FindObject(fillthisK))->Fill(ptProng[iprong]);
1462  }
1463 
1464 
1465  fillthisK="hd0KS_";
1466  fillthisK+=ptbin;
1467  ((TH1F*)listout->FindObject(fillthisK))->Fill(d0Prong[iprong]);
1468  if (recalcvtx){
1469  fillthisK="hd0vKS_";
1470  fillthisK+=ptbin;
1471  ((TH1F*)listout->FindObject(fillthisK))->Fill(d0[iprong]);
1472  }
1473  }
1474 
1475  if(fSys==0){
1476  fillthis="hcosthpointd0S_";
1477  fillthis+=ptbin;
1478  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[iprong]);
1479  }
1480  } //end loop on prongs
1481 
1482  fillthis="hdcaS_";
1483  fillthis+=ptbin;
1484  ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1485 
1486  fillthis="hcosthetapointS_";
1487  fillthis+=ptbin;
1488  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
1489 
1490  fillthis="hcosthetapointxyS_";
1491  fillthis+=ptbin;
1492  ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
1493 
1494 
1495  fillthis="hd0d0S_";
1496  fillthis+=ptbin;
1497  ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0());
1498 
1499  fillthis="hdeclS_";
1500  fillthis+=ptbin;
1501  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
1502 
1503  fillthis="hnormdeclS_";
1504  fillthis+=ptbin;
1505  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
1506 
1507  fillthis="hdeclxyS_";
1508  fillthis+=ptbin;
1509  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
1510 
1511  fillthis="hnormdeclxyS_";
1512  fillthis+=ptbin;
1513  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
1514 
1515  fillthis="hdeclxyd0d0S_";
1516  fillthis+=ptbin;
1517  ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
1518 
1519  fillthis="hnormdeclxyd0d0S_";
1520  fillthis+=ptbin;
1521  ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
1522 
1523  if(recalcvtx) {
1524  fillthis="hdeclvS_";
1525  fillthis+=ptbin;
1526  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1527 
1528  fillthis="hnormdeclvS_";
1529  fillthis+=ptbin;
1530  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1531 
1532  fillthis="hd0d0vS_";
1533  fillthis+=ptbin;
1534  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1535  }
1536 
1537  if(fSys==0){
1538  fillthis="hcosthetastarS_";
1539  fillthis+=ptbin;
1540  if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)) ((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1541  else {
1542  if (fReadMC || isSelectedPID>1)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
1543  if(isSelectedPID==1 || isSelectedPID==3)((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1544  }
1545 
1546  fillthis="hcosthpointd0d0S_";
1547  fillthis+=ptbin;
1548  ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1549  }
1550 
1551  if ((fReadMC && ((AliAODMCParticle*)arrMC->At(lab))->GetPdgCode() == 421)){
1552  for(Int_t it=0; it<2; it++){
1553  fillthis="hptD0S_";
1554  fillthis+=ptbin;
1555  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1556  fillthis="hphiD0S_";
1557  fillthis+=ptbin;
1558  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1559  Int_t nPointsITS = 0;
1560  for (Int_t il=0; il<6; il++){
1561  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
1562  }
1563  fillthis="hNITSpointsD0vsptS_";
1564  fillthis+=ptbin;
1565  ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(),nPointsITS);
1566  fillthis="hNSPDpointsD0S_";
1567  fillthis+=ptbin;
1568  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
1569  ((TH1I*)listout->FindObject(fillthis))->Fill(0);
1570  }
1571  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
1572  ((TH1I*)listout->FindObject(fillthis))->Fill(1);
1573  }
1574  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
1575  ((TH1I*)listout->FindObject(fillthis))->Fill(2);
1576  }
1577  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
1578  ((TH1I*)listout->FindObject(fillthis))->Fill(3);
1579  }
1580  fillthis="hNclsD0vsptS_";
1581  fillthis+=ptbin;
1582  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1583  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->GetTPCNcls();
1584  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1585  }
1586  }
1587  else {
1588  if (fReadMC || isSelectedPID>1){
1589  for(Int_t it=0; it<2; it++){
1590  fillthis="hptD0barS_";
1591  fillthis+=ptbin;
1592  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1593  fillthis="hphiD0barS_";
1594  fillthis+=ptbin;
1595  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1596  fillthis="hNclsD0barvsptS_";
1597  fillthis+=ptbin;
1598  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1599  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
1600  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1601  }
1602  }
1603  if(isSelectedPID==1 || isSelectedPID==3){
1604  for(Int_t it=0; it<2; it++){
1605  fillthis="hptD0S_";
1606  fillthis+=ptbin;
1607  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1608  fillthis="hphiD0S_";
1609  fillthis+=ptbin;
1610  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1611  Int_t nPointsITS = 0;
1612  for (Int_t il=0; il<6; il++){
1613  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
1614  }
1615  fillthis="hNITSpointsD0vsptS_";
1616  fillthis+=ptbin;
1617  ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(), nPointsITS);
1618  fillthis="hNSPDpointsD0S_";
1619  fillthis+=ptbin;
1620  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
1621  ((TH1I*)listout->FindObject(fillthis))->Fill(0);
1622  }
1623  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
1624  ((TH1I*)listout->FindObject(fillthis))->Fill(1);
1625  }
1626  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
1627  ((TH1I*)listout->FindObject(fillthis))->Fill(2);
1628  }
1629  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
1630  ((TH1I*)listout->FindObject(fillthis))->Fill(3);
1631  }
1632  fillthis="hNclsD0vsptS_";
1633  fillthis+=ptbin;
1634  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1635  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->GetTPCNcls();
1636  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1637  }
1638  }
1639  }
1640 
1641 
1642  } //end mass cut
1643 
1644  } else{ //Background or LS
1645  //if(!fReadMC){
1646  //cout<<"is background"<<endl;
1647 
1648  etaD = part->Eta();
1649  phiD = part->Phi();
1650 
1651  //no mass cut distributions: mass, ptbis
1652  fillthis="hMassB_";
1653  fillthis+=ptbin;
1654  fillthispt="histBkgPt";
1655 
1656  if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3))) {
1657  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0);
1658  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0,pt);
1659 
1660  fillthisetaphi="hetaphiD0candidateB_";
1661  fillthisetaphi+=ptbin;
1662  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1663 
1664  if(TMath::Abs(minvD0-mPDG)<0.05){
1665  fillthisetaphi="hetaphiD0candidatesignalregionB_";
1666  fillthisetaphi+=ptbin;
1667  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1668  }
1669  }
1670  if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) {
1671  ((TH1F*)listout->FindObject(fillthis))->Fill(minvD0bar);
1672  if(fFillPtHist && fReadMC) ((TH2F*)fOutputMassPt->FindObject(fillthispt))->Fill(minvD0bar,pt);
1673 
1674  fillthisetaphi="hetaphiD0barcandidateB_";
1675  fillthisetaphi+=ptbin;
1676  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1677 
1678  if(TMath::Abs(minvD0bar-mPDG)<0.05){
1679  fillthisetaphi="hetaphiD0barcandidatesignalregionB_";
1680  fillthisetaphi+=ptbin;
1681  ((TH2F*)listout->FindObject(fillthisetaphi))->Fill(etaD, phiD);
1682  }
1683 
1684  }
1685  if(fSys==0){
1686  fillthis="hptB1prongnoMcut_";
1687  fillthis+=ptbin;
1688 
1689  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1690 
1691  fillthis="hptB2prongsnoMcut_";
1692  fillthis+=ptbin;
1693  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1694  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt());
1695  }
1696 
1697 
1698  //apply cut on invariant mass on the pair
1699  if(TMath::Abs(minvD0-mPDG)<invmasscut || TMath::Abs(minvD0bar-mPDG)<invmasscut){
1700  if(fSys==0){
1701  ptProng[0]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt(); ptProng[1]=((AliAODTrack*)fDaughterTracks.UncheckedAt(0))->Pt();
1702  cosThetaStarD0 = part->CosThetaStarD0();
1703  cosThetaStarD0bar = part->CosThetaStarD0bar();
1704  }
1705 
1706  cosPointingAngle = part->CosPointingAngle();
1707  normalizedDecayLength2 = part->NormalizedDecayLength2();
1708  decayLength2 = part->DecayLength2();
1709  decayLengthxy = part->DecayLengthXY();
1710  normalizedDecayLengthxy=decayLengthxy/part->DecayLengthXYError();
1711  d0Prong[0]=part->Getd0Prong(0); d0Prong[1]=part->Getd0Prong(1);
1712 
1713 
1714  AliAODTrack *prongg=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
1715  if(!prongg) {
1716  if(fDebug>2) cout<<"No daughter found";
1717  return;
1718  }
1719  else{
1720  if(fArray==1){
1721  if(prongg->Charge()==1) {
1722  //fTotPosPairs[ptbin]++;
1723  ((TH1F*)fOutputMass->FindObject("hpospair"))->Fill(ptbin);
1724  } else {
1725  //fTotNegPairs[ptbin]++;
1726  ((TH1F*)fOutputMass->FindObject("hnegpair"))->Fill(ptbin);
1727  }
1728  }
1729  }
1730 
1731  //fill pt and phi distrib for prongs with M cut
1732 
1733  if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3))){
1734  for(Int_t it=0; it<2; it++){
1735  fillthis="hptD0B_";
1736  fillthis+=ptbin;
1737  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1738  fillthis="hphiD0B_";
1739  fillthis+=ptbin;
1740  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1741 
1742  Int_t nPointsITS = 0;
1743  for (Int_t il=0; il<6; il++){
1744  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(il)) nPointsITS++;
1745  }
1746  fillthis="hNITSpointsD0vsptB_";
1747  fillthis+=ptbin;
1748  ((TH2F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt(), nPointsITS);
1749  fillthis="hNSPDpointsD0B_";
1750  fillthis+=ptbin;
1751  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && !(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1))){ //no SPD points
1752  ((TH1I*)listout->FindObject(fillthis))->Fill(0);
1753  }
1754  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && !(((AliAODTrack*)(fDaughterTracks.UncheckedAt(it)))->HasPointOnITSLayer(1))){ //kOnlyFirst
1755  ((TH1I*)listout->FindObject(fillthis))->Fill(1);
1756  }
1757  if(!(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0)) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kOnlySecond
1758  ((TH1I*)listout->FindObject(fillthis))->Fill(2);
1759  }
1760  if(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(0) && ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->HasPointOnITSLayer(1)){ //kboth
1761  ((TH1I*)listout->FindObject(fillthis))->Fill(3);
1762  }
1763  fillthis="hNclsD0vsptB_";
1764  fillthis+=ptbin;
1765  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1766  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
1767  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1768  }
1769 
1770 
1771  }
1772 
1773  if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1)) {
1774  for(Int_t it=0; it<2; it++){
1775  fillthis="hptD0barB_";
1776  fillthis+=ptbin;
1777  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt());
1778  fillthis="hphiD0barB_";
1779  fillthis+=ptbin;
1780  ((TH1F*)listout->FindObject(fillthis))->Fill(((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Phi());
1781  fillthis="hNclsD0barvsptB_";
1782  fillthis+=ptbin;
1783  Float_t mom = ((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->Pt();
1784  Float_t ncls = (Float_t)((AliAODTrack*)fDaughterTracks.UncheckedAt(it))->GetTPCNcls();
1785  ((TH2F*)listout->FindObject(fillthis))->Fill(mom, ncls);
1786  }
1787  }
1788 
1789  fillthis="hd0B_";
1790  fillthis+=ptbin;
1791  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1792  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
1793 
1794  if(fReadMC){
1795  Int_t pdgMother[2]={0,0};
1796  Double_t factor[2]={1,1};
1797 
1798  for(Int_t iprong=0;iprong<2;iprong++){
1799  AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(iprong);
1800  lab=prong->GetLabel();
1801  if(lab>=0){
1802  AliAODMCParticle* mcprong=(AliAODMCParticle*)arrMC->At(lab);
1803  if(mcprong){
1804  Int_t labmom=mcprong->GetMother();
1805  if(labmom>=0){
1806  AliAODMCParticle* mcmother=(AliAODMCParticle*)arrMC->At(labmom);
1807  if(mcmother) pdgMother[iprong]=mcmother->GetPdgCode();
1808  }
1809  }
1810  }
1811 
1812  if(fSys==0){
1813 
1814  fillthis="hd0moresB_";
1815  fillthis+=ptbin;
1816 
1817  if(TMath::Abs(pdgMother[iprong])==310 || TMath::Abs(pdgMother[iprong])==130 || TMath::Abs(pdgMother[iprong])==321){ //K^0_S, K^0_L, K^+-
1818  if(ptProng[iprong]<=1)factor[iprong]=1./.7;
1819  else factor[iprong]=1./.6;
1820  fNentries->Fill(11);
1821  }
1822 
1823  if(TMath::Abs(pdgMother[iprong])==3122) { //Lambda
1824  factor[iprong]=1./0.25;
1825  fNentries->Fill(12);
1826  }
1827  fillthis="hd0moresB_";
1828  fillthis+=ptbin;
1829 
1830  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[iprong],factor[iprong]);
1831 
1832  if(recalcvtx){
1833  fillthis="hd0vmoresB_";
1834  fillthis+=ptbin;
1835  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[iprong],factor[iprong]);
1836  }
1837  }
1838  } //loop on prongs
1839 
1840  if(fSys==0){
1841  fillthis="hd0d0moresB_";
1842  fillthis+=ptbin;
1843  ((TH1F*)listout->FindObject(fillthis))->Fill(part->Prodd0d0(),factor[0]*factor[1]);
1844 
1845  fillthis="hcosthetapointmoresB_";
1846  fillthis+=ptbin;
1847  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,factor[0]*factor[1]);
1848 
1849  if(recalcvtx){
1850  fillthis="hd0d0vmoresB_";
1851  fillthis+=ptbin;
1852  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1],factor[0]*factor[1]);
1853  }
1854  }
1855  } //readMC
1856 
1857  if(fSys==0){
1858  //normalise pt distr to half afterwards
1859  fillthis="hptB_";
1860  fillthis+=ptbin;
1861  ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[0]);
1862  ((TH1F*)listout->FindObject(fillthis))->Fill(ptProng[1]);
1863 
1864  fillthis="hcosthetastarB_";
1865  fillthis+=ptbin;
1866  if (!fCutOnDistr || (fCutOnDistr && (fIsSelectedCandidate==1 || fIsSelectedCandidate==3)))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0);
1867  if (!fCutOnDistr || (fCutOnDistr && fIsSelectedCandidate>1))((TH1F*)listout->FindObject(fillthis))->Fill(cosThetaStarD0bar);
1868 
1869 
1870  fillthis="hd0p0B_";
1871  fillthis+=ptbin;
1872  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]);
1873  fillthis="hd0p1B_";
1874  fillthis+=ptbin;
1875  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[1]);
1876 
1877  fillthis="hcosthpointd0d0B_";
1878  fillthis+=ptbin;
1879  ((TH2F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,part->Prodd0d0());
1880 
1881  fillthis="hcosthpointd0B_";
1882  fillthis+=ptbin;
1883  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[0]);
1884  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle,d0Prong[1]);
1885 
1886 
1887  if(recalcvtx){
1888 
1889  fillthis="hd0vp0B_";
1890  fillthis+=ptbin;
1891  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1892  fillthis="hd0vp1B_";
1893  fillthis+=ptbin;
1894  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1895 
1896  fillthis="hd0vB_";
1897  fillthis+=ptbin;
1898  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]);
1899  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[1]);
1900 
1901  }
1902 
1903  }
1904 
1905  fillthis="hdcaB_";
1906  fillthis+=ptbin;
1907  ((TH1F*)listout->FindObject(fillthis))->Fill(part->GetDCA());
1908 
1909  fillthis="hd0d0B_";
1910  fillthis+=ptbin;
1911  ((TH1F*)listout->FindObject(fillthis))->Fill(d0Prong[0]*d0Prong[1]);
1912 
1913  if(recalcvtx){
1914  fillthis="hd0d0vB_";
1915  fillthis+=ptbin;
1916  ((TH1F*)listout->FindObject(fillthis))->Fill(d0[0]*d0[1]);
1917  }
1918 
1919  fillthis="hcosthetapointB_";
1920  fillthis+=ptbin;
1921  ((TH1F*)listout->FindObject(fillthis))->Fill(cosPointingAngle);
1922 
1923  fillthis="hcosthetapointxyB_";
1924  fillthis+=ptbin;
1925  ((TH1F*)listout->FindObject(fillthis))->Fill(part->CosPointingAngleXY());
1926 
1927  fillthis="hdeclB_";
1928  fillthis+=ptbin;
1929  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLength2);
1930 
1931  fillthis="hnormdeclB_";
1932  fillthis+=ptbin;
1933  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLength2);
1934 
1935  fillthis="hdeclxyB_";
1936  fillthis+=ptbin;
1937  ((TH1F*)listout->FindObject(fillthis))->Fill(decayLengthxy);
1938 
1939  fillthis="hnormdeclxyB_";
1940  fillthis+=ptbin;
1941  ((TH1F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy);
1942 
1943  fillthis="hdeclxyd0d0B_";
1944  fillthis+=ptbin;
1945  ((TH2F*)listout->FindObject(fillthis))->Fill(decayLengthxy,d0Prong[0]*d0Prong[1]);
1946 
1947  fillthis="hnormdeclxyd0d0B_";
1948  fillthis+=ptbin;
1949  ((TH2F*)listout->FindObject(fillthis))->Fill(normalizedDecayLengthxy,d0Prong[0]*d0Prong[1]);
1950 
1951 
1952  if(recalcvtx) {
1953 
1954  fillthis="hdeclvB_";
1955  fillthis+=ptbin;
1956  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[0]);
1957 
1958  fillthis="hnormdeclvB_";
1959  fillthis+=ptbin;
1960  ((TH1F*)listout->FindObject(fillthis))->Fill(decl[1]);
1961 
1962 
1963  }
1964  }//mass cut
1965  }//else (background)
1966 
1967  return;
1968 }
1969 
1970 //____________________________________________________________________________
1971 void AliAnalysisTaskSED0Mass::FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAODMCHeader *mcHeader, AliRDHFCutsD0toKpi* cuts, TList *listout){
1972  //
1974  //
1975 
1976 
1977  Double_t mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
1978 
1979  //cout<<"is selected = "<<fIsSelectedCandidate<<endl;
1980 
1981  // Fill candidate variable Tree (track selection, no candidate selection)
1982  if( fWriteVariableTree && !part->HasBadDaughters()
1983  && fCuts->AreDaughtersSelected(part) && fCuts->IsSelectedPID(part) ){
1984  fCandidateVariables[0] = part->InvMassD0();
1985  fCandidateVariables[1] = part->InvMassD0bar();
1986  fCandidateVariables[2] = part->Pt();
1987  fCandidateVariables[3] = part->GetDCA();
1988  Double_t ctsD0=0. ,ctsD0bar=0.; part->CosThetaStarD0(ctsD0,ctsD0bar);
1989  fCandidateVariables[4] = ctsD0;
1990  fCandidateVariables[5] = ctsD0bar;
1991  fCandidateVariables[6] = part->Pt2Prong(0);
1992  fCandidateVariables[7] = part->Pt2Prong(1);
1993  fCandidateVariables[8] = part->Getd0Prong(0);
1994  fCandidateVariables[9] = part->Getd0Prong(1);
1995  fCandidateVariables[10] = part->Prodd0d0();
1996  fCandidateVariables[11] = part->CosPointingAngle();
1997  fCandidateVariables[12] = part->CosPointingAngleXY();
1998  fCandidateVariables[13] = part->NormalizedDecayLengthXY();
1999  fCandidateVariables[14] = fCuts->IsSelectedSpecialCuts(part);
2000  fVariablesTree->Fill();
2001  }
2002 
2003  //cout<<"check cuts = "<<endl;
2004  //cuts->PrintAll();
2005  if (!fIsSelectedCandidate){
2006  //cout<<" cut " << cuts << " Rejected because "<<cuts->GetWhy()<<endl;
2007  return;
2008  }
2009 
2010 
2011  if(fDebug>2) cout<<"Candidate selected"<<endl;
2012 
2013  Double_t invmassD0 = part->InvMassD0(), invmassD0bar = part->InvMassD0bar();
2014  //printf("SELECTED\n");
2015  Int_t ptbin=cuts->PtBin(part->Pt());
2016  Double_t pt = part->Pt();
2017  Double_t y = part->YD0();
2018 
2019  Double_t impparXY=part->ImpParXY()*10000.;
2020  Double_t trueImpParXY=0.;
2021  Double_t arrayForSparse[3]={invmassD0,pt,impparXY};
2022  Double_t arrayForSparseTrue[3]={invmassD0,pt,trueImpParXY};
2023 
2024 
2025  // AliAODTrack *prong=(AliAODTrack*)fDaughterTracks.UncheckedAt(0);
2026  // if(!prong) {
2027  // AliDebug(2,"No daughter found");
2028  // return;
2029  // }
2030  // else{
2031  // if(prong->Charge()==1) {
2032  // ((TH1F*)fDistr->FindObject("hpospair"))->Fill(fCuts->GetNPtBins()+ptbin);
2033  // //fTotPosPairs[ptbin]++;
2034  // } else {
2035  // ((TH1F*)fDistr->FindObject("hnegpair"))->Fill(fCuts->GetNPtBins()+ptbin);
2036  // //fTotNegPairs[ptbin]++;
2037  // }
2038  // }
2039 
2040  // for(Int_t it=0;it<2;it++){
2041 
2042  // //request on spd points to be addes
2043  // if(/*nSPD==2 && */part->Pt() > 5. && (TMath::Abs(invmassD0-mPDG)<0.01 || TMath::Abs(invmassD0bar-mPDG)<0.01)){
2044  // FILE *f=fopen("4display.txt","a");
2045  // 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());
2046  // fclose(f);
2047  // //printf("PrimVtx NContributors: %d \n Prongs Rel Angle: %f \n \n",ncont,relangle);
2048  // }
2049  // }
2050 
2051  TString fillthis="", fillthispt="", fillthismasspt="", fillthismassy="";
2052  Int_t pdgDgD0toKpi[2]={321,211};
2053  Int_t labD0=-1;
2054  Bool_t isPrimary=kTRUE;
2055  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)
2056 
2057  //Define weights for Bayesian (if appropriate)
2058 
2059  Double_t weigD0=1.;
2060  Double_t weigD0bar=1.;
2061  if (fCuts->GetCombPID() && (fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeight || fCuts->GetBayesianStrategy() == AliRDHFCutsD0toKpi::kBayesWeightNoFilter)) {
2062  weigD0=fCuts->GetWeightsNegative()[AliPID::kKaon] * fCuts->GetWeightsPositive()[AliPID::kPion];
2063  weigD0bar=fCuts->GetWeightsPositive()[AliPID::kKaon] * fCuts->GetWeightsNegative()[AliPID::kPion];
2064  if (weigD0 > 1.0 || weigD0 < 0.) {weigD0 = 0.;}
2065  if (weigD0bar > 1.0 || weigD0bar < 0.) {weigD0bar = 0.;} //Prevents filling with weight > 1, or < 0
2066  }
2067 
2068  //count candidates selected by cuts
2069  fNentries->Fill(1);
2070  //count true D0 selected by cuts
2071  if (fReadMC && labD0>=0) fNentries->Fill(2);
2072 
2073  if ((fIsSelectedCandidate==1 || fIsSelectedCandidate==3) && fFillOnlyD0D0bar<2) { //D0
2074 
2075  arrayForSparse[0]=invmassD0; arrayForSparse[2]=impparXY;
2076 
2077  if(fReadMC){
2078  if(labD0>=0) {
2079  if(fArray==1) cout<<"LS signal ERROR"<<endl;
2080 
2081  AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
2082  Int_t pdgD0 = partD0->GetPdgCode();
2083  // cout<<"pdg = "<<pdgD0<<endl;
2084 
2085  if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE;
2086  if(!isPrimary)
2087  trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
2088  arrayForSparseTrue[0]=invmassD0; arrayForSparseTrue[2]=trueImpParXY;
2089 
2090  if (pdgD0==421){ //D0
2091  // cout<<"Fill S with D0"<<endl;
2092  fillthis="histSgn_";
2093  fillthis+=ptbin;
2094  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2095 
2096  if(fFillPtHist){
2097  fillthismasspt="histSgnPt";
2098  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2099  }
2100  if(fFillImpParHist){
2101  if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0);
2102  else {
2103  fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0);
2104  fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0);
2105  }
2106  }
2107 
2108  if(fFillYHist){
2109  fillthismassy="histSgnY_";
2110  fillthismassy+=ptbin;
2111  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2112  }
2113 
2114  if(fSys==0){
2115  if(TMath::Abs(invmassD0 - mPDG) < 0.027 && fFillVarHists){
2116  fillthis="histSgn27_";
2117  fillthis+=ptbin;
2118  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2119  }
2120  }
2121  } else{ //it was a D0bar
2122  fillthis="histRfl_";
2123  fillthis+=ptbin;
2124  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2125 
2126  if(fFillPtHist){
2127  fillthismasspt="histRflPt";
2128  // cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
2129  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2130  }
2131 
2132  if(fFillYHist){
2133  fillthismassy="histRflY_";
2134  fillthismassy+=ptbin;
2135  // cout << " Filling "<<fillthismassy<<" D0bar"<<endl;
2136  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2137  }
2138 
2139  }
2140  } else {//background
2141  fillthis="histBkg_";
2142  fillthis+=ptbin;
2143  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2144 
2145  if(fFillPtHist){
2146  fillthismasspt="histBkgPt";
2147  // cout << " Filling "<<fillthismasspt<<" D0bar"<<endl;
2148  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2149  }
2150  if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0);
2151 
2152  if(fFillYHist){
2153  fillthismassy="histBkgY_";
2154  fillthismassy+=ptbin;
2155  // cout << " Filling "<<fillthismassy<<" D0bar"<<endl;
2156  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2157  }
2158 
2159  }
2160 
2161  }else{
2162  fillthis="histMass_";
2163  fillthis+=ptbin;
2164  // cout<<"Filling "<<fillthis<<endl;
2165 
2166  // printf("Fill mass with D0");
2167  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0,weigD0);
2168 
2169 
2170  if(fFillPtHist){
2171  fillthismasspt="histMassPt";
2172  // cout<<"Filling "<<fillthismasspt<<endl;
2173  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0,pt,weigD0);
2174  }
2175  if(fFillImpParHist) {
2176  // cout << "Filling fHistMassPtImpParTC[0]"<<endl;
2177  fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0);
2178  }
2179 
2180  if(fFillYHist){
2181  fillthismassy="histMassY_";
2182  fillthismassy+=ptbin;
2183  // cout<<"Filling "<<fillthismassy<<endl;
2184  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0,y,weigD0);
2185  }
2186 
2187  }
2188 
2189  }
2190  if (fIsSelectedCandidate>1 && (fFillOnlyD0D0bar==0 || fFillOnlyD0D0bar==2)) { //D0bar
2191 
2192  arrayForSparse[0]=invmassD0bar; arrayForSparse[2]=impparXY;
2193 
2194  if(fReadMC){
2195  if(labD0>=0) {
2196  if(fArray==1) cout<<"LS signal ERROR"<<endl;
2197  AliAODMCParticle *partD0 = (AliAODMCParticle*)arrMC->At(labD0);
2198  Int_t pdgD0 = partD0->GetPdgCode();
2199  // cout<<" pdg = "<<pdgD0<<endl;
2200 
2201  if(CheckOrigin(arrMC,partD0)==5) isPrimary=kFALSE;
2202  if(!isPrimary)
2203  trueImpParXY=GetTrueImpactParameter(mcHeader,arrMC,partD0)*10000.;
2204  arrayForSparseTrue[0]=invmassD0bar; arrayForSparseTrue[2]=trueImpParXY;
2205 
2206  if (pdgD0==-421){ //D0bar
2207  fillthis="histSgn_";
2208  fillthis+=ptbin;
2209  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
2210  // if (TMath::Abs(invmassD0bar - mPDG) < 0.027){
2211  // fillthis="histSgn27_";
2212  // fillthis+=ptbin;
2213  // ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar);
2214  // }
2215 
2216  if(fFillPtHist){
2217  fillthismasspt="histSgnPt";
2218  // cout<<" Filling "<< fillthismasspt << endl;
2219  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2220  }
2221  if(fFillImpParHist){
2222  // cout << " Filling impact parameter thnsparse"<<endl;
2223  if(isPrimary) fHistMassPtImpParTC[1]->Fill(arrayForSparse,weigD0bar);
2224  else {
2225  fHistMassPtImpParTC[2]->Fill(arrayForSparse,weigD0bar);
2226  fHistMassPtImpParTC[3]->Fill(arrayForSparseTrue,weigD0bar);
2227  }
2228  }
2229 
2230  if(fFillYHist){
2231  fillthismassy="histSgnY_";
2232  fillthismassy+=ptbin;
2233  // cout<<" Filling "<< fillthismassy << endl;
2234  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2235  }
2236 
2237  } else{
2238  fillthis="histRfl_";
2239  fillthis+=ptbin;
2240  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
2241  if(fFillPtHist){
2242  fillthismasspt="histRflPt";
2243  // cout << " Filling "<<fillthismasspt<<endl;
2244  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2245  }
2246  if(fFillYHist){
2247  fillthismassy="histRflY_";
2248  fillthismassy+=ptbin;
2249  // cout << " Filling "<<fillthismassy<<endl;
2250  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2251  }
2252  }
2253  } else {//background or LS
2254  fillthis="histBkg_";
2255  fillthis+=ptbin;
2256  ((TH1F*)(listout->FindObject(fillthis)))->Fill(invmassD0bar,weigD0bar);
2257 
2258  if(fFillPtHist){
2259  fillthismasspt="histBkgPt";
2260  // cout<<" Filling "<< fillthismasspt << endl;
2261  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2262  }
2263  if(fFillImpParHist) fHistMassPtImpParTC[4]->Fill(arrayForSparse,weigD0bar);
2264  if(fFillYHist){
2265  fillthismassy="histBkgY_";
2266  fillthismassy+=ptbin;
2267  // cout<<" Filling "<< fillthismassy << endl;
2268  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2269  }
2270  }
2271  }else{
2272  fillthis="histMass_";
2273  fillthis+=ptbin;
2274  // printf("Fill mass with D0bar");
2275 
2276  ((TH1F*)listout->FindObject(fillthis))->Fill(invmassD0bar,weigD0bar);
2277 
2278 
2279  if(fFillPtHist){
2280  fillthismasspt="histMassPt";
2281  // cout<<" Filling "<< fillthismasspt << endl;
2282  ((TH2F*)(fOutputMassPt->FindObject(fillthismasspt)))->Fill(invmassD0bar,pt,weigD0bar);
2283  }
2284  if(fFillImpParHist) fHistMassPtImpParTC[0]->Fill(arrayForSparse,weigD0bar);
2285  if(fFillYHist){
2286  fillthismassy="histMassY_";
2287  fillthismassy+=ptbin;
2288  // cout<<" Filling "<< fillthismassy << endl;
2289  ((TH2F*)(fOutputMassY->FindObject(fillthismassy)))->Fill(invmassD0bar,y,weigD0bar);
2290  }
2291  }
2292  }
2293 
2294  return;
2295 }
2296 
2297 //__________________________________________________________________________
2298 AliAODVertex* AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped(AliAODEvent *aodev){
2300 
2301  Int_t skipped[2];
2302  Int_t nTrksToSkip=2;
2303  AliAODTrack *dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(0);
2304  if(!dgTrack){
2305  AliDebug(2,"no daughter found!");
2306  return 0x0;
2307  }
2308  skipped[0]=dgTrack->GetID();
2309  dgTrack = (AliAODTrack*)fDaughterTracks.UncheckedAt(1);
2310  if(!dgTrack){
2311  AliDebug(2,"no daughter found!");
2312  return 0x0;
2313  }
2314  skipped[1]=dgTrack->GetID();
2315 
2316  AliESDVertex *vertexESD=0x0;
2317  AliAODVertex *vertexAOD=0x0;
2318  AliVertexerTracks *vertexer = new AliVertexerTracks(aodev->GetMagneticField());
2319 
2320  //
2321  vertexer->SetSkipTracks(nTrksToSkip,skipped);
2322  vertexer->SetMinClusters(4);
2323  vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(aodev);
2324  if(!vertexESD) return vertexAOD;
2325  if(vertexESD->GetNContributors()<=0) {
2326  AliDebug(2,"vertexing failed");
2327  delete vertexESD; vertexESD=NULL;
2328  return vertexAOD;
2329  }
2330 
2331  delete vertexer; vertexer=NULL;
2332 
2333 
2334  // convert to AliAODVertex
2335  Double_t pos[3],cov[6],chi2perNDF;
2336  vertexESD->GetXYZ(pos); // position
2337  vertexESD->GetCovMatrix(cov); //covariance matrix
2338  chi2perNDF = vertexESD->GetChi2toNDF();
2339  delete vertexESD; vertexESD=NULL;
2340 
2341  vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
2342  return vertexAOD;
2343 
2344 }
2345 
2346 
2347 //________________________________________________________________________
2348 void AliAnalysisTaskSED0Mass::Terminate(Option_t */*option*/)
2349 {
2351  //
2352  if(fDebug > 1) printf("AnalysisTaskSED0Mass: Terminate() \n");
2353 
2354 
2355  fOutputMass = dynamic_cast<TList*> (GetOutputData(1));
2356  if (!fOutputMass) {
2357  printf("ERROR: fOutputMass not available\n");
2358  return;
2359  }
2360  fOutputMassPt = dynamic_cast<TList*> (GetOutputData(6));
2361  if ((fFillPtHist || fFillImpParHist) && !fOutputMassPt) {
2362  printf("ERROR: fOutputMass not available\n");
2363  return;
2364  }
2365 
2366  if(fFillVarHists){
2367  fDistr = dynamic_cast<TList*> (GetOutputData(2));
2368  if (!fDistr) {
2369  printf("ERROR: fDistr not available\n");
2370  return;
2371  }
2372  }
2373 
2374  fNentries = dynamic_cast<TH1F*>(GetOutputData(3));
2375 
2376  if(!fNentries){
2377  printf("ERROR: fNEntries not available\n");
2378  return;
2379  }
2380  fCuts = dynamic_cast<AliRDHFCutsD0toKpi*>(GetOutputData(4));
2381  if(!fCuts){
2382  printf("ERROR: fCuts not available\n");
2383  return;
2384  }
2385  fCounter = dynamic_cast<AliNormalizationCounter*>(GetOutputData(5));
2386  if (!fCounter) {
2387  printf("ERROR: fCounter not available\n");
2388  return;
2389  }
2390  if (fDrawDetSignal) {
2391  fDetSignal = dynamic_cast<TList*>(GetOutputData(8));
2392  if (!fDetSignal) {
2393  printf("ERROR: fDetSignal not available\n");
2394  return;
2395  }
2396  }
2397  if(fFillYHist){
2398  fOutputMassY = dynamic_cast<TList*> (GetOutputData(9));
2399  if (fFillYHist && !fOutputMassY) {
2400  printf("ERROR: fOutputMassY not available\n");
2401  return;
2402  }
2403  }
2404 
2405  Int_t nptbins=fCuts->GetNPtBins();
2406  for(Int_t ipt=0;ipt<nptbins;ipt++){
2407 
2408  if(fArray==1 && fFillVarHists){
2409  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
2410 
2411 
2412  if(fLsNormalization>1e-6) {
2413 
2414  TString massName="histMass_";
2415  massName+=ipt;
2416  ((TH1F*)fOutputMass->FindObject(massName))->Scale((1/fLsNormalization)*((TH1F*)fOutputMass->FindObject(massName))->GetEntries());
2417 
2418  }
2419 
2420 
2421  fLsNormalization = 2.*TMath::Sqrt(((TH1F*)fOutputMass->FindObject("hpospair"))->Integral(ipt+1,ipt+2)*((TH1F*)fOutputMass->FindObject("hnegpair"))->Integral(ipt+1,ipt+2));
2422  //fLsNormalization = 2.*TMath::Sqrt(fTotPosPairs[4]*fTotNegPairs[4]);
2423 
2424  if(fLsNormalization>1e-6) {
2425 
2426  TString nameDistr="hdcaB_";
2427  nameDistr+=ipt;
2428  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2429  nameDistr="hd0B_";
2430  nameDistr+=ipt;
2431  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2432  nameDistr="hd0d0B_";
2433  nameDistr+=ipt;
2434  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2435  nameDistr="hcosthetapointB_";
2436  nameDistr+=ipt;
2437  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2438  if(fSys==0){
2439  nameDistr="hptB_";
2440  nameDistr+=ipt;
2441  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2442  nameDistr="hcosthetastarB_";
2443  nameDistr+=ipt;
2444  ((TH1F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH1F*)fDistr->FindObject(nameDistr))->GetEntries());
2445  nameDistr="hcosthpointd0d0B_";
2446  nameDistr+=ipt;
2447  ((TH2F*)fDistr->FindObject(nameDistr))->Scale((1/fLsNormalization)*((TH2F*)fDistr->FindObject(nameDistr))->GetEntries());
2448  }
2449  }
2450  }
2451  }
2452  TString cvname,cstname;
2453 
2454  if (fArray==0){
2455  cvname="D0invmass";
2456  cstname="cstat0";
2457  } else {
2458  cvname="LSinvmass";
2459  cstname="cstat1";
2460  }
2461 
2462  TCanvas *cMass=new TCanvas(cvname,cvname);
2463  cMass->cd();
2464  ((TH1F*)fOutputMass->FindObject("histMass_3"))->Draw();
2465 
2466  TCanvas* cStat=new TCanvas(cstname,Form("Stat%s",fArray ? "LS" : "D0"));
2467  cStat->cd();
2468  cStat->SetGridy();
2469  fNentries->Draw("htext0");
2470 
2471  // TCanvas *ccheck=new TCanvas(Form("cc%d",fArray),Form("cc%d",fArray));
2472  // ccheck->cd();
2473 
2474  return;
2475 }
2476 
2477 
2478 //________________________________________________________________________
2479 void AliAnalysisTaskSED0Mass::CreateImpactParameterHistos(){
2481 
2482  Int_t nmassbins=200;
2483  Double_t fLowmasslimit=1.5648, fUpmasslimit=2.1648;
2484  Int_t fNImpParBins=400;
2485  Double_t fLowerImpPar=-2000., fHigherImpPar=2000.;
2486  Int_t nbins[3]={nmassbins,200,fNImpParBins};
2487  Double_t xmin[3]={fLowmasslimit,0.,fLowerImpPar};
2488  Double_t xmax[3]={fUpmasslimit,20.,fHigherImpPar};
2489 
2490 
2491  fHistMassPtImpParTC[0]=new THnSparseF("hMassPtImpParAll",
2492  "Mass vs. pt vs.imppar - All",
2493  3,nbins,xmin,xmax);
2494  fHistMassPtImpParTC[1]=new THnSparseF("hMassPtImpParPrompt",
2495  "Mass vs. pt vs.imppar - promptD",
2496  3,nbins,xmin,xmax);
2497  fHistMassPtImpParTC[2]=new THnSparseF("hMassPtImpParBfeed",
2498  "Mass vs. pt vs.imppar - DfromB",
2499  3,nbins,xmin,xmax);
2500  fHistMassPtImpParTC[3]=new THnSparseF("hMassPtImpParTrueBfeed",
2501  "Mass vs. pt vs.true imppar -DfromB",
2502  3,nbins,xmin,xmax);
2503  fHistMassPtImpParTC[4]=new THnSparseF("hMassPtImpParBkg",
2504  "Mass vs. pt vs.imppar - backgr.",
2505  3,nbins,xmin,xmax);
2506 
2507  for(Int_t i=0; i<5;i++){
2508  fOutputMassPt->Add(fHistMassPtImpParTC[i]);
2509  }
2510 }
2511 
2512 //_________________________________________________________________________________________________
2513 Float_t AliAnalysisTaskSED0Mass::GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray* arrayMC, AliAODMCParticle *partD0) const {
2515 
2516  printf(" AliAnalysisTaskSED0MassV1::GetTrueImpactParameter() \n");
2517 
2518  Double_t vtxTrue[3];
2519  mcHeader->GetVertex(vtxTrue);
2520  Double_t origD[3];
2521  partD0->XvYvZv(origD);
2522  Short_t charge=partD0->Charge();
2523  Double_t pXdauTrue[2],pYdauTrue[2],pZdauTrue[2];
2524  for(Int_t iDau=0; iDau<2; iDau++){
2525  pXdauTrue[iDau]=0.;
2526  pYdauTrue[iDau]=0.;
2527  pZdauTrue[iDau]=0.;
2528  }
2529 
2530  // Int_t nDau=partD0->GetNDaughters();
2531  Int_t labelFirstDau = partD0->GetDaughter(0);
2532 
2533  for(Int_t iDau=0; iDau<2; iDau++){
2534  Int_t ind = labelFirstDau+iDau;
2535  AliAODMCParticle* part = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ind));
2536  if(!part) continue;
2537  Int_t pdgCode = TMath::Abs( part->GetPdgCode() );
2538  if(!part){
2539  AliError("Daughter particle not found in MC array");
2540  return 99999.;
2541  }
2542  if(pdgCode==211 || pdgCode==321){
2543  pXdauTrue[iDau]=part->Px();
2544  pYdauTrue[iDau]=part->Py();
2545  pZdauTrue[iDau]=part->Pz();
2546  }
2547  }
2548 
2549  Double_t d0dummy[2]={0.,0.};
2550  AliAODRecoDecayHF aodDzeroMC(vtxTrue,origD,2,charge,pXdauTrue,pYdauTrue,pZdauTrue,d0dummy);
2551  return aodDzeroMC.ImpParXY();
2552 
2553 }
2554 
2555 //_________________________________________________________________________________________________
2556 Int_t AliAnalysisTaskSED0Mass::CheckOrigin(TClonesArray* arrayMC, AliAODMCParticle *mcPartCandidate) const {
2557  //
2559  //
2560  printf(" AliAnalysisTaskSED0Mass V1::CheckOrigin() \n");
2561 
2562  Int_t pdgGranma = 0;
2563  Int_t mother = 0;
2564  mother = mcPartCandidate->GetMother();
2565  Int_t istep = 0;
2566  Int_t abspdgGranma =0;
2567  Bool_t isFromB=kFALSE;
2568  Bool_t isQuarkFound=kFALSE;
2569  while (mother >0 ){
2570  istep++;
2571  AliAODMCParticle* mcGranma = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mother));
2572  if (mcGranma){
2573  pdgGranma = mcGranma->GetPdgCode();
2574  abspdgGranma = TMath::Abs(pdgGranma);
2575  if ((abspdgGranma > 500 && abspdgGranma < 600) || (abspdgGranma > 5000 && abspdgGranma < 6000)){
2576  isFromB=kTRUE;
2577  }
2578  if(abspdgGranma==4 || abspdgGranma==5) isQuarkFound=kTRUE;
2579  mother = mcGranma->GetMother();
2580  }else{
2581  AliError("Failed casting the mother particle!");
2582  break;
2583  }
2584  }
2585 
2586  if(isFromB) return 5;
2587  else return 4;
2588 }
charge
Int_t charge
Definition: MuonTrackingEfficiency.C:54
AliAODRecoDecayHF::NormalizedDecayLengthXY
Double_t NormalizedDecayLengthXY() const
Definition: AliAODRecoDecayHF.h:75
ClassImp
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Definition: AliAnalysisTaskTriggerRates.cxx:34
AliAnalysisTaskSED0Mass::fOutputMassPt
TList * fOutputMassPt
! list send on output slot 6
Definition: AliAnalysisTaskSED0Mass.h:94
AliAnalysisTaskSED0Mass::fDrawDetSignal
Bool_t fDrawDetSignal
flag to decide whether to fill "PID = x" bins in fNentrie
Definition: AliAnalysisTaskSED0Mass.h:122
AliAnalysisTaskSED0Mass::fFillYHist
Bool_t fFillYHist
flag to fill Pt and Impact Parameter Histograms
Definition: AliAnalysisTaskSED0Mass.h:114
AliNormalizationCounter::StoreCandidates
void StoreCandidates(AliVEvent *, Int_t nCand=0, Bool_t flagFilter=kTRUE)
Definition: AliNormalizationCounter.cxx:261
AliAnalysisTaskSED0Mass::FillVarHists
void FillVarHists(AliAODEvent *aodev, AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliRDHFCutsD0toKpi *cuts, TList *listout)
Definition: AliAnalysisTaskSED0Mass.cxx:1226
AliRDHFCutsD0toKpi::GetCombPID
Bool_t GetCombPID() const
Definition: AliRDHFCutsD0toKpi.h:103
AliAnalysisTaskSED0Mass::fCandidateVariables
Double_t * fCandidateVariables
! variables to be written to the tree
Definition: AliAnalysisTaskSED0Mass.h:120
AliAODRecoCascadeHF.h
AliAODRecoDecayHF::HasSelectionBit
Bool_t HasSelectionBit(Int_t i) const
Definition: AliAODRecoDecayHF.h:123
AliAnalysisVertexingHF.h
AliAODRecoDecayHF::NormalizedDecayLength2
Double_t NormalizedDecayLength2() const
Definition: AliAODRecoDecayHF.h:69
AliAODRecoDecayHF::ImpParXY
Double_t ImpParXY() const
Definition: AliAODRecoDecayHF.h:83
AliAnalysisTaskSED0Mass::fWriteVariableTree
Bool_t fWriteVariableTree
flag to check or not the selection bit
Definition: AliAnalysisTaskSED0Mass.h:118
AliAnalysisTaskSED0Mass::fDistr
TList * fDistr
! list send on output slot 2
Definition: AliAnalysisTaskSED0Mass.h:96
AliRDHFCuts::GetWhyRejection
Int_t GetWhyRejection() const
Definition: AliRDHFCuts.h:294
AliAnalysisTaskSED0Mass::fDetSignal
TList * fDetSignal
flag to decide whether to draw the TPC dE/dx and TOF signal before/after PID
Definition: AliAnalysisTaskSED0Mass.h:123
AliAODRecoDecayHF::CosPointingAngleXY
Double_t CosPointingAngleXY() const
Definition: AliAODRecoDecayHF.h:81
AliAnalysisTaskSED0Mass::fUseSelectionBit
Bool_t fUseSelectionBit
flag to fill Pt and Impact Parameter Histograms
Definition: AliAnalysisTaskSED0Mass.h:116
AliAnalysisVertexingHF::FillRecoCand
Bool_t FillRecoCand(AliVEvent *event, AliAODRecoDecayHF3Prong *rd3)
Definition: AliAnalysisVertexingHF.cxx:1669
AliAnalysisTaskSED0Mass::fSys
Int_t fSys
flag to enable filling variable histos
Definition: AliAnalysisTaskSED0Mass.h:111
AliAODRecoDecayHF::GetSelectionMap
ULong_t GetSelectionMap() const
Definition: AliAODRecoDecayHF.h:124
AliAnalysisTaskSED0Mass::fOutputMass
TList * fOutputMass
! list send on output slot 1
Definition: AliAnalysisTaskSED0Mass.h:93
AliRDHFCutsD0toKpi::IsInFiducialAcceptance
virtual Bool_t IsInFiducialAcceptance(Double_t pt, Double_t y) const
Definition: AliRDHFCutsD0toKpi.cxx:701
AliAnalysisTaskSED0Mass::fReadMC
Bool_t fReadMC
can be D0 or Like Sign candidates
Definition: AliAnalysisTaskSED0Mass.h:101
AliRDHFCutsD0toKpi::GetBayesianStrategy
Int_t GetBayesianStrategy() const
Definition: AliRDHFCutsD0toKpi.h:84
AliAODRecoDecayHF::HasBadDaughters
Bool_t HasBadDaughters() const
Definition: AliAODRecoDecayHF.h:241
AliAnalysisTaskSED0Mass::fHistMassPtImpParTC
THnSparseF * fHistMassPtImpParTC[5]
! histograms for impact paramter studies
Definition: AliAnalysisTaskSED0Mass.h:99
AliAnalysisTaskSED0Mass::fNentries
TH1F * fNentries
! histogram with number of events on output slot 3
Definition: AliAnalysisTaskSED0Mass.h:97
AliAnalysisTaskSED0Mass::UserCreateOutputObjects
virtual void UserCreateOutputObjects()
Implementation of interface methods.
Definition: AliAnalysisTaskSED0Mass.cxx:226
AliAnalysisTaskSED0Mass::fFillPtHist
Bool_t fFillPtHist
flag to reject events with SDD clusters
Definition: AliAnalysisTaskSED0Mass.h:113
AliAnalysisTaskSED0Mass::fFillImpParHist
Bool_t fFillImpParHist
flag to fill Y Histograms
Definition: AliAnalysisTaskSED0Mass.h:115
AliAODRecoDecayHF2Prong::CosThetaStarD0bar
Double_t CosThetaStarD0bar() const
angle of K
Definition: AliAODRecoDecayHF2Prong.h:44
AliAnalysisTaskSED0Mass::fIsRejectSDDClusters
Bool_t fIsRejectSDDClusters
fSys=0 -> p-p; fSys=1 ->PbPb (in this case fFillVarHists=kFALSE by default: set it to kTRUE after if ...
Definition: AliAnalysisTaskSED0Mass.h:112
AliRDHFCuts::AreDaughtersSelected
Bool_t AreDaughtersSelected(AliAODRecoDecayHF *rd) const
Definition: AliRDHFCuts.cxx:632
AliAnalysisTaskSED0Mass::CheckOrigin
Int_t CheckOrigin(TClonesArray *arrayMC, AliAODMCParticle *mcPartCandidate) const
Definition: AliAnalysisTaskSED0Mass.cxx:2556
AliAnalysisTaskSED0Mass::Terminate
virtual void Terminate(Option_t *option)
Definition: AliAnalysisTaskSED0Mass.cxx:2348
AliAnalysisTaskSED0Mass::GetTrueImpactParameter
Float_t GetTrueImpactParameter(AliAODMCHeader *mcHeader, TClonesArray *arrayMC, AliAODMCParticle *partD0) const
Definition: AliAnalysisTaskSED0Mass.cxx:2513
AliAnalysisTaskSED0Mass::DrawDetSignal
void DrawDetSignal(AliAODRecoDecayHF2Prong *part, TList *ListDetSignal)
Definition: AliAnalysisTaskSED0Mass.cxx:1188
AliAODRecoDecayHF2Prong.h
AliAnalysisTaskSED0Mass::fLsNormalization
Double_t fLsNormalization
number of pt bins
Definition: AliAnalysisTaskSED0Mass.h:106
AliAnalysisTaskSED0Mass::fUsePid4Distr
Bool_t fUsePid4Distr
flag to decide if apply cut also on distributions: 0 no cuts, 1 looser cuts, 2 tighter/ cuts ...
Definition: AliAnalysisTaskSED0Mass.h:103
AliAODRecoDecayHF::DecayLength2
Double_t DecayLength2() const
kinematics & topology
Definition: AliAODRecoDecayHF.h:61
Draw
void Draw(Int_t cell[], Int_t iBC, Int_t nBC, const Int_t cellref=151)
Definition: BadChannelAnalysis.C:53
Data
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)
Definition: charmCutsOptimization.C:371
AliAnalysisTaskSED0Mass::fCutOnDistr
Bool_t fCutOnDistr
flag for MC array: kTRUE = read it, kFALSE = do not read it
Definition: AliAnalysisTaskSED0Mass.h:102
AliAnalysisTaskSED0Mass::GetPrimaryVtxSkipped
AliAODVertex * GetPrimaryVtxSkipped(AliAODEvent *aodev)
Definition: AliAnalysisTaskSED0Mass.cxx:2298
AliRDHFCutsD0toKpi::IsSelectedSpecialCuts
Int_t IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
Definition: AliRDHFCutsD0toKpi.cxx:1136
AliAODRecoDecayHF::SetOwnPrimaryVtx
void SetOwnPrimaryVtx(const AliAODVertex *vtx)
Definition: AliAODRecoDecayHF.h:45
AliAODRecoDecayHF::DecayLengthXY
Double_t DecayLengthXY() const
Definition: AliAODRecoDecayHF.h:71
AliAnalysisTaskSED0Mass::FillMassHists
void FillMassHists(AliAODRecoDecayHF2Prong *part, TClonesArray *arrMC, AliAODMCHeader *mcHeader, AliRDHFCutsD0toKpi *cuts, TList *listout)
Definition: AliAnalysisTaskSED0Mass.cxx:1971
AliAnalysisTaskSED0Mass::fIsSelectedCandidate
Int_t fIsSelectedCandidate
keeps the daughter tracks
Definition: AliAnalysisTaskSED0Mass.h:109
AliRDHFCuts::GetIsPrimaryWithoutDaughters
Bool_t GetIsPrimaryWithoutDaughters() const
Definition: AliRDHFCuts.h:262
AliRDHFCutsD0toKpi::GetWeightsPositive
Double_t * GetWeightsPositive() const
Definition: AliRDHFCutsD0toKpi.h:76
AliRDHFCuts::IsEventSelected
Bool_t IsEventSelected(AliVEvent *event)
Definition: AliRDHFCuts.cxx:451
AliNormalizationCounter::StoreEvent
void StoreEvent(AliVEvent *, AliRDHFCuts *, Bool_t mc=kFALSE, Int_t multiplicity=-9999, Double_t spherocity=-99.)
Definition: AliNormalizationCounter.cxx:163
AliRDHFCutsD0toKpi::IsSelectedPID
virtual Int_t IsSelectedPID(AliAODRecoDecayHF *rd)
Definition: AliRDHFCutsD0toKpi.cxx:731
AliNormalizationCounter.h
AliAnalysisTaskSED0Mass::fOutputMassY
TList * fOutputMassY
! list send on output slot 9
Definition: AliAnalysisTaskSED0Mass.h:95
AliAnalysisTaskSED0Mass.h
AliAnalysisTaskSED0Mass::fVariablesTree
TTree * fVariablesTree
flag to decide whether to write the candidate variables on a tree variables
Definition: AliAnalysisTaskSED0Mass.h:119
AliRDHFCutsD0toKpi::IsSelected
virtual Int_t IsSelected(TObject *obj, Int_t selectionLevel)
Definition: AliRDHFCutsD0toKpi.h:39
AliRDHFCuts::GetNPtBins
Int_t GetNPtBins() const
Definition: AliRDHFCuts.h:233
AliAnalysisTaskSED0Mass::UserExec
virtual void UserExec(Option_t *option)
Definition: AliAnalysisTaskSED0Mass.cxx:979
AliAnalysisTaskSED0Mass::fDaughterTracks
TObjArray fDaughterTracks
flag to fill mass histogram with D0/D0bar only (0 = fill with both, 1 = fill with D0 only...
Definition: AliAnalysisTaskSED0Mass.h:108
nbins
const Int_t nbins
Definition: AverageDmesonRaa.C:63
AliAODRecoDecayHF::CosPointingAngle
Double_t CosPointingAngle() const
Definition: AliAODRecoDecayHF.h:79
AliAODRecoDecayHF::DecayLengthXYError
Double_t DecayLengthXYError() const
Definition: AliAODRecoDecayHF.h:73
AliAnalysisTaskSED0Mass::fFillOnlyD0D0bar
Int_t fFillOnlyD0D0bar
normalization
Definition: AliAnalysisTaskSED0Mass.h:107
AliRDHFCuts::PtBin
Int_t PtBin(Double_t pt) const
Definition: AliRDHFCuts.cxx:1140
nptbins
Int_t nptbins
Definition: CalculateAveragePt.C:56
AliAnalysisTaskSED0Mass::fFillVarHists
Bool_t fFillVarHists
selection outcome
Definition: AliAnalysisTaskSED0Mass.h:110
AliAnalysisTaskSED0Mass::fCounter
AliNormalizationCounter * fCounter
! AliNormalizationCounter on output slot 5
Definition: AliAnalysisTaskSED0Mass.h:104
AliRDHFCutsD0toKpi::GetWeightsNegative
Double_t * GetWeightsNegative() const
Definition: AliRDHFCutsD0toKpi.h:82