+* Fixed bug in LauIsobarDynamics where multiple incoherent amplitudes led to nonsensical fit fractions
+
17th October 2014 John Back
* Added the ability to calculate the transition amplitude matrix T in LauKMatrixPropagator,
as well as a few other minor speed-up changes and code checks. Example/PlotKMatrixTAmp.cc
can be used to check the T amplitude variation, phase shift and inelasticity, for a given
K matrix channel, as a function of the invariant mass squared variable s
15th October 2014 Thomas Latham
* Add methods to LauIsobarDynamics to make the integration binning more tunable by the user:
- setNarrowResonanceThreshold - modify the value below which a resonance is considered to be narrow (defaults to 0.02 GeV/c2)
- setIntegralBinningFactor - modify the factor by which the narrow resonance width is divided to obtain the bin size (defaults to 100)
* Print warning messages if the memory usage is likely to be very large
13th October 2014 Thomas Latham
* Modify Makefile to allow compilation with ROOT 6 (in addition to maintaining support for ROOT 5)
* Fix a few compilation errors on MacOSX 10.9
13th October 2014 Daniel Craik
* Update LauModIndPartWave to allow knots at kinematic limits to be modified
- Add new method setKnotAmp to modify existing knots (and the knot at the upper kinematic limit which is automatically added at initialisation)
* Update doxygen for LauIsobarDynamics::addIncoherentResonance to mention that incoherent resonances must be added last
10th October 2014 Thomas Latham
* Add new method to LauResonanceMaker to set whether the radius of a given Blatt-Weisskopf category should be fixed of floated
* Modify the methods of LauResonanceMaker to set the radius value and whether it should be floated so that they work before and after the resonances have been created
9th October 2014 John Back
* Corrected the eta-eta' and 4pi phase space factors in LauKMatrixPropagator,
which is used for the K-matrix amplitude:
- calcEtaEtaPRho() does not include the mass difference term m_eta - m_eta'
following the recommendation in hep-ph/0204328 and from advice from M Pennington
- calcFourPiRho() incorporates a better parameterisation of the double integral of Eq 4 in
hep-ph/0204328 which avoids the exponential increase for small values of s (~< 0.1)
- More detailed comments are provided in the above two functions to explain what is
going on and the reason for the choices made
6th October 2014 Thomas Latham
* Implement the mechanism for floating Blatt-Weisskopf barrier factor radius parameters
30th September 2014 Thomas Latham
* Fix issue in the checks on toy MC generation validity
- in the case of exceeding max iterations it was possible to enter an infinite loop
- the checks now detect all three possible states:
- aSqMaxSet_ is too low (generation is biased) => increase aSqMaxSet_ value
- aSqMaxSet_ is too high => reduce aSqMaxSet_ value to improve efficiency
- aSqMaxSet_ is high (causing low efficiency) but cannot be lowered without biasing => increase iterationsMax_ limit
* Update resonance parameter in LauResonanceMaker to match PDG 2014
* Modify behaviour when TTree objects are saved into files to avoid having multiple cycle numbers present
29th September 2014 Daniel Craik
* Add support for incoherent resonances in the signal model
- LauIsobarDynamics updated to include incoherent terms
- ABC for incoherent resonances, LauAbsIncohRes, added deriving from LauAbsResonance
- LauGaussIncohRes added to implement a Gaussian incoherent resonance, deriving from LauAbsIncohRes
- Small changes to various other classes to enable incoherent terms
* Fixed small bug in LauMagPhaseCoeffSet which could hang if phase is exactly pi or -pi
* Added charged version of the BelleNR resonances to LauResonanceMaker
* Updated parameters in LauConstants to match PDG 2014
14th July 2014 Thomas Latham
* Add intial support for fully-symmetric final states such as B0 -> KS KS KS
- Performs the symmetrisation of the signal model
- Background (and efficiency) histogram classes need some work if the user wants to provide folded histograms
8th July 2014 Daniel Craik
* Add class for model-independent partial wave
- Uses splines to produce a smooth amplitude from a set of magnitude and phase values at given invariant masses
- The individual magnitudes and phases can be floated in the fit
16th June 2014 Thomas Latham
* Allow floating of resonance parameters in simultaneous fits
13th June 2014 Thomas Latham
* Fix bug in LauResonanceInfo cloning method, where the width parameter was given a clone of the mass
10th June 2014 Thomas Latham
* Add new function to allow sharing of resonance parameters between components that are not charged conjugates, e.g. LASS_BW and LASS_NR
9th June 2014 Thomas Latham and Daniel Craik
* Fix bug in the new integration scheme
- Was not accounting for cases where several resonances share a floating parameter
- Meant that the integrals and caches for that resonance were not being updated
* Introduce a change in the implementation of the helicity flip for neutral parent decays
- Prior to this change the helicity was flipped for any neutral resonance in the decay of a neutral particle.
- Now the flip no longer occurs in flavour-specific decays (such as Bs0 -> D0bar K- pi+ or B0 -> K+ pi- pi0) since it is only required in flavour-conjugate modes (such as B0 -> KS pi+ pi-).
- This does not affect any physical results but it does contribute a pi phase flip to some odd-spin resonances (for example K*(892)0 in Bs0->D0barKpi).
- Therefore results are not necessarily comparable between fits run before and after this changeset.
// Calculate the initial fit fractions (for later comparison with Toy MC, if required)
this->updateCoeffs(coeffs);
this->calcExtraInfo(kTRUE);
for (UInt_t i = 0; i < nAmp_; i++) {
for (UInt_t j = i; j < nAmp_; j++) {
std::cout<<"INFO in LauIsobarDynamics::initialise : Initial fit fraction for amplitude ("<<i<<","<<j<<") = "<<fitFrac_[i][j].genValue()<<std::endl;
}
}
for (UInt_t i = 0; i < nIncohAmp_; i++) {
std::cout<<"INFO in LauIsobarDynamics::initialise : Initial fit fraction for incoherent amplitude ("<<i<<","<<i<<") = "<<fitFrac_[i+nAmp_][i+nAmp_].genValue()<<std::endl;
}
std::cout<<"INFO in LauIsobarDynamics::initialise : Initial efficiency = "<<meanDPEff_.initValue()<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::initialise : Initial DPRate = "<<DPRate_.initValue()<<std::endl;
}
void LauIsobarDynamics::initSummary()
{
UInt_t i(0);
TString nameP = daughters_->getNameParent();
TString name1 = daughters_->getNameDaug1();
TString name2 = daughters_->getNameDaug2();
TString name3 = daughters_->getNameDaug3();
std::cout<<"INFO in LauIsobarDynamics::initSummary : We are going to do a DP with "<<nameP<<" going to "<<name1<<" "<<name2<<" "<<name3<<std::endl;
std::cout<<" : For the following resonance combinations:"<<std::endl;
std::cout<<" : In tracks 2 and 3:"<<std::endl;
for (i = 0; i < nAmp_; i++) {
if (resPairAmp_[i] == 1) {
std::cout<<" : "<<resTypAmp_[i]<<" to "<<name2<<", "<< name3<<std::endl;
}
}
for (i = 0; i < nIncohAmp_; i++) {
if (incohResPairAmp_[i] == 1) {
std::cout<<" : "<<incohResTypAmp_[i]<<" (incoherent) to "<<name2<<", "<< name3<<std::endl;
}
}
std::cout<<std::endl;
std::cout<<" : In tracks 1 and 3:"<<std::endl;
for (i = 0; i < nAmp_; i++) {
if (resPairAmp_[i] == 2) {
std::cout<<" : "<<resTypAmp_[i]<<" to "<<name1<<", "<< name3<<std::endl;
}
}
for (i = 0; i < nIncohAmp_; i++) {
if (incohResPairAmp_[i] == 2) {
std::cout<<" : "<<incohResTypAmp_[i]<<" (incoherent) to "<<name1<<", "<< name3<<std::endl;
}
}
std::cout<<std::endl;
std::cout<<" : In tracks 1 and 2:"<<std::endl;
for (i = 0; i < nAmp_; i++) {
if (resPairAmp_[i] == 3) {
std::cout<<" : "<<resTypAmp_[i]<<" to "<<name1<<", "<< name2<<std::endl;
}
}
for (i = 0; i < nIncohAmp_; i++) {
if (incohResPairAmp_[i] == 3) {
std::cout<<" : "<<incohResTypAmp_[i]<<" (incoherent) to "<<name1<<", "<< name2<<std::endl;
}
}
std::cout<<std::endl;
for (i = 0; i < nAmp_; i++) {
if (resPairAmp_[i] == 0) {
std::cout<<" : and a non-resonant amplitude."<<std::endl;
// Increment the number of resonance amplitudes we have so far
++nIncohAmp_;
// Finally, add the resonance object to the internal array
sigIncohResonances_.push_back(theResonance);
std::cout<<"INFO in LauIsobarDynamics::addIncohResonance : Successfully added incoherent resonance. Total number of incoherent resonances so far = "<<nIncohAmp_<<std::endl;
LauAbsResonance* prodPole = new LauKMatrixProdPole(poleName, poleIndex, resPairAmpInt, thePropagator, daughters_);
resTypAmp_.push_back(poleName);
resPairAmp_.push_back(resPairAmpInt);
nAmp_++;
sigResonances_.push_back(prodPole);
// Also store the propName-poleName pair for calculating total fit fractions later on
// (avoiding the need to use dynamic casts to check which resonances are of the K-matrix type)
kMatrixPropSet_[poleName] = propName;
std::cout<<"INFO in LauIsobarDynamics::addKMatrixProdPole : Successfully added K-matrix production pole term. Total number of resonances so far = "<<nAmp_<<std::endl;
} else {
std::cerr<<"ERROR in LauIsobarDynamics::addKMatrixProdPole : The propagator of the name "<<propName
<<" could not be found for the production pole "<<poleName<<std::endl;
LauAbsResonance* prodSVP = new LauKMatrixProdSVP(SVPName, channelIndex, resPairAmpInt, thePropagator, daughters_);
resTypAmp_.push_back(SVPName);
resPairAmp_.push_back(resPairAmpInt);
++nAmp_;
sigResonances_.push_back(prodSVP);
// Also store the SVPName-propName pair for calculating total fit fractions later on
// (avoiding the need to use dynamic casts to check which resonances are of the K-matrix type)
kMatrixPropSet_[SVPName] = propName;
std::cout<<"INFO in LauIsobarDynamics::addKMatrixProdSVP : Successfully added K-matrix production slowly-varying (SVP) term. Total number of resonances so far = "<<nAmp_<<std::endl;
} else {
std::cerr<<"ERROR in LauIsobarDynamics::addKMatrixProdSVP : The propagator of the name "<<propName
<<" could not be found for the production slowly-varying part "<<SVPName<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Found narrow resonance: "<<name<<", mass = "<<mass<<", width = "<<width<<", pair int = "<<pair<<std::endl;
if ( pair == 1 ) {
if ( mass < minm23 || mass > maxm23 ){ continue; }
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Found narrow resonance: "<<name<<", mass = "<<mass<<", width = "<<width<<", pair int = "<<pair<<std::endl;
if ( pair == 1 ) {
if ( mass < minm23 || mass > maxm23 ){ continue; }
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : One or more narrow resonances found in m12, integrating over whole Dalitz plot with bin width of "<<m13BinWidth<<" GeV/c2..."<<std::endl;
// Divide the plot into 3 regions: the resonance band and
// the two areas either side.
Double_t mass = m13NarrowRes.begin()->second;
Double_t width = m13NarrowRes.begin()->first;
Double_t resMin = mass - 5.0*width;
Double_t resMax = mass + 5.0*width;
// if the resonance is close to threshold just go from
// threshold to resMax, otherwise treat threshold to resMin
// as a separate region
if ( resMin < (minm13+50.0*m13BinWidth_) ) {
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : One narrow resonance found in m13, close to threshold, dividing Dalitz plot into two regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : One narrow resonance found in m13, dividing Dalitz plot into three regions..."<<std::endl;
// if the resonance is close to threshold just go from
// threshold to resMax, otherwise treat threshold to resMin
// as a separate region
if ( res1Min < (minm13+50.0*m13BinWidth_) ) {
if ( res1Max > res2Min ) {
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Two narrow resonances found in m13, both close to threshold, dividing Dalitz plot into two regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Two narrow resonances found in m13, one close to threshold, dividing Dalitz plot into four regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Two narrow resonances found close together in m13, dividing Dalitz plot into three regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Two narrow resonances found in m13, dividing Dalitz plot into five regions..."<<std::endl;
// Divide the plot into 3 regions: the resonance band and
// the two areas either side.
Double_t mass = m23NarrowRes.begin()->second;
Double_t width = m23NarrowRes.begin()->first;
Double_t resMin = mass - 5.0*width;
Double_t resMax = mass + 5.0*width;
// if the resonance is close to threshold just go from
// threshold to resMax, otherwise treat threshold to resMin
// as a separate region
if ( resMin < (minm23+50.0*m23BinWidth_) ) {
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : One narrow resonance found in m23, close to threshold, dividing Dalitz plot into two regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : One narrow resonance found in m23, dividing Dalitz plot into three regions..."<<std::endl;
// if the resonance is close to threshold just go from
// threshold to resMax, otherwise treat threshold to resMin
// as a separate region
if ( res1Min < (minm23+50.0*m23BinWidth_) ) {
if ( res1Max > res2Min ) {
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Two narrow resonances found in m23, both close to threshold, dividing Dalitz plot into two regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Two narrow resonances found in m23, one close to threshold, dividing Dalitz plot into four regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Two narrow resonances found close together in m23, dividing Dalitz plot into three regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Two narrow resonances found in m23, dividing Dalitz plot into five regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : One narrow resonance found in m13 and one in m23, both close to threshold, dividing Dalitz plot into four regions..."<<std::endl;
} else if ( resMin13 < (minm13+50.0*m13BinWidth_) ) {
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : One narrow resonance found in m13, close to threshold, and one in m23, not close to threshold, dividing Dalitz plot into six regions..."<<std::endl;
} else if ( resMin23 < (minm23+50.0*m23BinWidth_) ) {
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : One narrow resonance found in m23, close to threshold, and one in m13, not close to threshold, dividing Dalitz plot into six regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : One narrow resonance found in both m13 and m23, neither close to threshold, dividing Dalitz plot into nine regions..."<<std::endl;
// We have multiple narrow resonances in m13 only.
// Divide the plot into 2 regions: threshold to the most
// massive of the narrow resonances, and the rest
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Multiple narrow resonances found in m13, dividing Dalitz plot into two regions..."<<std::endl;
Double_t mass = 0.0;
Double_t width = 0.0;
for ( std::map<Double_t,Double_t>::const_iterator iter = m13NarrowRes.begin(); iter != m13NarrowRes.end(); ++iter ) {
// We have multiple narrow resonances in m23 only.
// Divide the plot into 2 regions: threshold to the most
// massive of the narrow resonances, and the rest
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Multiple narrow resonances found in m23, dividing Dalitz plot into two regions..."<<std::endl;
Double_t mass = 0.0;
Double_t width = 0.0;
for ( std::map<Double_t,Double_t>::const_iterator iter = m23NarrowRes.begin(); iter != m23NarrowRes.end(); ++iter ) {
// We've got a single narrow resonance in m13 and multiple
// narrow resonances in m23. Divide the plot into 6 regions.
Double_t mass23 = 0.0;
Double_t width23 = 0.0;
for ( std::map<Double_t,Double_t>::const_iterator iter = m23NarrowRes.begin(); iter != m23NarrowRes.end(); ++iter ) {
if ( mass23 < iter->second ) {
mass23 = iter->second;
width23 = iter->first;
}
}
Double_t resMax23 = mass23 + 5.0*width23;
Double_t mass13 = m13NarrowRes.begin()->second;
Double_t width13 = m13NarrowRes.begin()->first;
Double_t resMin13 = mass13 - 5.0*width13;
Double_t resMax13 = mass13 + 5.0*width13;
// if the m13 resonance is close to threshold just go from
// threshold to resMax, otherwise treat threshold to resMin
// as a separate region
if ( resMin13 < (minm13+50.0*m13BinWidth_) ) {
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Multiple narrow resonances found in m23 and one in m13, close to threshold, dividing Dalitz plot into four regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Multiple narrow resonances found in m23 and one in m13, not close to threshold, dividing Dalitz plot into six regions..."<<std::endl;
// We've got a single narrow resonance in m23 and multiple
// narrow resonances in m13. Divide the plot into 6 regions.
Double_t mass13 = 0.0;
Double_t width13 = 0.0;
for ( std::map<Double_t,Double_t>::const_iterator iter = m13NarrowRes.begin(); iter != m13NarrowRes.end(); ++iter ) {
if ( mass13 < iter->second ) {
mass13 = iter->second;
width13 = iter->first;
}
}
Double_t resMax13 = mass13 + 5.0*width13;
Double_t mass23 = m23NarrowRes.begin()->second;
Double_t width23 = m23NarrowRes.begin()->first;
Double_t resMin23 = mass23 - 5.0*width23;
Double_t resMax23 = mass23 + 5.0*width23;
// if the m23 resonance is close to threshold just go from
// threshold to resMax, otherwise treat threshold to resMin
// as a separate region
if ( resMin23 < (minm23+50.0*m23BinWidth_) ) {
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Multiple narrow resonances found in m13 and one in m23, close to threshold, dividing Dalitz plot into four regions..."<<std::endl;
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Multiple narrow resonances found in m13 and one in m23, not close to threshold, dividing Dalitz plot into six regions..."<<std::endl;
// We've got multiple narrow resonances in both m13 and m23.
// Divide the plot into 4 regions.
std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : Multiple narrow resonances found in both m13 and m23, dividing Dalitz plot into four regions..."<<std::endl;
Double_t mass13 = 0.0;
Double_t width13 = 0.0;
for ( std::map<Double_t,Double_t>::const_iterator iter = m13NarrowRes.begin(); iter != m13NarrowRes.end(); ++iter ) {
if ( mass13 < iter->second ) {
mass13 = iter->second;
width13 = iter->first;
}
}
Double_t resMax13 = mass13 + 5.0*width13;
Double_t mass23 = 0.0;
Double_t width23 = 0.0;
for ( std::map<Double_t,Double_t>::const_iterator iter = m23NarrowRes.begin(); iter != m23NarrowRes.end(); ++iter ) {
// Exceeded maximum allowed iterations - the generation is too inefficient
if (printErrorMessages) {
std::cerr<<"WARNING in LauIsobarDynamics::checkToyMC : More than "<<iterationsMax_<<" iterations performed and no event accepted."<<std::endl;
}
if ( aSqMaxSet_ > 1.01 * aSqMaxVar_ ) {
if (printErrorMessages) {
std::cerr<<" : |A|^2 maximum was set to "<<aSqMaxSet_<<" but this appears to be too high."<<std::endl;
std::cerr<<" : Maximum value of |A|^2 found so far = "<<aSqMaxVar_<<std::endl;
std::cerr<<" : The value of |A|^2 maximum will be decreased and the generation restarted."<<std::endl;
}
aSqMaxSet_ = 1.01 * aSqMaxVar_;
std::cout<<"INFO in LauIsobarDynamics::checkToyMC : |A|^2 max reset to "<<aSqMaxSet_<<std::endl;
ok = MaxIterError;
} else {
if (printErrorMessages) {
std::cerr<<" : |A|^2 maximum was set to "<<aSqMaxSet_<<", which seems to be correct for the given model."<<std::endl;
std::cerr<<" : However, the generation is very inefficient - please check your model."<<std::endl;
std::cerr<<" : The maximum number of iterations will be increased and the generation restarted."<<std::endl;
}
iterationsMax_ *= 2;
std::cout<<"INFO in LauIsobarDynamics::checkToyMC : max number of iterations reset to "<<iterationsMax_<<std::endl;
ok = MaxIterError;
}
} else if (aSqMaxVar_ > aSqMaxSet_) {
// Found a value of ASq higher than the accept/reject ceiling - the generation is biased
if (printErrorMessages) {
std::cerr<<"WARNING in LauIsobarDynamics::checkToyMC : |A|^2 maximum was set to "<<aSqMaxSet_<<" but a value exceeding this was found: "<<aSqMaxVar_<<std::endl;
std::cerr<<" : Run was invalid, as any generated MC will be biased, according to the accept/reject method!"<<std::endl;
std::cerr<<" : The value of |A|^2 maximum be reset to be > "<<aSqMaxVar_<<" and the generation restarted."<<std::endl;
}
aSqMaxSet_ = 1.01 * aSqMaxVar_;
std::cout<<"INFO in LauIsobarDynamics::checkToyMC : |A|^2 max reset to "<<aSqMaxSet_<<std::endl;
ok = ASqMaxError;
} else if (printInfoMessages) {
std::cout<<"INFO in LauIsobarDynamics::checkToyMC : aSqMaxSet = "<<aSqMaxSet_<<" and aSqMaxVar = "<<aSqMaxVar_<<std::endl;
std::cerr<<"ERROR in LauIsobarDynamics::setDataEventNo : Event index too large: "<<iEvt<<" >= "<<data_.size()<<"."<<std::endl;
}
m13Sq_ = currentEvent_->retrievem13Sq();
m23Sq_ = currentEvent_->retrievem23Sq();
mPrime_ = currentEvent_->retrievemPrime();
thPrime_ = currentEvent_->retrievethPrime();
tagCat_ = currentEvent_->retrieveTagCat();
eff_ = currentEvent_->retrieveEff();
scfFraction_ = currentEvent_->retrieveScfFraction(); // These two are necessary, even though the dynamics don't actually use scfFraction_ or jacobian_,
jacobian_ = currentEvent_->retrieveJacobian(); // since this is at the heart of the caching mechanism.
std::cerr << "ERROR in LauIsobarDynamics::updateCoeffs : Expected " << this->getnTotAmp() << " but got " << coeffs.size() << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
// Now check if the coeffs have changed
Bool_t changed = (Amp_ != coeffs);
if (changed) {
// Copy the coeffs
Amp_ = coeffs;
}
// TODO should perhaps keep track of whether the resonance parameters have changed here and if none of those and none of the coeffs have changed then we don't need to update the norm
// Update the total normalisation for the signal likelihood