diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index c93d14b..83328d7 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,47 +1,46 @@ list(APPEND EXAMPLE_SOURCES B3piKMatrixMassProj B3piKMatrixPlots CalcChiSq GenFit3K GenFit3KS GenFit3pi GenFitBelleCPKpipi GenFitDs2KKpi GenFitEFKLLM GenFitIncoherent_Bs2KSKpi GenFitKpipi GenFitNoDP GenFitNoDPMultiDim GenFitTimeDep GenFitTimeDep_Bs2KSKpi KMatrixDto3pi KMatrixExample Master MergeDataFiles mixedSampleTest PlotKMatrixTAmp PlotResults point2PointTestSample QuasiFlatSqDalitz ResultsExtractor Slave SlaveRooFit Test_Dpipi - Test_Dpipi_efficiencyHist ) if(NOT LAURA_BUILD_ROOFIT_SLAVE) list(REMOVE_ITEM EXAMPLE_SOURCES SlaveRooFit) endif() foreach( _example ${EXAMPLE_SOURCES}) add_executable(${_example} ${_example}.cc) target_link_libraries(${_example} PRIVATE Laura++) 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/Test_Dpipi.cc b/examples/Test_Dpipi.cc index 369e120..1f3d47c 100644 --- a/examples/Test_Dpipi.cc +++ b/examples/Test_Dpipi.cc @@ -1,382 +1,398 @@ #include using std::cout; using std::cerr; using std::endl; #include #include #include #include "TFile.h" #include "TH2.h" #include "TRandom.h" #include "TString.h" #include "TSystem.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" void usage(const TString& progName) { cerr<<"Usage:"< vector_argv{argv, argc+argv}; enum efficiencyInputMethod {Spline, Hist}; efficiencyInputMethod effiMethod = efficiencyInputMethod::Spline; auto histTest = std::find(std::begin(vector_argv), std::end(vector_argv), "--hist"); if(histTest != std::end(vector_argv)) { - effiMethod = efficiencyInputMethod::Spline; + effiMethod = efficiencyInputMethod::Hist; vector_argv.erase(histTest); } const TString dtype = vector_argv[2]; Int_t port = 0; Int_t iFit = 0; Int_t firstExpt = 0; Int_t firstExptGen = 0; Int_t nExpt = 1; Int_t nExptGen = 1; LauTimeDepFitModel::CPEigenvalue eigenvalue = LauTimeDepFitModel::CPEven; Bool_t fixPhiMix(kFALSE); Bool_t useSinCos(kTRUE); // check the command line arguments if (argc<1) { usage(vector_argv[0]); return EXIT_FAILURE; } TString command = vector_argv[1]; if (command != "gen" && command != "fit") { usage(vector_argv[0]); return EXIT_FAILURE; } if (command == "fit") { if (argc>3) { port = atoi(vector_argv[3]); if (argc>4) { iFit = atoi(vector_argv[4]); if (argc>5) { firstExpt = atoi(vector_argv[5]); if (argc>6) { nExpt = atoi(vector_argv[6]); if (argc>7) { nExptGen = atoi(vector_argv[7]); } } } } } for (firstExptGen = 0; firstExptGen<(firstExpt+nExpt); firstExptGen+=nExptGen) { } firstExptGen -= nExptGen; if ( (nExpt > nExptGen) || (nExptGen%nExpt != 0) ) { cerr<<"Error, nExpt must be a factor of nExptGen."<3) { firstExptGen = atoi(vector_argv[3]); if (argc>4) { nExptGen = atoi(vector_argv[4]); if (argc>5) { Int_t eigval = atoi(vector_argv[5]); if ( eigval == 1 ) { eigenvalue = LauTimeDepFitModel::CPOdd; } else { eigenvalue = LauTimeDepFitModel::CPEven; } } } } } cout<<"dtype "<addMassVeto( 2, 2.00776, 2.01276 ); daughtersB0bar = new LauDaughters("B0_bar", "pi+", "pi-", "D0"); daughtersB0 = new LauDaughters("B0", "pi+", "pi-", "D0_bar"); // efficiency effModelB0bar = new LauEffModel(daughtersB0bar, vetoes); effModelB0 = new LauEffModel(daughtersB0, vetoes); LauIsobarDynamics* sigModelB0bar(0); LauIsobarDynamics* sigModelB0(0); LauTimeDepFitModel* fitModel(0); std::vector params; LauFlavTag* flavTag = new LauFlavTag(); flavTag->setTrueTagVarName("trueTag"); // Use alternative tagging calibration parameters? //flavTag->useAveOmegaDeltaOmega(); TFile* etafile(0); TH1* etahist(0); Lau1DHistPdf* etahistpdf(0); //if (command == "gen"){ etafile = TFile::Open("histogram.root"); etahist = dynamic_cast(etafile->Get("htemp")); etahistpdf = new Lau1DHistPdf("eta",etahist,0.0,0.54,kTRUE,kFALSE); //} // signal dynamics 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); 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 fitModel = new LauTimeDepFitModel(sigModelB0bar,sigModelB0,flavTag); fitModel->setASqMaxValue(1.45); 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( eigenvalue ); fitModel->setPhiMix( 2.0*LauConstants::beta, fixPhiMix, useSinCos ); fitModel->setAsymmetries(0.0,kTRUE); // Tag cat fractions, dilutions and Delta dilutions flavTag->addTagger("OSTagger","tagVal_OS","mistagVal_OS", etahistpdf,0.20,0.5,1.0,0.20,0.5,1.0); flavTag->addTagger("SSTagger","tagVal_SS","mistagVal_SS", etahistpdf,0.30,0.5,1.0,0.30,0.5,1.0); // 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 Int_t nGauss(3); std::vector scale(nGauss); scale[0] = kTRUE; scale[1] = kTRUE; scale[2] = kFALSE; std::vector dtPars(10); 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 * mean1 = new LauParameter(mean1Name, -1.27); LauParameter * mean2 = new LauParameter(mean2Name, 0.0); LauParameter * sigma1 = new LauParameter(sigma1Name, 3.0); LauParameter * sigma2 = new LauParameter(sigma2Name, 8.0); LauParameter * frac1 = new LauParameter(frac1Name, 0.0930); LauParameter * frac2 = new LauParameter(frac2Name, 0.0036); LauParameter * tau = new LauParameter(tauName, 1.520); LauParameter * freq = new LauParameter(freqName, 0.5064); TString mean0tagcat63Name("dt_"); mean0tagcat63Name += 63; mean0tagcat63Name += "_mean_0"; TString sigma0tagcat63Name("dt_"); sigma0tagcat63Name += 63; sigma0tagcat63Name += "_sigma_0"; LauParameter * mean0tagcat63 = new LauParameter(mean0tagcat63Name, -0.031); LauParameter * sigma0tagcat63 = new LauParameter(sigma0tagcat63Name, 0.972); //Decay time acceptance spline - same for all tag cats (though doesn't have to be) 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); //Lau1DCubicSpline* dtEffSpline = new Lau1DCubicSpline(dtvals,effvals); Lau1DCubicSpline* dtEffSpline = nullptr; TH1D* dtEffHist = nullptr; switch(effiMethod) { case efficiencyInputMethod::Spline: dtEffSpline = new Lau1DCubicSpline(dtvals,effvals,Lau1DCubicSpline::StandardSpline,Lau1DCubicSpline::Natural,Lau1DCubicSpline::Natural); break; case efficiencyInputMethod::Hist: const Int_t nBins = 9; const Double_t edges[nBins + 1] = {0,1,2,4,6,8,11,14,17,20}; const Double_t binFilling[nBins] = {0.0075,0.02,0.029,0.655,0.69,0.715,0.74,0.76,0.785}; dtEffHist = new TH1D("dtEffHist","Histogram of efficiencies", nBins, edges); for( Int_t i = 1 ; i <= nBins ; ++i ){ dtEffHist->SetBinContent(i, binFilling[i-1]); } break; } dtPars[0] = mean0tagcat63; dtPars[1] = mean1->createClone(); dtPars[2] = mean2->createClone(); dtPars[3] = sigma0tagcat63; dtPars[4] = sigma1->createClone(); dtPars[5] = sigma2->createClone(); dtPars[6] = frac1->createClone(); dtPars[7] = frac2->createClone(); dtPars[8] = tau->createClone(); dtPars[9] = freq->createClone(); if (dtype=="CPEven"){ LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "deltaTAvg", "deltaTAvgErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::SimFitSigBd, nGauss, scale, LauDecayTimePdf::DecayTime ); dtPdf->doSmearing(kFALSE); - dtPdf->setEffiSpline(dtEffSpline); - dtPdf->setEffiHist(dtEffHist); + switch(effiMethod) + { + case efficiencyInputMethod::Spline: + dtPdf->setEffiSpline(dtEffSpline); + break; + + case efficiencyInputMethod::Hist: + dtPdf->setEffiHist(dtEffHist); + break; + } fitModel->setSignalDtPdf( dtPdf ); } else { LauDecayTimePdf * dtPdf = new LauDecayTimePdf( "deltaTAvg", "deltaTAvgErr", dtPars, minDt, maxDt, minDtErr, maxDtErr, LauDecayTimePdf::SimFitNormBd, nGauss, scale, LauDecayTimePdf::DecayTime ); dtPdf->doSmearing(kFALSE); - dtPdf->setEffiSpline(dtEffSpline); - dtPdf->setEffiHist(dtEffHist); + switch(effiMethod) + { + case efficiencyInputMethod::Spline: + dtPdf->setEffiSpline(dtEffSpline); + break; + + case efficiencyInputMethod::Hist: + dtPdf->setEffiHist(dtEffHist); + break; + } fitModel->setSignalDtPdf( dtPdf ); } // set the number of signal events cout<<"nSigEvents = "<setNSigEvents(nSigPar); // set the number of experiments if (command == "fit") { fitModel->setNExpts(nExpt, firstExpt); } else { fitModel->setNExpts(nExptGen, firstExptGen); } fitModel->useAsymmFitErrors(kFALSE); //fitModel->useRandomInitFitPars(kTRUE); fitModel->useRandomInitFitPars(kFALSE); fitModel->doPoissonSmearing(kFALSE); fitModel->doEMLFit(kFALSE); fitModel->writeLatexTable(kFALSE); TString dataFile(""); TString treeName("fitTree"); TString rootFileName(""); TString tableFileName(""); TString fitToyFileName(""); TString splotFileName(""); if (command == "fit") { dataFile = "TEST-Dpipi_"+dtype; dataFile += "_expts"; dataFile += firstExptGen; dataFile += "-"; dataFile += firstExptGen+nExptGen-1; dataFile += "_CP"; if ( eigenvalue == LauTimeDepFitModel::CPEven ) { dataFile += "even"; } else { dataFile += "odd"; } dataFile += ".root"; rootFileName = "fits/fit"; rootFileName += iFit; rootFileName += "_expts"; rootFileName += firstExpt; rootFileName += "-"; rootFileName += firstExpt+nExpt-1; rootFileName += ".root"; tableFileName = "fitResults_"; tableFileName += iFit; tableFileName += "_expts"; tableFileName += firstExpt; tableFileName += "-"; tableFileName += firstExpt+nExpt-1; fitToyFileName = "fitToyMC_"+dtype; fitToyFileName += iFit; fitToyFileName += "_expts"; fitToyFileName += firstExpt; fitToyFileName += "-"; fitToyFileName += firstExpt+nExpt-1; fitToyFileName += ".root"; splotFileName = "splot_"; splotFileName += iFit; splotFileName += "_expts"; splotFileName += firstExpt; splotFileName += "-"; splotFileName += firstExpt+nExpt-1; splotFileName += ".root"; } else { dataFile = "TEST-Dpipi"; switch(effiMethod) { case efficiencyInputMethod::Spline: dataFile += "_Spline"; break; case efficiencyInputMethod::Hist: dataFile += "_Hist"; break; } dataFile += "_"+dtype+"_expts"; dataFile += firstExptGen; dataFile += "-"; dataFile += firstExptGen+nExptGen-1; dataFile += "_CP"; if ( eigenvalue == LauTimeDepFitModel::CPEven ) { dataFile += "even"; } else { dataFile += "odd"; } dataFile += ".root"; rootFileName = "dummy.root"; tableFileName = "genResults"; } // 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 if ( command == "fit" ){ fitModel->runSlave( dataFile, treeName, rootFileName, tableFileName, "localhost", port ); } else { fitModel->run( command, dataFile, treeName, rootFileName, tableFileName ); } return EXIT_SUCCESS; }