diff --git a/examples/Test_Dpipi.cc b/examples/Test_Dpipi.cc
index 5704276..d7a10c0 100644
--- a/examples/Test_Dpipi.cc
+++ b/examples/Test_Dpipi.cc
@@ -1,329 +1,330 @@
#include <iostream>
#include <vector>
#include <map>
#include <algorithm>
#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<TH1*>(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<Double_t,Double_t> tagEff {tagEffVal, tagEffVal};
// use a null calibration for the time being, so p0 = <eta> and p1 = 1
std::pair<Double_t,Double_t> calib0 {meanEta, meanEta};
std::pair<Double_t,Double_t> 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<LauAbsCoeffSet*> 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<LauAbsCoeffSet*>::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 Double_t maxDtErr(0.15);
const std::vector<Bool_t> 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 * 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);
+ 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<LauAbsRValue*> dtPars {
mean0,
mean1,
sigma0,
sigma1,
frac1,
tau,
freq
};
// Decay time error histogram
TFile* dteFile = TFile::Open("dte-hist.root");
TH1* dteHist = dynamic_cast<TH1*>(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<TH1*>(dtaFile->Get("dta_hist"));
// Create the spline knot positions and y-values
std::vector<Double_t> dtvals = {0.,0.5 ,1., 2., 3., 4., 7., 10., 15.};
// Starting Y values, to be fit to dtaHist
std::vector<Double_t> 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::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);
+ dtPdf -> getEffiPars().at(0) -> 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 = "<<nSigEvents<<std::endl;
LauParameter* nSigPar = new LauParameter("signalEvents", nSigEvents, -2.0*nSigEvents, 2.0*nSigEvents, kTRUE);
fitModel->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;
}