laura is hosted by Hepforge, IPPP Durham
Laura++  v3r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauDPPartialIntegralInfo.cc
Go to the documentation of this file.
1 
5 #include <iostream>
6 
8 #include "LauIntegrals.hh"
9 
11 
12 
13 LauDPPartialIntegralInfo::LauDPPartialIntegralInfo(const Double_t minm13, const Double_t maxm13,
14  const Double_t minm23, const Double_t maxm23,
15  const Double_t m13BinWidth, const Double_t m23BinWidth,
16  const Double_t precision,
17  const UInt_t nAmp,
18  const UInt_t nIncohAmp) :
19  minm13_(minm13),
20  maxm13_(maxm13),
21  minm23_(minm23),
22  maxm23_(maxm23),
23  m13BinWidth_(m13BinWidth),
24  m23BinWidth_(m23BinWidth),
25  nm13Points_(static_cast<UInt_t>((maxm13-minm13)/m13BinWidth)),
26  nm23Points_(static_cast<UInt_t>((maxm23-minm23)/m23BinWidth)),
27  nAmp_(nAmp),
28  nIncohAmp_(nIncohAmp)
29 {
30  const Double_t meanm13 = 0.5*(minm13 + maxm13);
31  const Double_t rangem13 = maxm13 - minm13;
32  const Double_t halfRangem13 = 0.5*rangem13;
33 
34  const Double_t meanm23 = 0.5*(minm23 + maxm23);
35  const Double_t rangem23 = maxm23 - minm23;
36  const Double_t halfRangem23 = 0.5*rangem23;
37 
38  const Double_t intFactor = halfRangem13*halfRangem23;
39 
40  // Avoid integral if we have no points in either x or y space
41  if (nm13Points_ == 0 || nm23Points_ == 0) {
42  std::cerr << "ERROR in LauDPPartialIntegralInfo constructor : Range has zero grid points in one or both of the dimensions!" << std::endl;
43  return;
44  }
45 
46  // Print a warning if we have a very large number of points
47  if ( (nm13Points_ * nm23Points_) > 8000000 ) {
48  std::cerr << "WARNING in LauDPPartialIntegralInfo constructor : The integration binning scheme has a very large number of bins, this could cause high memory consumption!" << std::endl;
49  std::cerr << " : In case of problems, consider using LauIsobarDynamics::setNarrowResonanceThreshold and/or LauIsobarDynamics::setIntegralBinningFactor to tune the binning behaviour." << std::endl;
50  }
51 
52  LauIntegrals dpIntegrals(precision);
53  dpIntegrals.calcGaussLegendreWeights(nm13Points_, m13Points_, m13Weights_);
54  dpIntegrals.calcGaussLegendreWeights(nm23Points_, m23Points_, m23Weights_);
55 
56  // Print out total weights for the integration
57  Double_t totm13Weight(0.0), totm23Weight(0.0);
58  for (UInt_t i = 0; i < nm13Points_; ++i) {
59  totm13Weight += m13Weights_[i];
60  }
61  for (UInt_t i = 0; i < nm23Points_; ++i) {
62  totm23Weight += m23Weights_[i];
63  }
64 
65  std::cout<<"INFO in LauDPPartialIntegralInfo constructor : nm13Points = "<<nm13Points_<<", nm23Points = "<<nm23Points_<<std::endl;
66  std::cout<<" : m13BinWidth = "<<m13BinWidth_<<", m23BinWidth = "<<m23BinWidth_<<std::endl;
67  std::cout<<" : Integrating over m13 = "<<minm13_<<" to "<<maxm13_<<", m23 = "<<minm23_<<" to "<<maxm23_<<std::endl;
68  std::cout<<" : totm13Weight = "<<totm13Weight<<", totm23Weight = "<<totm23Weight<<std::endl;
69 
70  // Calculate the m13 and m23 values at the grid points
71 
72  UInt_t midpoint = (nm13Points_ + 1)/2;
73  for (UInt_t i = 0; i < midpoint; ++i) {
74 
75  UInt_t ii = nm13Points_ - 1 - i; // symmetric i index
76 
77  Double_t dm13 = halfRangem13*m13Points_[i];
78  Double_t m13Val = meanm13 - dm13;
79  m13Points_[i] = m13Val;
80 
81  m13Val = meanm13 + dm13;
82  m13Points_[ii] = m13Val;
83 
84  }
85 
86  midpoint = (nm23Points_ +1)/2;
87  for (UInt_t i = 0; i < midpoint; i++) {
88 
89  UInt_t ii = nm23Points_ - 1 - i; // symmetric i index
90 
91  Double_t dm23 = halfRangem23*m23Points_[i];
92  Double_t m23Val = meanm23 - dm23;
93  m23Points_[i] = m23Val;
94 
95  m23Val = meanm23 + dm23;
96  m23Points_[ii] = m23Val;
97  }
98 
99  // Now compute the combined weights at each grid point
100  weights_.resize( nm13Points_ );
101  efficiencies_.resize( nm13Points_ );
102  amplitudes_.resize( nm13Points_ );
103  incohIntensities_.resize( nm13Points_ );
104  for (UInt_t i = 0; i < nm13Points_; ++i) {
105 
106  weights_[i].resize( nm23Points_ );
107  efficiencies_[i].resize( nm23Points_ );
108  amplitudes_[i].resize( nm23Points_ );
109  incohIntensities_[i].resize( nm23Points_ );
110 
111  for (UInt_t j = 0; j < nm23Points_; ++j) {
112 
113  Double_t weight = m13Weights_[i]*m23Weights_[j];
114  Double_t jacobian = 4.0*m13Points_[i]*m23Points_[j];
115  weight *= (jacobian*intFactor);
116 
117  weights_[i][j] = weight;
118 
119  amplitudes_[i][j].resize( nAmp_ );
120  incohIntensities_[i][j].resize( nIncohAmp_ );
121 
122  } // j weights loop
123  } // i weights loop
124 }
125 
127 {
128 }
129 
130 std::ostream& operator<<( std::ostream& stream, const LauDPPartialIntegralInfo& infoRecord )
131 {
132  stream << "minm13 = " << infoRecord.getMinm13() << ", ";
133  stream << "maxm13 = " << infoRecord.getMaxm13() << ", ";
134  stream << "minm23 = " << infoRecord.getMinm23() << ", ";
135  stream << "maxm23 = " << infoRecord.getMaxm23() << ", ";
136  stream << "m13BinWidth = " << infoRecord.getM13BinWidth() << ", ";
137  stream << "m23BinWidth = " << infoRecord.getM23BinWidth() << std::endl;
138  return stream;
139 }
140 
Double_t getMinm13() const
Retrieve the minm13 of DP.
virtual ~LauDPPartialIntegralInfo()
Destructor.
ClassImp(LauAbsCoeffSet)
Double_t getM23BinWidth() const
Retrieve the m23BinWidth of DP.
Class for defining (a section of) the Dalitz plot integration binning scheme.
Class for performing numerical integration routines.
Definition: LauIntegrals.hh:30
void calcGaussLegendreWeights(const Int_t numPoints, std::vector< Double_t > &abscissas, std::vector< Double_t > &weights)
Calculate the Gauss-Legendre weights.
Definition: LauIntegrals.cc:50
std::ostream & operator<<(std::ostream &os, const LauComplex &z)
Definition: LauComplex.cc:43
File containing declaration of LauDPPartialIntegralInfo class.
Double_t getMaxm13() const
Retrieve the maxm13 of DP.
Double_t getMinm23() const
Retrieve the minm23 of DP.
Double_t getM13BinWidth() const
Retrieve the m13BinWidth of DP.
Double_t getMaxm23() const
Retrieve the maxm23 of DP.
File containing declaration of LauIntegrals class.