diff --git a/examples/Test_Dpipi.cc b/examples/Test_Dpipi.cc index 4daf289..0e76fa0 100644 --- a/examples/Test_Dpipi.cc +++ b/examples/Test_Dpipi.cc @@ -1,435 +1,444 @@ #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}; flavTag->addTagger("OSTagger", "tagVal_OS", "mistagVal_OS", etaHistPDF, tagEff, calib0, calib1); 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.0, 0.01, 0.022, 0.035, 0.042, 0.05, 0.051, 0.052, 0.055}; + 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::AkimaSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); dtEffSpline->fitToTH1(dtaHist); - auto dtaModel = std::make_unique>( std::move(dtEffSpline) ); + auto dtaModel = std::make_unique>( "dteff_QFS", std::move(dtEffSpline) ); // 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 ); } - dtPdf->setSplineEfficiency( std::move(dtaModel) ); + 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::AkimaSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); + 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 dtaModel = std::make_unique( *dtaHist ); - dtPdf->setBinnedEfficiency( std::move(dtaModel) ); + 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; }