diff --git a/src/LauSimpleFitModel.cc b/src/LauSimpleFitModel.cc index 48bd16f..a84509e 100644 --- a/src/LauSimpleFitModel.cc +++ b/src/LauSimpleFitModel.cc @@ -1,2356 +1,2356 @@ /* Copyright 2004 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 LauSimpleFitModel.cc \brief File containing implementation of LauSimpleFitModel class. */ #include #include #include #include #include "TFile.h" #include "TH2.h" #include "TMinuit.h" #include "TRandom.h" #include "TSystem.h" #include "TVirtualFitter.h" #include "LauAbsBkgndDPModel.hh" #include "LauAbsCoeffSet.hh" #include "LauIsobarDynamics.hh" #include "LauAbsPdf.hh" #include "LauComplex.hh" #include "LauConstants.hh" #include "LauDaughters.hh" #include "LauEmbeddedData.hh" #include "LauEffModel.hh" #include "LauFitNtuple.hh" #include "LauGenNtuple.hh" #include "LauKinematics.hh" #include "LauPrint.hh" #include "LauRandom.hh" #include "LauScfMap.hh" #include "LauSimpleFitModel.hh" #include "TGraph2D.h" #include "TGraph.h" #include "TStyle.h" #include "TCanvas.h" ClassImp(LauSimpleFitModel) LauSimpleFitModel::LauSimpleFitModel(LauIsobarDynamics* sigModel) : LauAbsFitModel(), sigDPModel_(sigModel), kinematics_(sigModel ? sigModel->getKinematics() : 0), usingBkgnd_(kFALSE), nSigComp_(0), nSigDPPar_(0), nExtraPdfPar_(0), nNormPar_(0), meanEff_("meanEff",0.0,0.0,1.0), dpRate_("dpRate",0.0,0.0,100.0), //signalEvents_("signalEvents",1.0,0.0,1.0), signalEvents_(0), useSCF_(kFALSE), useSCFHist_(kFALSE), scfFrac_("scfFrac",0.0,0.0,1.0), scfFracHist_(0), scfMap_(0), compareFitData_(kFALSE), signalTree_(0), reuseSignal_(kFALSE), useReweighting_(kFALSE), sigDPLike_(0.0), scfDPLike_(0.0), sigExtraLike_(0.0), scfExtraLike_(0.0), sigTotalLike_(0.0), scfTotalLike_(0.0) { } LauSimpleFitModel::~LauSimpleFitModel() { delete signalTree_; for (LauBkgndEmbDataList::iterator iter = bkgndTree_.begin(); iter != bkgndTree_.end(); ++iter) { delete (*iter); } delete scfFracHist_; } void LauSimpleFitModel::setupBkgndVectors() { UInt_t nBkgnds = this->nBkgndClasses(); bkgndDPModels_.resize( nBkgnds ); bkgndPdfs_.resize( nBkgnds ); bkgndEvents_.resize( nBkgnds ); bkgndTree_.resize( nBkgnds ); reuseBkgnd_.resize( nBkgnds ); bkgndDPLike_.resize( nBkgnds ); bkgndExtraLike_.resize( nBkgnds ); bkgndTotalLike_.resize( nBkgnds ); } void LauSimpleFitModel::setNSigEvents(LauParameter* nSigEvents) { if ( nSigEvents == 0 ) { std::cerr << "ERROR in LauSimpleFitModel::setNSigEvents : The signal yield LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( signalEvents_ != 0 ) { std::cerr << "ERROR in LauSimpleFitModel::setNSigEvents : You are trying to overwrite the signal yield." << std::endl; return; } signalEvents_ = nSigEvents; TString name = signalEvents_->name(); if ( ! name.Contains("signalEvents") && !( name.BeginsWith("signal") && name.EndsWith("Events") ) ) { signalEvents_->name("signalEvents"); } Double_t value = nSigEvents->value(); signalEvents_->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0)); } void LauSimpleFitModel::setNBkgndEvents(LauAbsRValue* nBkgndEvents) { if ( nBkgndEvents == 0 ) { std::cerr << "ERROR in LauSimpleFitModel::setNBkgndEvents : The background yield LauParameter pointer is null." << std::endl; gSystem->Exit(EXIT_FAILURE); } if ( ! this->validBkgndClass( nBkgndEvents->name() ) ) { std::cerr << "ERROR in LauSimpleFitModel::setNBkgndEvents : Invalid background class \"" << nBkgndEvents->name() << "\"." << std::endl; std::cerr << " : Background class names must be provided in \"setBkgndClassNames\" before any other background-related actions can be performed." << std::endl; gSystem->Exit(EXIT_FAILURE); } UInt_t bkgndID = this->bkgndClassID( nBkgndEvents->name() ); if ( bkgndEvents_[bkgndID] != 0 ) { std::cerr << "ERROR in LauSimpleFitModel::setNBkgndEvents : You are trying to overwrite the background yield." << std::endl; return; } nBkgndEvents->name( nBkgndEvents->name()+"Events" ); if ( nBkgndEvents->isLValue() ) { Double_t value = nBkgndEvents->value(); LauParameter* yield = dynamic_cast( nBkgndEvents ); yield->range(-2.0*(TMath::Abs(value)+1.0), 2.0*(TMath::Abs(value)+1.0)); } bkgndEvents_[bkgndID] = nBkgndEvents; } void LauSimpleFitModel::splitSignalComponent( const TH2* dpHisto, const Bool_t upperHalf, const Bool_t fluctuateBins, LauScfMap* scfMap ) { if ( useSCF_ == kTRUE ) { std::cerr << "ERROR in LauSimpleFitModel::splitSignalComponent : Have already setup SCF." << std::endl; return; } if ( dpHisto == 0 ) { std::cerr << "ERROR in LauSimpleFitModel::splitSignalComponent : The histogram pointer is null." << std::endl; return; } const LauDaughters* daughters = sigDPModel_->getDaughters(); scfFracHist_ = new LauEffModel( daughters, 0 ); scfFracHist_->setEffHisto( dpHisto, kTRUE, fluctuateBins, 0.0, 0.0, upperHalf, daughters->squareDP() ); scfMap_ = scfMap; useSCF_ = kTRUE; useSCFHist_ = kTRUE; } void LauSimpleFitModel::splitSignalComponent( const Double_t scfFrac, const Bool_t fixed ) { if ( useSCF_ == kTRUE ) { std::cerr << "ERROR in LauSimpleFitModel::splitSignalComponent : Have already setup SCF." << std::endl; return; } scfFrac_.range( 0.0, 1.0 ); scfFrac_.value( scfFrac ); scfFrac_.initValue( scfFrac ); scfFrac_.genValue( scfFrac ); scfFrac_.fixed( fixed ); useSCF_ = kTRUE; useSCFHist_ = kFALSE; } void LauSimpleFitModel::setBkgndDPModel(const TString& bkgndClass, LauAbsBkgndDPModel* bkgndDPModel) { if (bkgndDPModel == 0) { std::cerr << "ERROR in LauSimpleFitModel::setBkgndDPModel : The model pointer is null." << std::endl; return; } // check that this background name is valid if ( ! this->validBkgndClass( bkgndClass ) ) { std::cerr << "ERROR in LauSimpleFitModel::setBkgndDPModel : 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; return; } UInt_t bkgndID = this->bkgndClassID( bkgndClass ); bkgndDPModels_[bkgndID] = bkgndDPModel; usingBkgnd_ = kTRUE; } void LauSimpleFitModel::setSignalPdf(LauAbsPdf* pdf) { if (pdf==0) { std::cerr << "ERROR in LauSimpleFitModel::setSignalPdf : The PDF pointer is null." << std::endl; return; } signalPdfs_.push_back(pdf); } void LauSimpleFitModel::setSCFPdf(LauAbsPdf* pdf) { if (pdf==0) { std::cerr << "ERROR in LauSimpleFitModel::setSCFPdf : The PDF pointer is null." << std::endl; return; } scfPdfs_.push_back(pdf); } void LauSimpleFitModel::setBkgndPdf(const TString& bkgndClass, LauAbsPdf* pdf) { if (pdf == 0) { std::cerr << "ERROR in LauSimpleFitModel::setBkgndPdf : The PDF pointer is null." << std::endl; return; } // check that this background name is valid if ( ! this->validBkgndClass( bkgndClass ) ) { std::cerr << "ERROR in LauSimpleFitModel::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; return; } UInt_t bkgndID = this->bkgndClassID( bkgndClass ); bkgndPdfs_[bkgndID].push_back(pdf); usingBkgnd_ = kTRUE; } void LauSimpleFitModel::setAmpCoeffSet(LauAbsCoeffSet* coeffSet) { // Resize the coeffPars vector if not already done if ( coeffPars_.empty() ) { coeffPars_.resize( sigDPModel_->getnTotAmp() ); for (std::vector::iterator iter = coeffPars_.begin(); iter != coeffPars_.end(); ++iter) { (*iter) = 0; } } // Is there a component called compName in the signal model? TString compName(coeffSet->name()); const Int_t index = sigDPModel_->resonanceIndex(compName); if ( index < 0 ) { std::cerr << "ERROR in LauSimpleFitModel::setAmpCoeffSet : Signal DP model doesn't contain component \"" << compName << "\"." << std::endl; return; } // Do we already have it in our list of names? if ( coeffPars_[index] != 0 && coeffPars_[index]->name() == compName) { std::cerr << "ERROR in LauSimpleFitModel::setAmpCoeffSet : Have already set coefficients for \"" << compName << "\"." << std::endl; return; } coeffSet->index(index); coeffPars_[index] = coeffSet; ++nSigComp_; std::cout << "INFO in LauSimpleFitModel::setAmpCoeffSet : Added coefficients for component A"<< index << ": \"" << compName << "\" to the fit model." << std::endl; coeffSet->printParValues(); } void LauSimpleFitModel::initialise() { // Initialisation if (!this->useDP() && signalPdfs_.empty()) { std::cerr << "ERROR in LauSimpleFitModel::initialise : Signal model doesn't exist for any variable." << std::endl; gSystem->Exit(EXIT_FAILURE); } if (this->useDP()) { // Check that we have all the Dalitz-plot models if ( sigDPModel_ == 0 ) { std::cerr << "ERROR in LauSimpleFitModel::initialise : The pointer to the signal DP model is null.\n"; std::cerr << " : Removing the Dalitz Plot from the model." << std::endl; this->useDP(kFALSE); } if ( usingBkgnd_ ) { for (LauBkgndDPModelList::const_iterator dpmodel_iter = bkgndDPModels_.begin(); dpmodel_iter != bkgndDPModels_.end(); ++dpmodel_iter ) { if ( (*dpmodel_iter) == 0 ) { std::cerr << "ERROR in LauSimpleFitModel::initialise : The pointer to one of the background DP models is null.\n"; std::cerr << " : Removing the Dalitz Plot from the model." << std::endl; this->useDP(kFALSE); break; } } } // Need to check that the number of components we have and that the dynamics has matches up UInt_t nAmp = sigDPModel_->getnTotAmp(); if (nAmp != nSigComp_) { std::cerr << "ERROR in LauSimpleFitModel::initialise : Number of signal DP components with magnitude and phase set not right." << std::endl; gSystem->Exit(EXIT_FAILURE); } // From the initial parameter values calculate the coefficients // so they can be passed to the signal model this->updateCoeffs(); // If all is well, go ahead and initialise them this->initialiseDPModels(); } // Next check that, if a given component is being used we've got the // right number of PDFs for all the variables involved // TODO - should probably check variable names and so on as well UInt_t nsigpdfvars(0); for ( LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter ) { std::vector varNames = (*pdf_iter)->varNames(); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { ++nsigpdfvars; } } } if (useSCF_) { UInt_t nscfpdfvars(0); for ( LauPdfList::const_iterator pdf_iter = scfPdfs_.begin(); pdf_iter != scfPdfs_.end(); ++pdf_iter ) { std::vector varNames = (*pdf_iter)->varNames(); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { ++nscfpdfvars; } } } if (nscfpdfvars != nsigpdfvars) { std::cerr << "ERROR in LauSimpleFitModel::initialise : There are " << nsigpdfvars << " TM signal PDF variables but " << nscfpdfvars << " SCF signal PDF variables." << std::endl; gSystem->Exit(EXIT_FAILURE); } } if (usingBkgnd_) { for (LauBkgndPdfsList::const_iterator bgclass_iter = bkgndPdfs_.begin(); bgclass_iter != bkgndPdfs_.end(); ++bgclass_iter) { UInt_t nbkgndpdfvars(0); const LauPdfList& pdfList = (*bgclass_iter); for ( LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter ) { std::vector varNames = (*pdf_iter)->varNames(); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { ++nbkgndpdfvars; } } } if (nbkgndpdfvars != nsigpdfvars) { std::cerr << "ERROR in LauSimpleFitModel::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 various extra PDFs this->setExtraPdfParameters(); // Set the initial bg and signal events this->setFitNEvents(); // Check that we have the expected number of fit variables const LauParameterPList& fitVars = this->fitPars(); if (fitVars.size() != (nSigDPPar_ + nExtraPdfPar_ + nNormPar_)) { std::cerr << "ERROR in LauSimpleFitModel::initialise : Number of fit parameters not of expected size." << std::endl; gSystem->Exit(EXIT_FAILURE); } this->setExtraNtupleVars(); } void LauSimpleFitModel::initialiseDPModels() { std::cout << "INFO in LauSimpleFitModel::initialiseDPModels : Initialising DP models" << std::endl; sigDPModel_->initialise(coeffs_); if (usingBkgnd_ == kTRUE) { for (LauBkgndDPModelList::iterator iter = bkgndDPModels_.begin(); iter != bkgndDPModels_.end(); ++iter) { (*iter)->initialise(); } } } void LauSimpleFitModel::recalculateNormalisation() { //std::cout << "INFO in LauSimpleFitModel::recalculateNormalizationInDPModels : Recalc Norm in DP model" << std::endl; sigDPModel_->recalculateNormalisation(); sigDPModel_->modifyDataTree(); } void LauSimpleFitModel::setSignalDPParameters() { // Set the fit parameters for the signal model. nSigDPPar_ = 0; if ( ! this->useDP() ) { return; } std::cout << "INFO in LauSimpleFitModel::setSignalDPParameters : Setting the initial fit parameters for the signal DP model." << std::endl; // Place isobar coefficient parameters in vector of fit variables LauParameterPList& fitVars = this->fitPars(); for (UInt_t i = 0; i < nSigComp_; i++) { LauParameterPList pars = coeffPars_[i]->getParameters(); for (LauParameterPList::iterator iter = pars.begin(); iter != pars.end(); ++iter) { if ( !(*iter)->clone() ) { fitVars.push_back(*iter); ++nSigDPPar_; } } } // Obtain the resonance parameters and place them in the vector of fit variables and in a separate vector LauParameterPSet& resVars = this->resPars(); resVars.clear(); LauParameterPList& sigDPPars = sigDPModel_->getFloatingParameters(); for ( LauParameterPList::iterator iter = sigDPPars.begin(); iter != sigDPPars.end(); ++iter ) { if ( resVars.insert(*iter).second ) { fitVars.push_back(*iter); ++nSigDPPar_; } } } void LauSimpleFitModel::setExtraPdfParameters() { // Include all the parameters of the various PDFs in the fit // NB all of them are passed to the fit, even though some have been fixed through parameter.fixed(kTRUE) // With the new "cloned parameter" scheme only "original" parameters are passed to the fit. // Their clones are updated automatically when the originals are updated. nExtraPdfPar_ = 0; nExtraPdfPar_ += this->addFitParameters(signalPdfs_); if (useSCF_ == kTRUE) { nExtraPdfPar_ += this->addFitParameters(scfPdfs_); } if (usingBkgnd_ == kTRUE) { for (LauBkgndPdfsList::iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) { nExtraPdfPar_ += this->addFitParameters(*iter); } } } void LauSimpleFitModel::setFitNEvents() { if ( signalEvents_ == 0 ) { std::cerr << "ERROR in LauSimpleFitModel::setFitNEvents : Signal yield not defined." << std::endl; return; } nNormPar_ = 0; // initialise the total number of events to be the sum of all the hypotheses Double_t nTotEvts = signalEvents_->value(); for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) { nTotEvts += (*iter)->value(); if ( (*iter) == 0 ) { std::cerr << "ERROR in LauSimpleFitModel::setFitNEvents : Background yield not defined." << std::endl; return; } } this->eventsPerExpt(TMath::FloorNint(nTotEvts)); LauParameterPList& fitVars = this->fitPars(); // if doing an extended ML fit add the number of signal events into the fit parameters if (this->doEMLFit()) { std::cout << "INFO in LauSimpleFitModel::setFitNEvents : Initialising number of events for signal and background components..." << std::endl; // add the signal events to the list of fit parameters fitVars.push_back(signalEvents_); ++nNormPar_; } else { std::cout << "INFO in LauSimpleFitModel::setFitNEvents : Initialising number of events for background components (and hence signal)..." << std::endl; } if (useSCF_ && !useSCFHist_) { fitVars.push_back(&scfFrac_); ++nNormPar_; } if (usingBkgnd_ == kTRUE) { for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) { std::vector parameters = (*iter)->getPars(); for ( LauParameter* parameter : parameters ) { if ( ! parameter->clone() ) { fitVars.push_back(parameter); ++nNormPar_; } } } } } void LauSimpleFitModel::setExtraNtupleVars() { // Set-up other parameters derived from the fit results, e.g. fit fractions. if (this->useDP() != kTRUE) { return; } // First clear the vectors so we start from scratch this->clearExtraVarVectors(); LauParameterList& extraVars = this->extraPars(); // Add a fit fraction for each signal component fitFrac_ = sigDPModel_->getFitFractions(); if (fitFrac_.size() != nSigComp_) { std::cerr << "ERROR in LauSimpleFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << fitFrac_.size() << std::endl; gSystem->Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } // Add the fit fraction that has not been corrected for the efficiency for each signal component fitFracEffUnCorr_ = sigDPModel_->getFitFractionsEfficiencyUncorrected(); if (fitFracEffUnCorr_.size() != nSigComp_) { std::cerr << "ERROR in LauSimpleFitModel::setExtraNtupleVars : Initial Fit Fraction array of unexpected dimension: " << fitFracEffUnCorr_.size() << std::endl; gSystem->Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } for (UInt_t i = 0; i < nSigComp_; i++) { for (UInt_t j = i; j < nSigComp_; j++) { extraVars.push_back(fitFrac_[i][j]); extraVars.push_back(fitFracEffUnCorr_[i][j]); } } // Store any extra parameters/quantities from the DP model (e.g. K-matrix total fit fractions) std::vector extraParams = sigDPModel_->getExtraParameters(); std::vector::iterator extraIter; for (extraIter = extraParams.begin(); extraIter != extraParams.end(); ++extraIter) { LauParameter extraParameter = (*extraIter); extraVars.push_back(extraParameter); } // Now add in the DP efficiency value Double_t initMeanEff = sigDPModel_->getMeanEff().initValue(); meanEff_.value(initMeanEff); meanEff_.initValue(initMeanEff); meanEff_.genValue(initMeanEff); extraVars.push_back(meanEff_); // Also add in the DP rate Double_t initDPRate = sigDPModel_->getDPRate().initValue(); dpRate_.value(initDPRate); dpRate_.initValue(initDPRate); dpRate_.genValue(initDPRate); extraVars.push_back(dpRate_); } void LauSimpleFitModel::finaliseFitResults(const TString& tablePrefixName) { // Retrieve parameters from the fit results for calculations and toy generation // and eventually store these in output root ntuples/text files // Now take the fit parameters and update them as necessary // e.g. to make mag > 0.0 and phase in the right range. // This function will also calculate any other values, such as the // fit fractions, using any errors provided by fitParErrors as appropriate. // Also obtain the pull values: (measured - generated)/(average error) if (this->useDP() == kTRUE) { for (UInt_t i = 0; i < nSigComp_; i++) { // Check whether we have mag > 0.0, and phase in the right range coeffPars_[i]->finaliseValues(); } } // update the pulls on the events if (this->doEMLFit()) { signalEvents_->updatePull(); } if (useSCF_ && !useSCFHist_) { scfFrac_.updatePull(); } if (usingBkgnd_ == kTRUE) { for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) { std::vector parameters = (*iter)->getPars(); for ( LauParameter* parameter : parameters ) { parameter->updatePull(); } } } // Update the pulls on all the extra PDFs' parameters this->updateFitParameters(signalPdfs_); if (useSCF_ == kTRUE) { this->updateFitParameters(scfPdfs_); } if (usingBkgnd_ == kTRUE) { for (LauBkgndPdfsList::iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) { this->updateFitParameters(*iter); } } // Fill the fit results to the ntuple for current experiment // update the coefficients and then calculate the fit fractions if (this->useDP() == kTRUE) { this->updateCoeffs(); sigDPModel_->updateCoeffs(coeffs_); sigDPModel_->calcExtraInfo(); LauParArray fitFrac = sigDPModel_->getFitFractions(); if (fitFrac.size() != nSigComp_) { std::cerr << "ERROR in LauSimpleFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << fitFrac.size() << std::endl; gSystem->Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } LauParArray fitFracEffUnCorr = sigDPModel_->getFitFractionsEfficiencyUncorrected(); if (fitFracEffUnCorr.size() != nSigComp_) { std::cerr << "ERROR in LauSimpleFitModel::finaliseFitResults : Fit Fraction array of unexpected dimension: " << fitFracEffUnCorr.size() << std::endl; gSystem->Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } for (UInt_t i(0); igetMeanEff().value()); dpRate_.value(sigDPModel_->getDPRate().value()); this->clearExtraVarVectors(); LauParameterList& extraVars = this->extraPars(); // Then store the final fit parameters, and any extra parameters for // the signal model (e.g. fit fractions) for (UInt_t i(0); i extraParams = sigDPModel_->getExtraParameters(); std::vector::iterator extraIter; for (extraIter = extraParams.begin(); extraIter != extraParams.end(); ++extraIter) { LauParameter extraParameter = (*extraIter); extraVars.push_back(extraParameter); } // Now add in the DP efficiency value extraVars.push_back(meanEff_); // Also add in the DP rate extraVars.push_back(dpRate_); this->printFitFractions(std::cout); } const LauParameterPList& fitVars = this->fitPars(); const LauParameterList& extraVars = this->extraPars(); LauFitNtuple* ntuple = this->fitNtuple(); ntuple->storeParsAndErrors(fitVars, extraVars); // find out the correlation matrix for the parameters ntuple->storeCorrMatrix(this->iExpt(), this->fitStatus(), this->covarianceMatrix()); // Fill the data into ntuple ntuple->updateFitNtuple(); // Print out the partial fit fractions, phases and the // averaged efficiency, reweighted by the dynamics (and anything else) if (this->writeLatexTable()) { TString sigOutFileName(tablePrefixName); sigOutFileName += "_"; sigOutFileName += this->iExpt(); sigOutFileName += "Expt.tex"; this->writeOutTable(sigOutFileName); } } void LauSimpleFitModel::printFitFractions(std::ostream& output) { // Print out Fit Fractions, total DP rate and mean efficiency for (UInt_t i = 0; i < nSigComp_; i++) { output << "FitFraction for component " << i << " (" << coeffPars_[i]->name() << ") = " << fitFrac_[i][i] << std::endl; } output << "Overall DP rate (integral of matrix element squared) = " << dpRate_ << std::endl; output << "Average efficiency weighted by whole DP dynamics = " << meanEff_ << std::endl; } void LauSimpleFitModel::writeOutTable(const TString& outputFile) { // Write out the results of the fit to a tex-readable table // TODO - need to include the yields in this table std::ofstream fout(outputFile); LauPrint print; std::cout << "INFO in LauSimpleFitModel::writeOutTable : Writing out results of the fit to the tex file " << outputFile << std::endl; if (this->useDP() == kTRUE) { // print the fit coefficients in one table coeffPars_.front()->printTableHeading(fout); for (UInt_t i = 0; i < nSigComp_; i++) { coeffPars_[i]->printTableRow(fout); } fout << "\\hline" << std::endl; fout << "\\end{tabular}" << std::endl << std::endl; // print the fit fractions in another fout << "\\begin{tabular}{|l|c|}" << std::endl; fout << "\\hline" << std::endl; fout << "Component & FitFraction \\\\" << std::endl; fout << "\\hline" << std::endl; Double_t fitFracSum(0.0); for (UInt_t i = 0; i < nSigComp_; i++) { TString resName = coeffPars_[i]->name(); resName = resName.ReplaceAll("_", "\\_"); fout << resName << " & $"; Double_t fitFrac = fitFrac_[i][i].value(); fitFracSum += fitFrac; print.printFormat(fout, fitFrac); fout << "$ \\\\" << std::endl; } fout << "\\hline" << std::endl; // Also print out sum of fit fractions fout << "Fit Fraction Sum & $"; print.printFormat(fout, fitFracSum); fout << "$ \\\\" << std::endl; fout << "\\hline" << std::endl; fout << "DP rate & $"; print.printFormat(fout, dpRate_.value()); fout << "$ \\\\" << std::endl; fout << "\\hline" << std::endl; fout << "$< \\varepsilon >$ & $"; print.printFormat(fout, meanEff_.value()); fout << "$ \\\\" << std::endl; fout << "\\hline" << std::endl; fout << "\\end{tabular}" << std::endl << std::endl; } if (!signalPdfs_.empty()) { fout << "\\begin{tabular}{|l|c|}" << std::endl; fout << "\\hline" << std::endl; if (useSCF_ == kTRUE) { fout << "\\Extra TM Signal PDFs' Parameters: & \\\\" << std::endl; } else { fout << "\\Extra Signal PDFs' Parameters: & \\\\" << std::endl; } this->printFitParameters(signalPdfs_, fout); if (useSCF_ == kTRUE && !scfPdfs_.empty()) { fout << "\\hline" << std::endl; fout << "\\Extra SCF Signal PDFs' Parameters: & \\\\" << std::endl; this->printFitParameters(scfPdfs_, fout); } if (usingBkgnd_ == kTRUE && !bkgndPdfs_.empty()) { fout << "\\hline" << std::endl; fout << "\\Extra Background PDFs' Parameters: & \\\\" << std::endl; for (LauBkgndPdfsList::const_iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) { this->printFitParameters(*iter, fout); } } fout << "\\hline \n\\end{tabular}" << std::endl << std::endl; } } void LauSimpleFitModel::checkInitFitParams() { // 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) { std::cout << "INFO in LauSimpleFitModel::checkInitFitParams : Setting random parameters for the signal DP model" << std::endl; this->randomiseInitFitPars(); } } void LauSimpleFitModel::randomiseInitFitPars() { // Only randomise those parameters that are not fixed! std::cout << "INFO in LauSimpleFitModel::randomiseInitFitPars : Randomising the initial fit magnitudes and phases of the resonances..." << std::endl; for (UInt_t i = 0; i < nSigComp_; i++) { coeffPars_[i]->randomiseInitValues(); } } std::pair LauSimpleFitModel::eventsToGenerate() { // Determine the number of events to generate for each hypothesis // If we're smearing then smear each one individually LauGenInfo nEvtsGen; // Keep track of whether any yield or asymmetry parameters are blinded Bool_t blind = kFALSE; // Signal if ( signalEvents_->blind() ) { blind = kTRUE; } Double_t evtWeight(1.0); Double_t nEvts = signalEvents_->genValue(); if ( nEvts < 0 ) { evtWeight = -1.0; nEvts = TMath::Abs( nEvts ); } Int_t nEvtsToGen { static_cast(nEvts) }; if (this->doPoissonSmearing()) { nEvtsToGen = LauRandom::randomFun()->Poisson(nEvts); } nEvtsGen["signal"] = std::make_pair( nEvtsToGen, evtWeight ); // Backgrounds const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { const LauAbsRValue* evtsPar = bkgndEvents_[bkgndID]; if ( evtsPar->blind() ) { blind = kTRUE; } evtWeight = 1.0; nEvts = evtsPar->genValue(); if ( nEvts < 0 ) { evtWeight = -1.0; nEvts = TMath::Abs( nEvts ); } nEvtsToGen = static_cast(nEvts); if (this->doPoissonSmearing()) { nEvtsToGen = LauRandom::randomFun()->Poisson(nEvts); } const TString& bkgndClass = this->bkgndClassName(bkgndID); - nEvtsGen[bkgndClass] = std::make_pair( nEvts, evtWeight ); + nEvtsGen[bkgndClass] = std::make_pair( nEvtsToGen, evtWeight ); } return std::make_pair( nEvtsGen, blind ); } Bool_t LauSimpleFitModel::genExpt() { // Routine to generate toy Monte Carlo events according to the various models we have defined. // Determine the number of events to generate for each hypothesis std::pair info = this->eventsToGenerate(); LauGenInfo nEvts = info.first; const Bool_t blind = info.second; Bool_t genOK(kTRUE); Int_t evtNum(0); const UInt_t nBkgnds = this->nBkgndClasses(); std::vector bkgndClassNames(nBkgnds); std::vector bkgndClassNamesGen(nBkgnds); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { TString name( this->bkgndClassName(iBkgnd) ); bkgndClassNames[iBkgnd] = name; bkgndClassNamesGen[iBkgnd] = "gen"+name; } const Bool_t storeSCFTruthInfo = ( useSCF_ || ( this->enableEmbedding() && signalTree_ != 0 && signalTree_->haveBranch("mcMatch") ) ); // Loop over the hypotheses and generate the requested number of events for each one for (LauGenInfo::const_iterator iter = nEvts.begin(); iter != nEvts.end(); ++iter) { const TString& type(iter->first); Int_t nEvtsGen( iter->second.first ); Double_t evtWeight( iter->second.second ); for (Int_t iEvt(0); iEvtsetGenNtupleDoubleBranchValue( "evtWeight", evtWeight ); // Add efficiency information this->setGenNtupleDoubleBranchValue( "efficiency", 1 ); if (type == "signal") { this->setGenNtupleIntegerBranchValue("genSig",1); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], 0 ); } genOK = this->generateSignalEvent(); this->setGenNtupleDoubleBranchValue( "efficiency", sigDPModel_->getEvtEff() ); } else { this->setGenNtupleIntegerBranchValue("genSig",0); if ( storeSCFTruthInfo ) { this->setGenNtupleIntegerBranchValue("genTMSig",0); this->setGenNtupleIntegerBranchValue("genSCFSig",0); } UInt_t bkgndID(0); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { Int_t gen(0); if ( bkgndClassNames[iBkgnd] == type ) { gen = 1; bkgndID = iBkgnd; } this->setGenNtupleIntegerBranchValue( bkgndClassNamesGen[iBkgnd], gen ); } genOK = this->generateBkgndEvent(bkgndID); } if (!genOK) { // If there was a problem with the generation then break out and return. // The problem model will have adjusted itself so that all should be OK next time. break; } if (this->useDP() == kTRUE) { this->setDPBranchValues(); } // Store the event number (within this experiment) this->setGenNtupleIntegerBranchValue("iEvtWithinExpt",evtNum); // and then increment it ++evtNum; this->fillGenNtupleBranches(); if ( !blind && (iEvt%500 == 0) ) { std::cout << "INFO in LauSimpleFitModel::genExpt : Generated event number " << iEvt << " out of " << nEvtsGen << " " << type << " events." << std::endl; } } if (!genOK) { break; } } if (this->useDP() && genOK) { sigDPModel_->checkToyMC(kTRUE,kTRUE); // Get the fit fractions if they're to be written into the latex table if (this->writeLatexTable()) { LauParArray fitFrac = sigDPModel_->getFitFractions(); if (fitFrac.size() != nSigComp_) { std::cerr << "ERROR in LauSimpleFitModel::genExpt : Fit Fraction array of unexpected dimension: " << fitFrac.size() << std::endl; gSystem->Exit(EXIT_FAILURE); } for (UInt_t i(0); iExit(EXIT_FAILURE); } } for (UInt_t i(0); igetMeanEff().value()); dpRate_.value(sigDPModel_->getDPRate().value()); } } // If we're reusing embedded events or if the generation is being // reset then clear the lists of used events if (reuseSignal_ || !genOK) { if (signalTree_) { signalTree_->clearUsedList(); } } for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { LauEmbeddedData* data = bkgndTree_[bkgndID]; if (reuseBkgnd_[bkgndID] || !genOK) { if (data) { data->clearUsedList(); } } } return genOK; } Bool_t LauSimpleFitModel::generateSignalEvent() { // Generate signal event Bool_t genOK(kTRUE); Bool_t genSCF(kFALSE); if (this->useDP()) { if (this->enableEmbedding() && signalTree_) { if (useReweighting_) { // Select a (random) event from the generated data. Then store the // reconstructed DP co-ords, together with other pdf information, // as the embedded data. genOK = signalTree_->getReweightedEvent(sigDPModel_); } else { signalTree_->getEmbeddedEvent(kinematics_); } if (signalTree_->haveBranch("mcMatch")) { Int_t match = static_cast(signalTree_->getValue("mcMatch")); if (match) { this->setGenNtupleIntegerBranchValue("genTMSig",1); this->setGenNtupleIntegerBranchValue("genSCFSig",0); genSCF = kFALSE; } else { this->setGenNtupleIntegerBranchValue("genTMSig",0); this->setGenNtupleIntegerBranchValue("genSCFSig",1); genSCF = kTRUE; } } } else { genOK = sigDPModel_->generate(); if ( genOK && useSCF_ ) { Double_t frac(0.0); if ( useSCFHist_ ) { frac = scfFracHist_->calcEfficiency( kinematics_ ); } else { frac = scfFrac_.genValue(); } if ( frac < LauRandom::randomFun()->Rndm() ) { this->setGenNtupleIntegerBranchValue("genTMSig",1); this->setGenNtupleIntegerBranchValue("genSCFSig",0); genSCF = kFALSE; } else { this->setGenNtupleIntegerBranchValue("genTMSig",0); this->setGenNtupleIntegerBranchValue("genSCFSig",1); genSCF = kTRUE; // Optionally smear the DP position // of the SCF event if ( scfMap_ != 0 ) { Double_t xCoord(0.0), yCoord(0.0); if ( kinematics_->squareDP() ) { xCoord = kinematics_->getmPrime(); yCoord = kinematics_->getThetaPrime(); } else { xCoord = kinematics_->getm13Sq(); yCoord = kinematics_->getm23Sq(); } // Find the bin number where this event is generated Int_t binNo = scfMap_->binNumber( xCoord, yCoord ); // Retrieve the migration histogram TH2* histo = scfMap_->trueHist( binNo ); const LauAbsEffModel * effModel = sigDPModel_->getEffModel(); do { // Get a random point from the histogram histo->GetRandom2( xCoord, yCoord ); // Update the kinematics if ( kinematics_->squareDP() ) { kinematics_->updateSqDPKinematics( xCoord, yCoord ); } else { kinematics_->updateKinematics( xCoord, yCoord ); } } while ( ! effModel->passVeto( kinematics_ ) ); } } } } } else { if (this->enableEmbedding() && signalTree_) { signalTree_->getEmbeddedEvent(0); if (signalTree_->haveBranch("mcMatch")) { Int_t match = static_cast(signalTree_->getValue("mcMatch")); if (match) { this->setGenNtupleIntegerBranchValue("genTMSig",1); this->setGenNtupleIntegerBranchValue("genSCFSig",0); genSCF = kFALSE; } else { this->setGenNtupleIntegerBranchValue("genTMSig",0); this->setGenNtupleIntegerBranchValue("genSCFSig",1); genSCF = kTRUE; } } } else if ( useSCF_ ) { Double_t frac = scfFrac_.genValue(); if ( frac < LauRandom::randomFun()->Rndm() ) { this->setGenNtupleIntegerBranchValue("genTMSig",1); this->setGenNtupleIntegerBranchValue("genSCFSig",0); genSCF = kFALSE; } else { this->setGenNtupleIntegerBranchValue("genTMSig",0); this->setGenNtupleIntegerBranchValue("genSCFSig",1); genSCF = kTRUE; } } } if (genOK) { if ( useSCF_ ) { if ( genSCF ) { this->generateExtraPdfValues(&scfPdfs_, signalTree_); } else { this->generateExtraPdfValues(&signalPdfs_, signalTree_); } } else { this->generateExtraPdfValues(&signalPdfs_, signalTree_); } } // Check for problems with the embedding if (this->enableEmbedding() && signalTree_ && (signalTree_->nEvents() == signalTree_->nUsedEvents())) { std::cerr << "WARNING in LauSimpleFitModel::generateSignalEvent : Source of embedded signal events used up, clearing the list of used events." << std::endl; signalTree_->clearUsedList(); } return genOK; } Bool_t LauSimpleFitModel::generateBkgndEvent(UInt_t bkgndID) { // Generate background event Bool_t genOK(kTRUE); LauAbsBkgndDPModel* model = bkgndDPModels_[bkgndID]; LauEmbeddedData* embeddedData(0); if (this->enableEmbedding()) { embeddedData = bkgndTree_[bkgndID]; } LauPdfList* extraPdfs = &bkgndPdfs_[bkgndID]; if (this->useDP()) { if (embeddedData) { embeddedData->getEmbeddedEvent(kinematics_); } else { if (model == 0) { std::cerr << "ERROR in LauSimpleFitModel::generateBkgndEvent : Can't find the DP model for background class \"" << this->bkgndClassName(bkgndID) << "\"." << std::endl; gSystem->Exit(EXIT_FAILURE); } genOK = model->generate(); } } else { if (embeddedData) { embeddedData->getEmbeddedEvent(0); } } if (genOK) { this->generateExtraPdfValues(extraPdfs, embeddedData); } // Check for problems with the embedding if (embeddedData && (embeddedData->nEvents() == embeddedData->nUsedEvents())) { std::cerr << "WARNING in LauSimpleFitModel::generateBkgndEvent : Source of embedded " << this->bkgndClassName(bkgndID) << " events used up, clearing the list of used events." << std::endl; embeddedData->clearUsedList(); } return genOK; } void LauSimpleFitModel::setupGenNtupleBranches() { // Setup the required ntuple branches this->addGenNtupleDoubleBranch("evtWeight"); this->addGenNtupleIntegerBranch("genSig"); this->addGenNtupleDoubleBranch("efficiency"); if ( useSCF_ || ( this->enableEmbedding() && signalTree_ != 0 && signalTree_->haveBranch("mcMatch") ) ) { this->addGenNtupleIntegerBranch("genTMSig"); this->addGenNtupleIntegerBranch("genSCFSig"); } const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { TString name( this->bkgndClassName(iBkgnd) ); name.Prepend("gen"); this->addGenNtupleIntegerBranch(name); } if (this->useDP() == kTRUE) { this->addGenNtupleDoubleBranch("m12"); this->addGenNtupleDoubleBranch("m23"); this->addGenNtupleDoubleBranch("m13"); this->addGenNtupleDoubleBranch("m12Sq"); this->addGenNtupleDoubleBranch("m23Sq"); this->addGenNtupleDoubleBranch("m13Sq"); this->addGenNtupleDoubleBranch("cosHel12"); this->addGenNtupleDoubleBranch("cosHel23"); this->addGenNtupleDoubleBranch("cosHel13"); if (kinematics_->squareDP()) { this->addGenNtupleDoubleBranch("mPrime"); this->addGenNtupleDoubleBranch("thPrime"); } } for (LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter) { std::vector varNames = (*pdf_iter)->varNames(); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { this->addGenNtupleDoubleBranch( (*var_iter) ); } } } } void LauSimpleFitModel::setDPBranchValues() { // Store all the DP information this->setGenNtupleDoubleBranchValue("m12", kinematics_->getm12()); this->setGenNtupleDoubleBranchValue("m23", kinematics_->getm23()); this->setGenNtupleDoubleBranchValue("m13", kinematics_->getm13()); this->setGenNtupleDoubleBranchValue("m12Sq", kinematics_->getm12Sq()); this->setGenNtupleDoubleBranchValue("m23Sq", kinematics_->getm23Sq()); this->setGenNtupleDoubleBranchValue("m13Sq", kinematics_->getm13Sq()); this->setGenNtupleDoubleBranchValue("cosHel12", kinematics_->getc12()); this->setGenNtupleDoubleBranchValue("cosHel23", kinematics_->getc23()); this->setGenNtupleDoubleBranchValue("cosHel13", kinematics_->getc13()); if (kinematics_->squareDP()) { this->setGenNtupleDoubleBranchValue("mPrime", kinematics_->getmPrime()); this->setGenNtupleDoubleBranchValue("thPrime", kinematics_->getThetaPrime()); } } void LauSimpleFitModel::generateExtraPdfValues(LauPdfList* extraPdfs, LauEmbeddedData* embeddedData) { // Generate from the extra PDFs if (extraPdfs) { for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) { LauFitData genValues; if (embeddedData) { genValues = embeddedData->getValues( (*pdf_iter)->varNames() ); } else { genValues = (*pdf_iter)->generate(kinematics_); } for ( LauFitData::const_iterator var_iter = genValues.begin(); var_iter != genValues.end(); ++var_iter ) { TString varName = var_iter->first; if ( varName != "m13Sq" && varName != "m23Sq" ) { Double_t value = var_iter->second; this->setGenNtupleDoubleBranchValue(varName,value); } } } } } void LauSimpleFitModel::propagateParUpdates() { // Update the total normalisation for the signal likelihood if (this->useDP() == kTRUE) { this->updateCoeffs(); sigDPModel_->updateCoeffs(coeffs_); } // Update the signal events from the background events if not doing an extended fit if ( !this->doEMLFit() && !signalEvents_->fixed() ) { this->updateSigEvents(); } } void LauSimpleFitModel::updateSigEvents() { // The background parameters will have been set from Minuit. // We need to update the signal events using these. Double_t nTotEvts = this->eventsPerExpt(); signalEvents_->range(-2.0*nTotEvts,2.0*nTotEvts); for (LauBkgndYieldList::iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) { LauAbsRValue* nBkgndEvents = (*iter); if ( nBkgndEvents->isLValue() ) { LauParameter* yield = dynamic_cast( nBkgndEvents ); yield->range(-2.0*nTotEvts,2.0*nTotEvts); } } if ( signalEvents_->fixed() ) { return; } // Subtract background events (if any) from signal. Double_t signalEvents = nTotEvts; if ( usingBkgnd_ == kTRUE ) { for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) { signalEvents -= (*iter)->value(); } } signalEvents_->value(signalEvents); } void LauSimpleFitModel::cacheInputFitVars() { // Fill the internal data trees of the signal and backgrond models. LauFitDataTree* inputFitData = this->fitData(); // First the Dalitz plot variables (m_ij^2) if (this->useDP() == kTRUE) { // need to append SCF smearing bins before caching DP amplitudes if ( scfMap_ != 0 ) { this->appendBinCentres( inputFitData ); } sigDPModel_->fillDataTree(*inputFitData); if (usingBkgnd_ == kTRUE) { for (LauBkgndDPModelList::iterator iter = bkgndDPModels_.begin(); iter != bkgndDPModels_.end(); ++iter) { (*iter)->fillDataTree(*inputFitData); } } } // ...and then the extra PDFs this->cacheInfo(signalPdfs_, *inputFitData); this->cacheInfo(scfPdfs_, *inputFitData); for (LauBkgndPdfsList::iterator iter = bkgndPdfs_.begin(); iter != bkgndPdfs_.end(); ++iter) { this->cacheInfo((*iter), *inputFitData); } // and finally the SCF fractions and jacobians if ( useSCF_ && useSCFHist_ ) { if ( !inputFitData->haveBranch( "m13Sq" ) || !inputFitData->haveBranch( "m23Sq" ) ) { std::cerr << "ERROR in LauSimpleFitModel::cacheInputFitVars : Input data does not contain DP branches and so can't cache the SCF fraction." << std::endl; gSystem->Exit(EXIT_FAILURE); } UInt_t nEvents = inputFitData->nEvents(); recoSCFFracs_.clear(); recoSCFFracs_.reserve( nEvents ); if ( kinematics_->squareDP() ) { recoJacobians_.clear(); recoJacobians_.reserve( nEvents ); } for (UInt_t iEvt = 0; iEvt < nEvents; iEvt++) { const LauFitData& dataValues = inputFitData->getData(iEvt); LauFitData::const_iterator m13_iter = dataValues.find("m13Sq"); LauFitData::const_iterator m23_iter = dataValues.find("m23Sq"); kinematics_->updateKinematics( m13_iter->second, m23_iter->second ); Double_t scfFrac = scfFracHist_->calcEfficiency( kinematics_ ); recoSCFFracs_.push_back( scfFrac ); if ( kinematics_->squareDP() ) { recoJacobians_.push_back( kinematics_->calcSqDPJacobian() ); } } } } void LauSimpleFitModel::appendBinCentres( LauFitDataTree* inputData ) { // We'll be caching the DP amplitudes and efficiencies of the centres of the true bins. // To do so, we attach some fake points at the end of inputData, the number of the entry // minus the total number of events corresponding to the number of the histogram for that // given true bin in the LauScfMap object. (What this means is that when Laura is provided with // the LauScfMap object by the user, it's the latter who has to make sure that it contains the // right number of histograms and in exactly the right order!) // Get the x and y co-ordinates of the bin centres std::vector binCentresXCoords; std::vector binCentresYCoords; scfMap_->listBinCentres(binCentresXCoords, binCentresYCoords); // The SCF histograms could be in square Dalitz plot histograms. // The dynamics takes normal Dalitz plot coords, so we might have to convert them back. Bool_t sqDP = kinematics_->squareDP(); UInt_t nBins = binCentresXCoords.size(); fakeSCFFracs_.clear(); fakeSCFFracs_.reserve( nBins ); if ( sqDP ) { fakeJacobians_.clear(); fakeJacobians_.reserve( nBins ); } for (UInt_t iBin = 0; iBin < nBins; ++iBin) { if ( sqDP ) { kinematics_->updateSqDPKinematics(binCentresXCoords[iBin],binCentresYCoords[iBin]); binCentresXCoords[iBin] = kinematics_->getm13Sq(); binCentresYCoords[iBin] = kinematics_->getm23Sq(); fakeJacobians_.push_back( kinematics_->calcSqDPJacobian() ); } else { kinematics_->updateKinematics(binCentresXCoords[iBin],binCentresYCoords[iBin]); } fakeSCFFracs_.push_back( scfFracHist_->calcEfficiency( kinematics_ ) ); } // Set up inputFitVars_ object to hold the fake events inputData->appendFakePoints(binCentresXCoords,binCentresYCoords); } Double_t LauSimpleFitModel::getTotEvtLikelihood(UInt_t iEvt) { // Get the DP likelihood for signal and backgrounds this->getEvtDPLikelihood(iEvt); // Get the combined extra PDFs likelihood for signal and backgrounds this->getEvtExtraLikelihoods(iEvt); // If appropriate, combine the TM and SCF likelihoods Double_t sigLike = sigDPLike_ * sigExtraLike_; if ( useSCF_ ) { Double_t scfFrac(0.0); if (useSCFHist_) { scfFrac = recoSCFFracs_[iEvt]; } else { scfFrac = scfFrac_.unblindValue(); } sigLike *= (1.0 - scfFrac); if ( (scfMap_ != 0) && (this->useDP() == kTRUE) ) { // if we're smearing the SCF DP PDF then the SCF frac // is already included in the SCF DP likelihood sigLike += (scfDPLike_ * scfExtraLike_); } else { sigLike += (scfFrac * scfDPLike_ * scfExtraLike_); } } // Construct the total event likelihood Double_t likelihood = signalEvents_->unblindValue() * sigLike; const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { likelihood += (bkgndEvents_[bkgndID]->unblindValue() * bkgndDPLike_[bkgndID] * bkgndExtraLike_[bkgndID]); } return likelihood; } Double_t LauSimpleFitModel::getEventSum() const { Double_t eventSum(0.0); eventSum += signalEvents_->unblindValue(); for (LauBkgndYieldList::const_iterator iter = bkgndEvents_.begin(); iter != bkgndEvents_.end(); ++iter) { eventSum += (*iter)->unblindValue(); } return eventSum; } void LauSimpleFitModel::getEvtDPLikelihood(UInt_t iEvt) { // Function to return the signal and background likelihoods for the // Dalitz plot for the given event evtNo. if (this->useDP() == kTRUE) { sigDPModel_->calcLikelihoodInfo(iEvt); sigDPLike_ = sigDPModel_->getEvtLikelihood(); if ( useSCF_ == kTRUE ) { if ( scfMap_ == 0 ) { // we're not smearing the SCF DP position // so the likelihood is the same as the TM scfDPLike_ = sigDPLike_; } else { // calculate the smeared SCF DP likelihood scfDPLike_ = this->getEvtSCFDPLikelihood(iEvt); } } const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { if (usingBkgnd_ == kTRUE) { bkgndDPLike_[bkgndID] = bkgndDPModels_[bkgndID]->getLikelihood(iEvt); } else { bkgndDPLike_[bkgndID] = 0.0; } } } else { // There's always going to be a term in the likelihood for the // signal, so we'd better not zero it. sigDPLike_ = 1.0; scfDPLike_ = 1.0; const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { if (usingBkgnd_ == kTRUE) { bkgndDPLike_[bkgndID] = 1.0; } else { bkgndDPLike_[bkgndID] = 0.0; } } } } Double_t LauSimpleFitModel::getEvtSCFDPLikelihood(UInt_t iEvt) { Double_t scfDPLike(0.0); Double_t recoJacobian(1.0); Double_t xCoord(0.0); Double_t yCoord(0.0); Bool_t squareDP = kinematics_->squareDP(); if ( squareDP ) { xCoord = sigDPModel_->getEvtmPrime(); yCoord = sigDPModel_->getEvtthPrime(); recoJacobian = recoJacobians_[iEvt]; } else { xCoord = sigDPModel_->getEvtm13Sq(); yCoord = sigDPModel_->getEvtm23Sq(); } // Find the bin that our reco event falls in Int_t recoBin = scfMap_->binNumber( xCoord, yCoord ); // Find out which true Bins contribute to the given reco bin const std::vector* trueBins = scfMap_->trueBins(recoBin); Int_t nDataEvents = this->eventsPerExpt(); // Loop over the true bins for (std::vector::const_iterator iter = trueBins->begin(); iter != trueBins->end(); ++iter) { Int_t trueBin = (*iter); // prob of a true event in the given true bin migrating to the reco bin Double_t pRecoGivenTrue = scfMap_->prob( recoBin, trueBin ); // We've cached the DP amplitudes and the efficiency for the // true bin centres, just after the data points sigDPModel_->calcLikelihoodInfo( nDataEvents + trueBin ); Double_t pTrue = sigDPModel_->getEvtLikelihood(); // Get the cached SCF fraction (and jacobian if we're using the square DP) Double_t scfFraction = fakeSCFFracs_[ trueBin ]; Double_t jacobian(1.0); if ( squareDP ) { jacobian = fakeJacobians_[ trueBin ]; } scfDPLike += pTrue * jacobian * scfFraction * pRecoGivenTrue; } // Divide by the reco jacobian scfDPLike /= recoJacobian; return scfDPLike; } void LauSimpleFitModel::getEvtExtraLikelihoods(UInt_t iEvt) { // Function to return the signal and background likelihoods for the // extra variables for the given event evtNo. sigExtraLike_ = 1.0; for (LauPdfList::iterator iter = signalPdfs_.begin(); iter != signalPdfs_.end(); ++iter) { (*iter)->calcLikelihoodInfo(iEvt); sigExtraLike_ *= (*iter)->getLikelihood(); } if (useSCF_) { scfExtraLike_ = 1.0; for (LauPdfList::iterator iter = scfPdfs_.begin(); iter != scfPdfs_.end(); ++iter) { (*iter)->calcLikelihoodInfo(iEvt); scfExtraLike_ *= (*iter)->getLikelihood(); } } const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t bkgndID(0); bkgndID < nBkgnds; ++bkgndID ) { if (usingBkgnd_) { bkgndExtraLike_[bkgndID] = 1.0; LauPdfList& pdfList = bkgndPdfs_[bkgndID]; for (LauPdfList::iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) { (*pdf_iter)->calcLikelihoodInfo(iEvt); bkgndExtraLike_[bkgndID] *= (*pdf_iter)->getLikelihood(); } } else { bkgndExtraLike_[bkgndID] = 0.0; } } } void LauSimpleFitModel::updateCoeffs() { coeffs_.clear(); coeffs_.reserve(nSigComp_); for (UInt_t i = 0; i < nSigComp_; i++) { coeffs_.push_back(coeffPars_[i]->particleCoeff()); } } void LauSimpleFitModel::setupSPlotNtupleBranches() { // add branches for storing the experiment number and the number of // the event within the current experiment this->addSPlotNtupleIntegerBranch("iExpt"); this->addSPlotNtupleIntegerBranch("iEvtWithinExpt"); // Store the efficiency of the event (for inclusive BF calculations). if (this->storeDPEff()) { this->addSPlotNtupleDoubleBranch("efficiency"); if ( sigDPModel_->usingScfModel() ) { this->addSPlotNtupleDoubleBranch("scffraction"); } } // Store the total event likelihood for each species. if (useSCF_) { this->addSPlotNtupleDoubleBranch("sigTMTotalLike"); this->addSPlotNtupleDoubleBranch("sigSCFTotalLike"); this->addSPlotNtupleDoubleBranch("sigSCFFrac"); } else { this->addSPlotNtupleDoubleBranch("sigTotalLike"); } if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { TString name( this->bkgndClassName(iBkgnd) ); name += "TotalLike"; this->addSPlotNtupleDoubleBranch(name); } } // Store the DP likelihoods if (this->useDP()) { if (useSCF_) { this->addSPlotNtupleDoubleBranch("sigTMDPLike"); this->addSPlotNtupleDoubleBranch("sigSCFDPLike"); } else { this->addSPlotNtupleDoubleBranch("sigDPLike"); } if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { TString name( this->bkgndClassName(iBkgnd) ); name += "DPLike"; this->addSPlotNtupleDoubleBranch(name); } } } // Store the likelihoods for each extra PDF if (useSCF_) { this->addSPlotNtupleBranches(&signalPdfs_, "sigTM"); this->addSPlotNtupleBranches(&scfPdfs_, "sigSCF"); } else { this->addSPlotNtupleBranches(&signalPdfs_, "sig"); } if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { const TString& bkgndClass = this->bkgndClassName(iBkgnd); const LauPdfList* pdfList = &(bkgndPdfs_[iBkgnd]); this->addSPlotNtupleBranches(pdfList, bkgndClass); } } } void LauSimpleFitModel::addSPlotNtupleBranches(const LauPdfList* extraPdfs, const TString& prefix) { if (extraPdfs) { // Loop through each of the PDFs for (LauPdfList::const_iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) { // Count the number of input variables that are not // DP variables (used in the case where there is DP // dependence for e.g. MVA) UInt_t nVars(0); std::vector varNames = (*pdf_iter)->varNames(); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { ++nVars; } } if ( nVars == 1 ) { // If the PDF only has one variable then // simply add one branch for that variable TString varName = (*pdf_iter)->varName(); TString name(prefix); name += varName; name += "Like"; this->addSPlotNtupleDoubleBranch(name); } else if ( nVars == 2 ) { // If the PDF has two variables then we // need a branch for them both together and // branches for each TString allVars(""); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { allVars += (*var_iter); TString name(prefix); name += (*var_iter); name += "Like"; this->addSPlotNtupleDoubleBranch(name); } TString name(prefix); name += allVars; name += "Like"; this->addSPlotNtupleDoubleBranch(name); } else { std::cerr << "WARNING in LauSimpleFitModel::addSPlotNtupleBranches : Can't yet deal with 3D PDFs." << std::endl; } } } } Double_t LauSimpleFitModel::setSPlotNtupleBranchValues(LauPdfList* extraPdfs, const TString& prefix, UInt_t iEvt) { // Store the PDF value for each variable in the list Double_t totalLike(1.0); Double_t extraLike(0.0); if (extraPdfs) { for (LauPdfList::iterator pdf_iter = extraPdfs->begin(); pdf_iter != extraPdfs->end(); ++pdf_iter) { // calculate the likelihood for this event (*pdf_iter)->calcLikelihoodInfo(iEvt); extraLike = (*pdf_iter)->getLikelihood(); totalLike *= extraLike; // Count the number of input variables that are not // DP variables (used in the case where there is DP // dependence for e.g. MVA) UInt_t nVars(0); std::vector varNames = (*pdf_iter)->varNames(); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { ++nVars; } } if ( nVars == 1 ) { // If the PDF only has one variable then // simply store the value for that variable TString varName = (*pdf_iter)->varName(); TString name(prefix); name += varName; name += "Like"; this->setSPlotNtupleDoubleBranchValue(name, extraLike); } else if ( nVars == 2 ) { // If the PDF has two variables then we // store the value for them both together // and for each on their own TString allVars(""); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { allVars += (*var_iter); TString name(prefix); name += (*var_iter); name += "Like"; Double_t indivLike = (*pdf_iter)->getLikelihood( (*var_iter) ); this->setSPlotNtupleDoubleBranchValue(name, indivLike); } TString name(prefix); name += allVars; name += "Like"; this->setSPlotNtupleDoubleBranchValue(name, extraLike); } else { std::cerr << "WARNING in LauSimpleFitModel::setSPlotNtupleBranchValues : Can't yet deal with 3D PDFs." << std::endl; } } } return totalLike; } LauSPlot::NameSet LauSimpleFitModel::variableNames() const { LauSPlot::NameSet nameSet; if (this->useDP()) { nameSet.insert("DP"); } // Loop through all the signal PDFs for (LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter) { // Loop over the variables involved in each PDF std::vector varNames = (*pdf_iter)->varNames(); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { // If they are not DP coordinates then add them if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { nameSet.insert( (*var_iter) ); } } } return nameSet; } LauSPlot::NumbMap LauSimpleFitModel::freeSpeciesNames() const { LauSPlot::NumbMap numbMap; if (!signalEvents_->fixed() && this->doEMLFit()) { numbMap["sig"] = signalEvents_->genValue(); } if ( usingBkgnd_ ) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { const TString& bkgndClass = this->bkgndClassName(iBkgnd); const LauAbsRValue* par = bkgndEvents_[iBkgnd]; if (!par->fixed()) { numbMap[bkgndClass] = par->genValue(); if ( ! par->isLValue() ) { std::cerr << "WARNING in LauSimpleFitModel::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; } } } } return numbMap; } LauSPlot::NumbMap LauSimpleFitModel::fixdSpeciesNames() const { LauSPlot::NumbMap numbMap; if (signalEvents_->fixed() && this->doEMLFit()) { numbMap["sig"] = signalEvents_->genValue(); } if ( usingBkgnd_ ) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { const TString& bkgndClass = this->bkgndClassName(iBkgnd); const LauAbsRValue* par = bkgndEvents_[iBkgnd]; if (par->fixed()) { numbMap[bkgndClass] = par->genValue(); } } } return numbMap; } LauSPlot::TwoDMap LauSimpleFitModel::twodimPDFs() const { LauSPlot::TwoDMap twodimMap; for (LauPdfList::const_iterator pdf_iter = signalPdfs_.begin(); pdf_iter != signalPdfs_.end(); ++pdf_iter) { // Count the number of input variables that are not DP variables UInt_t nVars(0); std::vector varNames = (*pdf_iter)->varNames(); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { ++nVars; } } if ( nVars == 2 ) { if (useSCF_) { twodimMap.insert( std::make_pair( "sigTM", std::make_pair( varNames[0], varNames[1] ) ) ); } else { twodimMap.insert( std::make_pair( "sig", std::make_pair( varNames[0], varNames[1] ) ) ); } } } if ( useSCF_ ) { for (LauPdfList::const_iterator pdf_iter = scfPdfs_.begin(); pdf_iter != scfPdfs_.end(); ++pdf_iter) { // Count the number of input variables that are not DP variables UInt_t nVars(0); std::vector varNames = (*pdf_iter)->varNames(); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { ++nVars; } } if ( nVars == 2 ) { twodimMap.insert( std::make_pair( "sigSCF", std::make_pair( varNames[0], varNames[1] ) ) ); } } } if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { const TString& bkgndClass = this->bkgndClassName(iBkgnd); const LauPdfList& pdfList = bkgndPdfs_[iBkgnd]; for (LauPdfList::const_iterator pdf_iter = pdfList.begin(); pdf_iter != pdfList.end(); ++pdf_iter) { // Count the number of input variables that are not DP variables UInt_t nVars(0); std::vector varNames = (*pdf_iter)->varNames(); for ( std::vector::const_iterator var_iter = varNames.begin(); var_iter != varNames.end(); ++var_iter ) { if ( (*var_iter) != "m13Sq" && (*var_iter) != "m23Sq" ) { ++nVars; } } if ( nVars == 2 ) { twodimMap.insert( std::make_pair( bkgndClass, std::make_pair( varNames[0], varNames[1] ) ) ); } } } } return twodimMap; } void LauSimpleFitModel::storePerEvtLlhds() { std::cout << "INFO in LauSimpleFitModel::storePerEvtLlhds : Storing per-event likelihood values..." << std::endl; // if we've not been using the DP model then we need to cache all // the info here so that we can get the efficiency from it LauFitDataTree* inputFitData = this->fitData(); if (!this->useDP() && this->storeDPEff()) { sigDPModel_->initialise(coeffs_); sigDPModel_->fillDataTree(*inputFitData); } UInt_t evtsPerExpt(this->eventsPerExpt()); for (UInt_t iEvt = 0; iEvt < evtsPerExpt; ++iEvt) { this->setSPlotNtupleIntegerBranchValue("iExpt",this->iExpt()); this->setSPlotNtupleIntegerBranchValue("iEvtWithinExpt",iEvt); // the DP information this->getEvtDPLikelihood(iEvt); if (this->storeDPEff()) { if (!this->useDP()) { sigDPModel_->calcLikelihoodInfo(iEvt); } this->setSPlotNtupleDoubleBranchValue("efficiency",sigDPModel_->getEvtEff()); if ( sigDPModel_->usingScfModel() ) { this->setSPlotNtupleDoubleBranchValue("scffraction",sigDPModel_->getEvtScfFraction()); } } if (this->useDP()) { sigTotalLike_ = sigDPLike_; if (useSCF_) { scfTotalLike_ = scfDPLike_; this->setSPlotNtupleDoubleBranchValue("sigTMDPLike",sigDPLike_); this->setSPlotNtupleDoubleBranchValue("sigSCFDPLike",scfDPLike_); } else { this->setSPlotNtupleDoubleBranchValue("sigDPLike",sigDPLike_); } if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { TString name = this->bkgndClassName(iBkgnd); name += "DPLike"; this->setSPlotNtupleDoubleBranchValue(name,bkgndDPLike_[iBkgnd]); } } } else { sigTotalLike_ = 1.0; if (useSCF_) { scfTotalLike_ = 1.0; } if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { bkgndTotalLike_[iBkgnd] = 1.0; } } } // the signal PDF values if ( useSCF_ ) { sigTotalLike_ *= this->setSPlotNtupleBranchValues(&signalPdfs_, "sigTM", iEvt); scfTotalLike_ *= this->setSPlotNtupleBranchValues(&scfPdfs_, "sigSCF", iEvt); } else { sigTotalLike_ *= this->setSPlotNtupleBranchValues(&signalPdfs_, "sig", iEvt); } // the background PDF values if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { const TString& bkgndClass = this->bkgndClassName(iBkgnd); LauPdfList& pdfs = bkgndPdfs_[iBkgnd]; bkgndTotalLike_[iBkgnd] *= this->setSPlotNtupleBranchValues(&(pdfs), bkgndClass, iEvt); } } // the total likelihoods if (useSCF_) { Double_t scfFrac(0.0); if ( useSCFHist_ ) { scfFrac = recoSCFFracs_[iEvt]; } else { scfFrac = scfFrac_.unblindValue(); } this->setSPlotNtupleDoubleBranchValue("sigSCFFrac",scfFrac); sigTotalLike_ *= ( 1.0 - scfFrac ); if ( scfMap_ == 0 ) { scfTotalLike_ *= scfFrac; } this->setSPlotNtupleDoubleBranchValue("sigTMTotalLike",sigTotalLike_); this->setSPlotNtupleDoubleBranchValue("sigSCFTotalLike",scfTotalLike_); } else { this->setSPlotNtupleDoubleBranchValue("sigTotalLike",sigTotalLike_); } if (usingBkgnd_) { const UInt_t nBkgnds = this->nBkgndClasses(); for ( UInt_t iBkgnd(0); iBkgnd < nBkgnds; ++iBkgnd ) { TString name = this->bkgndClassName(iBkgnd); name += "TotalLike"; this->setSPlotNtupleDoubleBranchValue(name,bkgndTotalLike_[iBkgnd]); } } // fill the tree this->fillSPlotNtupleBranches(); } std::cout << "INFO in LauSimpleFitModel::storePerEvtLlhds : Finished storing per-event likelihood values." << std::endl; } void LauSimpleFitModel::embedSignal(const TString& fileName, const TString& treeName, Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment, Bool_t useReweighting) { if (signalTree_) { std::cerr << "ERROR in LauSimpleFitModel::embedSignal : Already embedding signal from a file." << std::endl; return; } signalTree_ = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment); Bool_t dataOK = signalTree_->findBranches(); if (!dataOK) { delete signalTree_; signalTree_ = 0; std::cerr << "ERROR in LauSimpleFitModel::embedSignal : Problem creating data tree for embedding." << std::endl; return; } reuseSignal_ = reuseEventsWithinEnsemble; useReweighting_ = useReweighting; if (this->enableEmbedding() == kFALSE) { this->enableEmbedding(kTRUE); } } void LauSimpleFitModel::embedBkgnd(const TString& bkgndClass, const TString& fileName, const TString& treeName, Bool_t reuseEventsWithinEnsemble, Bool_t reuseEventsWithinExperiment) { if ( ! this->validBkgndClass( bkgndClass ) ) { 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; return; } UInt_t bkgndID = this->bkgndClassID( bkgndClass ); if (bkgndTree_[bkgndID]) { std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Already embedding background from a file." << std::endl; return; } bkgndTree_[bkgndID] = new LauEmbeddedData(fileName,treeName,reuseEventsWithinExperiment); Bool_t dataOK = bkgndTree_[bkgndID]->findBranches(); if (!dataOK) { delete bkgndTree_[bkgndID]; bkgndTree_[bkgndID] = 0; std::cerr << "ERROR in LauSimpleFitModel::embedBkgnd : Problem creating data tree for embedding." << std::endl; return; } reuseBkgnd_[bkgndID] = reuseEventsWithinEnsemble; if (this->enableEmbedding() == kFALSE) { this->enableEmbedding(kTRUE); } } void LauSimpleFitModel::weightEvents( const TString& dataFileName, const TString& dataTreeName ) { // Routine to provide weights for events that are uniformly distributed // in the DP (or square DP) so as to reproduce the given DP model if ( kinematics_->squareDP() ) { std::cout << "INFO in LauSimpleFitModel::weightEvents : will create weights assuming events were generated flat in the square DP" << std::endl; } else { std::cout << "INFO in LauSimpleFitModel::weightEvents : will create weights assuming events were generated flat in phase space" << std::endl; } // This reads in the given dataFile and creates an input // fit data tree that stores them for all events and experiments. Bool_t dataOK = this->verifyFitData(dataFileName,dataTreeName); if (!dataOK) { std::cerr << "ERROR in LauSimpleFitModel::weightEvents : Problem caching the data." << std::endl; return; } LauFitDataTree* inputFitData = this->fitData(); if ( ! inputFitData->haveBranch( "m13Sq_MC" ) || ! inputFitData->haveBranch( "m23Sq_MC" ) ) { std::cerr << "WARNING in LauSimpleFitModel::weightEvents : Cannot find MC truth DP coordinate branches in supplied data, aborting." << std::endl; return; } // Create the ntuple to hold the DP weights TString weightsFileName( dataFileName ); Ssiz_t index = weightsFileName.Last('.'); weightsFileName.Insert( index, "_DPweights" ); LauGenNtuple * weightsTuple = new LauGenNtuple( weightsFileName, dataTreeName ); weightsTuple->addIntegerBranch("iExpt"); weightsTuple->addIntegerBranch("iEvtWithinExpt"); weightsTuple->addDoubleBranch("dpModelWeight"); UInt_t nExpmt = this->nExpt(); UInt_t firstExpmt = this->firstExpt(); for (UInt_t iExpmt = firstExpmt; iExpmt < (firstExpmt+nExpmt); ++iExpmt) { inputFitData->readExperimentData(iExpmt); UInt_t nEvents = inputFitData->nEvents(); if (nEvents < 1) { std::cerr << "WARNING in LauSimpleFitModel::weightEvents : Zero events in experiment " << iExpmt << ", skipping..." << std::endl; continue; } weightsTuple->setIntegerBranchValue( "iExpt", iExpmt ); // Calculate and store the weights for the events in this experiment for ( UInt_t iEvent(0); iEvent < nEvents; ++iEvent ) { weightsTuple->setIntegerBranchValue( "iEvtWithinExpt", iEvent ); const LauFitData& evtData = inputFitData->getData( iEvent ); Double_t m13Sq_MC = evtData.find("m13Sq_MC")->second; Double_t m23Sq_MC = evtData.find("m23Sq_MC")->second; Double_t dpModelWeight(0.0); if ( kinematics_->withinDPLimits( m13Sq_MC, m23Sq_MC ) ) { kinematics_->updateKinematics( m13Sq_MC, m23Sq_MC ); dpModelWeight = sigDPModel_->getEventWeight(); if ( kinematics_->squareDP() ) { dpModelWeight *= kinematics_->calcSqDPJacobian(); } dpModelWeight /= sigDPModel_->getDPNorm(); } weightsTuple->setDoubleBranchValue( "dpModelWeight", dpModelWeight ); weightsTuple->fillBranches(); } } weightsTuple->buildIndex( "iExpt", "iEvtWithinExpt" ); weightsTuple->addFriendTree(dataFileName, dataTreeName); weightsTuple->writeOutGenResults(); delete weightsTuple; } void LauSimpleFitModel::savePDFPlots(const TString& label) { savePDFPlotsWave(label, 0); savePDFPlotsWave(label, 1); savePDFPlotsWave(label, 2); std::cout << "LauCPFitModel::plot" << std::endl; // ((LauIsobarDynamics*)sigDPModel_)->plot(); //Double_t minm13 = sigDPModel_->getKinematics()->getm13Min(); Double_t minm13 = 0.0; Double_t maxm13 = sigDPModel_->getKinematics()->getm13Max(); //Double_t minm23 = sigDPModel_->getKinematics()->getm23Min(); Double_t minm23 = 0.0; Double_t maxm23 = sigDPModel_->getKinematics()->getm23Max(); Double_t mins13 = minm13*minm13; Double_t maxs13 = maxm13*maxm13; Double_t mins23 = minm23*minm23; Double_t maxs23 = maxm23*maxm23; Double_t s13, s23, chPdf; Int_t n13=200.00, n23=200.00; Double_t delta13, delta23; delta13 = (maxs13 - mins13)/n13; delta23 = (maxs23 - mins23)/n23; UInt_t nAmp = sigDPModel_->getnCohAmp(); for (UInt_t resID = 0; resID <= nAmp; ++resID) { TGraph2D *dt = new TGraph2D(); TString resName = "TotalAmp"; if (resID != nAmp){ TString tStrResID = Form("%d", resID); const LauIsobarDynamics* model = sigDPModel_; const LauAbsResonance* resonance = model->getResonance(resID); resName = resonance->getResonanceName(); std::cout << "resName = " << resName << std::endl; } resName.ReplaceAll("(", ""); resName.ReplaceAll(")", ""); resName.ReplaceAll("*", "Star"); TCanvas *c = new TCanvas("c"+resName+label,resName+" ("+label+")",0,0,600,400); dt->SetName(resName+label); dt->SetTitle(resName+" ("+label+")"); Int_t count=0; for (Int_t i=0; igetKinematics()->withinDPLimits2(s23, s13)) { //if (s13 > 4) continue; sigDPModel_->calcLikelihoodInfo(s13, s23); LauComplex chAmp = sigDPModel_->getEvtDPAmp(); if (resID != nAmp){ chAmp = sigDPModel_->getFullAmplitude(resID); } chPdf = chAmp.abs2(); //if ((z > 0.04)||(z < -0.03)) continue; //z = TMath::Log(z); if (sigDPModel_->getDaughters()->gotSymmetricalDP()){ Double_t sLow = s13; Double_t sHigh = s23; if (sLow>sHigh) { continue; //sLow = s23; //sHigh = s13; } //if (sLow > 3.5) continue; //if (i<10) std::cout << "SymmetricalDP" << std::endl; //if (z>0.02) z = 0.02; //if (z<-0.02) z = -0.02; dt->SetPoint(count,sHigh,sLow,chPdf); count++; } else { dt->SetPoint(count,s13,s23,chPdf); count++; } } } } gStyle->SetPalette(1); dt->GetXaxis()->SetTitle("m_{K#pi}(low)"); dt->GetYaxis()->SetTitle("m_{K#pi}(high)"); dt->Draw("SURF1"); dt->GetXaxis()->SetTitle("m_{K#pi}(low)"); dt->GetYaxis()->SetTitle("m_{K#pi}(high)"); c->SaveAs("plot_2D_"+resName + "_"+label+".C"); } } void LauSimpleFitModel::savePDFPlotsWave(const TString& label, const Int_t& spin) { std::cout << "label = "<< label << ", spin = "<< spin << std::endl; TString tStrResID = "S_Wave"; if (spin == 1) tStrResID = "P_Wave"; if (spin == 2) tStrResID = "D_Wave"; std::cout << "LauSimpleFitModel::savePDFPlotsWave: "<< tStrResID << std::endl; TCanvas *c = new TCanvas("c"+tStrResID+label,tStrResID+" ("+label+")",0,0,600,400); Double_t minm13 = 0.0; Double_t maxm13 = sigDPModel_->getKinematics()->getm13Max(); Double_t minm23 = 0.0; Double_t maxm23 = sigDPModel_->getKinematics()->getm23Max(); Double_t mins13 = minm13*minm13; Double_t maxs13 = maxm13*maxm13; Double_t mins23 = minm23*minm23; Double_t maxs23 = maxm23*maxm23; Double_t s13, s23, chPdf; TGraph2D *dt = new TGraph2D(); dt->SetName(tStrResID+label); dt->SetTitle(tStrResID+" ("+label+")"); Int_t n13=200.00, n23=200.00; Double_t delta13, delta23; delta13 = (maxs13 - mins13)/n13; delta23 = (maxs23 - mins23)/n23; UInt_t nAmp = sigDPModel_->getnCohAmp(); Int_t count=0; for (Int_t i=0; igetKinematics()->withinDPLimits2(s23, s13)) { //if (s13 > 4) continue; LauComplex chAmp(0,0); Bool_t noWaveRes = kTRUE; sigDPModel_->calcLikelihoodInfo(s13, s23); for (UInt_t resID = 0; resID < nAmp; ++resID) { const LauIsobarDynamics* model = sigDPModel_; const LauAbsResonance* resonance = model->getResonance(resID); Int_t spin_res = resonance->getSpin(); if (spin != spin_res) continue; noWaveRes = kFALSE; chAmp += sigDPModel_->getFullAmplitude(resID); } if (noWaveRes) return; chPdf = chAmp.abs2(); if (sigDPModel_->getDaughters()->gotSymmetricalDP()){ Double_t sLow = s13; Double_t sHigh = s23; if (sLow>sHigh) { continue; //sLow = s23; //sHigh = s13; } dt->SetPoint(count,sHigh,sLow,chPdf); count++; } else { dt->SetPoint(count,s13,s23,chPdf); count++; } } } } gStyle->SetPalette(1); dt->GetXaxis()->SetTitle("pipi"); dt->GetYaxis()->SetTitle("pipi"); dt->Draw("SURF1"); dt->GetXaxis()->SetTitle("pipi"); dt->GetYaxis()->SetTitle("pipi"); c->SaveAs("plot_2D_"+tStrResID+"_"+label+".C"); }