laura is hosted by Hepforge, IPPP Durham
Laura++  v2r1p1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauScfMap.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2006 - 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 <cstdlib>
16 #include <iostream>
17 #include <vector>
18 #include <map>
19 using std::cout;
20 using std::cerr;
21 using std::flush;
22 using std::endl;
23 
24 #include "TAxis.h"
25 #include "TH2.h"
26 #include "TSystem.h"
27 
28 #include "LauScfMap.hh"
29 
31 
32 
34  nBinsX_(-1),
35  nBinsY_(-1),
36  numberOfBins_(-1)
37 {
38 }
39 
41 {
42  for (std::vector<TH2*>::iterator iter = histos_.begin(); iter != histos_.end(); ++iter) {
43  delete (*iter);
44  }
45 }
46 
47 void LauScfMap::setHistos(const std::vector<const TH2*>& histos)
48 {
49  // Check that we have enough histograms
50  Int_t nBinsX = histos.front()->GetNbinsX();
51  Int_t nBinsY = histos.front()->GetNbinsY();
52  Int_t nHist = histos.size();
53 
54  if (nHist != nBinsX*nBinsY) {
55  cerr<<"ERROR in LauScfMap::setHistos : There should be as many histograms supplied as there are bins in each of them."<<endl;
56  cerr<<" : There are "<<nHist<<" histograms and "<<nBinsX*nBinsY<<" bins."<<endl;
57  return;
58  }
59 
60  // number of (tru, and reco as well) bins in a 2D histogram
61  // Later we will create as many fake points to cache DP amplitudes
62  nBinsX_ = nBinsX;
63  nBinsY_ = nBinsY;
65 
66  // Clear out any old histos we might have lying around
67  for (std::vector<TH2*>::iterator iter = histos_.begin(); iter != histos_.end(); ++iter) {
68  delete (*iter);
69  }
70  histos_.clear();
71 
72  // Make sure we have enough space for the new histograms
73  histos_.reserve(nHist);
74 
75  // Loop through the supplied vector
76  // Clone each histogram and store it
77  for (std::vector<const TH2*>::const_iterator iter = histos.begin(); iter != histos.end(); ++iter) {
78  TH2* histo = dynamic_cast<TH2*>((*iter)->Clone());
79  if ( histo==0 ) {
80  cerr<<"ERROR in LauScfMap::setHistos : Problem cloning one of the histograms."<<endl;
81  gSystem->Exit(EXIT_FAILURE);
82  }
83  histos_.push_back(histo);
84  }
85 
86  if ( histos_.size() != histos.size() ) {
87  cerr<<"ERROR in LauScfMap::setHistos : Problem cloning the histograms."<<endl;
88  gSystem->Exit(EXIT_FAILURE);
89  }
90 
91  // Now we need to setup the map
92 
93  // First clear it
94  contribs_.clear();
95 
96  // Fill in the map that links, to each reco bin, a vector with the numbers of the tru bins
97  // that have non-zero contributions in that reco bin
98  // Loop over each histogram bin
99  for (Int_t j(0); j<nBinsY_; ++j) {
100  for (Int_t i(0); i<nBinsX_; ++i) {
101  Int_t binNo = histos_.front()->GetBin(i+1,j+1);
102  std::vector<Int_t>& truBins = contribs_[binNo];
103  truBins.reserve(nHist);
104  for (Int_t k(0); k<nHist; ++k) {
105  Double_t content = histos_[k]->GetBinContent(binNo);
106  if (content>0.0) {
107  truBins.push_back(k);
108  }
109  }
110  }
111  }
112 }
113 
114 void LauScfMap::listBinCentres( std::vector<Double_t>& xCoords, std::vector<Double_t>& yCoords) const
115 {
116  // Create the list of fake points, located at the centres of the
117  // true bins, that we will add to the data points, so as to cache
118  // the DP amplitudes, scfFractions and jacobians at those points.
119  // that have non-zero contributions in that reco bin
120 
121  const TAxis* xAxis = histos_.front()->GetXaxis();
122  Double_t xMax = xAxis->GetXmax();
123  Double_t xMin = xAxis->GetXmin();
124  Double_t xIncrement = (Double_t)((xMax - xMin)/nBinsX_);
125 
126  const TAxis* yAxis = histos_.front()->GetYaxis();
127  Double_t yMax = yAxis->GetXmax();
128  Double_t yMin = yAxis->GetXmin();
129  Double_t yIncrement = (Double_t)((yMax - yMin)/nBinsY_);
130 
131  // Loop over each histogram bin
132  Double_t yCoord = yMin + yIncrement/2;
133  for (Int_t j(1); j<=nBinsY_; ++j) {
134  Double_t xCoord = xMin + xIncrement/2;
135  for (Int_t i(1); i<=nBinsX_; ++i) {
136  xCoords.push_back(xCoord);
137  yCoords.push_back(yCoord);
138  xCoord += xIncrement;
139  }
140  yCoord += yIncrement;
141  }
142 }
143 
144 Int_t LauScfMap::binNumber(Double_t xCoord, Double_t yCoord) const
145 {
146  const TH2* histo = histos_.front();
147  if (histo != 0) {
148  return histo->FindFixBin(xCoord,yCoord);
149  } else {
150  cerr<<"ERROR in LauScfMap::recoBin : No valid histograms found."<<endl;
151  return -1;
152  }
153 }
154 
155 const std::vector<Int_t>* LauScfMap::trueBins(Int_t recoBin) const
156 {
157  std::map< Int_t, std::vector<Int_t> >::const_iterator iter = contribs_.find(recoBin);
158  if (iter != contribs_.end()) {
159  return &(iter->second);
160  } else {
161  return 0;
162  }
163 }
164 
165 Double_t LauScfMap::prob(Int_t recoBin, Int_t trueBin) const
166 {
167  const TH2* histo = histos_[trueBin];
168  if (histo != 0) {
169  return histo->GetBinContent(recoBin);
170  } else {
171  cerr<<"ERROR in LauScfMap::prob : No histogram found for true bin "<<trueBin<<endl;
172  return 0.0;
173  }
174 }
175 
176 TH2* LauScfMap::trueHist(Int_t trueBin)
177 {
178  // Need to turn the histogram bin number into the vector element
179  Int_t x(0), y(0), z(0);
180  histos_.front()->GetBinXYZ( trueBin, x, y, z );
181 
182  Int_t theBin = nBinsX_*(y-1) + (x-1);
183 
184  if ( theBin < 0 || static_cast<UInt_t>(theBin) >= histos_.size() ) {
185  cerr<<"ERROR in LauScfMap::trueHist : No histogram found for true bin "<<trueBin<<", which corresponds to x="<<x<<", y="<<y<<", or entry "<<theBin<<" in the vector."<<endl;
186  gSystem->Exit(EXIT_FAILURE);
187  }
188 
189  TH2* histo = histos_[theBin];
190  if ( histo == 0 ) {
191  cerr<<"ERROR in LauScfMap::trueHist : Null histogram pointer found for true bin "<<trueBin<<", which corresponds to x="<<x<<", y="<<y<<", or entry "<<theBin<<" in the vector."<<endl;
192  gSystem->Exit(EXIT_FAILURE);
193  }
194 
195  return histo;
196 }
197 
ClassImp(LauAbsCoeffSet)
Int_t numberOfBins_
Number of bins in a 2D histogram (nBinsX_*nBinsY_)
Definition: LauScfMap.hh:116
TH2 * trueHist(Int_t trueBin)
Retrieve the migration histogram for trueBin.
Definition: LauScfMap.cc:176
void listBinCentres(std::vector< Double_t > &xCoords, std::vector< Double_t > &yCoords) const
Create lists of the co-ordinates of the centres of true bins.
Definition: LauScfMap.cc:114
const std::vector< Int_t > * trueBins(Int_t recoBin) const
Find which true bins contribute to the given reco bin.
Definition: LauScfMap.cc:155
File containing declaration of LauScfMap class.
void setHistos(const std::vector< const TH2 * > &histos)
Construct the smearing matrix from the supplied histograms.
Definition: LauScfMap.cc:47
Int_t binNumber(Double_t xCoord, Double_t yCoord) const
Find the global bin number for the given co-ordinates.
Definition: LauScfMap.cc:144
std::vector< TH2 * > histos_
The vector of two-dimensional histograms.
Definition: LauScfMap.hh:106
Int_t nBinsY_
Number of bins in the y-axis.
Definition: LauScfMap.hh:114
Int_t nBinsX_
Number of bins in the x-axis.
Definition: LauScfMap.hh:112
virtual ~LauScfMap()
Destructor.
Definition: LauScfMap.cc:40
Class for representing the 4D smearing matrix for mis-reconstructed signal (self cross feed) ...
Definition: LauScfMap.hh:29
std::map< Int_t, std::vector< Int_t > > contribs_
Map that links each reco global bin number and the vector of true bin indices into LauScfMap::histos_...
Definition: LauScfMap.hh:109
Double_t prob(Int_t recoBin, Int_t trueBin) const
Probability of a true event in the given true bin migrating to the reco bin.
Definition: LauScfMap.cc:165