diff --git a/examples/ProgOpts/Test_Dpipi_ProgOpts.cc b/examples/ProgOpts/Test_Dpipi_ProgOpts.cc index dc786cf..496494a 100644 --- a/examples/ProgOpts/Test_Dpipi_ProgOpts.cc +++ b/examples/ProgOpts/Test_Dpipi_ProgOpts.cc @@ -1,165 +1,167 @@ #include #include #include "boost/program_options.hpp" #include "boost/tokenizer.hpp" #include "Test_Dpipi_ProgOpts.hh" namespace po = boost::program_options; void validate( boost::any& v, const std::vector& values, Command* /*target_type*/, int) { // Make sure no previous assignment to 'v' has been made po::validators::check_first_occurrence(v); // Extract the first string from 'values'. If there is more than // one string, it's an error, and exception will be thrown. const std::string& s { po::validators::get_single_string(values) }; if ( s == "gen" ) { v = boost::any( Command::Generate ); } else if ( s == "fit" ) { v = boost::any( Command::Fit ); } else if ( s == "simfit" ) { v = boost::any( Command::SimFit ); } else { throw po::validation_error(po::validation_error::invalid_option_value, "command", s, 3); } } void validate( boost::any& v, const std::vector& values, LauTimeDepFitModel::CPEigenvalue* /*target_type*/, int) { // Make sure no previous assignment to 'v' has been made po::validators::check_first_occurrence(v); // Extract the first string from 'values'. If there is more than // one string, it's an error, and exception will be thrown. const std::string& s { po::validators::get_single_string(values) }; if ( s == "QFS" ) { v = boost::any( LauTimeDepFitModel::CPEigenvalue::QFS ); } else if ( s == "CPEven" ) { v = boost::any( LauTimeDepFitModel::CPEigenvalue::CPEven ); } else if ( s == "CPOdd" ) { v = boost::any( LauTimeDepFitModel::CPEigenvalue::CPOdd ); } else { throw po::validation_error(po::validation_error::invalid_option_value); } } namespace LauDecayTime { void validate( boost::any& v, const std::vector& values, EfficiencyMethod* /*target_type*/, int) { // Make sure no previous assignment to 'v' has been made po::validators::check_first_occurrence(v); // Extract the first string from 'values'. If there is more than // one string, it's an error, and exception will be thrown. const std::string& s { po::validators::get_single_string(values) }; if ( s == "uniform" || s == "flat" ) { v = boost::any( EfficiencyMethod::Uniform ); } else if ( s == "binned" || s == "hist" ) { v = boost::any( EfficiencyMethod::Binned ); } else if ( s == "spline" ) { v = boost::any( EfficiencyMethod::Spline ); } else { throw po::validation_error(po::validation_error::invalid_option_value); } } }; TestDpipi_ProgramSettings::TestDpipi_ProgramSettings(const int argc, const char ** argv) { po::options_description main_desc{"Main options"}; main_desc.add_options() ("command", po::value(&command)->required(), "main command: gen, fit, or simfit") ; po::positional_options_description p; p.add("command", 1); po::options_description common_desc{"Common options"}; common_desc.add_options() ("help", "produce help message") ("dtype", po::value(&dType)->default_value(LauTimeDepFitModel::CPEigenvalue::QFS,"QFS"), "type of D decay: QFS, CPOdd, or CPEven") ("dta-model", po::value(&timeEffModel)->default_value(LauDecayTime::EfficiencyMethod::Uniform,"uniform"), "decay-time acceptance model: uniform/flat, binned/hist, spline") ("dtr", po::value(&timeResolution)->default_value(kTRUE), "enable/disable decay-time resolution") ("dtr-perevent", po::value(&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(&RNGseed)->default_value(0), "set the seed for the RNG; if not set, the time is used to generate a seed") ("dir", po::value(&directory)->default_value("",""), "set the directory used to find the nTuples within the root file. Defaults to no directory") ("bkgndList", po::value(&bkgndListStr)->default_value("comb,Bd2DKpi,Bs2DKpi,Lb2Dppi"), "the comma-separated list of background components to include") ("fitBackInput", po::value(&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") ; po::options_description gen_desc{"Generation options"}; gen_desc.add_options() ("firstExptGen", po::value(&firstExptGen)->default_value(0), "first experiment to generate") ("nExptGen", po::value(&nExptGen)->default_value(1), "number of experiments to generate") ; po::options_description fit_desc{"Fitting options"}; fit_desc.add_options() ("firstExpt", po::value(&firstExptFit)->default_value(0), "first experiment to fit") ("nExpt", po::value(&nExptFit)->default_value(1), "number of experiments to fit") ("iFit", po::value(&iFit)->default_value(0), "fit ID - used to distinguish multiple fits to a given dataset (used to disentangle multiple minima)") ("fixTau", po::value(&fixLifetime)->default_value(kTRUE), "enable/disable floating of B lifetime parameter") ("fixDm", po::value(&fixDeltaM)->default_value(kTRUE), "enable/disable floating of B mixing frequency parameter") ("fixPhiMix", po::value(&fixPhiMix)->default_value(kFALSE), "enable/disable floating of B mixing phase parameter(s)") ("fixSplineKnots", po::value(&fixSplineKnots)->default_value(kFALSE), "enable/disable floating of the decay-time-acceptance spline") ("useSinCos", po::value(&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(&useAveDeltaCalibVals)->default_value(kTRUE), "enable/disable fitting calib values as their average +/- delta values rather than having separate values for B and Bbar") ("addTaggers", po::value(&addTaggers)->default_value(kTRUE), "add/omit taggers") ("floatCalibPars", po::value(&floatCalibPars)->default_value(kTRUE), "enable/disable floating of the FT calibration parameters") ("floatYields", po::value(&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(&fitBack)->default_value(kFALSE), "enable/disable refitting some generated data from a model") ("blindFit", po::value(&blindFit)->default_value(kFALSE), "enable/disable blinding of the physics parameters") + ("dataFit", po::value(&dataFit)->default_value(kFALSE), "fit to data (not MC or toys)") + ("genToy", po::value(&genToy)->default_value(kFALSE), "enable/disable generation of a toy to a fit result") ; po::options_description simfit_desc{"Simultaneous fitting options"}; simfit_desc.add_options() ("port", po::value(&port)->default_value(0), "port number on which sim fit coordinator is listening") ; po::options_description all_desc; all_desc.add(main_desc).add(common_desc).add(gen_desc).add(fit_desc).add(simfit_desc); po::variables_map vm; try { po::store(po::command_line_parser(argc, argv). options(all_desc). positional(p). run(), vm); po::notify(vm); } catch ( boost::wrapexcept& e ) { std::cout << argv[0] << " requires a command, one of 'gen', 'fit', or 'simfit'\n" << "Append --help to those commands to see help on the related options" << std::endl; parsedOK = kFALSE; return; } catch ( po::validation_error& e ) { std::cerr << e.what() << std::endl; parsedOK = kFALSE; return; } if ( vm.count("help") ) { po::options_description help_desc; help_desc.add(common_desc).add(gen_desc).add(fit_desc).add(simfit_desc); std::cout << help_desc << std::endl; helpRequested = kTRUE; return; } if ( ! bkgndListStr.empty() ) { const boost::char_separator sep{","}; const boost::tokenizer> tokens{bkgndListStr, sep}; for (const auto& t : tokens) { bkgndList.insert( t ); } } } diff --git a/examples/ProgOpts/Test_Dpipi_ProgOpts.hh b/examples/ProgOpts/Test_Dpipi_ProgOpts.hh index e2b0bcf..8a79ae9 100644 --- a/examples/ProgOpts/Test_Dpipi_ProgOpts.hh +++ b/examples/ProgOpts/Test_Dpipi_ProgOpts.hh @@ -1,53 +1,55 @@ #ifndef TEST_DPIPI_PROGOPTS #define TEST_DPIPI_PROGOPTS #include "Rtypes.h" #include "LauDecayTime.hh" #include "LauTimeDepFitModel.hh" #include "Command.hh" #include #include struct TestDpipi_ProgramSettings { TestDpipi_ProgramSettings(const int argc, const char ** argv); Bool_t parsedOK{kTRUE}; Bool_t helpRequested{kFALSE}; Command command{Command::Generate}; LauTimeDepFitModel::CPEigenvalue dType{LauTimeDepFitModel::CPEigenvalue::QFS}; LauDecayTime::EfficiencyMethod timeEffModel{LauDecayTime::EfficiencyMethod::Uniform}; UInt_t firstExptGen{0}; UInt_t nExptGen{1}; UInt_t firstExptFit{0}; UInt_t nExptFit{1}; UInt_t iFit{0}; UInt_t port{0}; UInt_t RNGseed{0}; Bool_t timeResolution{kTRUE}; Bool_t perEventTimeErr{kFALSE}; Bool_t fixLifetime{kTRUE}; Bool_t fixDeltaM{kTRUE}; Bool_t fixPhiMix{kFALSE}; Bool_t fixSplineKnots{kFALSE}; Bool_t useSinCos{kTRUE}; Bool_t useAveDeltaCalibVals{kTRUE}; Bool_t addTaggers{kTRUE}; Bool_t floatCalibPars{kTRUE}; Bool_t floatYields{kFALSE}; Bool_t fitBack{kFALSE}; Bool_t blindFit{kFALSE}; + Bool_t dataFit{kFALSE}; + Bool_t genToy{kFALSE}; std::string bkgndListStr{""}; std::string directory{""}; std::string fitBackInput{""}; std::set bkgndList; }; #endif diff --git a/src/Lau1DHistPdf.cc b/src/Lau1DHistPdf.cc index cc54568..609ae85 100644 --- a/src/Lau1DHistPdf.cc +++ b/src/Lau1DHistPdf.cc @@ -1,268 +1,270 @@ /* Copyright 2006 University of Warwick Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. */ /* Laura++ package authors: John Back Paul Harrison Thomas Latham */ /*! \file Lau1DHistPdf.cc \brief File containing implementation of Lau1DHistPdf class. */ #include #include #include #include "TAxis.h" #include "TH1.h" #include "TRandom.h" #include "TSystem.h" #include "Lau1DHistPdf.hh" #include "LauRandom.hh" class LauParameter; ClassImp(Lau1DHistPdf) Lau1DHistPdf::Lau1DHistPdf(const TString& theVarName, const TH1* hist, Double_t minAbscissa, Double_t maxAbscissa, Bool_t useInterpolation, Bool_t fluctuateBins) : LauAbsPdf(theVarName, std::vector(), minAbscissa, maxAbscissa), hist_(hist ? dynamic_cast(hist->Clone()) : 0), useInterpolation_(useInterpolation), fluctuateBins_(fluctuateBins), nBins_(0), axisMin_(0.0), axisMax_(0.0), axisRange_(0.0) { // Constructor // Set the directory for the histogram hist_->SetDirectory(0); // Save various attributes of the histogram nBins_ = hist_->GetNbinsX(); TAxis* xAxis = hist_->GetXaxis(); axisMin_ = xAxis->GetXmin(); axisMax_ = xAxis->GetXmax(); axisRange_ = axisMax_ - axisMin_; // Check that axis range corresponds to range of abscissa if (TMath::Abs(this->getMinAbscissa() - axisMin_)>1e-6) { std::cerr << "ERROR in Lau1DHistPdf::Lau1DHistPdf : Histogram axis minimum: " << axisMin_ << + " for histogram : " << hist_->GetName() << " does not correspond to abscissa minimum: " << this->getMinAbscissa() << "." << std::endl; gSystem->Exit(EXIT_FAILURE); } if (TMath::Abs(this->getMaxAbscissa() - axisMax_)>1e-6) { std::cerr << "ERROR in Lau1DHistPdf::Lau1DHistPdf : Histogram axis maximum: " << axisMax_ << + " for histogram : " << hist_->GetName() << " does not correspond to abscissa maximum: " << this->getMaxAbscissa() << "." << std::endl; gSystem->Exit(EXIT_FAILURE); } // If the bins are to be fluctuated then do so now before // calculating anything that depends on the bin content. if (fluctuateBins) { this->doBinFluctuation(); } // Calculate the PDF normalisation. this->calcNorm(); // And check it is OK. this->checkNormalisation(); } Lau1DHistPdf::~Lau1DHistPdf() { // Destructor delete hist_; hist_ = 0; } void Lau1DHistPdf::calcPDFHeight( const LauKinematics* /*kinematics*/ ) { if (this->heightUpToDate()) { return; } // Get the maximum height of the histogram Int_t maxBin = hist_->GetMaximumBin(); Double_t height = hist_->GetBinContent(maxBin); this->setMaxHeight(height); } void Lau1DHistPdf::calcNorm() { // Calculate the histogram normalisation. // Loop over the range to get the total area. // Just sum the contributions up using 1e-3 increments of the range. // Multiply the end result by dx. Double_t dx(1e-3*axisRange_); Double_t area(0.0); Double_t x(axisMin_ + dx/2.0); while (x > axisMin_ && x < axisMax_) { area += this->interpolate(x); x += dx; } Double_t norm = area*dx; this->setNorm(norm); } void Lau1DHistPdf::checkNormalisation() { Double_t dx(1e-3*axisRange_); Double_t area(0.0); Double_t areaNoNorm(0.0); const Bool_t useInterpolationOrig { useInterpolation_ }; Double_t x(axisMin_ + dx/2.0); while (x > axisMin_ && x < axisMax_) { area += this->interpolateNorm(x); useInterpolation_ = kFALSE; areaNoNorm += this->interpolate(x); useInterpolation_ = useInterpolationOrig; x += dx; } Double_t norm = area*dx; std::cout << "INFO in Lau1DHistPdf::checkNormalisation : Area = " << area << ", dx = " << dx << std::endl; std::cout << " : Area with no norm = " << areaNoNorm << "*dx = " << areaNoNorm*dx << std::endl; std::cout << " : The total area of the normalised histogram PDF is " << norm << std::endl; } Double_t Lau1DHistPdf::getBinHistValue(Int_t bin) const { // Check that bin is in range [1 , nBins_] if ((bin < 1) || (bin > nBins_)) { return 0.0; } Double_t value = static_cast(hist_->GetBinContent(bin)); // protect against negative values if ( value < 0.0 ) { std::cerr << "WARNING in Lau1DHistPdf::getBinHistValue : Negative bin content set to zero!" << (hist_->GetName() ? Form(" Culprit histogram: %s",hist_->GetName()) : "") << std::endl; value = 0.0; } return value; } Double_t Lau1DHistPdf::interpolateNorm(Double_t x) const { // Get the normalised interpolated value. Double_t value = this->interpolate(x); Double_t norm = this->getNorm(); return value/norm; } void Lau1DHistPdf::calcLikelihoodInfo(const LauAbscissas& abscissas) { // Check that the given abscissa is within the allowed range if (!this->checkRange(abscissas)) { gSystem->Exit(EXIT_FAILURE); } // Get our abscissa Double_t abscissa = abscissas[0]; // Calculate the interpolated value Double_t value = this->interpolate(abscissa); this->setUnNormPDFVal(value); } Double_t Lau1DHistPdf::interpolate(Double_t x) const { // This function returns the interpolated value of the histogram function // for the given value of x by finding the adjacent bins and extrapolating // using weights based on the inverse distance of the point from the adajcent // bin centres. // Find the histogram bin Int_t bin = hist_->FindFixBin(x); // Ask whether we want to do the interpolation, since this method is // not reliable for low statistics histograms. if (useInterpolation_ == kFALSE) { return this->getBinHistValue(bin); } // Find the bin centres (actual co-ordinate positions, not histogram indices) Double_t cbinx = hist_->GetBinCenter(bin); // Find the adjacent bins Double_t deltax = x - cbinx; Int_t bin_adj(0); if (deltax > 0.0) { bin_adj = bin + 1; } else { bin_adj = bin - 1; } Bool_t isBoundary(kFALSE); if ( bin_adj > nBins_ || bin_adj < 1 ) { isBoundary = kTRUE; } // At the edges, do no interpolation, use entry in bin. if (isBoundary == kTRUE) { return this->getBinHistValue(bin); } // Linear interpolation using inverse distance as weights. // Find the adjacent bin centre Double_t cbinx_adj = hist_->GetBinCenter(bin_adj); Double_t deltax_adj = cbinx_adj - x; Double_t dx0 = TMath::Abs(deltax); Double_t dx1 = TMath::Abs(deltax_adj); Double_t denom = dx0 + dx1; Double_t value0 = this->getBinHistValue(bin); Double_t value1 = this->getBinHistValue(bin_adj); Double_t value = value0*dx1 + value1*dx0; value /= denom; return value; } void Lau1DHistPdf::doBinFluctuation() { TRandom* random = LauRandom::randomFun(); for (Int_t bin(0); binGetBinContent(bin+1); Double_t currentError = hist_->GetBinError(bin+1); Double_t newContent = random->Gaus(currentContent,currentError); if (newContent<0.0) { hist_->SetBinContent(bin+1,0.0); } else { hist_->SetBinContent(bin+1,newContent); } } }