laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0p1
A maximum likelihood fitting package for performing Dalitz-plot analysis.
Lau2DAbsHistDP.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2004 - 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 
17 #include "TAxis.h"
18 #include "TH2.h"
19 #include "TRandom.h"
20 #include "TSystem.h"
21 
22 #include "Lau2DAbsHistDP.hh"
23 #include "LauBifurcatedGaussPdf.hh"
24 #include "LauDaughters.hh"
25 #include "LauKinematics.hh"
26 #include "LauRandom.hh"
27 
29 
30 
31 Lau2DAbsHistDP::Lau2DAbsHistDP(const LauDaughters* daughters, Bool_t useUpperHalfOnly, Bool_t squareDP) :
32  kinematics_( (daughters!=0) ? daughters->getKinematics() : 0 ),
33  upperHalf_(useUpperHalfOnly),
34  squareDP_(squareDP)
35 {
36 }
37 
39 {
40 }
41 
43 {
44  TRandom* random = LauRandom::randomFun();
45 
46  Int_t nBinsX = static_cast<Int_t>(hist->GetNbinsX());
47  Int_t nBinsY = static_cast<Int_t>(hist->GetNbinsY());
48 
49  for (Int_t i(0); i<nBinsX; ++i) {
50  for (Int_t j(0); j<nBinsY; ++j) {
51  Double_t currentContent = hist->GetBinContent(i+1,j+1);
52  Double_t currentError = hist->GetBinError(i+1,j+1);
53  Double_t newContent = random->Gaus(currentContent,currentError);
54  if (newContent<0.0) {
55  hist->SetBinContent(i+1,j+1,0.0);
56  } else {
57  hist->SetBinContent(i+1,j+1,newContent);
58  }
59  }
60  }
61 }
62 
63 void Lau2DAbsHistDP::doBinFluctuation(TH2* hist, const TH2* errorHi, const TH2* errorLo)
64 {
65  LauParameter* mean = new LauParameter("mean", 0.5, 0.0, 1.0, kFALSE);
66  LauParameter* sigL = new LauParameter("sigmaL", 0.5, 0.0, 1.0, kFALSE);
67  LauParameter* sigR = new LauParameter("sigmaR", 0.5, 0.0, 1.0, kFALSE);
68 
69  std::vector<LauAbsRValue*> pars(3);
70  pars[0] = mean;
71  pars[1] = sigL;
72  pars[2] = sigR;
73 
74  const TString varName("tmp");
75 
76  Int_t nBinsX = static_cast<Int_t>(hist->GetNbinsX());
77  Int_t nBinsY = static_cast<Int_t>(hist->GetNbinsY());
78 
79  LauFitData genData;
80 
81  for (Int_t i(0); i<nBinsX; ++i) {
82  for (Int_t j(0); j<nBinsY; ++j) {
83  const Double_t currentContent = hist->GetBinContent(i+1,j+1);
84  const Double_t currentErrorLo = errorLo->GetBinContent(i+1,j+1);
85  const Double_t currentErrorHi = errorHi->GetBinContent(i+1,j+1);
86 
87  mean->value( currentContent );
88  sigL->value( currentErrorLo );
89  sigR->value( currentErrorHi );
90 
91  const Double_t minVal = TMath::Max( 0.0, currentContent-5.0*currentErrorLo );
92  const Double_t maxVal = TMath::Min( 1.0, currentContent+5.0*currentErrorHi );
93 
94  LauBifurcatedGaussPdf bfgaus(varName, pars, minVal, maxVal);
95  bfgaus.heightUpToDate(kFALSE);
96  genData = bfgaus.generate(0);
97 
98  const Double_t newContent = genData[varName];
99 
100  hist->SetBinContent(i+1,j+1,newContent);
101  }
102  }
103 
104  delete pars[0];
105  delete pars[1];
106  delete pars[2];
107  pars.clear();
108 }
109 
110 void Lau2DAbsHistDP::raiseOrLowerBins(TH2* hist, const Double_t avEff, const Double_t avEffError)
111 {
112  TRandom* random = LauRandom::randomFun();
113 
114  Double_t curAvg = this->computeAverageContents(hist);
115  Double_t newAvg = random->Gaus(avEff,avEffError);
116 
117  hist->Scale( newAvg / curAvg );
118 }
119 
120 Double_t Lau2DAbsHistDP::computeAverageContents(const TH2* hist) const
121 {
122  Double_t totalContent(0.0);
123  Int_t binsWithinDPBoundary(0);
124 
125  Int_t nBinsX = static_cast<Int_t>(hist->GetNbinsX());
126  Int_t nBinsY = static_cast<Int_t>(hist->GetNbinsY());
127 
128  // Loop through the bins and include any that have their centre or any
129  // of the four corners within the kinematic boundary
130  for ( Int_t i(0); i<nBinsX; ++i ) {
131  Double_t binXCentre = hist->GetXaxis()->GetBinCenter(i+1);
132  Double_t binXLowerEdge = hist->GetXaxis()->GetBinLowEdge(i+1);
133  Double_t binXUpperEdge = hist->GetXaxis()->GetBinUpEdge(i+1);
134 
135  for ( Int_t j(0); j<nBinsY; ++j ) {
136  Double_t binYCentre = hist->GetYaxis()->GetBinCenter(i+1);
137  Double_t binYLowerEdge = hist->GetYaxis()->GetBinLowEdge(i+1);
138  Double_t binYUpperEdge = hist->GetYaxis()->GetBinUpEdge(i+1);
139 
140  if ( withinDPBoundaries( binXCentre, binYCentre ) ||
141  withinDPBoundaries( binXLowerEdge, binYLowerEdge ) ||
142  withinDPBoundaries( binXUpperEdge, binYUpperEdge ) ||
143  withinDPBoundaries( binXLowerEdge, binYUpperEdge ) ||
144  withinDPBoundaries( binXUpperEdge, binYLowerEdge ) ) {
145 
146  totalContent += hist->GetBinContent(i+1, j+1);
147  ++binsWithinDPBoundary;
148  }
149  }
150  }
151 
152  return totalContent/binsWithinDPBoundary;
153 }
154 
155 Bool_t Lau2DAbsHistDP::withinDPBoundaries(Double_t x, Double_t y) const
156 {
158 }
159 
160 void Lau2DAbsHistDP::getUpperHalf(Double_t& x, Double_t& y) const
161 {
162  if ( upperHalf_ == kTRUE ) {
163  if ( squareDP_ == kFALSE && x > y ) {
164  Double_t temp = y;
165  y = x;
166  x = temp;
167  } else if ( squareDP_ == kTRUE && y > 0.5 ) {
168  y = 1.0 - y;
169  }
170  }
171 }
TRandom * randomFun()
Access the singleton random number generator with a particular seed.
Definition: LauRandom.cc:20
virtual Bool_t heightUpToDate() const
Check if the maximum height of the PDF is up to date.
Definition: LauAbsPdf.hh:264
ClassImp(LauAbsCoeffSet)
void raiseOrLowerBins(TH2 *hist, const Double_t avEff, const Double_t avEffError)
Rescale the histogram bin contents based on the desired average efficiency and its uncertainty...
void doBinFluctuation(TH2 *hist)
Fluctuate the contents of each histogram bin independently, in accordance with their errors...
LauParameter()
Default constructor.
Definition: LauParameter.cc:30
Class that defines the particular 3-body decay under study.
Definition: LauDaughters.hh:33
File containing declaration of LauDaughters class.
virtual LauFitData generate(const LauKinematics *kinematics)
Generate an event from the PDF.
Definition: LauAbsPdf.cc:298
std::map< TString, Double_t > LauFitData
Type for holding event data.
File containing declaration of LauKinematics class.
File containing declaration of LauBifurcatedGaussPdf class.
Abstract base class for defining a variation across a 2D DP based on a histogram. ...
Bool_t withinSqDPLimits(Double_t mPrime, Double_t thetaPrime) const
Check whether a given (m&#39;,theta&#39;) point is within the kinematic limits of the Dalitz plot...
Bool_t withinDPBoundaries(Double_t x, Double_t y) const
Check whether the given co-ordinates are within the kinematic boundary.
virtual ~Lau2DAbsHistDP()
Copy constructor.
Class for defining the fit parameter objects.
Definition: LauParameter.hh:34
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...
File containing LauRandom namespace.
Class for defining a bifurcated Gaussian PDF.
Bool_t squareDP_
Boolean for using square DP variables.
const LauKinematics * kinematics_
Kinematics used to check events are within DP boundary.
File containing declaration of Lau2DAbsHistDP class.
Bool_t upperHalf_
Boolean for using the upper half of DP.
void getUpperHalf(Double_t &x, Double_t &y) const
If only using the upper half of the (symmetric) DP then transform into the correct half...
Double_t computeAverageContents(const TH2 *hist) const
Compute the average bin content for bins within the kinematic boundary.
Double_t value() const
The value of the parameter.