diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 521aac2..66531ac 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,49 +1,50 @@ add_subdirectory(ProgOpts) list(APPEND EXAMPLE_SOURCES B3piKMatrixMassProj B3piKMatrixPlots CalcChiSq GenFit3K GenFit3KS GenFit3pi GenFitBelleCPKpipi GenFitDpipi GenFitDs2KKpi GenFitEFKLLM GenFitIncoherent_Bs2KSKpi GenFitKpipi GenFitNoDP GenFitNoDPMultiDim GenFitTimeDep GenFitTimeDep_Bs2KSKpi + GenFitTimeDep_Bd2Dpipi KMatrixDto3pi KMatrixExample MergeDataFiles mixedSampleTest PlotKMatrixTAmp PlotResults point2PointTestSample QuasiFlatSqDalitz QuasiFlatSqDalitz-CustomMasses ResultsExtractor SimFitCoordinator SimFitTask SimFitTaskRooFit ) if(NOT LAURA_BUILD_ROOFIT_TASK) list(REMOVE_ITEM EXAMPLE_SOURCES SimFitTaskRooFit) endif() foreach( _example ${EXAMPLE_SOURCES}) add_executable(${_example} ${_example}.cc) target_link_libraries(${_example} PRIVATE Laura++ ProgramOptions) install(TARGETS ${_example} RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) endforeach() # Also install the python script version of GenFit3pi configure_file(GenFit3pi.py.in GenFit3pi.py @ONLY) install(PROGRAMS ${CMAKE_CURRENT_BINARY_DIR}/GenFit3pi.py DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/examples/GenFitTimeDep_Bd2Dpipi.cc b/examples/GenFitTimeDep_Bd2Dpipi.cc new file mode 100644 index 0000000..649ad21 --- /dev/null +++ b/examples/GenFitTimeDep_Bd2Dpipi.cc @@ -0,0 +1,1082 @@ + +#include +#include +#include +#include +#include +#include + +#include "TFile.h" +#include "TH2.h" +#include "TRandom.h" +#include "TString.h" +#include "TSystem.h" +#include "TF1.h" +#include "TCanvas.h" + +#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 "TimeDep_Dpipi_ProgOpts.hh" + +// 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 controlSpline = Lau1DCubicSpline::readFromJson( jsonFileName, controlSplineName ); + LauCubicSplineDecayTimeEfficiencyPtr controlEff { std::make_unique>( controlSplineName, std::move( controlSpline ) ) }; + + LauSixthSplineDecayTimeEfficiencyPtr signalEff; + if ( signalSplineName != "" ) { + auto signalSpline = Lau1DCubicSpline::readFromJson( jsonFileName, signalSplineName ); + signalEff = std::make_unique>( signalSplineName, *controlEff, std::move( signalSpline ) ); + } + + std::vector bkgndEffs; + bkgndEffs.reserve( bkgndSplineNames.size() ); + for ( auto& bkgndSplineName : bkgndSplineNames ) { + auto bkgndSpline = Lau1DCubicSpline::readFromJson( jsonFileName, bkgndSplineName ); + bkgndEffs.push_back( std::make_unique>( bkgndSplineName, *controlEff, std::move( bkgndSpline ) ) ); + } + + 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} + {} + + 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; + 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 + LauVetoes* vetoes = new LauVetoes(); + LauEffModel* effModelB0bar = new LauEffModel(daughtersB0bar, vetoes); + LauEffModel* effModelB0 = new LauEffModel(daughtersB0, vetoes); + + 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 ); + } + + // 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) ) ); + } + + for(auto& [name,hists] : BkgndEtas) + { + if (settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS){ + if(name == "comb") + { + const BgTagDecayFlvAsym bgt = + runNo == 1 ? BgTagDecayFlvAsym(535.,576.,456.,588.) : BgTagDecayFlvAsym(1514.,1497.,1220.,1575.); + if(useAveDeltaCalibVals) + { + flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), + hists.first , bgt.aveDeltaB0(), + hists.second , bgt.aveDeltaB0bar() ); + } + else + { + flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), + hists.first , bgt.decayB0(), + hists.second , bgt.decayB0bar()); + } + } + else{flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), hists.first, tagEff, hists.second, tagEff );} + //flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), hist, name == "comb" ? std::make_pair(0.5,0.038) : tagEff ); + } else if (settings.dType == LauTimeDepFitModel::CPEigenvalue::CPEven){ + if(name == "comb") + { + BgTagDecayFlvAsym* bgt = nullptr; + if(settings.directory == "B2Dhh_KK") + { + bgt = runNo == 1 ? new BgTagDecayFlvAsym(344.,442.,344.,442.) : new BgTagDecayFlvAsym(710.,814.,710.,814.); + } + else + { + bgt = runNo == 1 ? new BgTagDecayFlvAsym(55.,68.,55.,68.) : new BgTagDecayFlvAsym(210.,241.,210.,241.); + } + + std::pair effs{bgt->decayB0()}; + if(useAveDeltaCalibVals){effs = bgt->aveDeltaB0();} + flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), hists.first, effs); + delete bgt; + } + else{flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), hists.first, tagEff);} +// else{flavTag->setBkgndParams(name, Form("IFT_Run%u",runNo), hists.first, tagEff, hists.second, tagEff );} FIXME once we have hists split by decay flv, uncomment this line and remove the above... actually don't FIXME, this is fine as above since the data isn't flavour specific in the CP modes. Keep it combined + } + } + + // 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); + 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); + 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.200103); + //const Double_t maxDt(14.9998); + const Double_t minDt(0.2); + const Double_t maxDt(15.0); + const Double_t minDtErr(0.0); + const Double_t maxDtErr(0.15); + + LauParameter * tauB0 = new LauParameter("dt_tau", 1.519, 0.5, 5.0, settings.fixLifetime); + LauParameter * deltaMd = new LauParameter("dt_deltaM", 0.5065, 0.0, 1.0, settings.fixDeltaM); + + std::vector dtPhysPars { + tauB0, + deltaMd + }; + + auto dtPhysModel = std::make_unique( LauDecayTime::FuncType::ExpTrig, dtPhysPars ); + + const std::vector scale { + false, false, false + // settings.perEventTimeErr && true, + // settings.perEventTimeErr && true, + }; + const std::size_t nGauss{scale.size()}; + + LauParameter * mean0 = new LauParameter("dt_mean_2", runNo == 1 ? -0.00174035 : -0.00218543 , -0.01, 0.01, kTRUE ); + LauParameter * mean1 = new LauParameter("dt_mean_1", runNo == 1 ? 0.043133 : 0.046862 , -0.10, 0.10, kTRUE ); + LauParameter * mean2 = new LauParameter("dt_mean_0", runNo == 1 ? -0.00108159 : -0.00159348 , -0.01, 0.01, kTRUE ); + LauParameter * sigma0 = new LauParameter("dt_sigma_2", runNo == 1 ? 0.046272 : 0.046064 , 0.0, 2.0, kTRUE ); + LauParameter * sigma1 = new LauParameter("dt_sigma_1", runNo == 1 ? 0.11815 : 0.11870 , 0.0, 2.5, kTRUE ); + LauParameter * sigma2 = new LauParameter("dt_sigma_0", runNo == 1 ? 0.025127 : 0.025077 , 0.0, 2.5, kTRUE ); + LauParameter * frac0 = new LauParameter("dt_frac_2", runNo == 1 ? 0.44359 : 0.43327 , 0.0, 1.0, kTRUE ); + LauParameter * frac1 = new LauParameter("dt_frac_1", runNo == 1 ? 0.045545 : 0.045310 , 0.0, 1.0, kTRUE ); + + std::vector dtResoPars { + mean0, + mean1, + mean2, + sigma0, + sigma1, + sigma2, + frac0, + frac1 + }; + + auto dtResoModel = std::make_unique( nGauss, dtResoPars, scale ); + + //Bs backgrounds + LauParameter * tauBs = new LauParameter("dt_tau_Bs", 1.527, 1.526, 1.528, kTRUE); + LauParameter * deltaMs = new LauParameter("dt_deltaM_Bs", 17.765, 17.60, 17.80, kTRUE); + LauParameter * deltaGs = new LauParameter("dt_deltaGamma_Bs", 0.084, 0.080, 0.090, kTRUE); + + //Lb backgrounds + LauParameter * tauLb = new LauParameter("dt_tau_Lb", 1.464, 1.450, 1.470, kTRUE); + + // 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 ) { + 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}; + + std::unique_ptr dtEffSpline = nullptr; + Lau1DCubicSpline* dtEffSplinePtr = nullptr; + std::unique_ptr dtCorrSpline = nullptr; + Lau1DCubicSpline* dtCorrSplinePtr = nullptr; + + LauSplineDecayTimeEfficiency* sigDTA = nullptr; + + std::unique_ptr> dtaModel; + + 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); + + // 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(); + } else { + dtaModel->floatKnots(); + dtaModel->fixKnot( 0, true ); + dtaModel->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) ); + } 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) ); + } + + 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); + } + + for (auto& [bg, _] : BkgndTypes) { + + if ( bkgndDecayTimeTypes[bg] == LauDecayTime::FuncType::Hist ) { + + TString histname = bg + DTAnameSuffix; + if(histname.Contains("comb")) + { + histname = "Sideband" + DTAnameSuffix; + 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 { + + std::vector bgPhysParams; + if (bg == "Bd2DKpi") { + bgPhysParams = + { + tauB0->createClone(), + deltaMd->createClone() + }; + } + else if (bg == "Bs2DKpi") { + bgPhysParams = + { + tauBs, + deltaGs, + deltaMs + }; + } + else if (bg == "Lb2Dppi") { + bgPhysParams = + { + tauLb + }; + } + auto bgPhysModel = std::make_unique( bkgndDecayTimeTypes[bg], bgPhysParams ); + + std::vector bgResoParams { + mean0->createClone(), + mean1->createClone(), + mean2->createClone(), + sigma0->createClone(), + sigma1->createClone(), + sigma2->createClone(), + frac0->createClone(), + frac1->createClone() + }; + auto bgResoModel = std::make_unique( nGauss, bgResoParams, scale ); + + 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 ) { + 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)); + 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; +} diff --git a/examples/ProgOpts/CMakeLists.txt b/examples/ProgOpts/CMakeLists.txt index 24de35c..30fa487 100644 --- a/examples/ProgOpts/CMakeLists.txt +++ b/examples/ProgOpts/CMakeLists.txt @@ -1,7 +1,7 @@ find_package(Boost REQUIRED COMPONENTS program_options) -add_library(ProgramOptions STATIC Test_Dpipi_ProgOpts.cc Test_Dpipi_ProgOpts.hh Command.hh) +add_library(ProgramOptions STATIC TimeDep_Dpipi_ProgOpts.cc TimeDep_Dpipi_ProgOpts.hh Command.hh) target_include_directories(ProgramOptions PUBLIC ${CMAKE_CURRENT_LIST_DIR}) target_link_libraries(ProgramOptions PUBLIC Laura++ Boost::boost Boost::program_options) diff --git a/examples/ProgOpts/Test_Dpipi_ProgOpts.cc b/examples/ProgOpts/TimeDep_Dpipi_ProgOpts.cc similarity index 99% rename from examples/ProgOpts/Test_Dpipi_ProgOpts.cc rename to examples/ProgOpts/TimeDep_Dpipi_ProgOpts.cc index a30c73c..9ac5d3f 100644 --- a/examples/ProgOpts/Test_Dpipi_ProgOpts.cc +++ b/examples/ProgOpts/TimeDep_Dpipi_ProgOpts.cc @@ -1,174 +1,174 @@ #include #include #include "boost/program_options.hpp" #include "boost/tokenizer.hpp" -#include "Test_Dpipi_ProgOpts.hh" +#include "TimeDep_Dpipi_ProgOpts.hh" namespace po = boost::program_options; void validate( boost::any& v, const std::vector& values, Command* /*target_type*/, int) { // Make sure no previous assignment to 'v' has been made po::validators::check_first_occurrence(v); // Extract the first string from 'values'. If there is more than // one string, it's an error, and exception will be thrown. const std::string& s { po::validators::get_single_string(values) }; if ( s == "gen" ) { v = boost::any( Command::Generate ); } else if ( s == "fit" ) { v = boost::any( Command::Fit ); } else if ( s == "simfit" ) { v = boost::any( Command::SimFit ); } else { throw po::validation_error(po::validation_error::invalid_option_value, "command", s, 3); } } void validate( boost::any& v, const std::vector& values, LauTimeDepFitModel::CPEigenvalue* /*target_type*/, int) { // Make sure no previous assignment to 'v' has been made po::validators::check_first_occurrence(v); // Extract the first string from 'values'. If there is more than // one string, it's an error, and exception will be thrown. const std::string& s { po::validators::get_single_string(values) }; if ( s == "QFS" ) { v = boost::any( LauTimeDepFitModel::CPEigenvalue::QFS ); } else if ( s == "CPEven" ) { v = boost::any( LauTimeDepFitModel::CPEigenvalue::CPEven ); } else if ( s == "CPOdd" ) { v = boost::any( LauTimeDepFitModel::CPEigenvalue::CPOdd ); } else { throw po::validation_error(po::validation_error::invalid_option_value); } } namespace LauDecayTime { void validate( boost::any& v, const std::vector& values, EfficiencyMethod* /*target_type*/, int) { // Make sure no previous assignment to 'v' has been made po::validators::check_first_occurrence(v); // Extract the first string from 'values'. If there is more than // one string, it's an error, and exception will be thrown. const std::string& s { po::validators::get_single_string(values) }; if ( s == "uniform" || s == "flat" ) { v = boost::any( EfficiencyMethod::Uniform ); } else if ( s == "binned" || s == "hist" ) { v = boost::any( EfficiencyMethod::Binned ); } else if ( s == "spline" ) { v = boost::any( EfficiencyMethod::Spline ); } else { throw po::validation_error(po::validation_error::invalid_option_value); } } }; TestDpipi_ProgramSettings::TestDpipi_ProgramSettings(const int argc, const char ** argv) { po::options_description main_desc{"Main options"}; main_desc.add_options() ("command", po::value(&command)->required(), "main command: gen, fit, or simfit") ; po::positional_options_description p; p.add("command", 1); po::options_description common_desc{"Common options"}; common_desc.add_options() ("help", "produce help message") ("dtype", po::value(&dType)->default_value(LauTimeDepFitModel::CPEigenvalue::QFS,"QFS"), "type of D decay: QFS, CPOdd, or CPEven") ("dta-model", po::value(&timeEffModel)->default_value(LauDecayTime::EfficiencyMethod::Uniform,"uniform"), "decay-time acceptance model: uniform/flat, binned/hist, spline") ("dtr", po::value(&timeResolution)->default_value(kTRUE), "enable/disable decay-time resolution") ("dtr-perevent", po::value(&perEventTimeErr)->default_value(kFALSE), "enable/disable use of per-event decay-time error (requires decay-time resolution to be enabled to take effect)") ("scaleYields", po::value(&scaleYields)->default_value(kTRUE), "toggle whether (kTRUE) the yields are scaled to account for the neglected components (the default strategy) or (kFALSE) the neglected components' yields are added to the combinatorial yield") ("setNegToZero", po::value(&setNegToZero)->default_value(kFALSE), "toggle whether (kTRUE) the negative yields from the mass fits for Run 1 KK are set to 0") ("constrainedYields", po::value(&constrainedYields)->default_value(kTRUE), "toggle whether (kTRUE) the yields in the mass fits for Run 1 KK are constrained based on the other samples") ("seed", po::value(&RNGseed)->default_value(0), "set the seed for the RNG; if not set, the time is used to generate a seed") ("run", po::value(&run)->default_value(2), "set the Run number (1 or 2)") ("dir", po::value(&directory)->default_value(""), "set the directory used to find the nTuples within the root file. Defaults to no directory") ("bkgndList", po::value(&bkgndListStr)->default_value("comb,Bd2DKpi,Bs2DKpi,Lb2Dppi"), "the comma-separated list of background components to include") ("fitBackInput", po::value(&fitBackInput)->default_value("fit0_ToyMC_QFS_expts0-0_expt0.root"), "set the file used as input to the fit: only used if fitBack is also selected") ("eosRoot", po::value(&eosRoot)->default_value("root://eoslhcb.cern.ch//eos/lhcb/wg/b2oc/B2D0pipi"), "set the base EOS area from which all the input data/histograms/json files/etc. should be read") ("eosConfigDir", po::value(&eosConfigDir)->default_value("timedep_fit_histos"), "set the directory name from which the configuration (json files / histogram files) should be read") ("eosDataDir", po::value(&eosDataDir)->default_value("Tuples22/Data"), "set the directory name from which the data should be read") ; po::options_description gen_desc{"Generation options"}; gen_desc.add_options() ("firstExptGen", po::value(&firstExptGen)->default_value(0), "first experiment to generate") ("nExptGen", po::value(&nExptGen)->default_value(1), "number of experiments to generate") ; po::options_description fit_desc{"Fitting options"}; fit_desc.add_options() ("firstExpt", po::value(&firstExptFit)->default_value(0), "first experiment to fit") ("nExpt", po::value(&nExptFit)->default_value(1), "number of experiments to fit") ("iFit", po::value(&iFit)->default_value(0), "fit ID - used to distinguish multiple fits to a given dataset (used to disentangle multiple minima)") ("fixTau", po::value(&fixLifetime)->default_value(kTRUE), "enable/disable floating of B lifetime parameter") ("fixDm", po::value(&fixDeltaM)->default_value(kTRUE), "enable/disable floating of B mixing frequency parameter") ("fixPhiMix", po::value(&fixPhiMix)->default_value(kFALSE), "enable/disable floating of B mixing phase parameter(s)") ("fixSplineKnots", po::value(&fixSplineKnots)->default_value(kFALSE), "enable/disable floating of the decay-time-acceptance spline") ("useSinCos", po::value(&useSinCos)->default_value(kTRUE), "enable/disable using sine and cosine of B mixing phase as the floating parameters rather than the phase itself") ("useAveDeltaCalibVals", po::value(&useAveDeltaCalibVals)->default_value(kTRUE), "enable/disable fitting calib values as their average +/- delta values rather than having separate values for B and Bbar") ("addTaggers", po::value(&addTaggers)->default_value(kTRUE), "add/omit taggers") ("floatCalibPars", po::value(&floatCalibPars)->default_value(kTRUE), "enable/disable floating of the FT calibration parameters") ("floatYields", po::value(&floatYields)->default_value(kFALSE), "enable/disable floating of the yields of each component - if enabled/disabled implies inclusion/exclusion of non-DP PDFs in/from the likelihood") ("fitBack", po::value(&fitBack)->default_value(kFALSE), "enable/disable refitting some generated data from a model") ("blindFit", po::value(&blindFit)->default_value(kFALSE), "enable/disable blinding of the physics parameters") ("dataFit", po::value(&dataFit)->default_value(kFALSE), "fit to data (not MC or toys)") ("genToy", po::value(&genToy)->default_value(kFALSE), "enable/disable generation of a toy to a fit result") ; po::options_description simfit_desc{"Simultaneous fitting options"}; simfit_desc.add_options() ("port", po::value(&port)->default_value(0), "port number on which sim fit coordinator is listening") ; po::options_description all_desc; all_desc.add(main_desc).add(common_desc).add(gen_desc).add(fit_desc).add(simfit_desc); po::variables_map vm; try { po::store(po::command_line_parser(argc, argv). options(all_desc). positional(p). run(), vm); po::notify(vm); } catch ( boost::wrapexcept& e ) { std::cout << argv[0] << " requires a command, one of 'gen', 'fit', or 'simfit'\n" << "Append --help to those commands to see help on the related options" << std::endl; parsedOK = kFALSE; return; } catch ( po::validation_error& e ) { std::cerr << e.what() << std::endl; parsedOK = kFALSE; return; } if ( vm.count("help") ) { po::options_description help_desc; help_desc.add(common_desc).add(gen_desc).add(fit_desc).add(simfit_desc); std::cout << help_desc << std::endl; helpRequested = kTRUE; return; } if ( ! bkgndListStr.empty() ) { const boost::char_separator sep{","}; const boost::tokenizer> tokens{bkgndListStr, sep}; for (const auto& t : tokens) { bkgndList.insert( t ); } } } diff --git a/examples/ProgOpts/Test_Dpipi_ProgOpts.hh b/examples/ProgOpts/TimeDep_Dpipi_ProgOpts.hh similarity index 100% rename from examples/ProgOpts/Test_Dpipi_ProgOpts.hh rename to examples/ProgOpts/TimeDep_Dpipi_ProgOpts.hh diff --git a/examples/Test_Dpipi.cc b/examples/Test_Dpipi.cc deleted file mode 100644 index 5a98ab0..0000000 --- a/examples/Test_Dpipi.cc +++ /dev/null @@ -1,448 +0,0 @@ - -#include - -#include -#include -#include - -#include "TFile.h" -#include "TH2.h" -#include "TRandom.h" -#include "TString.h" -#include "TSystem.h" -#include "TF1.h" -#include "TCanvas.h" - -#include "Lau1DHistPdf.hh" -#include "Lau1DCubicSpline.hh" -#include "LauBinnedDecayTimeEfficiency.hh" -#include "LauBkgndDPModel.hh" -#include "LauDaughters.hh" -#include "LauDecayTime.hh" -#include "LauDecayTimePdf.hh" -#include "LauDecayTimePhysicsModel.hh" -#include "LauDecayTimeResolution.hh" -#include "LauEffModel.hh" -#include "LauFlavTag.hh" -#include "LauIsobarDynamics.hh" -#include "LauMagPhaseCoeffSet.hh" -#include "LauRandom.hh" -#include "LauRealImagCoeffSet.hh" -#include "LauSplineDecayTimeEfficiency.hh" -#include "LauTimeDepFitModel.hh" -#include "LauVetoes.hh" - -#include "Test_Dpipi_ProgOpts.hh" - -int main(const int argc, const char ** argv) -{ - const TestDpipi_ProgramSettings settings{argc,argv}; - if ( settings.helpRequested ) { - return EXIT_SUCCESS; - } - if ( ! settings.parsedOK ) { - return EXIT_FAILURE; - } - - LauRandom::setSeed(settings.RNGseed); - - LauDaughters* daughtersB0bar = new LauDaughters("B0_bar", "pi+", "pi-", "D0"); - LauDaughters* daughtersB0 = new LauDaughters("B0", "pi+", "pi-", "D0_bar"); - - // efficiency - LauVetoes* vetoes = new LauVetoes(); - LauEffModel* effModelB0bar = new LauEffModel(daughtersB0bar, vetoes); - LauEffModel* effModelB0 = new LauEffModel(daughtersB0, vetoes); - - // background types - /* - std::map bkgndInfo { - {"Bkgnd1",LauFlavTag::BkgndType::Combinatorial}, - {"Bkgnd2",LauFlavTag::BkgndType::SignalLike} - }; - */ - - // 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"); - } - flavTag->setTrueTagVarName("trueTag"); - - TFile* etaFile = TFile::Open("ft-eta-hist.root"); - TH1* etaHist = dynamic_cast(etaFile->Get("ft_eta_hist")); - - // 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 - const Double_t meanEta { etaHist->GetMean() }; - const Bool_t interpolateEtaHist { meanEta > etaHist->GetBinWidth(1) }; - Lau1DHistPdf* etaHistPDF = new Lau1DHistPdf( "eta", etaHist, 0.0, 0.5, interpolateEtaHist, kFALSE ); - - const Double_t tagEffVal { (interpolateEtaHist) ? 0.5 : 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("OSTagger", "tagVal_OS", "mistagVal_OS", etaHistPDF, tagEff, calib0, calib1); - if (settings.floatCalibPars) { - flavTag->floatAllCalibPars(); - } - } - - // flavour tagging setup for backgrounds - /* - std::map BkgndEtas; - TFile* bkgndEtaFile = TFile::Open("ft-bkgnd-eta-hists.root"); - for(auto& [name, _] : BkgndTypes) - { - TH1* bkgndEtaHist = dynamic_cast(bkgndEtaFile->Get( name+"_eta_hist" )); - Lau1DHistPdf* bkgndEtaHistPDF = new Lau1DHistPdf("eta",bkgndEtaHist,0.0,0.5,kTRUE,kFALSE); - BkgndEtas.emplace( std::make_pair(name, bkgndEtaHistPDF) ); - } - for(auto& [name,hist] : BkgndEtas) - {flavTag->setBkgndParams(name, "IFT", hist, tagEff );} - */ - - // signal dynamics - LauIsobarDynamics* sigModelB0bar = new LauIsobarDynamics(daughtersB0bar, effModelB0bar); - sigModelB0bar->setIntFileName("integ_B0bar.dat"); - sigModelB0bar->addResonance("D*+_2", 2, LauAbsResonance::RelBW); - sigModelB0bar->addResonance("D*+_0", 2, LauAbsResonance::RelBW); - sigModelB0bar->addResonance("rho0(770)", 3, LauAbsResonance::RelBW); - sigModelB0bar->addResonance("f_0(980)", 3, LauAbsResonance::RelBW); - sigModelB0bar->addResonance("f_2(1270)", 3, LauAbsResonance::RelBW); - - LauIsobarDynamics* sigModelB0 = new LauIsobarDynamics(daughtersB0, effModelB0); - sigModelB0->setIntFileName("integ_B0.dat"); - sigModelB0->addResonance("D*-_2", 1, LauAbsResonance::RelBW); - sigModelB0->addResonance("D*-_0", 1, LauAbsResonance::RelBW); - sigModelB0->addResonance("rho0(770)", 3, LauAbsResonance::RelBW); - sigModelB0->addResonance("f_0(980)", 3, LauAbsResonance::RelBW); - sigModelB0->addResonance("f_2(1270)", 3, LauAbsResonance::RelBW); - - // fit model - LauTimeDepFitModel* fitModel = new LauTimeDepFitModel(sigModelB0bar,sigModelB0,flavTag); - - std::vector coeffset; - coeffset.push_back( new LauRealImagCoeffSet("D*+_2", 1.00, 0.00, kTRUE, kTRUE) ); - coeffset.push_back( new LauRealImagCoeffSet("D*+_0", 0.53*TMath::Cos( 3.00), 0.53*TMath::Sin( 3.00), kFALSE, kFALSE) ); - coeffset.push_back( new LauRealImagCoeffSet("rho0(770)", 1.22*TMath::Cos( 2.25), 1.22*TMath::Sin( 2.25), kFALSE, kFALSE) ); - coeffset.push_back( new LauRealImagCoeffSet("f_0(980)", 0.19*TMath::Cos(-2.48), 0.19*TMath::Sin(-2.48), kFALSE, kFALSE) ); - coeffset.push_back( new LauRealImagCoeffSet("f_2(1270)", 0.75*TMath::Cos( 2.97), 0.75*TMath::Sin( 2.97), kFALSE, kFALSE) ); - - for (std::vector::iterator iter=coeffset.begin(); iter!=coeffset.end(); ++iter) { - fitModel->setAmpCoeffSet(*iter); - } - - // background DP models - /* - LauBkgndDPModel* Bkgnd1Model = new LauBkgndDPModel(daughtersB0, vetoes); - LauBkgndDPModel* Bkgnd1Modelbar = new LauBkgndDPModel(daughtersB0bar, vetoes); - LauBkgndDPModel* Bkgnd2Model = new LauBkgndDPModel(daughtersB0, vetoes); - LauBkgndDPModel* Bkgnd2Modelbar = new LauBkgndDPModel(daughtersB0bar, vetoes); - fitModel->setBkgndDPModels( "Bkgnd1", Bkgnd1Model, Bkgnd1Modelbar ); - fitModel->setBkgndDPModels( "Bkgnd2", Bkgnd2Model, Bkgnd2Modelbar ); - */ - - // 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 ); - - // 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.0); - const Double_t maxDt(15.0); - const Double_t minDtErr(0.0); - const Double_t maxDtErr(0.15); - - LauParameter * tau = new LauParameter("dt_tau", 1.519, 0.5, 5.0, settings.fixLifetime); - LauParameter * freq = new LauParameter("dt_deltaM", 0.5064, 0.0, 1.0, settings.fixDeltaM); - - std::vector dtPhysPars { - tau, - freq - }; - - auto dtPhysModel = std::make_unique( LauDecayTime::FuncType::ExpTrig, dtPhysPars ); - - const std::vector scale { - settings.perEventTimeErr && true, - settings.perEventTimeErr && true, - }; - const std::size_t nGauss{scale.size()}; - - LauParameter * mean0 = new LauParameter("dt_mean_0", scale[0] ? -2.01290e-03 : -2.25084e-03, -0.01, 0.01, kTRUE ); - LauParameter * mean1 = new LauParameter("dt_mean_1", scale[1] ? -2.01290e-03 : -5.04275e-03, -0.01, 0.01, kTRUE ); - LauParameter * sigma0 = new LauParameter("dt_sigma_0", scale[0] ? 9.95145e-01 : 3.03923e-02, 0.0, 2.0, kTRUE ); - LauParameter * sigma1 = new LauParameter("dt_sigma_1", scale[1] ? 1.81715e+00 : 6.22376e-02, 0.0, 2.5, kTRUE ); - LauParameter * frac1 = new LauParameter("dt_frac_1", scale[0] && scale[1] ? 1.-9.35940e-01 : 1.-7.69603e-01, 0.0, 1.0, kTRUE); - - std::vector dtResoPars { - mean0, - mean1, - sigma0, - sigma1, - frac1 - }; - - auto dtResoModel = std::make_unique( nGauss, dtResoPars, scale ); - - // 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("dte-hist.root"); - TH1* dtErrHist = dynamic_cast(dtErrFile->Get("dte_hist")); - - LauDecayTimePdf * dtPdf{nullptr}; - if ( settings.timeResolution ) { - if ( settings.perEventTimeErr ) { - 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("dta-hist.root"); - TH1* dtaHist = dynamic_cast(dtaFile->Get("dta_hist")); - - // 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}; - - switch(settings.timeEffModel) - { - case LauDecayTime::EfficiencyMethod::Spline: - { - fitModel->setASqMaxValue(0.06); - - auto dtEffSpline = std::make_unique( dtvals, effvals, Lau1DCubicSpline::SplineType::AkimaSpline ); - - auto dtaModel = std::make_unique>( "dteff_QFS", std::move(dtEffSpline) ); - dtaModel->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(); - } else { - dtaModel->floatKnots(); - dtaModel->fixKnot( 0, true ); - dtaModel->fixKnot( 3, true ); - } - - if ( settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS ) { - // For the QFS mode we just use the cubic model as it is - dtPdf->setSplineEfficiency( std::move(dtaModel) ); - } else { - // For the CP modes we modify it using a corrective spline - auto dtCorrSpline = std::make_unique( dtvals, corvals, Lau1DCubicSpline::SplineType::AkimaSpline ); - auto dtaCPModel = std::make_unique>( "dteff_CP", *dtaModel, std::move( dtCorrSpline ) ); - dtPdf->setSplineEfficiency( std::move(dtaCPModel) ); - } - - 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 ); - - //Background dt part - /* - TFile* background_dt = TFile::Open("Lifetimes_PV_WL.root"); - TH1* bkgnd1DtHist = dynamic_cast( background_dt->Get("Bkgnd1") ); - TH1* bkgnd2DtHist = dynamic_cast( background_dt->Get("Bkgnd2") ); - TH1* bkgnd1DtErrHist = dynamic_cast( background_dt->Get("Bkgnd1_Err") ); - TH1* bkgnd2DtErrHist = dynamic_cast( background_dt->Get("Bkgnd2_Err") ); - LauDecayTimePdf* bkgnd1DtPdf{nullptr}; - LauDecayTimePdf* bkgnd2DtPdf{nullptr}; - if ( settings.timeResolution and settings.perEventTimeErr ) { - bkgnd1DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd1DtHist, bkgnd1DtErrHist ); - bkgnd2DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd2DtHist, bkgnd2DtErrHist ); - } else { - bkgnd1DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd1DtHist ); - bkgnd2DtPdf= new LauDecayTimePdf( "decayTime", minDt, maxDt, bkgnd2DtHist ); - } - fitModel->setBkgndDtPdf("Bkgnd1",bkgnd1DtPdf); - fitModel->setBkgndDtPdf("Bkgnd2",bkgnd2DtPdf); - */ - - // set the signal yield - Double_t nSigEvents{0}; - switch (settings.dType) { - case LauTimeDepFitModel::CPEigenvalue::CPEven : - nSigEvents = 15000; - break; - case LauTimeDepFitModel::CPEigenvalue::CPOdd : - nSigEvents = 5000; - break; - case LauTimeDepFitModel::CPEigenvalue::QFS : - nSigEvents = 50000; - break; - } - - // set the background yields - /* - const Double_t nBkgnd1(100), nBkgnd2(200); - LauParameter* Bkgnd1Yield = new LauParameter("Bkgnd1",nBkgnd1,-1.0*nBkgnd1,2.0*nBkgnd1,kFALSE); - LauParameter* Bkgnd2Yield = new LauParameter("Bkgnd2",nBkgnd2,-1.0*nBkgnd2,2.0*nBkgnd2,kFALSE); - fitModel->setNBkgndEvents(Bkgnd1Yield); - fitModel->setNBkgndEvents(Bkgnd2Yield); - */ - - std::cout<<"nSigEvents = "<setNSigEvents(nSigPar); - - // 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; - - 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 += ".root"; - - if (settings.command == Command::Generate) { - rootFileName = "dummy.root"; - tableFileName = "genResults"; - } else { - rootFileName = "fit"; rootFileName += settings.iFit; - rootFileName += "_Results_"; rootFileName += dTypeStr; - rootFileName += "_expts"; rootFileName += settings.firstExptFit; rootFileName += "-"; rootFileName += settings.firstExptFit+settings.nExptFit-1; - rootFileName += ".root"; - - tableFileName = "fit"; tableFileName += settings.iFit; - tableFileName += "_Results_"; tableFileName += dTypeStr; - tableFileName += "_expts"; tableFileName += settings.firstExptFit; tableFileName += "-"; tableFileName += settings.firstExptFit+settings.nExptFit-1; - - fitToyFileName = "fit"; fitToyFileName += settings.iFit; - fitToyFileName += "_ToyMC_"; fitToyFileName += dTypeStr; - fitToyFileName += "_expts"; fitToyFileName += settings.firstExptFit; fitToyFileName += "-"; fitToyFileName += settings.firstExptFit+settings.nExptFit-1; - fitToyFileName += ".root"; - - splotFileName = "fit"; splotFileName += settings.iFit; - splotFileName += "_sPlot_"; splotFileName += dTypeStr; - splotFileName += "_expts"; splotFileName += settings.firstExptFit; splotFileName += "-"; splotFileName += settings.firstExptFit+settings.nExptFit-1; - splotFileName += ".root"; - } - - // 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; -} diff --git a/examples/runTimeDepTest.sh b/examples/runTimeDepTest.sh index 508c0ee..99d2fdd 100755 --- a/examples/runTimeDepTest.sh +++ b/examples/runTimeDepTest.sh @@ -1,242 +1,242 @@ #!/bin/bash #SBATCH --partition=epp #SBATCH --ntasks=1 #SBATCH --cpus-per-task=6 #SBATCH --mem-per-cpu=3997 #SBATCH --time=02:00:00 #Run with sbatch runTimeDepTest.sh # Copyright 2020 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 processes=() # trap ^C and call ctrl_c() trap ctrl_c INT function ctrl_c() { echo " ** Caught ^C event!" for process in ${processes[@]} do if ps -u $USER | grep -q "$process" then echo "Killing : $process" kill $process fi done exit 1 } # Check that we have a kerberos token, since nothing will work without one! klist -s > /dev/null 2>&1 if [ $? -ne 0 ] then echo "You don't have a valid Kerberos ticket!" exit 1 fi iFit=0 dta_model="spline" dtr=1 dtr_perevent=0 fixTau=1 fixDeltaM=1 seed=140279 # $(pwd | xargs basename) dataFit=1 runNo=12 # Modify these if you want to run more experiments nExpt=1 firstExpt=0 logName_genQFS_run1=Logs/gen-QFS-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log logName_genQFS_run2=Logs/gen-QFS-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log logName_genCPE_KK_run1=Logs/gen-CPEvenKK-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log logName_genCPE_KK_run2=Logs/gen-CPEvenKK-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log logName_genCPE_pipi_run1=Logs/gen-CPEvenpipi-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log logName_genCPE_pipi_run2=Logs/gen-CPEvenpipi-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent.log logName_coord=Logs/coordinator-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log logName_taskQFS_run1=Logs/task-QFS-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log logName_taskQFS_run2=Logs/task-QFS-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log logName_taskCPE_KK_run1=Logs/task-CPEvenKK-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log logName_taskCPE_KK_run2=Logs/task-CPEvenKK-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log logName_taskCPE_pipi_run1=Logs/task-CPEvenpipi-run1-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log logName_taskCPE_pipi_run2=Logs/task-CPEvenpipi-run2-dtamodel_$dta_model-dtr_$dtr-dtrperevent_$dtr_perevent-$iFit.log mkdir -p Logs # make sure we have the most up-to-date json files eosArea="root://eoslhcb.cern.ch//eos/lhcb/wg/b2oc/B2D0pipi/timedep_fit_histos" rm -f run1DTAsplines.json run2DTAsplines.json Bd2D0pipi_DP_Model_Coeffs.json xrdcp $eosArea/run1DTAsplines.json . if [ $? -ne 0 ] then echo "Problem downloading json files from EOS, exiting" exit 1 fi xrdcp $eosArea/run2DTAsplines.json . if [ $? -ne 0 ] then echo "Problem downloading json files from EOS, exiting" exit 1 fi xrdcp $eosArea/Bd2D0pipi_DP_Model_Coeffs.json . if [ $? -ne 0 ] then echo "Problem downloading json files from EOS, exiting" exit 1 fi if [ $iFit == 0 -a $dataFit == 0 ] then #run generation of the modes in parallel (put the tasks in the bg) echo "Generating samples..." if [ $runNo -eq 1 -o $runNo -eq 12 ] then - Test_Dpipi gen --dtype QFS --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_Kpi > $logName_genQFS_run1 2>&1 & + GenFitTimeDep_Bd2Dpipi gen --dtype QFS --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_Kpi > $logName_genQFS_run1 2>&1 & task=$! processes+=( $task ) - Test_Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_KK > $logName_genCPE_KK_run1 2>&1 & + GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_KK > $logName_genCPE_KK_run1 2>&1 & task=$! processes+=( $task ) - Test_Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_pipi > $logName_genCPE_pipi_run1 2>&1 & + GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_pipi > $logName_genCPE_pipi_run1 2>&1 & task=$! processes+=( $task ) fi if [ $runNo -eq 2 -o $runNo -eq 12 ] then - Test_Dpipi gen --dtype QFS --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_Kpi > $logName_genQFS_run2 2>&1 & + GenFitTimeDep_Bd2Dpipi gen --dtype QFS --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_Kpi > $logName_genQFS_run2 2>&1 & task=$! processes+=( $task ) - Test_Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_KK > $logName_genCPE_KK_run2 2>&1 & + GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_KK > $logName_genCPE_KK_run2 2>&1 & task=$! processes+=( $task ) - Test_Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_pipi > $logName_genCPE_pipi_run2 2>&1 & + GenFitTimeDep_Bd2Dpipi gen --dtype CPEven --sigOnly $sigOnly --fixTau $fixTau --seed $seed --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_pipi > $logName_genCPE_pipi_run2 2>&1 & task=$! processes+=( $task ) fi wait #wait for the generation to complete fi echo "Running fit $iFit" if [ $runNo -eq 12 ] then nTasks=6 else nTasks=3 fi SimFitCoordinator $iFit $nExpt $firstExpt $nTasks > $logName_coord 2>&1 & task=$! echo "SimFitCoordinator process : $task" processes+=( $task ) echo "Initialised coordinator process for fit $iFit" sleep 5 port="" while [ "x$port" == "x" ] do sleep 2 port=`tail $logName_coord | grep "Waiting for connection" | awk '{print $NF}'` done echo "Coordinator is listening on port $port" echo "Initialising tasks..." if [ $runNo -eq 1 -o $runNo -eq 12 ] then - Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype QFS --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_Kpi --fixDm $fixDeltaM > $logName_taskQFS_run1 2>&1 & + GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype QFS --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_Kpi --fixDm $fixDeltaM > $logName_taskQFS_run1 2>&1 & task=$! echo "QFS Run 1 process : $task" processes+=( $task ) sleep 5 - Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_KK --fixDm $fixDeltaM > $logName_taskCPE_KK_run1 2>&1 & + GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_KK --fixDm $fixDeltaM > $logName_taskCPE_KK_run1 2>&1 & task=$! echo "CP KK Run 1 process : $task" processes+=( $task ) sleep 5 - Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_pipi --fixDm $fixDeltaM > $logName_taskCPE_pipi_run1 2>&1 & + GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 1 --dir B2Dhh_pipi --fixDm $fixDeltaM > $logName_taskCPE_pipi_run1 2>&1 & task=$! echo "CP pipi Run 1 process : $task" processes+=( $task ) fi if [ $runNo -eq 2 -o $runNo -eq 12 ] then - Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype QFS --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_Kpi --fixDm $fixDeltaM > $logName_taskQFS_run2 2>&1 & + GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype QFS --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_Kpi --fixDm $fixDeltaM > $logName_taskQFS_run2 2>&1 & task=$! echo "QFS Run 2 process : $task" processes+=( $task ) sleep 5 - Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_KK --fixDm $fixDeltaM > $logName_taskCPE_KK_run2 2>&1 & + GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_KK --fixDm $fixDeltaM > $logName_taskCPE_KK_run2 2>&1 & task=$! echo "CP KK Run 2 process : $task" processes+=( $task ) sleep 5 - Test_Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_pipi --fixDm $fixDeltaM > $logName_taskCPE_pipi_run2 2>&1 & + GenFitTimeDep_Bd2Dpipi simfit --blindFit $dataFit --dataFit $dataFit --dtype CPEven --seed $seed --port $port --iFit $iFit --fixTau $fixTau --dta-model $dta_model --dtr $dtr --dtr-perevent $dtr_perevent --firstExpt $firstExpt --nExpt $nExpt --firstExptGen $firstExpt --nExptGen $nExpt --run 2 --dir B2Dhh_pipi --fixDm $fixDeltaM > $logName_taskCPE_pipi_run2 2>&1 & task=$! echo "CP pipi Run 2 process : $task" processes+=( $task ) fi wait #wait for the fits and coordinator to terminate echo "Fit completed" echo "Warnings/Errors from coordinator:" grep -e ERROR -e WARNING -e Error -e Warning $logName_coord echo "Warnings/Errors from CPeven KK tasks:" if [ $runNo -eq 1 -o $runNo -eq 12 ] then grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_KK_run1 fi if [ $runNo -eq 2 -o $runNo -eq 12 ] then grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_KK_run2 fi echo "Warnings/Errors from CPeven pipi tasks:" if [ $runNo -eq 1 -o $runNo -eq 12 ] then grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_pipi_run1 fi if [ $runNo -eq 2 -o $runNo -eq 12 ] then grep -e ERROR -e WARNING -e Error -e Warning $logName_taskCPE_pipi_run2 fi echo "Warnings/Errors from QFS tasks:" if [ $runNo -eq 1 -o $runNo -eq 12 ] then grep -e ERROR -e WARNING -e Error -e Warning $logName_taskQFS_run1 fi if [ $runNo -eq 2 -o $runNo -eq 12 ] then grep -e ERROR -e WARNING -e Error -e Warning $logName_taskQFS_run2 fi