// 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 required."<<std::endl;
- std::cerr<<" : Try to decrease the maximum allowed value of the total amplitude squared using the "
- <<"LauIsobarDynamics::setASqMaxValue(Double_t) function and re-run."<<std::endl;
- std::cerr<<" : This needs to be, perhaps significantly, less than "<<aSqMaxSet_<<std::endl;
- std::cerr<<" : Maximum value of ASq so far = "<<aSqMaxVar_<<std::endl;
+ 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;
}
- aSqMaxSet_ = 1.01 * aSqMaxVar_;
- std::cout<<"INFO in LauIsobarDynamics::checkToyMC : |A|^2 max reset to "<<aSqMaxSet_<<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 : aSqMaxSet_ was set to "<<aSqMaxSet_<<" but actual aSqMax was "<<aSqMaxVar_<<std::endl;
+ 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<<" : Please set aSqMaxSet >= "<<aSqMaxVar_<<" using the LauIsobarDynamics::setASqMaxValue(Double_t) function and re-run."<<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