("dtype", po::value<LauTimeDepFitModel::CPEigenvalue>(&dType)->default_value(LauTimeDepFitModel::CPEigenvalue::QFS,"QFS"), "type of D decay: QFS, CPOdd, or CPEven")
("dtr-perevent", po::value<Bool_t>(&perEventTimeErr)->default_value(kFALSE), "enable/disable use of per-event decay-time error (requires decay-time resolution to be enabled to take effect)")
("seed", po::value<UInt_t>(&RNGseed)->default_value(0), "set the seed for the RNG; if not set, the time is used to generate a seed")
("dir", po::value<std::string>(&directory)->default_value("","<none>"), "set the directory used to find the nTuples within the root file. Defaults to no directory")
("bkgndList", po::value<std::string>(&bkgndListStr)->default_value("comb,Bd2DKpi,Bs2DKpi,Lb2Dppi"), "the comma-separated list of background components to include")
("fitBackInput", po::value<std::string>(&fitBackInput)->default_value("","fit0_ToyMC_QFS_expts0-0_expt0.root"), "set the file used as input to the fit: only used if fitBack is also selected")
("firstExpt", po::value<UInt_t>(&firstExptFit)->default_value(0), "first experiment to fit")
("nExpt", po::value<UInt_t>(&nExptFit)->default_value(1), "number of experiments to fit")
("iFit", po::value<UInt_t>(&iFit)->default_value(0), "fit ID - used to distinguish multiple fits to a given dataset (used to disentangle multiple minima)")
("fixTau", po::value<Bool_t>(&fixLifetime)->default_value(kTRUE), "enable/disable floating of B lifetime parameter")
("fixDm", po::value<Bool_t>(&fixDeltaM)->default_value(kTRUE), "enable/disable floating of B mixing frequency parameter")
("fixPhiMix", po::value<Bool_t>(&fixPhiMix)->default_value(kFALSE), "enable/disable floating of B mixing phase parameter(s)")
("fixSplineKnots", po::value<Bool_t>(&fixSplineKnots)->default_value(kFALSE), "enable/disable floating of the decay-time-acceptance spline")
("useSinCos", po::value<Bool_t>(&useSinCos)->default_value(kTRUE), "enable/disable using sine and cosine of B mixing phase as the floating parameters rather than the phase itself")
("useAveDeltaCalibVals", po::value<Bool_t>(&useAveDeltaCalibVals)->default_value(kTRUE), "enable/disable fitting calib values as their average +/- delta values rather than having separate values for B and Bbar")
("floatCalibPars", po::value<Bool_t>(&floatCalibPars)->default_value(kTRUE), "enable/disable floating of the FT calibration parameters")
("floatYields", po::value<Bool_t>(&floatYields)->default_value(kFALSE), "enable/disable floating of the yields of each component - if enabled/disabled implies inclusion/exclusion of non-DP PDFs in/from the likelihood")
("fitBack", po::value<Bool_t>(&fitBack)->default_value(kFALSE), "enable/disable refitting some generated data from a model")
+ ("blindFit", po::value<Bool_t>(&blindFit)->default_value(kFALSE), "enable/disable blinding of the physics parameters")
std::cerr << "WARNING in LauParameter::blindParameter : blinding has already been set up for a clone of this parameter - it will be replaced!" << std::endl;
// first check that min is less than max (or they are the same - this is allowed)
if (minVal > maxVal) {
std::cerr<<"ERROR in LauParameter::checkRange : minValue: "<<minVal<<" greater than maxValue: "<<maxVal<<std::endl;
if (minValue_ > maxValue_) {
minValue_ = maxValue_;
std::cerr<<" : Setting both to "<<maxValue_<<"."<<std::endl;
} else {
std::cerr<<" : Not changing anything."<<std::endl;
}
return;
}
minValue_ = minVal;
maxValue_ = maxVal;
// now check that the value is still within the range
if ((val < minVal) || (val > maxVal)) {
if (name_ != "") {
std::cerr<<"ERROR in LauParameter::checkRange : value: "<<val<<" not in allowed range ["<<minVal<<" , "<<maxVal<<"] for parameter \""<<name_<<"\"."<<std::endl;
} else {
std::cerr<<"ERROR in LauParameter::checkRange : value: "<<val<<" not in allowed range ["<<minVal<<" , "<<maxVal<<"]."<<std::endl;
}
if (val < minVal) {
std::cerr<<" : Setting value to minValue: "<<minVal<<std::endl;
value_ = minVal;
} else {
std::cerr<<" : Setting value to maxValue: "<<maxVal<<std::endl;
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The LauParameter pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( signalEvents_ != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl;
return;
}
if ( signalAsym_ != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl;
return;
}
signalEvents_ = nSigEvents;
TString yieldName = signalEvents_->name();
if ( ! yieldName.BeginsWith("signalEvents") ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The signal yield parameter name should start with \"signalEvents\" (plus any optional suffix)." << std::endl;
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The event LauParameter pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( sigAsym == 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The asym LauParameter pointer is null." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( signalEvents_ != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl;
return;
}
if ( signalAsym_ != 0 ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : You are trying to overwrite the signal asymmetry." << std::endl;
return;
}
signalEvents_ = nSigEvents;
TString yieldName = signalEvents_->name();
if ( ! yieldName.BeginsWith("signalEvents") ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNSigEvents : The signal yield parameter name should start with \"signalEvents\" (plus any optional suffix)." << std::endl;
if ( ! this->validBkgndClass( bkgndClassName ) ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : Invalid background class \"" << bkgndClassName << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
if ( ! this->validBkgndClass( bkgndClassName ) ) {
std::cerr << "ERROR in LauTimeDepFitModel::setNBkgndEvents : Invalid background class \"" << bkgndClassName << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
// TODO If these are all histograms shouldn't need to add much more code in other functions
if (pdf==0) {
std::cerr<<"ERROR in LauTimeDepFitModel::setBkgndDtPdf : The PDF pointer is null, not adding it."<<std::endl;
return;
}
// check that this background name is valid
if ( ! this->validBkgndClass( bkgndClass) ) {
std::cerr << "ERROR in LauTimeDepFitModel::setBkgndDtPdf : Invalid background class \"" << bkgndClass << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
std::cerr << "ERROR in LauTimeDepFitModel::setBkgndDPModels : the model pointer is null for the particle model." << std::endl;
return;
}
// check that this background name is valid
if ( ! this->validBkgndClass( bkgndClass) ) {
std::cerr << "ERROR in LauTimeDepFitModel::setBkgndDPModels : Invalid background class \"" << bkgndClass << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
std::cout << "INFO in LauTimeDepFitModel::setBkgndDPModels : the model pointer is null for the anti-particle model. Using only the particle model." << std::endl;
std::cerr << "ERROR in LauTimeDepFitModel::setBkgndPdf : PDF pointer is null." << std::endl;
return;
}
// check that this background name is valid
if ( ! this->validBkgndClass( bkgndClass ) ) {
std::cerr << "ERROR in LauTimeDepFitModel::setBkgndPdf : Invalid background class \"" << bkgndClass << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
+ std::cerr << "WARNING in LauTimeDepFitModel::blindPhiMix : not using sine/cosine parameters, so using the sine blinding string to blind the phase itself." << std::endl;
// std::cerr << "ERROR in LauTimeDepFitModel::initialise : There are " << nsigpdfvars << " signal PDF variables but " << nbkgndpdfvars << " bkgnd PDF variables." << std::endl;
// gSystem->Exit(EXIT_FAILURE);
// }
// }
//}
// Clear the vectors of parameter information so we can start from scratch
this->clearFitParVectors();
// Set the fit parameters for signal and background models
this->setSignalDPParameters();
// Set the fit parameters for the decay time models
this->setDecayTimeParameters();
// Set the fit parameters for the extra PDFs
this->setExtraPdfParameters();
// Set the initial bg and signal events
this->setFitNEvents();
// Handle flavour-tagging calibration parameters
this->setCalibParams();
// Add tagging efficiency parameters
this->setTagEffParams();
//Asymmetry terms AProd and in setAsymmetries()?
this->setAsymParams();
// Check that we have the expected number of fit variables
std::cerr << "ERROR in LauTimeDepFitModel::initialiseDPModels : Unequal number of signal DP components in the particle and anti-particle models: " << nAmpB0bar << " != " << nAmpB0 << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
if ( nAmpB0bar != nSigComp_ ) {
std::cerr << "ERROR in LauTimeDepFitModel::initialiseDPModels : Number of signal DP components in the model (" << nAmpB0bar << ") not equal to number of coefficients supplied (" << nSigComp_ << ")." << std::endl;
gSystem->Exit(EXIT_FAILURE);
}
std::cout<<"INFO in LauTimeDepFitModel::initialiseDPModels : Initialising signal DP model"<<std::endl;
// TODO should check (first time) that they match in terms of number of entries in the vectors and that each entry has the same number of points, ranges, weights etc.
std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0bar signal DP model doesn't contain component \""<<compName<<"\"."<<std::endl;
return;
}
std::cerr<<"WARNING in LauTimeDepFitModel::setAmpCoeffSet : B0bar signal DP model doesn't contain component \""<<compName<<"\" but does contain the conjugate \""<<conjName<<"\", resetting name to use the conjugate."<<std::endl;
TString tmp = compName;
compName = conjName;
conjName = tmp;
coeffSet->name( compName );
}
if ( conjugate ) {
if ( ! sigModelB0_->hasResonance(conjName) ) {
std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<<conjName<<"\"."<<std::endl;
return;
}
} else {
if ( ! sigModelB0_->hasResonance(compName) ) {
std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : B0 signal DP model doesn't contain component \""<<compName<<"\"."<<std::endl;
return;
}
}
// Do we already have it in our list of names?
for (auto& coeffset : coeffPars_){
if (coeffset->name() == compName) {
std::cerr<<"ERROR in LauTimeDepFitModel::setAmpCoeffSet : Have already set coefficients for \""<<compName<<"\"."<<std::endl;
std::cout << "INFO in LauTimeDepFitModel::setCalibParams : Setting the initial fit parameters of the flavour tagging calibration parameters." << std::endl;
std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0bar_.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0bar_[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0bar_[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j = i; j < nSigComp_; j++) {
TString name = fitFracB0bar_[i][j].name();
name.Insert( name.Index("FitFrac"), "B0bar" );
fitFracB0bar_[i][j].name(name);
extraVars.push_back(fitFracB0bar_[i][j]);
}
}
fitFracB0_ = sigModelB0_->getFitFractions();
if (fitFracB0_.size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0_.size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
for (UInt_t i(0); i<nSigComp_; ++i) {
if (fitFracB0_[i].size() != nSigComp_) {
std::cerr<<"ERROR in LauTimeDepFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: "<<fitFracB0_[i].size()<<std::endl;
gSystem->Exit(EXIT_FAILURE);
}
}
for (UInt_t i(0); i<nSigComp_; ++i) {
for (UInt_t j = i; j < nSigComp_; j++) {
TString name = fitFracB0_[i][j].name();
name.Insert( name.Index("FitFrac"), "B0" );
fitFracB0_[i][j].name(name);
extraVars.push_back(fitFracB0_[i][j]);
}
}
// Calculate the ACPs and FitFrac asymmetries
this->calcAsymmetries(kTRUE);
// Add the Fit Fraction asymmetry for each signal component
for (UInt_t i = 0; i < nSigComp_; i++) {
extraVars.push_back(fitFracAsymm_[i]);
}
// Add the calculated CP asymmetry for each signal component
std::cerr << "ERROR in LauTimeDepFitModel::setBkgndAsymmetries : Invalid background class \"" << bkgndClass << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
// Update the number of signal events to be total-sum(background events)
this->updateSigEvents();
// Check whether we want to have randomised initial fit parameters for the signal model
if (this->useRandomInitFitPars() == kTRUE) {
this->randomiseInitFitPars();
}
}
void LauTimeDepFitModel::randomiseInitFitPars()
{
// Only randomise those parameters that are not fixed!
std::cout<<"INFO in LauTimeDepFitModel::randomiseInitFitPars : Randomising the initial values of the coefficients of the DP components (and phiMix)..."<<std::endl;
if (iEvt%1000 == 0) {std::cout<<"INFO in LauTimeDepFitModel::genExpt : Generated event number "<<iEvt<<" out of "<<nEvtsGen<<" "<<evtCategory<<" events."<<std::endl;}
}
auto t2 = std::chrono::high_resolution_clock::now();
if (signalTree_ && (signalTree_->nEvents() == signalTree_->nUsedEvents())) {
std::cerr<<"WARNING in LauTimeDepFitModel::generateSignalEvent : Source of embedded signal events used up, clearing the list of used events."<<std::endl;
// std::cerr << "WARNING in LauCPFitModel::generateBkgndEvent : Source of embedded " << bkgndClass << " events used up, clearing the list of used events." << std::endl;
// embeddedData->clearUsedList();
//}
//
LauKinematics* kinematics{nullptr};
switch ( BkgndTypes_[bkgndID] ) {
case LauFlavTag::BkgndType::Combinatorial:
{
// First choose the true tag, accounting for the production asymmetry
// CONVENTION WARNING regarding meaning of sign of AProd
// NB the true tag doesn't really mean anything for combinatorial background
Double_t random = LauRandom::randomFun()->Rndm();
if ( random <= 0.5 * ( 1.0 - AProdBkgnd_[bkgndID]->unblindValue() ) ) {
curEvtTrueTagFlv_ = LauFlavTag::Flavour::B;
} else {
curEvtTrueTagFlv_ = LauFlavTag::Flavour::Bbar;
}
if ( cpEigenValue_ == CPEigenvalue::QFS ) {
if ( BkgndDPModelsBbar_[bkgndID] != nullptr ) {
// generate the true decay flavour and the corresponding DP position
// (the supply of two DP models indicates a possible asymmetry)
std::cerr<<"ERROR in LauTimeDepFitModel::setupGenNtupleBranches : Decay flavour variable not set for QFS decay, see LauFlavTag::setDecayFlvVarName()."<<std::endl;
// TODO - any other decay time function types that make sense for combinatorial?
// - should also have a set of checks in initialise that we have everything we need for the backgrounds and that the various settings make sense
default :
// TODO as per comment just above, once the above mentioned checks are implemented this error message can be removed
std::cerr << "WARNING in LauTimeDepFitModel::getEvtDPDtLikelihood : bkgnd types other than Hist and Exp don't make sense for combinatorial!" << std::endl;
break;
}
// ...include flavour tagging
for (std::size_t taggerID{0}; taggerID<nTaggers; ++taggerID){
std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Invalid background class \"" << bkgndClass << "\"." << std::endl;
std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl;
std::cerr << "WARNING in LauTimeDepFitModel::freeSpeciesNames : \"" << par->name() << "\" is a LauFormulaPar, which implies it is perhaps not entirely free to float in the fit, so the sWeight calculation may not be reliable" << std::endl;