laura is hosted by Hepforge, IPPP Durham
Laura++  v2r1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauEmbeddedData.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2007 - 2013.
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 
6 // Authors:
7 // Thomas Latham
8 // John Back
9 // Paul Harrison
10 
15 #include <iostream>
16 #include <vector>
17 using std::cerr;
18 using std::endl;
19 
20 #include "TRandom.h"
21 
22 #include "LauAbsDPDynamics.hh"
23 #include "LauEmbeddedData.hh"
24 #include "LauKinematics.hh"
25 #include "LauRandom.hh"
26 
28 
29 
30 LauEmbeddedData::LauEmbeddedData(const TString& fileName, const TString& treeName, Bool_t allowReuseOfEvents) :
31  theDataTree_(new LauFitDataTree(fileName,treeName)),
32  reuseEvents_(allowReuseOfEvents)
33 {
34 }
35 
37 {
38  delete theDataTree_;
39 }
40 
42 {
43  if ( theDataTree_ == 0 ) {
44  cerr<<"ERROR in LauEmbeddedData::findBranches : Invalid pointer to the data tree object."<<endl;
45  return kFALSE;
46  }
47 
48  Bool_t branchesReadOK = theDataTree_->findBranches();
49  if ( ! branchesReadOK ) {
50  cerr<<"ERROR in LauEmbeddedData::findBranches : Problem finding branches."<<endl;
51  return kFALSE;
52  }
53 
55 
56  return kTRUE;
57 }
58 
60 {
61  if (!theDataTree_) {
62  cerr<<"ERROR in LauEmbeddedData::getReweightedEvent : Invalid pointer to the data tree object."<<endl;
63  return kFALSE;
64  }
65 
66  if (!dynamics) {
67  cerr<<"ERROR in LauEmbeddedData::getReweightedEvent : Amplitude model is null."<<endl;
68  return kFALSE;
69  }
70 
71  UInt_t numEvents = this->nEvents();
72  UInt_t iEvt(0);
73 
74  // First select a random event within the embedded data sample
75  if (this->reuseEvents()) {
76  iEvt = LauRandom::randomFun()->Integer(numEvents);
77  } else {
78  if (this->nUsedEvents() == numEvents) {
79  cerr<<"ERROR in LauEmbeddedData::getReweightedEvent : Have already used all events in the tree."<<endl;
80  return kFALSE;
81  }
82  Bool_t ok(kFALSE);
83  while (!ok) {
84  iEvt = LauRandom::randomFun()->Integer(numEvents);
85  ok = usedEvents_.insert(iEvt).second;
86  }
87  }
88 
89  // Retrieve the data for the selected event
90  theData_ = theDataTree_->getData(iEvt);
91 
92  LauKinematics* kinematics = dynamics->getKinematics();
93 
94  if (kinematics != 0) {
95 
96  // Get the true m13Sq, m23Sq variables.
97  Double_t m13Sq_MC = this->getValue("m13Sq_MC");
98  Double_t m23Sq_MC = this->getValue("m23Sq_MC");
99 
100  if (kinematics->withinDPLimits(m13Sq_MC,m23Sq_MC)) {
101 
102  // Update the kinematics with the true m13Sq, m23Sq variables
103  kinematics->updateKinematics(m13Sq_MC, m23Sq_MC);
104 
105  // Use this event or not, according to whether the amplitude model selects it.
106  // The LauIsobarDynamics uses an accept/reject method based on the
107  // ratio of the current to maximum ASq value. Other models may use alternative methods,
108  // but they must provide a yes/no answer if the event is accepted.
109  Bool_t gotReweightedEvent = dynamics->gotReweightedEvent();
110  if (gotReweightedEvent == kTRUE) {
111 
112  // Event is accepted.
113  // Update the kinematics to use the reco variables.
114  Double_t m13Sq = this->getValue("m13Sq");
115  Double_t m23Sq = this->getValue("m23Sq");
116  kinematics->updateKinematics(m13Sq, m23Sq);
117 
118  if (LauAbsDPDynamics::GenOK != dynamics->checkToyMC(kFALSE,kFALSE)) {
119  return kFALSE;
120  }
121 
122  } else {
123 
124  // Recursively run this function until we get an accepted event
125  return this->getReweightedEvent(dynamics);
126 
127  }
128 
129  } else {
130 
131  cerr<<"WARNING in LauEmbeddedData::getReweightedEvent : Skipping event "<<iEvt<<", which isn't within the DP boundary."<<endl;
132  return this->getReweightedEvent(dynamics);
133 
134  }
135  }
136 
137  return kTRUE;
138 }
139 
141 {
142  if (!theDataTree_) {
143  cerr<<"ERROR in LauEmbeddedData::getEmbeddedEvent : Invalid pointer to the data tree object."<<endl;
144  return;
145  }
146  UInt_t numEvents = this->nEvents();
147  UInt_t iEvt(0);
148  if (this->reuseEvents()) {
149  iEvt = LauRandom::randomFun()->Integer(numEvents);
150  } else {
151  if (this->nUsedEvents() == numEvents) {
152  cerr<<"ERROR in LauEmbeddedData::getEmbeddedEvent : Have already used all events in the tree."<<endl;
153  return;
154  }
155  Bool_t ok(kFALSE);
156  while (!ok) {
157  iEvt = LauRandom::randomFun()->Integer(numEvents);
158  ok = usedEvents_.insert(iEvt).second;
159  }
160  }
161  theData_ = theDataTree_->getData(iEvt);
162 
163  if (kinematics!=0) {
164  Double_t m13Sq = this->getValue("m13Sq");
165  Double_t m23Sq = this->getValue("m23Sq");
166  if (kinematics->withinDPLimits(m13Sq,m23Sq)) {
167  kinematics->updateKinematics(m13Sq,m23Sq);
168  } else {
169  cerr<<"WARNING in LauEmbeddedData::getEmbeddedEvent : Skipping event that isn't within the DP boundary."<<endl;
170  this->getEmbeddedEvent(kinematics);
171  }
172  }
173 }
174 
175 Double_t LauEmbeddedData::getValue(const TString& name) const
176 {
177  LauFitData::const_iterator iter = theData_.find(name);
178  if (iter == theData_.end()) {
179  cerr<<"ERROR in LauEmbeddedData::getValue : Could not find branch \""<<name<<"\" in embedded event."<<endl;
180  return 0.0;
181  } else {
182  return iter->second;
183  }
184 }
185 
186 LauFitData LauEmbeddedData::getValues(const std::vector<TString>& names) const
187 {
188  LauFitData retVal;
189 
190  for ( std::vector<TString>::const_iterator iter = names.begin(); iter != names.end(); ++iter ) {
191  const TString& name = (*iter);
192  LauFitData::const_iterator data_iter = theData_.find(name);
193  if ( data_iter == theData_.end() ) {
194  cerr<<"WARNING in LauEmbeddedData::getValues : Could not find branch \""<<name<<"\" in embedded event."<<endl;
195  } else {
196  retVal.insert( *data_iter );
197  }
198  }
199 
200  return retVal;
201 }
202 
Bool_t findBranches()
Find all of the branches in the tree.
Class for defining the abstract interface for signal Dalitz plot dynamics.
LauKinematics * getKinematics()
Retrieve the Dalitz plot kinematics.
TRandom * randomFun()
Access the singleton random number generator with a particular seed.
Definition: LauRandom.cc:20
UInt_t nEvents() const
Retrieve the number of events.
ClassImp(LauAbsCoeffSet)
const TString & name() const
The parameter name.
void readAllData()
Read all events from the tree.
Bool_t reuseEvents() const
Boolean determining whether events should be reused.
LauFitDataTree * theDataTree_
The structure containing the data.
std::map< TString, Double_t > LauFitData
Type for holding event data.
std::set< UInt_t > usedEvents_
Used events.
Double_t getValue(const TString &name) const
Get the value of a specified branch.
File containing declaration of LauKinematics class.
File containing declaration of LauEmbeddedData class.
void updateKinematics(Double_t m13Sq, Double_t m23Sq)
Update all kinematic quantities based on the DP co-ordinates m13Sq and m23Sq.
LauFitData getValues(const std::vector< TString > &names) const
Get values of specified branches.
Class to store the data for embedding in toy experiments.
void getEmbeddedEvent(LauKinematics *kinematics)
Retrieve an event from the data sample.
Bool_t withinDPLimits(Double_t m13Sq, Double_t m23Sq) const
Check whether a given (m13Sq,m23Sq) point is within the kinematic limits of the Dalitz plot...
virtual ToyMCStatus checkToyMC(Bool_t printErrorMessages=kTRUE, Bool_t printInfoMessages=kFALSE)=0
Check the status of the toy MC generation.
const LauFitData & getData(UInt_t iEvt) const
Retrieve the data for a given event.
File containing LauRandom namespace.
File containing declaration of LauAbsDPDynamics class.
UInt_t nUsedEvents() const
Retrieve the number of events that have already been sampled.
virtual ~LauEmbeddedData()
Destructor.
Bool_t getReweightedEvent(LauAbsDPDynamics *dynamics)
Retrieve an event from the data sample, applying an accept/reject based on the given DP model...
LauFitData theData_
The data for the currently retrieved event.
Class for calculating 3-body kinematic quantities.
Bool_t findBranches()
Find and read the branches in data tree.
Class to store the input fit variables.
virtual Bool_t gotReweightedEvent()=0
Calculates whether an event with the current kinematics should be accepted in order to produce a dist...