\param [in] effModel the model to describe the efficiency across the Dalitz plot
\param [in] scfFractionModel the models to describe the fraction of poorly constructed events (the self cross feed fraction) across the Dalitz plot for various tagging categories
//! Calculate the fit fractions, mean efficiency and total DP rate
/*!
\param [in] init whether the calculated values should be stored as the initial/generated values or the fitted values
*/
void calcExtraInfo(const Bool_t init = kFALSE);
//! Calculates whether an event with the current kinematics should be accepted in order to produce a distribution of events that matches the model e.g. when reweighting embedded data
/*!
Uses the accept/reject method.
\return kTRUE if the event has been accepted, kFALSE otherwise
*/
Bool_t gotReweightedEvent();
//! Calculate the acceptance rate, for events with the current kinematics, when generating events according to the model
/*!
\return the weight for the current kinematics
*/
Double_t getEventWeight();
//! Retrieve the total amplitude for the current event
- //! Calculate the integration scheme for narrow m23Regions that are integrated according to m23Binnings grid size, and the overlap according to m13Binnings x m23Binnings
- /*!
- \return vector of pointers to LauDPPartialIntegralInfo objects that contain the individual regions
// 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<<" : A"<<i<<": "<<resTypAmp_[i]<<" to "<<name2<<", "<< name3<<std::endl;
}
}
for (i = 0; i < nIncohAmp_; i++) {
if (incohResPairAmp_[i] == 1) {
std::cout<<" : A"<<nAmp_+i<<": "<<incohResTypAmp_[i]<<" (incoherent) to "<<name2<<", "<< name3<<std::endl;
}
}
std::cout<<" : In tracks 1 and 3:"<<std::endl;
for (i = 0; i < nAmp_; i++) {
if (resPairAmp_[i] == 2) {
std::cout<<" : A"<<i<<": "<<resTypAmp_[i]<<" to "<<name1<<", "<< name3<<std::endl;
}
}
for (i = 0; i < nIncohAmp_; i++) {
if (incohResPairAmp_[i] == 2) {
std::cout<<" : A"<<nAmp_+i<<": "<<incohResTypAmp_[i]<<" (incoherent) to "<<name1<<", "<< name3<<std::endl;
}
}
std::cout<<" : In tracks 1 and 2:"<<std::endl;
for (i = 0; i < nAmp_; i++) {
if (resPairAmp_[i] == 3) {
std::cout<<" : A"<<i<<": "<<resTypAmp_[i]<<" to "<<name1<<", "<< name2<<std::endl;
}
}
for (i = 0; i < nIncohAmp_; i++) {
if (incohResPairAmp_[i] == 3) {
std::cout<<" : A"<<nAmp_+i<<": "<<incohResTypAmp_[i]<<" (incoherent) to "<<name1<<", "<< name2<<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_, useProdAdler);
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_, useProdAdler);
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 ){
+ std::cout << std::string(53, ' ') << ": But its pole is outside the kinematically allowed range, so will not consider it narrow for the purposes of integration" << std::endl;
+ std::cout << std::string(53, ' ') << ": But its pole is outside the kinematically allowed range, so will not consider it narrow for the purposes of integration" << std::endl;
+ std::cout << std::string(53, ' ') << ": But its pole is outside the kinematically allowed range, so will not consider it narrow for the purposes of integration" << std::endl;
+ const TString& name = (*iter)->getResonanceName();
Int_t pair = (*iter)->getPairInt();
+ Double_t mass = (*iter)->getMass();
+ Double_t width = (*iter)->getWidth();
- TString name = (*iter)->getResonanceName();
+ if ( width > narrowWidth_ || width == 0.0 ) {
+ 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 ){
std::cout << std::string(53, ' ') << ": But its pole is outside the kinematically allowed range, so will not consider it narrow for the purposes of integration" << std::endl;
std::cout << std::string(53, ' ') << ": But its pole is outside the kinematically allowed range, so will not consider it narrow for the purposes of integration" << std::endl;
std::cout << std::string(53, ' ') << ": But its pole is outside the kinematically allowed range, so will not consider it narrow for the purposes of integration" << std::endl;
+ // Depending on how many narrow resonances we have and where they are
+ // we adopt different approaches
+ if ( ! m12NarrowRes.empty() ) {
- } else {
+ // We have at least one narrow resonance in m12
+ // Switch to using the square DP for the integration
+ // TODO - for the time being just use a single, reasonably fine by default and tunable, grid
+ // - can later consider whether there's a need to split up the mPrime axis into regions around particularly narrow resonances in m12
+ // - but it seems that this isn't really needed since even the default tune gives a good resolution for most narrow resonances such as phi / chi_c0
+ std::cout<<"INFO in LauIsobarDynamics::calcDPNormalisationScheme : One or more narrow resonances found in m12, integrating over whole square Dalitz plot with bin widths of "<<mPrimeBinWidth_<<" in mPrime and "<<thPrimeBinWidth_<<" in thetaPrime..."<<std::endl;
+ for ( std::vector<std::pair<Double_t, Double_t> >::const_iterator iter = m13NarrowRes.begin(); iter != m13NarrowRes.end(); ++iter ) {
+ Double_t mass = iter->first;
+ Double_t width = iter->second;
- for ( std::vector<std::pair<Double_t, Double_t> >::const_iterator iter = m23NarrowRes.begin(); iter != m23NarrowRes.end(); ++iter ) {
- Double_t width = iter->first;
- Double_t mass = iter->second;
+ Double_t regionBegin = mass - 5.0 * width;
+ Double_t regionEnd = mass + 5.0 * width;
+ Double_t binning = width / binningFactor_;
- Double_t regionBegin = mass - 5.0 * width;
- Double_t regionEnd = mass + 5.0 * width;
+ // check if we ought to extend the region to the edge of the phase space (in either direction)
+ if ( regionBegin < (minm13+50.0*m13BinWidth_) ) {
+ std::cout << "INFO in LauIsobarDynamics::calcDPNormalisationScheme : Resonance at m13 = " << mass << " is close to threshold, extending integration region" << std::endl;
+ regionBegin = minm13;
+ }
+ if ( regionEnd > (maxm13-50.0*m13BinWidth_) ) {
+ std::cout << "INFO in LauIsobarDynamics::calcDPNormalisationScheme : Resonance at m13 = " << mass << " is close to upper edge of phase space, extending integration region" << std::endl;
- // Get coarse regions by calculating the gaps between the narrow resonances and using the same functions to calculate the integrate scheme for each
+ // check if we ought to extend the region to the edge of the phase space (in either direction)
+ if ( regionBegin < (minm23+50.0*m23BinWidth_) ) {
+ std::cout << "INFO in LauIsobarDynamics::calcDPNormalisationScheme : Resonance at m23 = " << mass << " is close to threshold, extending integration region" << std::endl;
+ regionBegin = minm23;
+ }
+ if ( regionEnd > (maxm23-50.0*m23BinWidth_) ) {
+ std::cout << "INFO in LauIsobarDynamics::calcDPNormalisationScheme : Resonance at m23 = " << mass << " is close to upper edge of phase space, extending integration region" << std::endl;
// 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