diff --git a/examples/GenFitTimeDep_Bd2Dpipi.cc b/examples/GenFitTimeDep_Bd2Dpipi.cc index c1c0592..8ee7b34 100644 --- a/examples/GenFitTimeDep_Bd2Dpipi.cc +++ b/examples/GenFitTimeDep_Bd2Dpipi.cc @@ -1,1023 +1,978 @@ #include "TimeDep_Dpipi_ProgOpts.hh" #include "Lau1DHistPdf.hh" #include "Lau1DCubicSpline.hh" #include "LauBinnedDecayTimeEfficiency.hh" #include "LauBkgndDPModel.hh" #include "LauCrystalBallPdf.hh" #include "LauDaughters.hh" #include "LauDecayTime.hh" #include "LauDecayTimePdf.hh" #include "LauDecayTimePhysicsModel.hh" #include "LauDecayTimeResolution.hh" #include "LauEffModel.hh" #include "LauExponentialPdf.hh" #include "LauFlavTag.hh" #include "LauGounarisSakuraiRes.hh" #include "LauIsobarDynamics.hh" #include "LauLHCbNR.hh" #include "LauMagPhaseCoeffSet.hh" #include "LauRandom.hh" #include "LauRealImagCoeffSet.hh" #include "LauResonanceMaker.hh" #include "LauSplineDecayTimeEfficiency.hh" #include "LauSumPdf.hh" #include "LauTimeDepFitModel.hh" #include "LauVetoes.hh" #include "TFile.h" #include "TH2.h" #include "TRandom.h" #include "TString.h" #include "TSystem.h" #include "TF1.h" #include "TCanvas.h" #include "nlohmann/json.hpp" #include #include #include #include #include #include // Helper function for creating Double-Crystal-Ball PDFs LauAbsPdf* makeDCBPDF( const TString& varName, const Double_t varMin, const Double_t varMax, const TString& componentName, const std::map paramVals ); // Helper function for creating sets of DTA splines // TODO - not yet used using LauCubicSplineDecayTimeEfficiencyPtr = std::unique_ptr>; using LauSixthSplineDecayTimeEfficiencyPtr = std::unique_ptr>; -auto createSplineDTASet( const TString& jsonFileName, const TString& controlSplineName, const TString& signalSplineName, const std::vector& bkgndSplineNames ) +auto createSplineDTASet( const TString& jsonFileName, const TString& dModeStr, const TString& runStr, const TString& suffix, const TString& signalName, const std::vector& bkgndNames ) { + const TString controlSplineName { Form("%s_Kpi_%s%s",signalName.Data(),runStr.Data(),suffix.Data()) }; + const TString controlEffName { Form("dteff_QFS_%s",runStr.Data()) }; auto controlSpline = Lau1DCubicSpline::readFromJson( jsonFileName, controlSplineName ); - LauCubicSplineDecayTimeEfficiencyPtr controlEff { std::make_unique>( controlSplineName, std::move( controlSpline ) ) }; + LauCubicSplineDecayTimeEfficiencyPtr controlEff { std::make_unique>( controlEffName, std::move( controlSpline ) ) }; LauSixthSplineDecayTimeEfficiencyPtr signalEff; - if ( signalSplineName != "" ) { + if ( dModeStr != "Kpi" ) { + const TString signalSplineName { Form("%s_%s_%s%s",signalName.Data(),dModeStr.Data(),runStr.Data(),suffix.Data()) }; + const TString signalEffName { Form("dteff_CP%s_%s",dModeStr.Data(),runStr.Data()) }; auto signalSpline = Lau1DCubicSpline::readFromJson( jsonFileName, signalSplineName ); - signalEff = std::make_unique>( signalSplineName, *controlEff, std::move( signalSpline ) ); + signalEff = std::make_unique>( signalEffName, *controlEff, std::move( signalSpline ) ); } - std::vector bkgndEffs; - bkgndEffs.reserve( bkgndSplineNames.size() ); - for ( auto& bkgndSplineName : bkgndSplineNames ) { + std::map bkgndEffs; + for ( auto& name : bkgndNames ) { + const TString bkgndSplineName { Form("%s_%s_%s%s",name.Data(),dModeStr.Data(),runStr.Data(),suffix.Data()) }; + const TString bkgndEffName { Form("%s_DTA_ratio_%s_%s",name.Data(),runStr.Data(),dModeStr.Data()) }; auto bkgndSpline = Lau1DCubicSpline::readFromJson( jsonFileName, bkgndSplineName ); - bkgndEffs.push_back( std::make_unique>( bkgndSplineName, *controlEff, std::move( bkgndSpline ) ) ); + auto bkgndEff = std::make_unique>( bkgndEffName, *controlEff, std::move( bkgndSpline ) ); + bkgndEffs.emplace( std::make_pair( name, std::move(bkgndEff) ) ); } return std::make_tuple( std::move(controlEff), std::move(signalEff), std::move(bkgndEffs) ); } class BgTagDecayFlvAsym { private: const Double_t eff_decayB0_taggedB0 {0.}; const Double_t eff_decayB0_taggedB0bar {0.}; const Double_t eff_decayB0bar_taggedB0 {0.}; const Double_t eff_decayB0bar_taggedB0bar {0.}; const Double_t eff_decayB0_ave {0.}; const Double_t eff_decayB0bar_ave {0.}; const Double_t eff_decayB0_delta {0.}; const Double_t eff_decayB0bar_delta {0.}; public: BgTagDecayFlvAsym() = default; BgTagDecayFlvAsym( //This constructor assumes a tagging efficincy of 100% const Double_t n_decayB0_taggedB0, const Double_t n_decayB0_taggedB0bar, const Double_t n_decayB0bar_taggedB0, const Double_t n_decayB0bar_taggedB0bar ): eff_decayB0_taggedB0 {n_decayB0_taggedB0 /(n_decayB0_taggedB0+n_decayB0_taggedB0bar)}, eff_decayB0_taggedB0bar {n_decayB0_taggedB0bar/(n_decayB0_taggedB0+n_decayB0_taggedB0bar)}, eff_decayB0bar_taggedB0 {n_decayB0bar_taggedB0 /(n_decayB0bar_taggedB0+n_decayB0bar_taggedB0bar)}, eff_decayB0bar_taggedB0bar{n_decayB0bar_taggedB0bar/(n_decayB0bar_taggedB0+n_decayB0bar_taggedB0bar)}, eff_decayB0_ave {0.5*(eff_decayB0_taggedB0 + eff_decayB0_taggedB0bar)}, eff_decayB0bar_ave {0.5*(eff_decayB0bar_taggedB0 + eff_decayB0bar_taggedB0bar)}, eff_decayB0_delta {eff_decayB0_taggedB0 - eff_decayB0_taggedB0bar}, eff_decayB0bar_delta {eff_decayB0bar_taggedB0 - eff_decayB0bar_taggedB0bar} {} BgTagDecayFlvAsym( const nlohmann::json& j ) : BgTagDecayFlvAsym( j[0], j[1], j[2], j[3] ) {} std::pair decayB0() const {return std::make_pair(eff_decayB0_taggedB0 ,eff_decayB0_taggedB0bar );} std::pair decayB0bar() const {return std::make_pair(eff_decayB0bar_taggedB0,eff_decayB0bar_taggedB0bar);} std::pair aveDeltaB0() const {return std::make_pair(eff_decayB0_ave , eff_decayB0_delta) ;} std::pair aveDeltaB0bar() const {return std::make_pair(eff_decayB0bar_ave, eff_decayB0bar_delta);} //TODO add constructor for if tagging eff isn't 100%? }; int main(const int argc, const char ** argv) { const TestDpipi_ProgramSettings settings{argc,argv}; if ( ! settings.parsedOK ) { return EXIT_FAILURE; } if ( settings.helpRequested ) { return EXIT_SUCCESS; } if(settings.dataFit and settings.command == Command::Generate) { std::cerr << "I can't generate data! Wait for Run 3" << std::endl; return EXIT_FAILURE; } if(settings.dataFit and not settings.blindFit and not settings.fixPhiMix) { std::cerr << "We don't have permission to do unblinded fits yet!" << std::endl; return EXIT_FAILURE; } if(settings.run == 0 or settings.run > 2) { std::cerr << "The Run number must be either 1 or 2" << std::endl; return EXIT_FAILURE; } const UInt_t runNo = settings.run; const std::string runStr = Form("Run%u", runNo); + const TString dModeStr { TString{settings.directory}.ReplaceAll("B2Dhh_","") }; + LauRandom::setSeed(settings.RNGseed); LauDaughters* daughtersB0bar = new LauDaughters("B0_bar", "pi+", "pi-", "D0", kTRUE); LauDaughters* daughtersB0 = new LauDaughters("B0", "pi+", "pi-", "D0_bar",kTRUE); // efficiency auto vetoes { LauVetoes::readFromJson( "Bd2D0pipi_DP_Vetoes.json", "Bd2D0pipi_Vetoes" ) }; LauEffModel* effModelB0bar = new LauEffModel(daughtersB0bar, vetoes.get()); LauEffModel* effModelB0 = new LauEffModel(daughtersB0, vetoes.get()); TFile* effFile(nullptr); TH2* effHist(nullptr); if (settings.directory == "B2Dhh_Kpi"){ effFile = TFile::Open(Form("%s/%s/eff_map_mprime_thprime_model_d2kpi_run%u.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str(),runNo)); effHist = dynamic_cast(effFile->Get(Form("eff_map_mprime_thprime_model_d2kpi_run%u",runNo))); } else if (settings.directory == "B2Dhh_KK"){ effFile = TFile::Open(Form("%s/%s/eff_map_mprime_thprime_model_d2kk_run%u.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str(),runNo)); effHist = dynamic_cast(effFile->Get(Form("eff_map_mprime_thprime_model_d2kk_run%u",runNo))); } else if (settings.directory == "B2Dhh_pipi"){ effFile = TFile::Open(Form("%s/%s/eff_map_mprime_thprime_model_d2pipi_run%u.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str(),runNo)); effHist = dynamic_cast(effFile->Get(Form("eff_map_mprime_thprime_model_d2pipi_run%u",runNo))); } else {std::cerr << "bad dir given!" << std::endl; return EXIT_FAILURE;} effModelB0->setEffHisto(effHist, kTRUE, kFALSE, 0.0, 0.0, kFALSE, kTRUE); effModelB0bar->setEffHisto(effHist, kTRUE, kFALSE, 0.0, 0.0, kFALSE, kTRUE); // background types std::map BkgndTypes { {"comb", LauFlavTag::BkgndType::Combinatorial}, {"Bs2DKpi", LauFlavTag::BkgndType::FlavourSpecific}, {"Bd2DKpi", LauFlavTag::BkgndType::FlavourSpecific}, {"Lb2Dppi", LauFlavTag::BkgndType::FlavourSpecific} }; auto notInBkgndList = [&settings](const std::map::value_type& e) { return ( settings.bkgndList.find( e.first.Data() ) == settings.bkgndList.end() ); }; std::vector toRemove; for ( const auto& elem : BkgndTypes ) { if ( notInBkgndList(elem) ) { toRemove.push_back( elem.first ); } } for ( const auto& name : toRemove ) { BkgndTypes.erase( name ); } + std::map bkgndDecayTimeTypes = { + {"comb", LauDecayTime::FuncType::Hist}, + {"Bd2DKpi", LauDecayTime::FuncType::ExpTrig}, + {"Bs2DKpi", LauDecayTime::FuncType::ExpHypTrig}, + {"Lb2Dppi", LauDecayTime::FuncType::Exp} + }; + // setup flavour tagging const Bool_t useAveDeltaCalibVals { settings.useAveDeltaCalibVals }; const Bool_t useEtaPrime { kFALSE }; //LauFlavTag* flavTag = new LauFlavTag(useAveDeltaCalibVals,useEtaPrime); LauFlavTag* flavTag = new LauFlavTag(useAveDeltaCalibVals,useEtaPrime,BkgndTypes); if (settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS) { flavTag->setDecayFlvVarName("decayFlv"); } if(not settings.dataFit) { flavTag->setTrueTagVarName("trueTag"); } TFile* etaFile = TFile::Open(Form("%s/%s/ETAhists.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str())); TH1* etaHist{nullptr}; TH1* etaHistB0{nullptr}; TH1* etaHistB0bar{nullptr}; if (settings.directory == "B2Dhh_Kpi"){ etaHistB0 = dynamic_cast(etaFile->Get(Form("signalMC_Kpi_Run%u_etaHist_dcyB",runNo))); etaHistB0bar = dynamic_cast(etaFile->Get(Form("signalMC_Kpi_Run%u_etaHist_dcyBbar",runNo))); etaHist = dynamic_cast(etaHistB0->Clone(Form("signalMC_Kpi_Run%u_etaHist",runNo))); etaHist->Add(etaHistB0bar); } else if (settings.directory == "B2Dhh_KK"){ etaHist = dynamic_cast(etaFile->Get(Form("signalMC_KK_Run%u_etaHist",runNo))); } else if (settings.directory == "B2Dhh_pipi"){ etaHist = dynamic_cast(etaFile->Get(Form("signalMC_pipi_Run%u_etaHist",runNo))); } else {std::cerr << "bad dir given!" << std::endl; return EXIT_FAILURE;} // Crude check as to whether we're doing perfect vs realistic mis-tag // - in the former case all entries should be in the first bin // If the tagging is perfect then don't interpolate the eta histogram // and also make it perfectly efficient, otherwise do interpolate and // make it 50% efficient Lau1DHistPdf* etaHistPdf = new Lau1DHistPdf( "eta", etaHist, 0.0, 0.5, kTRUE, kFALSE ); const Double_t meanEta { etaHistPdf->getMean() }; const Double_t tagEffVal{1.0}; std::pair tagEff {tagEffVal, useAveDeltaCalibVals ? 0.0 : tagEffVal}; // use a null calibration for the time being, so p0 = and p1 = 1 std::pair calib0 {meanEta, useAveDeltaCalibVals ? 0.0 : meanEta}; std::pair calib1 {1.0, useAveDeltaCalibVals ? 0.0 : 1.0}; if(settings.addTaggers) { flavTag->addTagger(Form("IFT_Run%u",runNo), "IFT_TAGDEC", "IFT_TAGETA", etaHistPdf, tagEff, calib0, calib1); if (settings.floatCalibPars) { flavTag->floatAllCalibPars(); } } // flavour tagging setup for backgrounds std::map< TString, std::pair > BkgndEtas; TH1* bkgndEtaHist1{nullptr}; Lau1DHistPdf* bkgndEtaHistPDF1{nullptr}; TH1* bkgndEtaHistMinus1{nullptr}; Lau1DHistPdf* bkgndEtaHistPDFMinus1{nullptr}; for(auto& [name, _] : BkgndTypes) { const char* decayName = name == "comb" ? "sideband" : name.Data(); if (settings.directory == "B2Dhh_Kpi"){ bkgndEtaHist1 = dynamic_cast(etaFile->Get( Form("%s_Kpi_Run%u_etaHist_dcyB" ,decayName,runNo) )); bkgndEtaHistPDF1 = new Lau1DHistPdf("eta",bkgndEtaHist1 ,0.0,0.5,kTRUE,kFALSE); bkgndEtaHistMinus1 = dynamic_cast(etaFile->Get( Form("%s_Kpi_Run%u_etaHist_dcyBbar",decayName,runNo) )); bkgndEtaHistPDFMinus1 = new Lau1DHistPdf("eta",bkgndEtaHistMinus1,0.0,0.5,kTRUE,kFALSE); } else if (settings.directory == "B2Dhh_KK"){ bkgndEtaHist1 = dynamic_cast(etaFile->Get( Form("%s_KK_Run%u_etaHist" ,decayName,runNo) )); bkgndEtaHistPDF1 = new Lau1DHistPdf("eta",bkgndEtaHist1 ,0.0,0.5,kTRUE,kFALSE); } else if (settings.directory == "B2Dhh_pipi"){ bkgndEtaHist1 = dynamic_cast(etaFile->Get( Form("%s_pipi_Run%u_etaHist" ,decayName,runNo) )); bkgndEtaHistPDF1 = new Lau1DHistPdf("eta",bkgndEtaHist1 ,0.0,0.5,kTRUE,kFALSE); } else {std::cerr << "bad dir given!" << std::endl; return EXIT_FAILURE;} BkgndEtas.emplace( std::make_pair(name, std::make_pair(bkgndEtaHistPDF1,bkgndEtaHistPDFMinus1) ) ); } // NB deliberately not using uniform initialisation here because of this issue: // https://json.nlohmann.me/home/faq/#brace-initialization-yields-arrays const TString combTagAsymJsonFileName { "Bd2D0pipi_CombTagAsym.json" }; if ( not std::filesystem::exists(combTagAsymJsonFileName.Data()) ) { std::cerr << "Warning : couldn't find the json file of combinatorial tagging asymmetry parameters; attempting to fetch from eos ..." << std::endl; int r = system( Form("xrdcp %s/%s/%s .", settings.eosRoot.c_str(), settings.eosConfigDir.c_str(), combTagAsymJsonFileName.Data()) ); if (r != 0) { std::cerr << "FATAL : couldn't get the json file from remote!" << std::endl; return 1;} } const nlohmann::json combTagAsymJson = LauJsonTools::readJsonFile( combTagAsymJsonFileName.Data(), "", LauJsonTools::JsonType::Object ); if ( combTagAsymJson.is_null() ) { std::cerr << "ERROR : unable to retrieve JSON object from root element of file \"" << combTagAsymJsonFileName << "\"" << std::endl; return 1; } for(auto& [name,hists] : BkgndEtas) { if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) { if ( name == "comb" ) { const BgTagDecayFlvAsym bgt { combTagAsymJson.at(runStr).at(settings.directory) }; const std::pair effsB0 { useAveDeltaCalibVals ? bgt.aveDeltaB0() : bgt.decayB0() }; const std::pair effsB0bar { useAveDeltaCalibVals ? bgt.aveDeltaB0bar() : bgt.decayB0bar() }; flavTag->setBkgndParams( name, Form("IFT_Run%u",runNo), hists.first, effsB0, hists.second, effsB0bar ); } else { flavTag->setBkgndParams( name, Form("IFT_Run%u",runNo), hists.first, tagEff, hists.second, tagEff ); } } else if (settings.dType == LauTimeDepFitModel::CPEigenvalue::CPEven){ if ( name == "comb" ) { const BgTagDecayFlvAsym bgt { combTagAsymJson.at(runStr).at(settings.directory) }; const std::pair effs { useAveDeltaCalibVals ? bgt.aveDeltaB0() : bgt.decayB0() }; flavTag->setBkgndParams( name, Form("IFT_Run%u",runNo), hists.first, effs); } else { flavTag->setBkgndParams( name, Form("IFT_Run%u",runNo), hists.first, tagEff ); } } } // signal dynamics if ( not std::filesystem::exists("Bd2D0pipi_DP_Model_Coeffs.json") ) { std::cerr << "Warning : couldn't find the json file of DP model coeffs; attempting to fetch from eos ..." << std::endl; int r = system( Form("xrdcp %s/%s/Bd2D0pipi_DP_Model_Coeffs.json .", settings.eosRoot.c_str(), settings.eosConfigDir.c_str()) ); if (r != 0) { std::cerr << "FATAL : couldn't get the json file from remote!" << std::endl; return 1;} } LauIsobarDynamics* sigModelB0bar = new LauIsobarDynamics(daughtersB0bar, effModelB0bar); if (settings.directory == "B2Dhh_Kpi"){ sigModelB0bar->setIntFileName(Form("integ_Kpi_Run%u_B0bar.dat",runNo)); } else if (settings.directory == "B2Dhh_KK"){ sigModelB0bar->setIntFileName(Form("integ_KK_Run%u_B0bar.dat",runNo)); } else if (settings.directory == "B2Dhh_pipi"){ sigModelB0bar->setIntFileName(Form("integ_pipi_Run%u_B0bar.dat",runNo)); } sigModelB0bar->constructModelFromJson( "Bd2D0pipi_DP_Model_Coeffs.json", "Bd2D0pipi_B0bar_Model" ); LauIsobarDynamics* sigModelB0 = new LauIsobarDynamics(daughtersB0, effModelB0); if (settings.directory == "B2Dhh_Kpi"){ sigModelB0->setIntFileName(Form("integ_Kpi_Run%u_B0.dat",runNo)); } else if (settings.directory == "B2Dhh_KK"){ sigModelB0->setIntFileName(Form("integ_KK_Run%u_B0.dat",runNo)); } else if (settings.directory == "B2Dhh_pipi"){ sigModelB0->setIntFileName(Form("integ_pipi_Run%u_B0.dat",runNo)); } sigModelB0->constructModelFromJson( "Bd2D0pipi_DP_Model_Coeffs.json", "Bd2D0pipi_B0_Model" ); // fit model LauTimeDepFitModel* fitModel = new LauTimeDepFitModel(sigModelB0bar,sigModelB0,flavTag); fitModel->breakOnBadLikelihood(kFALSE); fitModel->sendFixedYieldsToFitter(kTRUE); std::vector> coeffset { LauAbsCoeffSet::readFromJson( "Bd2D0pipi_DP_Model_Coeffs.json", "Bd2D0pipi_Coeffs" ) }; for ( auto& coeff : coeffset ) { fitModel->setAmpCoeffSet( std::move(coeff) ); } // background DP models LauBkgndDPModel* BkgndModel{nullptr}; LauBkgndDPModel* BkgndModelbar{nullptr}; TFile* BkgndHistFile = TFile::Open(Form("%s/%s/BG_SDPs.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str())); TFile* BkgndHistFile2 = TFile::Open(Form("%s/%s/sidebandSDP.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str())); TH2* BkgndHistB0{nullptr}; TH2* BkgndHistB0bar{nullptr}; for(auto& [name, _] : BkgndTypes) { //Only have Kpi histos for now for misid if( name=="comb"){ if( settings.directory=="B2Dhh_Kpi"){ BkgndHistB0 = dynamic_cast(BkgndHistFile2-> Get(Form("sideband_Run%u_Kpi_SDP_decayFlv1",runNo))); BkgndHistB0bar = dynamic_cast(BkgndHistFile2->Get(Form("sideband_Run%u_Kpi_SDP_decayFlvMinus1",runNo))); } else if( settings.directory=="B2Dhh_KK"){ BkgndHistB0 = dynamic_cast(BkgndHistFile2->Get(Form("sideband_Run%u_KK_SDP",runNo))); } else if( settings.directory=="B2Dhh_pipi"){ BkgndHistB0 = dynamic_cast(BkgndHistFile2->Get(Form("sideband_Run%u_pipi_SDP",runNo))); } } else { TString mode = settings.directory; mode.ReplaceAll("B2Dhh_",""); TString histName = Form("%s_%s_Run%u_SDP_decFlv",name.Data(),mode.Data(),runNo); BkgndHistB0 = dynamic_cast(BkgndHistFile->Get(histName+"1")); BkgndHistB0bar = dynamic_cast(BkgndHistFile->Get(histName+"Minus1")); } BkgndModel = new LauBkgndDPModel(daughtersB0, vetoes.get()); BkgndHistB0->Smooth(); BkgndModel->setBkgndSpline(BkgndHistB0, kFALSE, kFALSE, kTRUE ); //For CP comb we don't split based on decay flavour since it's impossible if(name == "comb" and settings.directory!="B2Dhh_Kpi") {BkgndModelbar=nullptr;} else { BkgndModelbar = new LauBkgndDPModel(daughtersB0bar, vetoes.get()); BkgndHistB0bar->Smooth(); BkgndModelbar->setBkgndSpline(BkgndHistB0bar, kFALSE, kFALSE, kTRUE ); } fitModel->setBkgndDPModels( name, BkgndModel, BkgndModelbar ); } // decay type and mixing parameter const Bool_t fixPhiMix{ settings.fixPhiMix || settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS }; const Bool_t useSinCos{ settings.useSinCos }; fitModel->setCPEigenvalue( settings.dType ); fitModel->setPhiMix( 2.0*LauConstants::beta, fixPhiMix, useSinCos ); if ( useSinCos ) { if ( settings.blindFit ) { //Blinding strings generated from: https://www.thewordfinder.com/random-sentence-generator/ const TString sinBlindString{settings.dataFit ? "Tom realized he could be making a big mistake." : "Don't Look a Gift Horse In The Mouth"}; const TString cosBlindString{settings.dataFit ? "He admitted defeat." : "Long In The Tooth"}; fitModel->blindPhiMix( sinBlindString, cosBlindString ); } std::cout << "Is sin(phiMix) blinded? " << std::boolalpha << fitModel->sinPhiMixBlinded() << std::endl; std::cout << "Is cos(phiMix) blinded? " << std::boolalpha << fitModel->cosPhiMixBlinded() << std::endl; } else { if ( settings.blindFit ) { const TString blindString{"A Cut Above"}; fitModel->blindPhiMix( blindString ); } std::cout << "Is phiMix blinded? " << std::boolalpha << fitModel->phiMixBlinded() << std::endl; } // production asymmetries fitModel->setAsymmetries( 0.0, kTRUE ); for(auto& [name, _] : BkgndTypes) { fitModel->setBkgndAsymmetries( name, 0.0, kTRUE ); } // decay time PDFs const Double_t minDt{0.2}; const Double_t maxDt{15.0}; const Double_t minDtErr{0.0}; const Double_t maxDtErr{0.15}; // physics models auto dtPhysicsModels { LauDecayTimePhysicsModel::readFromJson( "Bd2D0pipi_DecayTime.json", "PhysicsModels" ) }; auto& dtPhysModel { dtPhysicsModels.at("signal") }; auto lifetimePar = dynamic_cast( dtPhysModel->findParameter("dt_sig_tau_Bd") ); if ( lifetimePar ) { lifetimePar->fixed( settings.fixLifetime ); } auto deltaMPar = dynamic_cast( dtPhysModel->findParameter("dt_sig_deltaM_Bd") ); if ( deltaMPar ) { deltaMPar->fixed( settings.fixDeltaM ); } // resolution models auto dtResolutionModels { LauDecayTimeResolution::readFromJson( "Bd2D0pipi_DecayTime.json", Form("ResolutionModels_Run%d", runNo) ) }; auto& dtResoModel { dtResolutionModels.at("signal") }; // Decay time error histogram // (always set this so that it gets generated properly, whether we're using it in the PDF or not) TFile* dtErrFile = TFile::Open(Form("%s/%s/dte-hist.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str())); TH1* dtErrHist = dynamic_cast(dtErrFile->Get("dte_hist")); LauDecayTimePdf * dtPdf{nullptr}; if ( settings.timeResolution ) { if ( settings.perEventTimeErr ) { std::cerr << "Warning : the resolution model has not been constructed with per-event error scaling, but you have requested it" << std::endl; dtPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, "decayTimeErr", minDtErr, maxDtErr, std::move(dtPhysModel), std::move(dtResoModel), dtErrHist ); } else { dtPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, std::move(dtPhysModel), std::move(dtResoModel) ); } } else { dtPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, std::move(dtPhysModel) ); } // Decay time acceptance histogram TFile* dtaFile = TFile::Open(Form("%s/%s/DTAhists.root",settings.eosRoot.c_str(),settings.eosConfigDir.c_str())); TH1* dtaHist{nullptr}; if (settings.directory == "B2Dhh_Kpi"){ dtaHist = dynamic_cast(dtaFile->Get(Form("signalMC_Kpi_Run%u_DTAhist",runNo))); } else if (settings.directory == "B2Dhh_KK"){ dtaHist = dynamic_cast(dtaFile->Get(Form("signalMC_KK_Run%u_DTAhist",runNo))); } else if (settings.directory == "B2Dhh_pipi"){ dtaHist = dynamic_cast(dtaFile->Get(Form("signalMC_pipi_Run%u_DTAhist",runNo))); } else {std::cerr << "bad dir given!" << std::endl; return EXIT_FAILURE;} if( not std::filesystem::exists("run1DTAsplines.json") and runNo == 1 ) { std::cerr << "Warning : couldn't find the json file of splines in Run 1; attempting to fetch from eos ..." << std::endl; int r = system( Form("xrdcp %s/%s/run1DTAsplines.json .", settings.eosRoot.c_str(), settings.eosConfigDir.c_str()) ); if(r != 0){ std::cerr << "FATAL : couldn't get the json file from remote!" << std::endl; return 1;} } if( not std::filesystem::exists("run2DTAsplines.json") and runNo == 2 ) { std::cerr << "Warning : couldn't find the json file of splines in Run 2; attempting to fetch from eos ..." << std::endl; int r = system( Form("xrdcp %s/%s/run2DTAsplines.json .", settings.eosRoot.c_str(), settings.eosConfigDir.c_str()) ); if(r != 0){ std::cerr << "FATAL : couldn't get the json file from remote!" << std::endl; return 1;} } - //The below is replaced with the new Lau1DCubicSpline::readFromJson - // Create the spline knot positions and - // starting Y values, to be fit to dtaHist - //const std::vector dtvals { 0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 7.0, 10.0, 15.0}; - //const std::vector effvals {0.000, 0.010, 0.022, 0.035, 0.042, 0.050, 0.051, 0.052, 0.055}; - //const std::vector corvals { 1.10, 1.08, 1.06, 1.04, 1.02, 1.00, 0.98, 0.96, 0.94}; + const TString jsonFileName { Form("run%uDTAsplines.json",runNo) }; + const TString suffix { "_DTAhist" }; - std::unique_ptr dtEffSpline = nullptr; - Lau1DCubicSpline* dtEffSplinePtr = nullptr; - std::unique_ptr dtCorrSpline = nullptr; - Lau1DCubicSpline* dtCorrSplinePtr = nullptr; - - LauSplineDecayTimeEfficiency* sigDTA = nullptr; + const TString signalName { "signalMC" }; + std::vector bkgndNames; + bkgndNames.reserve( bkgndDecayTimeTypes.size() ); + for ( auto& [name, type] : bkgndDecayTimeTypes ) { + if ( type != LauDecayTime::FuncType::Hist ) { + bkgndNames.push_back( name ); + } + } - std::unique_ptr> dtaModel; + auto [ qfsDTAModel, cpDTAModel, bkgndDTAModels ] = createSplineDTASet( jsonFileName, dModeStr, runStr, suffix, signalName, bkgndNames ); switch(settings.timeEffModel) { case LauDecayTime::EfficiencyMethod::Spline: { fitModel->setASqMaxValue(0.01); if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) { fitModel->setASqMaxValue(0.009); } - //auto dtEffSpline = std::make_unique( dtvals, effvals, Lau1DCubicSpline::SplineType::AkimaSpline ); - - TString jsonFilename = Form("run%uDTAsplines.json",runNo); - dtEffSpline = Lau1DCubicSpline::readFromJson(jsonFilename, Form("signalMC_Kpi_Run%u_DTAhist",runNo)); - dtEffSplinePtr = dtEffSpline.get(); - - dtaModel = std::make_unique>( Form("dteff_QFS_Run%u",runNo), std::move(dtEffSpline) ); - dynamic_cast(dtaModel->getParameters().front())->maxValue(1e-1); - //dtaModel->fitToTH1(dtaHist); + dynamic_cast(qfsDTAModel->getParameters().front())->maxValue(1e-1); + //qfsDTAModel->fitToTH1(dtaHist); // Set which knots to float and which to fix (at least 1 knot must be fixed, not the first one) // Knots should only be floating if requested AND the B lifetime is fixed! if ( settings.fixSplineKnots or not settings.fixLifetime ) { - dtaModel->fixKnots(); + qfsDTAModel->fixKnots(); } else { - dtaModel->floatKnots(); - dtaModel->fixKnot( 0, true ); - dtaModel->fixKnot( 3, true ); + qfsDTAModel->floatKnots(); + qfsDTAModel->fixKnot( 0, true ); + qfsDTAModel->fixKnot( 3, true ); } - sigDTA = dtaModel.get(); + if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) { - // For the QFS mode we just use the cubic model as it is - sigDTA = dtaModel.get(); - dtPdf->setSplineEfficiency( std::move(dtaModel) ); + dtPdf->setSplineEfficiency( std::move(qfsDTAModel) ); } else { - // For the CP modes we modify it using a corrective spline - // Json available? - //auto dtCorrSpline = std::make_unique( dtvals, corvals, Lau1DCubicSpline::SplineType::AkimaSpline ); - std::unique_ptr> dtaCPModel; - if (settings.directory == "B2Dhh_KK"){ - dtCorrSpline = Lau1DCubicSpline::readFromJson(jsonFilename, Form("signalMC_KK_Run%u_DTAhist",runNo)); - dtCorrSplinePtr = dtCorrSpline.get(); - dtaCPModel = std::make_unique>( Form("dteff_CPKK_Run%u",runNo), *dtaModel, std::move( dtCorrSpline ) ); - } else if (settings.directory == "B2Dhh_pipi"){ - dtCorrSpline = Lau1DCubicSpline::readFromJson(jsonFilename, Form("signalMC_pipi_Run%u_DTAhist",runNo)); - dtCorrSplinePtr = dtCorrSpline.get(); - dtaCPModel = std::make_unique>( Form("dteff_CPpipi_Run%u",runNo), *dtaModel, std::move( dtCorrSpline ) ); - } - dtPdf->setSplineEfficiency( std::move(dtaCPModel) ); + dtPdf->setSplineEfficiency( std::move(cpDTAModel) ); } break; } case LauDecayTime::EfficiencyMethod::Binned: { fitModel->setASqMaxValue(0.06); auto dtaBinnedModel = std::make_unique( *dtaHist ); dtPdf->setBinnedEfficiency( std::move(dtaBinnedModel) ); break; } case LauDecayTime::EfficiencyMethod::Uniform: { fitModel->setASqMaxValue(4.45); break; } } fitModel->setSignalDtPdf( dtPdf ); - std::map bkgndDecayTimeTypes = { - {"comb", LauDecayTime::FuncType::Hist}, - {"Bd2DKpi", LauDecayTime::FuncType::ExpTrig}, - {"Bs2DKpi", LauDecayTime::FuncType::ExpHypTrig}, - {"Lb2Dppi", LauDecayTime::FuncType::Exp} - }; - - TH1* bgDTA{nullptr}; - std::map > decayFlvYields = { - {"comb" , runNo == 1 ? std::make_pair( 1111, 1044) : std::make_pair( 3011, 2795)}, - {"Bd2DKpi", runNo == 1 ? std::make_pair( 9552, 8933) : std::make_pair(13842,12036)}, - {"Bs2DKpi", runNo == 1 ? std::make_pair(16989,15634) : std::make_pair(18294,17005)}, - {"Lb2Dppi", runNo == 1 ? std::make_pair( 3695, 3431) : std::make_pair( 2389, 2501)} - }; - TString DTAnameSuffix{""}; - - if (settings.directory == "B2Dhh_Kpi"){ - DTAnameSuffix = Form("_Kpi_Run%u_DTAhist",runNo); - } else if (settings.directory == "B2Dhh_KK"){ - DTAnameSuffix = Form("_KK_Run%u_DTAhist",runNo); - } else if (settings.directory == "B2Dhh_pipi"){ - DTAnameSuffix = Form("_pipi_Run%u_DTAhist",runNo); - } + const TString dtaHistSuffix { Form("_%s_%s_DTAhist",dModeStr.Data(),runStr.c_str()) }; for (auto& [bg, _] : BkgndTypes) { if ( bkgndDecayTimeTypes[bg] == LauDecayTime::FuncType::Hist ) { - TString histname = bg + DTAnameSuffix; + TString histname = bg + dtaHistSuffix; if(histname.Contains("comb")) { - histname = "Sideband" + DTAnameSuffix; + histname = "Sideband" + dtaHistSuffix; histname.ReplaceAll("DTAhist","DThist"); } TH1* dt_hist = dynamic_cast(dtaFile->Get(histname.Data())); TH1* dt_err_hist = nullptr; // TODO get these! if ( not dt_hist ) { std::cerr << "PROBLEM FINDING THE LIFETIME HIST:" << histname << std::endl; return EXIT_FAILURE; } if ( settings.perEventTimeErr and not dt_err_hist ) { std::cerr << "PROBLEM FINDING THE LIFETIME ERROR HIST:" << bg << std::endl; return EXIT_FAILURE; } LauDecayTimePdf* backgroundPdf {nullptr}; if ( settings.perEventTimeErr ) { backgroundPdf = new LauDecayTimePdf("decayTime", minDt, maxDt, "decayTimeErr", minDtErr, maxDtErr, dt_hist, dt_err_hist ); } else { backgroundPdf = new LauDecayTimePdf("decayTime", minDt, maxDt, dt_hist ); } fitModel->setBkgndDtPdf(bg,backgroundPdf); } else { auto& bgPhysModel { dtPhysicsModels.at(bg) }; auto& bgResoModel { dtResolutionModels.at(bg) }; + auto& bgDTAModel { bkgndDTAModels.at(bg) }; TH1* dt_err_hist = nullptr; // TODO get these! if ( settings.perEventTimeErr and not dt_err_hist ) { std::cerr << "PROBLEM FINDING THE LIFETIME ERROR HIST:" << bg << std::endl; return EXIT_FAILURE; } LauDecayTimePdf* backgroundPdf {nullptr}; if ( settings.timeResolution ) { if ( settings.perEventTimeErr ) { std::cerr << "Warning : the resolution model has not been constructed with per-event error scaling, but you have requested it" << std::endl; backgroundPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, "decayTimeErr", minDtErr, maxDtErr, std::move(bgPhysModel), std::move(bgResoModel), dt_err_hist ); } else { backgroundPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, std::move(bgPhysModel), std::move(bgResoModel) ); } } else { backgroundPdf = new LauDecayTimePdf( "decayTime", minDt, maxDt, std::move(bgPhysModel) ); } // background decay time acceptance histogram - TString histname = bg + DTAnameSuffix; - bgDTA = dynamic_cast(dtaFile->Get(histname)); + TString histname = bg + dtaHistSuffix; + auto bgDTA = dynamic_cast(dtaFile->Get(histname)); if ( not bgDTA and bg!="comb" ) { std::cerr << "Couldn't find DTA histogram for background: " << bg << std::endl; return EXIT_FAILURE; } //XXX if we decide to rebin, do it here bgDTA->Divide(dtaHist); for (Int_t bin = 1; bin <= bgDTA->GetNbinsX(); ++bin) { if( bgDTA->GetBinContent(bin)==0. ){bgDTA->SetBinContent(bin,1.);} } - TString jsonFilename = Form("run%iDTAsplines.json",runNo); - auto dtBkgndCorrSpline = Lau1DCubicSpline::readFromJson(jsonFilename,histname); - - std::unique_ptr> bgDTAModel; - if (sigDTA==nullptr){ - std::cout << "Null pointer" << std::endl; - } - if (settings.directory == "B2Dhh_Kpi"){ - bgDTAModel = std::make_unique>( Form("%s_DTA_ratio_Run%u_Kpi",bg.Data(), runNo), *sigDTA, std::move( dtBkgndCorrSpline ) ); - } else if (settings.directory == "B2Dhh_KK"){ - bgDTAModel = std::make_unique>( Form("%s_DTA_ratio_Run%u_KK",bg.Data(), runNo), *sigDTA, std::move( dtBkgndCorrSpline ) ); - } else if (settings.directory == "B2Dhh_pipi"){ - bgDTAModel = std::make_unique>( Form("%s_DTA_ratio_Run%u_pipi",bg.Data(),runNo), *sigDTA, std::move( dtBkgndCorrSpline ) ); - } - //bgDTAModel->fitToTH1(bgDTA, LauOutputLevel::Quiet, true); backgroundPdf->setSplineEfficiency( std::move(bgDTAModel) ); fitModel->setBkgndDtPdf(bg,backgroundPdf); } } // Signal and bkgnd yields Double_t nSigEvents{0}; Double_t nCombEvents{0}; std::map bkgndYields; if (settings.directory == "B2Dhh_Kpi") { nSigEvents = (runNo == 1) ? 28362 : 115032; nCombEvents = (runNo == 1) ? 2415 : 7637; if ( ! settings.scaleYields ) { nCombEvents += (runNo == 1) ? (303 + 111 + 123 + 143 + 2) : (974 + 390 + 256 + 389 + 4); } bkgndYields = {{"comb", nCombEvents }, {"Bs2DKpi", runNo == 1 ? 796 : 1713 }, {"Bd2DKpi", runNo == 1 ? 254 : 639 }, {"Lb2Dppi", runNo == 1 ? 245 : 304 }}; } else if (settings.directory == "B2Dhh_KK"){ if (settings.setNegToZero and runNo == 1) { // Mass fit with mis-ID yields fixed to 0 nSigEvents = 3098; nCombEvents = 800; if( ! settings.scaleYields){ nCombEvents += 20 + 22 + 53 + 39 + 0; } bkgndYields = {{"comb", nCombEvents }, {"Bs2DKpi", 0 }, {"Bd2DKpi", 0 }, {"Lb2Dppi", 0 }}; } else if (settings.constrainedYields and runNo == 1) { // Mass fit with mis-ID yields constrained based on // Run1/Run2 in Kpi and pipi and Run2 yield of KK nSigEvents = 3093; nCombEvents = 777; if( ! settings.scaleYields){ nCombEvents += 19 + 21 + 52 + 41 + 0; } bkgndYields = {{"comb", nCombEvents }, {"Bs2DKpi", 19 }, {"Bd2DKpi", 6 }, {"Lb2Dppi", 7 }}; } else { // Original mass fit nSigEvents = runNo == 1 ? 3140 : 12099; nCombEvents = (runNo == 1) ? 1016 : 1501; if ( ! settings.scaleYields ) { nCombEvents += (runNo == 1) ? (63 + 23 + 26 + 13 + 1) : (107 + 43 + 28 + 5 + 1); } bkgndYields = {{"comb", nCombEvents }, {"Bs2DKpi", runNo == 1 ? -172 : 167 }, {"Bd2DKpi", runNo == 1 ? -64 : 73 }, {"Lb2Dppi", runNo == 1 ? -61 : 33 }}; } } else if (settings.directory == "B2Dhh_pipi"){ nSigEvents = runNo == 1 ? 1175 : 4566; nCombEvents = (runNo == 1) ? 115 : 301; if ( ! settings.scaleYields ) { nCombEvents += (runNo == 1) ? (7 + 3 + 3 + 79 + 0) : (36 + 15 + 10 + 138 + 1); } bkgndYields = {{"comb", nCombEvents }, {"Bs2DKpi", runNo == 1 ? 57 : 99 }, {"Bd2DKpi", runNo == 1 ? 18 : 41 }, {"Lb2Dppi", runNo == 1 ? 21 : 19 }}; } else {std::cerr << "bad dir given!" << std::endl; return EXIT_FAILURE;} if(settings.command == Command::Generate) { for(auto& [mode, yield] : bkgndYields) { if(yield < 0.) { std::cerr << "WARNING : setting negative yield for mode " << mode << " to 0 for generation" << std::endl; yield = 0.; } } } const Bool_t fixYields { settings.bkgndList.empty() or not settings.floatYields }; TString yieldName; for (auto& [bg, _] : BkgndTypes) { const Double_t nBkg = bkgndYields[bg]; std::cout<setNBkgndEvents( nBkgndEvents ); } std::cout<<"nSigEvents = "<setNSigEvents(nSigPar); if ( settings.command == Command::Generate or not fixYields ) { // mB PDFs // PDFs for B+ and B- are set to be identical here but could be different in general const TString mBName { "mB" }; const Double_t mBMin { 5.200 }; const Double_t mBMax { 5.700 }; // Signal PDF is a double Crystal Ball function const std::map sigParVals { {"mB_sig_mean", 5.2802}, {"mB_sig_sigma", 0.01125}, {"mB_sig_alphaL", 1.226}, {"mB_sig_alphaR", -1.722}, {"mB_sig_orderL", 1.817}, {"mB_sig_orderR", 2.861}, {"mB_sig_frac", 0.480} }; fitModel->setSignalPdfs( makeDCBPDF( mBName, mBMin, mBMax, "sig", sigParVals ) ); // Combinatorial PDF is an exponential if ( BkgndTypes.find("comb") != BkgndTypes.end() ) { LauParameter* mB_comb_slope = new LauParameter("mB_comb_slope", -0.00419); LauAbsPdf* mB_comb_pdf = new LauExponentialPdf(mBName, {mB_comb_slope}, mBMin, mBMax); fitModel->setBkgndPdf("comb", mB_comb_pdf); } // Bd2DKpi PDF is double Crystal Ball function if ( BkgndTypes.find("Bd2DKpi") != BkgndTypes.end() ) { const std::map Bd2DKpiParVals { {"mB_Bd2DKpi_mean", 5.2401}, {"mB_Bd2DKpi_sigma", 0.01580}, {"mB_Bd2DKpi_alphaL", 0.560}, {"mB_Bd2DKpi_alphaR", -2.212}, {"mB_Bd2DKpi_orderL", 1.076}, {"mB_Bd2DKpi_orderR", 2.453}, {"mB_Bd2DKpi_frac", 0.544} }; fitModel->setBkgndPdf("Bd2DKpi", makeDCBPDF( mBName, mBMin, mBMax, "Bd2DKpi", Bd2DKpiParVals ) ); } // Bs2DKpi PDF is double Crystal Ball function if ( BkgndTypes.find("Bs2DKpi") != BkgndTypes.end() ) { const std::map Bs2DKpiParVals { {"mB_Bs2DKpi_mean", 5.3244}, {"mB_Bs2DKpi_sigma", 0.01827}, {"mB_Bs2DKpi_alphaL", 0.865}, {"mB_Bs2DKpi_alphaR", -0.001}, {"mB_Bs2DKpi_orderL", 10.000}, {"mB_Bs2DKpi_orderR", 10.000}, {"mB_Bs2DKpi_frac", 0.954} }; fitModel->setBkgndPdf("Bs2DKpi", makeDCBPDF( mBName, mBMin, mBMax, "Bs2DKpi", Bs2DKpiParVals ) ); } // Lb2Dppi PDF is double Crystal Ball function if ( BkgndTypes.find("Lb2Dppi") != BkgndTypes.end() ) { const std::map Lb2DppiParVals { {"mB_Lb2Dppi_mean", 5.4884}, {"mB_Lb2Dppi_sigma", 0.02561}, {"mB_Lb2Dppi_alphaL", 0.091}, {"mB_Lb2Dppi_alphaR", -18.053}, {"mB_Lb2Dppi_orderL", 2.906}, {"mB_Lb2Dppi_orderR", 6.087}, {"mB_Lb2Dppi_frac", 0.968} }; fitModel->setBkgndPdf("Lb2Dppi", makeDCBPDF( mBName, mBMin, mBMax, "Lb2Dppi", Lb2DppiParVals ) ); } } // set the number of experiments if (settings.command == Command::Generate) { fitModel->setNExpts(settings.nExptGen, settings.firstExptGen); } else { fitModel->setNExpts(settings.nExptFit, settings.firstExptFit); } fitModel->useAsymmFitErrors(kFALSE); fitModel->useRandomInitFitPars(kFALSE); fitModel->writeLatexTable(kFALSE); const Bool_t haveBkgnds = ( fitModel->nBkgndClasses() > 0 ); fitModel->doPoissonSmearing(haveBkgnds); fitModel->doEMLFit(haveBkgnds); TString dTypeStr; switch (settings.dType) { case LauTimeDepFitModel::CPEigenvalue::CPEven : dTypeStr = "CPEven"; break; case LauTimeDepFitModel::CPEigenvalue::CPOdd : dTypeStr = "CPOdd"; break; case LauTimeDepFitModel::CPEigenvalue::QFS : dTypeStr = "QFS"; break; } TString dataFile(""); TString treeName("fitTree"); TString rootFileName(""); TString tableFileName(""); TString fitToyFileName(""); TString splotFileName(""); dataFile = "TEST-Dpipi"; dataFile += "_"+dTypeStr+"_"+settings.directory; switch(settings.timeEffModel) { case LauDecayTime::EfficiencyMethod::Spline: dataFile += "_Spline"; break; case LauDecayTime::EfficiencyMethod::Binned: dataFile += "_Binned"; break; case LauDecayTime::EfficiencyMethod::Uniform: dataFile += "_Uniform"; break; } if (settings.timeResolution) { if (settings.perEventTimeErr) { dataFile += "_DTRperevt"; } else { dataFile += "_DTRavg"; } } else { dataFile += "_DTRoff"; } dataFile += "_expts"; dataFile += settings.firstExptGen; dataFile += "-"; dataFile += settings.firstExptGen+settings.nExptGen-1; dataFile += Form("-Run%u",runNo); dataFile += ".root"; if (settings.dataFit) { if(runNo == 2) { dataFile = Form("%s/%s/Run2/B2Dhh-CollisionCombined-MagCombined-Stripping24r2-withDNN-withPIDcorr-withMVA-withPIDMVA-WithDpiMatching-withSwappedMassHypotheses-withKKvetoes-withSelection-LauraPrepped.root",settings.eosRoot.c_str(),settings.eosDataDir.c_str()); } else { dataFile = Form("%s/%s/Run1/B2Dhh-CollisionCombined-MagCombined-Stripping21r1p2-withDNN-withPIDcorr-withMVA-withPIDMVA-WithDpiMatching-withSwappedMassHypotheses-withKKvetoes-withSelection-LauraPrepped.root",settings.eosRoot.c_str(),settings.eosDataDir.c_str()); } treeName = Form("%s/B2DhhReco",settings.directory.data()); } if (settings.command == Command::Generate) { rootFileName = "dummy.root"; tableFileName = "genResults"; } else { rootFileName = "fit"; rootFileName += settings.iFit; rootFileName += "_Results_"; rootFileName += dTypeStr; rootFileName += "_"; rootFileName += settings.directory; rootFileName += "_expts"; rootFileName += settings.firstExptFit; rootFileName += "-"; rootFileName += settings.firstExptFit+settings.nExptFit-1; rootFileName += Form("_Run%u",runNo); rootFileName += ".root"; tableFileName = "fit"; tableFileName += settings.iFit; tableFileName += "_Results_"; tableFileName += dTypeStr; tableFileName += "_"; tableFileName += settings.directory; tableFileName += "_expts"; tableFileName += settings.firstExptFit; tableFileName += "-"; tableFileName += settings.firstExptFit+settings.nExptFit-1; tableFileName += Form("_Run%u",runNo); fitToyFileName = "fit"; fitToyFileName += settings.iFit; fitToyFileName += "_ToyMC_"; fitToyFileName += dTypeStr; fitToyFileName += "_"; fitToyFileName += settings.directory; fitToyFileName += "_expts"; fitToyFileName += settings.firstExptFit; fitToyFileName += "-"; fitToyFileName += settings.firstExptFit+settings.nExptFit-1; fitToyFileName += Form("_Run%u",runNo); fitToyFileName += ".root"; splotFileName = "fit"; splotFileName += settings.iFit; splotFileName += "_sPlot_"; splotFileName += dTypeStr; splotFileName += "_"; splotFileName += settings.directory; splotFileName += "_expts"; splotFileName += settings.firstExptFit; splotFileName += "-"; splotFileName += settings.firstExptFit+settings.nExptFit-1; splotFileName += Form("_Run%u",runNo); splotFileName += ".root"; } if( settings.genToy and not settings.blindFit ) { fitModel->compareFitData(50, fitToyFileName); } // Generate toy from the fitted parameters //fitModel->compareFitData(1, fitToyFileName); // Write out per-event likelihoods and sWeights //fitModel->writeSPlotData(splotFileName, "splot", kFALSE); // Execute the generation/fit switch (settings.command) { case Command::Generate : fitModel->run( "gen", dataFile, treeName, rootFileName, tableFileName ); break; case Command::Fit : fitModel->run( "fit", dataFile, treeName, rootFileName, tableFileName ); break; case Command::SimFit : fitModel->runTask( dataFile, treeName, rootFileName, tableFileName, "localhost", settings.port ); break; } return EXIT_SUCCESS; } LauAbsPdf* makeDCBPDF( const TString& varName, const Double_t varMin, const Double_t varMax, const TString& componentName, const std::map paramVals ) { const TString meanName { Form( "%s_%s_mean", varName.Data(), componentName.Data() ) }; const TString sigmaName { Form( "%s_%s_sigma", varName.Data(), componentName.Data() ) }; const TString alphaLName { Form( "%s_%s_alphaL", varName.Data(), componentName.Data() ) }; const TString alphaRName { Form( "%s_%s_alphaR", varName.Data(), componentName.Data() ) }; const TString orderLName { Form( "%s_%s_orderL", varName.Data(), componentName.Data() ) }; const TString orderRName { Form( "%s_%s_orderR", varName.Data(), componentName.Data() ) }; const TString fracName { Form( "%s_%s_frac", varName.Data(), componentName.Data() ) }; LauParameter* mean = new LauParameter( meanName, paramVals.at( meanName ) ); LauParameter* sigma = new LauParameter( sigmaName, paramVals.at( sigmaName ) ); LauParameter* alphaL = new LauParameter( alphaLName, paramVals.at( alphaLName ) ); LauParameter* alphaR = new LauParameter( alphaRName, paramVals.at( alphaRName ) ); LauParameter* orderL = new LauParameter( orderLName, paramVals.at( orderLName ) ); LauParameter* orderR = new LauParameter( orderRName, paramVals.at( orderRName ) ); LauParameter* frac = new LauParameter( fracName, paramVals.at( fracName ) ); std::vector pars; pars.reserve(4); pars.clear(); pars.push_back(mean); pars.push_back(sigma); pars.push_back(alphaL); pars.push_back(orderL); LauAbsPdf* pdf_L = new LauCrystalBallPdf(varName, pars, varMin, varMax); pars.clear(); pars.push_back(mean); pars.push_back(sigma); pars.push_back(alphaR); pars.push_back(orderR); LauAbsPdf* pdf_R = new LauCrystalBallPdf(varName, pars, varMin, varMax); LauAbsPdf* pdf = new LauSumPdf(pdf_L, pdf_R, frac); return pdf; }