laura is hosted by Hepforge, IPPP Durham
Laura++  v3r1
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 #include "LauKinematics.hh"
10 
12 
13 
14 LauDPPartialIntegralInfo::LauDPPartialIntegralInfo(const Double_t minm13, const Double_t maxm13,
15  const Double_t minm23, const Double_t maxm23,
16  const Double_t m13BinWidth, const Double_t m23BinWidth,
17  const Double_t precision,
18  const UInt_t nAmp,
19  const UInt_t nIncohAmp,
20  const Bool_t squareDP,
21  const LauKinematics* kinematics) :
22  minm13_(minm13),
23  maxm13_(maxm13),
24  minm23_(minm23),
25  maxm23_(maxm23),
26  m13BinWidth_(m13BinWidth),
27  m23BinWidth_(m23BinWidth),
28  nm13Points_(static_cast<UInt_t>((maxm13-minm13)/m13BinWidth)),
29  nm23Points_(static_cast<UInt_t>((maxm23-minm23)/m23BinWidth)),
30  nAmp_(nAmp),
31  nIncohAmp_(nIncohAmp),
32  squareDP_(squareDP)
33 {
34  const Double_t meanm13 = 0.5*(minm13 + maxm13);
35  const Double_t rangem13 = maxm13 - minm13;
36  const Double_t halfRangem13 = 0.5*rangem13;
37 
38  const Double_t meanm23 = 0.5*(minm23 + maxm23);
39  const Double_t rangem23 = maxm23 - minm23;
40  const Double_t halfRangem23 = 0.5*rangem23;
41 
42  const Double_t intFactor = halfRangem13*halfRangem23;
43 
44  // Raise error if squareDP is true but the kinematics object is not provided
45  if ( squareDP_ && kinematics == 0 ) {
46  std::cerr << "ERROR in LauDPPartialIntegralInfo constructor : Integration in the square DP has been specified but no valid kinematics object has been provided!" << std::endl;
47  return;
48  }
49 
50  // Avoid integral if we have no points in either x or y space
51  if (nm13Points_ == 0 || nm23Points_ == 0) {
52  std::cerr << "ERROR in LauDPPartialIntegralInfo constructor : Range has zero grid points in one or both of the dimensions!" << std::endl;
53  return;
54  }
55 
56  // Print a warning if we have a very large number of points
57  if ( (nm13Points_ * nm23Points_) > 8000000 ) {
58  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;
59  std::cerr << " : In case of problems, consider using LauIsobarDynamics::setNarrowResonanceThreshold and/or LauIsobarDynamics::setIntegralBinningFactor to tune the binning behaviour." << std::endl;
60  }
61 
62  LauIntegrals dpIntegrals(precision);
63  dpIntegrals.calcGaussLegendreWeights(nm13Points_, m13Points_, m13Weights_);
64  dpIntegrals.calcGaussLegendreWeights(nm23Points_, m23Points_, m23Weights_);
65 
66  // Print out total weights for the integration
67  Double_t totm13Weight(0.0), totm23Weight(0.0);
68  for (UInt_t i = 0; i < nm13Points_; ++i) {
69  totm13Weight += m13Weights_[i];
70  }
71  for (UInt_t i = 0; i < nm23Points_; ++i) {
72  totm23Weight += m23Weights_[i];
73  }
74 
75  if ( squareDP_ ) {
76  std::cout<<"INFO in LauDPPartialIntegralInfo constructor : nmPrimePoints = "<<nm13Points_<<", nthPrimePoints = "<<nm23Points_<<std::endl;
77  std::cout<<" : mPrimeBinWidth = "<<m13BinWidth_<<", thPrimeBinWidth = "<<m23BinWidth_<<std::endl;
78  std::cout<<" : Integrating over mPrime = "<<minm13_<<" to "<<maxm13_<<", thPrime = "<<minm23_<<" to "<<maxm23_<<std::endl;
79  std::cout<<" : totmPrimeWeight = "<<totm13Weight<<", totthPrimeWeight = "<<totm23Weight<<std::endl;
80  } else {
81  std::cout<<"INFO in LauDPPartialIntegralInfo constructor : nm13Points = "<<nm13Points_<<", nm23Points = "<<nm23Points_<<std::endl;
82  std::cout<<" : m13BinWidth = "<<m13BinWidth_<<", m23BinWidth = "<<m23BinWidth_<<std::endl;
83  std::cout<<" : Integrating over m13 = "<<minm13_<<" to "<<maxm13_<<", m23 = "<<minm23_<<" to "<<maxm23_<<std::endl;
84  std::cout<<" : totm13Weight = "<<totm13Weight<<", totm23Weight = "<<totm23Weight<<std::endl;
85  }
86 
87  // Calculate the m13 and m23 values at the grid points
88 
89  UInt_t midpoint = (nm13Points_ + 1)/2;
90  for (UInt_t i = 0; i < midpoint; ++i) {
91 
92  UInt_t ii = nm13Points_ - 1 - i; // symmetric i index
93 
94  Double_t dm13 = halfRangem13*m13Points_[i];
95  Double_t m13Val = meanm13 - dm13;
96  m13Points_[i] = m13Val;
97 
98  m13Val = meanm13 + dm13;
99  m13Points_[ii] = m13Val;
100 
101  }
102 
103  midpoint = (nm23Points_ +1)/2;
104  for (UInt_t i = 0; i < midpoint; i++) {
105 
106  UInt_t ii = nm23Points_ - 1 - i; // symmetric i index
107 
108  Double_t dm23 = halfRangem23*m23Points_[i];
109  Double_t m23Val = meanm23 - dm23;
110  m23Points_[i] = m23Val;
111 
112  m23Val = meanm23 + dm23;
113  m23Points_[ii] = m23Val;
114  }
115 
116  // Now compute the combined weights at each grid point
117  weights_.resize( nm13Points_ );
118  efficiencies_.resize( nm13Points_ );
119  amplitudes_.resize( nm13Points_ );
120  incohIntensities_.resize( nm13Points_ );
121  for (UInt_t i = 0; i < nm13Points_; ++i) {
122 
123  weights_[i].resize( nm23Points_ );
124  efficiencies_[i].resize( nm23Points_ );
125  amplitudes_[i].resize( nm23Points_ );
126  incohIntensities_[i].resize( nm23Points_ );
127 
128  for (UInt_t j = 0; j < nm23Points_; ++j) {
129 
130  Double_t weight = m13Weights_[i]*m23Weights_[j];
131 
132  Double_t jacobian(0.0);
133  if ( squareDP_ ) {
134  jacobian = kinematics->calcSqDPJacobian( m13Points_[i], m23Points_[j] );
135  } else {
136  jacobian = 4.0*m13Points_[i]*m23Points_[j];
137  }
138  weight *= (jacobian*intFactor);
139 
140  weights_[i][j] = weight;
141 
142  amplitudes_[i][j].resize( nAmp_ );
143  incohIntensities_[i][j].resize( nIncohAmp_ );
144 
145  } // j weights loop
146  } // i weights loop
147 }
148 
150 {
151 }
152 
153 std::ostream& operator<<( std::ostream& stream, const LauDPPartialIntegralInfo& infoRecord )
154 {
155  stream << "minm13 = " << infoRecord.getMinm13() << ", ";
156  stream << "maxm13 = " << infoRecord.getMaxm13() << ", ";
157  stream << "minm23 = " << infoRecord.getMinm23() << ", ";
158  stream << "maxm23 = " << infoRecord.getMaxm23() << ", ";
159  stream << "m13BinWidth = " << infoRecord.getM13BinWidth() << ", ";
160  stream << "m23BinWidth = " << infoRecord.getM23BinWidth() << std::endl;
161  return stream;
162 }
163 
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
File containing declaration of LauKinematics class.
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.
Class for calculating 3-body kinematic quantities.