#ifndef PSEUDO_EXPERIMENT_INL_
#define PSEUDO_EXPERIMENT_INL_
#include <iostream>
#include <assert.h>
#include <time.h>
#include <chrono>
#include <random>
#include <future>
#include <random>
#include <tclap/CmdLine.h>
#include <Minuit2/FunctionMinimum.h>
#include <Minuit2/MnUserParameterState.h>
#include <Minuit2/MnPrint.h>
#include <Minuit2/MnMigrad.h>
#include <Minuit2/MnMinimize.h>
#include <Minuit2/MnMinos.h>
#include <Minuit2/MnContours.h>
#include <Minuit2/CombinedMinimizer.h>
#include <Minuit2/MnPlot.h>
#include <Minuit2/MinosError.h>
#include <Minuit2/ContoursError.h>
#include <Minuit2/VariableMetricMinimizer.h>
#ifdef _ROOT_AVAILABLE_
#include <TROOT.h>
#include <TH1D.h>
#include <TApplication.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TStyle.h>
#endif //_ROOT_AVAILABLE_
{
size_t nstudies = 0;
bool verbose = false;
try {
TCLAP::CmdLine cmd("Command line arguments for ", '=');
TCLAP::ValueArg<size_t>
EArg(
"n",
"number-of-events",
"Number of events",
true, 10e6,
"size_t");
TCLAP::ValueArg<size_t> DArg("m", "number-of-studies","Number of MC studies", true, 1000, "size_t");
cmd.add(DArg);
TCLAP::SwitchArg VArg("v", "verbose","Print results of fits and covariance matrices");
cmd.add(VArg);
nstudies = DArg.getValue();
verbose = VArg.getValue();
}
catch (TCLAP::ArgException &e) {
std::cerr << "error: " << e.error() << " for arg " << e.argId()
<< std::endl;
}
double data_min = 0.0;
double data_max = 10.0;
auto splot_model =
hydra::add_pdfs( {N_Gaussian, N_Exponential}, Gaussian_PDF, Exponential_PDF);
splot_model.SetExtended(1);
double obs_min = 0.0;
double obs_max = 15.0;
BreitWigner_PDF, ChiSquare_PDF);
full_model.SetExtended(1);
{
auto seed = seeder();
auto gaussian_handler = std::async(std::launch::async,
std::normal_distribution<> dist(
mean,
sigma);
auto first = dataset.begin<XVar>();
auto last = dataset.begin<XVar>() +
int(
nentries/2);
auto it = first;
do
{
double x= dist(gen);
if((x > data_min) && (x < data_max) )
{
*it=x;
++it;
}
}while(it != last);
} );
seed = seeder();
auto exp_handler = std::async(std::launch::async,
[seed, data_min, data_max,
tau,
nentries, &dataset]( ){
std::exponential_distribution<> dist(-
tau);
auto first = dataset.begin<XVar>() +
int(
nentries/2);
auto last = dataset.end<XVar>();
auto it = first;
do
{
double x= dist(gen);
if((x > data_min) && (x < data_max) )
{
*it=x;
++it;
}
}while(it != last);
} );
seed = seeder();
auto breit_wigner_handler = std::async(std::launch::async,
std::uniform_real_distribution<> dist(0.0, 1.0);
auto first = dataset.begin<YVar>() ;
auto last = dataset.begin<YVar>() +
int(
nentries/2);
auto breit_wigner_dist = [](
double mean,
double width,
double rnd){
};
auto it = first;
do
{
double rnd = dist(gen);
double x =breit_wigner_dist(
mass,
width, rnd);
if((x > obs_min) && (x < obs_max) )
{
*it=x;
++it;
}
}while(it != last);
} );
seed = seeder();
auto noise_handler = std::async(std::launch::async,
std::chi_squared_distribution<> dist(
ndof);
auto first = dataset.begin<YVar>() +
int(
nentries/2);
auto last = dataset.end<YVar>();
auto it = first;
do
{
double x= dist(gen);
if((x > obs_min) && (x < obs_max) )
{
*it=x;
++it;
}
}while(it != last);
} );
gaussian_handler.wait();
exp_handler.wait();
breit_wigner_handler.wait();
noise_handler.wait();
std::random_shuffle(dataset.begin(), dataset.end());
if(verbose){
std::cout << " Dataset " << " size: " << dataset.size() << std::endl;
for(int i=0; i<100; i++){
std::cout << i << ") "<< dataset[i] << std::endl;
}
}
}
#ifdef _ROOT_AVAILABLE_
TH1::SetDefaultBufferSize(nstudies);
TH1D hist_data_sfit("data_sfit", "Control variable", 100, data_min, data_max);
TH1D hist_data_observable("data_observable", "Observable", 100, obs_min, obs_max);
TH1D hist_N_Gaussian("N_Gaussian", "N_Gaussian", 100, 1.0, 0.0);
TH1D hist_N_Exponential("N_Exponential", "N_Exponential", 100, 1.0, 0.0);
TH1D hist_mean("mean", "mean", 100, 1.0, 0.0);
TH1D hist_sigma("sigma", "sigma", 100, 1.0, 0.0 );
TH1D hist_tau("tau", "tau", 100, -0.25, -0.15 );
TH1D hist_mass("mass", "mass", 100, 1.0, 0.0 );
TH1D hist_mass_error("mass_error", "mass error", 100, 1.0, 0.0 );
TH1D hist_width("width", "width", 100, 1.0, 0.0 );
TH1D hist_width_error("width_error", "width error", 100, 1.0, 0.0 );
TH1D hist_full_mass("mass_full", "mass (complete model)", 100, 1.0, 0.0 );
TH1D hist_full_mass_error("mass_error_full", "mass error (complete model)", 100, 1.0, 0.0 );
TH1D hist_full_width("width_full", "width (complete model)", 100, 1.0, 0.0 );
TH1D hist_full_width_error("width_error_full", "width error (complete model)", 100, 1.0, 0.0 );
for(auto x: dataset){
hist_data_sfit.Fill( hydra::get<0>(x) );
hist_data_observable.Fill( hydra::get<1>(x) );
}
#endif //_ROOT_AVAILABLE_
{
for(size_t study=0; study<nstudies; study++)
{
bs_range.begin() + dataset.size());
if(verbose){
ROOT::Minuit2::MnPrint::SetGlobalLevel(3);
}
else{
ROOT::Minuit2::MnPrint::SetGlobalLevel(-1);
}
MnMigrad migrad_splot(splot_fcn, splot_fcn.GetParameters().GetMnState(),
strategy);
auto start = std::chrono::high_resolution_clock::now();
FunctionMinimum minimum_splot = FunctionMinimum(migrad_splot(5000, 1));
auto end = std::chrono::high_resolution_clock::now();
if(verbose){
std::cout<<"SFit minimum: "<< minimum_splot << std::endl;
std::cout << "-----------------------------------------"<<std::endl;
std::cout <<
"| [Fit Time] (ms) = " <<
elapsed.count() <<std::endl;
std::cout << "-----------------------------------------"<<std::endl;
}
auto covar = sweigts.GetCovMatrix();
if(verbose){
std::cout << "SFit covariance matrix: "
<< std::endl
<< covar
<< std::endl;
std::cout << std::endl
<< "sWeights ( Gaussian, Exponential):"
<< std::endl;
for(size_t i = 0; i< 10; i++)
std::cout << i << ") : "
<< sweigts[i]
<< std::endl;
}
start = std::chrono::high_resolution_clock::now();
FunctionMinimum minimum = FunctionMinimum(migrad(1000, 5));
end = std::chrono::high_resolution_clock::now();
if(verbose){
std::cout << std::endl <<"Fit background subtracted minimum: "
<< minimum
<< std::endl
<< std::endl;
std::cout << "-----------------------------------------"<<std::endl;
std::cout <<
"| Fit time (ms) = " <<
elapsed.count() <<std::endl;
std::cout << "-----------------------------------------"<<std::endl;
}
MnMigrad full_migrad(full_fcn, full_fcn.GetParameters().GetMnState(),
strategy);
start = std::chrono::high_resolution_clock::now();
FunctionMinimum full_minimum = FunctionMinimum(full_migrad(5000, 1));
end = std::chrono::high_resolution_clock::now();
if(1){
std::cout << std::endl <<"Full fit minimum: "
<< full_minimum
<< std::endl
<< std::endl;
std::cout << "-----------------------------------------"<<std::endl;
std::cout <<
"| Full fit time (ms) = " <<
elapsed.count() <<std::endl;
std::cout << "-----------------------------------------"<<std::endl;
std::cout << "----------------------------------------" << std::endl;
std::cout << "Logging result #" << study << " " << std::endl;
std::cout << "----------------------------------------" << std::endl;
std::cout << "N_Gaussian ........... " << minimum_splot.UserParameters().Value("N_Gaussian") <<std::endl;
std::cout << "N_Exponential ........ " << minimum_splot.UserParameters().Value("N_Exponential") <<std::endl;
std::cout << "mean ................. " << minimum_splot.UserParameters().Value("mean") <<std::endl;
std::cout << "sigma ................ " << minimum_splot.UserParameters().Value("sigma")<<std::endl ;
std::cout << "tau .................. " << minimum_splot.UserParameters().Value("tau") <<std::endl;
std::cout << "mass ................. " << minimum.UserParameters().Value("mass") <<std::endl;
std::cout << "mass error ........... " << minimum.UserParameters().Error("mass") <<std::endl;
std::cout << "width ................ " << minimum.UserParameters().Value("width") <<std::endl;
std::cout << "width error .......... " << minimum.UserParameters().Error("width") <<std::endl;
std::cout << "mass (full) .......... " << full_minimum.UserParameters().Value("mass") <<std::endl;
std::cout << "mass error (full) .... " << full_minimum.UserParameters().Error("mass") <<std::endl;
std::cout << "width (full) ......... " << full_minimum.UserParameters().Value("width") <<std::endl;
std::cout << "width error (full) ... " << full_minimum.UserParameters().Error("width") <<std::endl;
std::cout << "----------------------------------------" << std::endl;
}
minimum_splot.UserParameters().Value("N_Gaussian"),
minimum_splot.UserParameters().Value("N_Exponential"),
minimum_splot.UserParameters().Value("mean"),
minimum_splot.UserParameters().Value("sigma"),
minimum_splot.UserParameters().Value("tau"),
minimum.UserParameters().Value("mass"),
minimum.UserParameters().Error("mass"),
minimum.UserParameters().Value("width"),
minimum.UserParameters().Error("width"),
full_minimum.UserParameters().Value("mass"),
full_minimum.UserParameters().Error("mass"),
full_minimum.UserParameters().Value("width"),
full_minimum.UserParameters().Error("width")
));
}
if(verbose){
std::cout << std::endl << std::endl << std::endl;
std::cout << "----------------------------------------" << std::endl;
std::cout << "Dumping logged studies: " << variable_log.size() << std::endl;
std::cout << "----------------------------------------" << std::endl;
}
for(auto x:variable_log)
{
if(verbose) std::cout << x << std::endl;
hist_N_Gaussian.Fill(hydra::get<0>(x));
hist_N_Exponential.Fill(hydra::get<1>(x));
hist_mean.Fill(hydra::get<2>(x));
hist_sigma.Fill(hydra::get<3>(x));
hist_tau.Fill(hydra::get<4>(x));
hist_mass.Fill(hydra::get<5>(x));
hist_mass_error.Fill(hydra::get<6>(x));
hist_width.Fill(hydra::get<7>(x));
hist_width_error.Fill(hydra::get<8>(x));
hist_full_mass.Fill(hydra::get<9>(x));
hist_full_mass_error.Fill(hydra::get<10>(x));
hist_full_width.Fill(hydra::get<11>(x));
hist_full_width_error.Fill(hydra::get<12>(x));
}
}
#ifdef _ROOT_AVAILABLE_
TApplication *myapp=new TApplication("myapp",0,0);
TCanvas canvas_13("canvas_13" ,"Dataset", 1000, 500);
canvas_13.Divide(2);
canvas_13.cd(1);
hist_data_sfit.Draw("hist");
canvas_13.cd(2);
hist_data_observable.Draw("hist");
gStyle->SetOptFit(111);
TCanvas canvas_0("canvas_0" ,"", 500, 500);
hist_N_Gaussian.SetLineWidth(2);
hist_N_Gaussian.Draw("E");
hist_N_Gaussian.Fit("gaus","ML");
canvas_0.Update();
TCanvas canvas_1("canvas_1" ,"", 500, 500);
hist_N_Exponential.SetLineWidth(2);
hist_N_Exponential.Draw("E");
hist_N_Exponential.Fit("gaus","ML");
canvas_1.Update();
TCanvas canvas_2("canvas_2" ,"", 500, 500);
hist_mean.SetLineWidth(2);
hist_mean.Draw("E");
hist_mean.Fit("gaus","ML");
canvas_2.Update();
TCanvas canvas_3("canvas_3" ,"", 500, 500);
hist_sigma.SetLineWidth(2);
hist_sigma.Draw("E");
hist_sigma.Fit("gaus","ML");
canvas_3.Update();
TCanvas canvas_4("canvas_4" ,"", 500, 500);
hist_tau.SetLineWidth(2);
hist_tau.Draw("E");
hist_tau.Fit("gaus","ML");
canvas_4.Update();
TCanvas canvas_5("canvas_5" ,"", 500, 500);
hist_mass.SetLineWidth(2);
hist_mass.Draw("E");
hist_mass.Fit("gaus","ML");
canvas_5.Update();
TCanvas canvas_6("canvas_6" ,"", 500, 500);
hist_mass_error.SetLineWidth(2);
hist_mass_error.Draw("E");
hist_mass_error.Fit("gaus","ML");
canvas_6.Update();
TCanvas canvas_7("canvas_7" ,"", 500, 500);
hist_width.SetLineWidth(2);
hist_width.Draw("E");
hist_width.Fit("gaus","ML");
canvas_7.Update();
TCanvas canvas_8("canvas_8" ,"", 500, 500);
hist_width_error.Draw("hist");
hist_width_error.SetLineWidth(2);
hist_width_error.Draw("E");
hist_width_error.Fit("gaus","ML");
canvas_8.Update();
TCanvas canvas_9("canvas_9" ,"", 500, 500);
hist_full_mass.SetLineWidth(2);
hist_full_mass.Draw("E");
hist_full_mass.Fit("gaus","ML");
canvas_9.Update();
TCanvas canvas_10("canvas_10" ,"", 500, 500);
hist_full_mass_error.SetLineWidth(2);
hist_full_mass_error.Draw("E");
hist_full_mass_error.Fit("gaus","ML");
canvas_10.Update();
TCanvas canvas_11("canvas_11" ,"", 500, 500);
hist_full_width.SetLineWidth(2);
hist_full_width.Draw("E");
hist_full_width.Fit("gaus","ML");
canvas_11.Update();
TCanvas canvas_12("canvas_12" ,"", 500, 500);
hist_full_width_error.Draw("hist");
hist_full_width_error.SetLineWidth(2);
hist_full_width_error.Draw("E");
hist_full_width_error.Fit("gaus","ML");
canvas_12.Update();
myapp->Run();
#endif //_ROOT_AVAILABLE_
return 0;
}
#endif