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