Hydra  4.0.1
A header-only templated C++ framework to perform data analysis on massively parallel platforms.
fit_gaussian.C
/*----------------------------------------------------------------------------
*
* Copyright (C) 2016 - 2023 Antonio Augusto Alves Junior
*
* This file is part of Hydra Data Analysis Framework.
*
* Hydra is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* Hydra is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with Hydra. If not, see <http://www.gnu.org/licenses/>.
*
*---------------------------------------------------------------------------*/
/*
* fit_gaussian.C
*
* Created on: 23/03/2018
* Author: Antonio Augusto Alves Junior
*/
/*!
* \example fit_gaussian.C
*/
#include <iostream>
#include <assert.h>
#include <time.h>
#include <chrono>
/**
*
*/
#ifndef HYDRA_HOST_SYSTEM
#define HYDRA_HOST_SYSTEM CPP
#endif
#ifndef HYDRA_DEVICE_SYSTEM
#define HYDRA_DEVICE_SYSTEM TBB
#endif
//this lib
#include <hydra/Function.h>
#include <hydra/Pdf.h>
#include <hydra/Zip.h>
//Minuit2
#include "Minuit2/FunctionMinimum.h"
#include "Minuit2/MnUserParameterState.h"
#include "Minuit2/MnPrint.h"
#include "Minuit2/MnMigrad.h"
#include "Minuit2/MnMinimize.h"
/*-------------------------------------
* Include classes from ROOT to fill
* and draw histograms and plots.
*-------------------------------------
*/
#include <TROOT.h>
#include <TH1D.h>
#include <TApplication.h>
#include <TCanvas.h>
using namespace ROOT::Minuit2;
using namespace hydra::arguments;
declarg(xvar, double)
/**
*
* This macro generates data and fit a Gaussian distribution in parallel.
*
* \note Call this macro in ROOT compile mode:
* ```
* root [0] .L fit_gaussian.C+
* root [1] fit_gaussian(50000000)
* ```
*
* \param nentries Number of points to generate.
*/
void fit_gaussian(size_t nentries=1000000 )
{
//-----------------
// some definitions
double min = -6.0;
double max = 6.0;
//Parameters for X direction
auto xmean = hydra::Parameter::Create("X-mean" ).Value(0.0).Error(0.0001).Limits(-1.0, 1.0);
auto xsigma = hydra::Parameter::Create("X-sigma").Value(1.0).Error(0.0001).Limits(0.01, 1.5);
//Gaussian distribution for X direction
//Model for X direction
//------------------------
TH1D* hist_xvar = new TH1D("hist_xvar", "X-axis", 100, min, max);
TH1D* hist_fit_xvar = new TH1D("hist_fit_xvar", "X-axis", 100, min, max);
//begin raii scope
{
//1D device buffer
//-------------------------------------------------------
//gaussian range
std::cout<< std::endl<< "Generated data:"<< std::endl;
for(size_t i=0; i<10; i++)
std::cout << "[" << i << "] :" << dataset[i] << std::endl;
auto xfcn = hydra::make_loglikehood_fcn(xmodel, dataset);
//-------------------------------------------------------
//fit
ROOT::Minuit2::MnPrint::SetGlobalLevel(3);
//minimization strategy
MnStrategy strategy(2);
//create Migrad minimizer
MnMigrad xmigrad(xfcn, xfcn.GetParameters().GetMnState() , strategy);
//print parameters before fitting
std::cout << xfcn.GetParameters().GetMnState() << std::endl;
//Minimize and profile the time
auto start = std::chrono::high_resolution_clock::now();
FunctionMinimum xminimum = FunctionMinimum( xmigrad(std::numeric_limits<unsigned int>::max(), 5));
auto stop = std::chrono::high_resolution_clock::now();
std::chrono::duration<double, std::milli> elapsed = stop - start;
//print minuit result
std::cout << " minimum: " << xminimum << std::endl;
//time
std::cout << "-----------------------------------------"<<std::endl;
std::cout << "| Time (ms) ="<< elapsed.count() <<std::endl;
std::cout << "-----------------------------------------"<<std::endl;
Hist_Data.Fill( dataset );
//draw data
for(size_t i=0; i<100; i++){
hist_xvar->SetBinContent(i+1, Hist_Data.GetBinContent(i) );
}
//draw fitted function
hist_fit_xvar->Sumw2();
for (size_t i=0 ; i<=100 ; i++) {
double x = hist_fit_xvar->GetBinCenter(i);
hist_fit_xvar->SetBinContent(i, xfcn.GetPDF()(x) );
}
hist_fit_xvar->Scale(hist_xvar->Integral()/hist_fit_xvar->Integral() );
}//end raii scope
//draw histograms
TCanvas* canvas = new TCanvas("canvas_d" ,"Distributions - Device", 500, 500);
hist_xvar->Draw("hist");
hist_fit_xvar->Draw("histsameC");
hist_fit_xvar->SetLineColor(2);
}