diff --git a/examples/Test_Dpipi.cc b/examples/Test_Dpipi.cc index 0d9efc0..086c509 100644 --- a/examples/Test_Dpipi.cc +++ b/examples/Test_Dpipi.cc @@ -1,346 +1,336 @@ #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.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; } 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"); TFile* etafile = TFile::Open("histogram.root"); TH1* etahist = dynamic_cast(etafile->Get("htemp")); Lau1DHistPdf* etahistpdf = new Lau1DHistPdf("eta",etahist,0.0,0.54,kTRUE,kFALSE); std::pair tagEff {0.5, 0.5}; std::pair calib0 {0.4, 0.4}; 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(20.0); const Double_t minDtErr(0.0); const Double_t maxDtErr(2.5); const std::vector scale { settings.perEventTimeErr && kTRUE, settings.perEventTimeErr && kTRUE, - settings.perEventTimeErr && kFALSE }; const UInt_t nGauss(scale.size()); - TString mean0Name("dt_"); mean0Name += "_mean_0"; - TString mean1Name("dt_"); mean1Name += "_mean_1"; - TString mean2Name("dt_"); mean2Name += "_mean_2"; - TString sigma0Name("dt_"); sigma0Name += "_sigma_0"; - TString sigma1Name("dt_"); sigma1Name += "_sigma_1"; - TString sigma2Name("dt_"); sigma2Name += "_sigma_2"; - TString frac1Name("dt_"); frac1Name += "_frac_1"; - TString frac2Name("dt_"); frac2Name += "_frac_2"; - TString tauName("dt_"); tauName += "_tau"; - TString freqName("dt_"); freqName += "_deltaM"; - - LauParameter * mean0 = new LauParameter(mean0Name, 0.0); - LauParameter * mean1 = new LauParameter(mean1Name, 0.0); - LauParameter * mean2 = new LauParameter(mean2Name, 0.0); - LauParameter * sigma0 = new LauParameter(sigma0Name, 1.0); - LauParameter * sigma1 = new LauParameter(sigma1Name, 2.0); - LauParameter * sigma2 = new LauParameter(sigma2Name, 8.0); - LauParameter * frac1 = new LauParameter(frac1Name, 0.095); - LauParameter * frac2 = new LauParameter(frac2Name, 0.005); - LauParameter * tau = new LauParameter(tauName, 1.520); - LauParameter * freq = new LauParameter(freqName, 0.5064); + TString mean0Name("dt_mean_0"); + TString mean1Name("dt_mean_1"); + TString sigma0Name("dt_sigma_0"); + TString sigma1Name("dt_sigma_1"); + TString frac1Name("dt_frac_1"); + TString tauName("dt_tau"); + TString freqName("dt_deltaM"); + + LauParameter * mean0 = new LauParameter(mean0Name, scale[0] ? -1.63e-3 : -1.84e-03, -0.01, 0.01, kTRUE ); + LauParameter * mean1 = new LauParameter(mean1Name, scale[1] ? -1.63e-3 : -3.62e-03, -0.01, 0.01, kTRUE ); + LauParameter * sigma0 = new LauParameter(sigma0Name,scale[0] ? 0.991 : 3.05e-02, 0.5, 2.5, kTRUE ); + LauParameter * sigma1 = new LauParameter(sigma1Name,scale[1] ? 1.80 : 6.22e-02, 0.5, 2.5, kTRUE ); + LauParameter * frac1 = new LauParameter(frac1Name, scale[0] && scale[1] ? 0.065 : 0.761 ,0.,1.,kTRUE); + LauParameter * tau = new LauParameter(tauName, 1.520,1.3,20.,kTRUE); + LauParameter * freq = new LauParameter(freqName, 0.5064, 0.3, 0.7,kTRUE); std::vector dtPars { mean0, mean1, - mean2, sigma0, sigma1, - sigma2, frac1, - frac2, tau, freq }; //Decay time acceptance spline knot x- and y-values std::vector dtvals; dtvals.push_back(0.0); dtvals.push_back(1.0); dtvals.push_back(2.0); dtvals.push_back(4.0); dtvals.push_back(6.0); dtvals.push_back(8.0); dtvals.push_back(11.0); dtvals.push_back(14.0); dtvals.push_back(17.0); dtvals.push_back(20.0); std::vector effvals; effvals.push_back(0.0); effvals.push_back(0.15); effvals.push_back(0.25); effvals.push_back(0.33); effvals.push_back(0.38); effvals.push_back(0.4); effvals.push_back(0.43); effvals.push_back(0.45); effvals.push_back(0.47); effvals.push_back(0.50); LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "decayTime", "decayTimeErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::ExpTrig, nGauss, scale, LauDecayTimePdf::DecayTime, settings.timeEffModel ); dtPdf->doSmearing(settings.timeResolution); switch(settings.timeEffModel) { case LauDecayTimePdf::EfficiencyMethod::Spline: { fitModel->setASqMaxValue(0.34); Lau1DCubicSpline* dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,Lau1DCubicSpline::AkimaSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); dtPdf->setEffiSpline(dtEffSpline); break; } case LauDecayTimePdf::EfficiencyMethod::Binned: { fitModel->setASqMaxValue(0.34); const UInt_t upscale = 5; //This value must be >= 1 this linearly interpolates a histogram between the points of the spline const Int_t nBins = upscale * (dtvals.size()-1) ; Double_t edges[nBins + 1]; for(Int_t i = 0; i < nBins + 1; ++i) { Double_t binWidth = dtvals[i/upscale + 1] - dtvals[i/upscale]; //This is the width of the bin if upscale is 1, it is the distance between the current 2 knots edges[i] = dtvals[i/upscale] + (i%upscale)*binWidth/upscale; } std::cout << "Bins: "; for (int i = 0; i < nBins + 1; ++i){std::cout << edges[i] << ", ";} std::cout << std::endl; Double_t binFilling[nBins]; TH1D* dtEffHist = new TH1D("dtEffHist","Histogram of efficiencies", nBins, edges); for(Int_t i = 0; i < nBins; ++i) { Double_t lknot = effvals[i/upscale]; Double_t rknot = effvals[i/upscale + 1]; Double_t pos = 0.5*(edges[i]+edges[i+1]); //edges[i] + (edges[i+1]-edges[i]) * (1.0*(i%upscale)/upscale); Double_t weight = (pos - dtvals[i/upscale])/(dtvals[i/upscale + 1] - dtvals[i/upscale]); Double_t filling= lknot*(1-weight) + rknot*weight; std::cout << "lknot: " << lknot << "; rknot: " << rknot << "; pos: " << pos << "; filling: " << filling << std::endl; binFilling[i] = filling; dtEffHist->SetBinContent(i+1, binFilling[i]); } dtPdf->setEffiHist(dtEffHist); break; } case LauDecayTimePdf::EfficiencyMethod::Flat: { fitModel->setASqMaxValue(3.4); 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(kTRUE); 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; }