diff --git a/examples/Test_Dpipi.cc b/examples/Test_Dpipi.cc index 128a88b..5704276 100644 --- a/examples/Test_Dpipi.cc +++ b/examples/Test_Dpipi.cc @@ -1,329 +1,329 @@ #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 "LauDaughters.hh" #include "LauDecayTimePdf.hh" #include "LauEffModel.hh" #include "LauIsobarDynamics.hh" #include "LauMagPhaseCoeffSet.hh" #include "LauRandom.hh" #include "LauRealImagCoeffSet.hh" #include "LauTimeDepFitModel.hh" #include "LauVetoes.hh" #include "LauFlavTag.hh" #include "Lau1DHistPdf.hh" #include "Lau1DCubicSpline.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; } const Bool_t fixPhiMix{ settings.fixPhiMix || settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS }; const Bool_t useSinCos{kTRUE}; 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; } 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(); //vetoes->addMassVeto( 2, 2.00776, 2.01276 ); LauEffModel* effModelB0bar = new LauEffModel(daughtersB0bar, vetoes); LauEffModel* effModelB0 = new LauEffModel(daughtersB0, vetoes); //Args for flavTag: useAveDelta - kFALSE and useEtaPrime - kFALSE LauFlavTag* flavTag = new LauFlavTag(kFALSE,kFALSE); flavTag->setTrueTagVarName("trueTag"); if (settings.dType == LauTimeDepFitModel::CPEigenvalue::QFS) { flavTag->setDecayFlvVarName("decayFlv"); } TFile* etaFile = TFile::Open("ft-eta-hist.root"); TH1* etaHist = dynamic_cast(etaFile->Get("ft_eta_hist")); Lau1DHistPdf* etaHistPDF = new Lau1DHistPdf("eta",etaHist,0.0,0.5,kTRUE,kFALSE); const Double_t meanEta { etaHistPDF->getMean() }; // if the tagging is perfect then also make it perfectly efficient, otherwise 50% efficient const Double_t tagEffVal { (meanEta == 0.0) ? 1.0 : 0.5 }; std::pair tagEff {tagEffVal, tagEffVal}; // use a null calibration for the time being, so p0 = and p1 = 1 std::pair calib0 {meanEta, meanEta}; std::pair calib1 {1.0, 1.0}; flavTag->addTagger("OSTagger", "tagVal_OS", "mistagVal_OS", etaHistPDF, tagEff, calib0, calib1); // 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); } fitModel->setCPEigenvalue( settings.dType ); fitModel->setPhiMix( 2.0*LauConstants::beta, fixPhiMix, useSinCos ); // production asymmetry fitModel->setAsymmetries(0.0,kTRUE); // Delta t PDFs const Double_t minDt(0.0); const Double_t maxDt(15.0); const Double_t minDtErr(0.0); const Double_t maxDtErr(0.215); const std::vector scale { settings.perEventTimeErr && kTRUE, settings.perEventTimeErr && kTRUE, }; const UInt_t nGauss(scale.size()); LauParameter * mean0 = new LauParameter("dt_mean_0", scale[0] ? -1.63e-3 : -1.84e-03, -0.01, 0.01, kTRUE ); LauParameter * mean1 = new LauParameter("dt_mean_1", scale[1] ? -1.63e-3 : -3.62e-03, -0.01, 0.01, kTRUE ); LauParameter * sigma0 = new LauParameter("dt_sigma_0", scale[0] ? 0.991 : 3.05e-02, 0.0, 2.0, kTRUE ); LauParameter * sigma1 = new LauParameter("dt_sigma_1", scale[1] ? 1.80 : 6.22e-02, 0.0, 2.5, kTRUE ); LauParameter * frac1 = new LauParameter("dt_frac_1", scale[0] && scale[1] ? 0.065 : 0.761, 0.0, 1.0, kTRUE); LauParameter * tau = new LauParameter("dt_tau", 1.520, 0.5, 5.0, settings.fixLifetime); LauParameter * freq = new LauParameter("dt_deltaM", 0.5064, 0.0, 1.0, settings.fixDeltaM); std::vector dtPars { mean0, mean1, sigma0, sigma1, frac1, tau, freq }; - - // 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 y-values - std::vector dtvals = {0.,1., 2., 4., 7., 10., 15.}; - // Starting Y values, to be fit to dtaHist - std::vector effvals = {0.,0.022,0.035,0.05,0.051,0.052,0.055}; - + // Decay time error histogram TFile* dteFile = TFile::Open("dte-hist.root"); TH1* dteHist = dynamic_cast(dteFile->Get("dte_hist")); LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "decayTime", "decayTimeErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::ExpTrig, nGauss, scale, LauDecayTimePdf::DecayTime, settings.timeEffModel ); dtPdf->doSmearing(settings.timeResolution); if ( settings.perEventTimeErr ) { dtPdf->setErrorHisto( dteHist ); } + // 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 y-values + std::vector dtvals = {0.,0.5 ,1., 2., 3., 4., 7., 10., 15.}; + // Starting Y values, to be fit to dtaHist + std::vector effvals = {0.,0.01 ,0.022,0.035,0.042,0.05,0.051,0.052,0.055}; + Lau1DCubicSpline* dtEffSpline = nullptr; switch(settings.timeEffModel) { case LauDecayTimePdf::EfficiencyMethod::Spline: { fitModel->setASqMaxValue(0.06); - dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,Lau1DCubicSpline::StandardSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); + dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,Lau1DCubicSpline::AkimaSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); dtEffSpline->fitToTH1(dtaHist); dtPdf->setEffiSpline(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 both DeltaM and the B lifetime are fixed! for(auto& par : dtPdf -> getEffiPars()) {par -> fixed( not (settings.fixDeltaM and settings.fixLifetime) );} dtPdf -> getEffiPars().at(3) -> fixed(kTRUE); break; } case LauDecayTimePdf::EfficiencyMethod::Binned: { fitModel->setASqMaxValue(0.06); dtPdf->setEffiHist(dtaHist); break; } case LauDecayTimePdf::EfficiencyMethod::Flat: { fitModel->setASqMaxValue(4.1); break; } } fitModel->setSignalDtPdf( dtPdf ); // set the number of signal events 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->doPoissonSmearing(kFALSE); fitModel->doEMLFit(kFALSE); fitModel->writeLatexTable(kFALSE); 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 LauDecayTimePdf::EfficiencyMethod::Spline: dataFile += "_Spline"; break; case LauDecayTimePdf::EfficiencyMethod::Binned: dataFile += "_Hist"; break; case LauDecayTimePdf::EfficiencyMethod::Flat: dataFile += "_Flat"; 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/test/TestFitSplineToTH1.cc b/test/TestFitSplineToTH1.cc index 83528bb..732bb7b 100644 --- a/test/TestFitSplineToTH1.cc +++ b/test/TestFitSplineToTH1.cc @@ -1,35 +1,35 @@ #include #include "TFile.h" #include "TH1.h" #include "TF1.h" #include "TCanvas.h" #include "Lau1DCubicSpline.hh" //int main(const int argc, const char ** argv) int main() { TFile* dtaFile = TFile::Open("dta-hist.root"); TH1* dtaHist = dynamic_cast(dtaFile->Get("dta_hist")); - std::vector dtvals = {0.,1., 2., 4., 7., 10., 15.}; - std::vector effvals = {0.,0.022,0.035,0.05,0.051,0.052,0.055}; + std::vector dtvals = {0.,0.5 ,1., 2., 3., 4., 7., 10., 15.}; + std::vector effvals = {0.,0.01 ,0.022,0.035,0.042,0.05,0.051,0.052,0.055}; Lau1DCubicSpline* dtEffSpline = nullptr; dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,Lau1DCubicSpline::AkimaSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); dtEffSpline->fitToTH1(dtaHist); TF1* tf1 = dtEffSpline->makeTF1(); TCanvas canvas("canvas","canvas",800,450); dtaHist->Draw(); tf1->Draw("SAME"); canvas.SaveAs("TF1fitTest.pdf"); delete dtEffSpline; delete dtaHist; dtaFile->Close(); delete dtaFile; return 0; }