//// Study the fit with simple model, giving negligible bias

#include <RooRealVar.h>
#include <RooGaussian.h>
#include <RooPolynomial.h>
#include <RooArgusBG.h>
#include <RooFFTConvPdf.h>
#include <RooAddPdf.h>
#include <RooDataSet.h>
#include <RooPlot.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TTree.h>
#include <RooRandom.h>
#include "RooNumIntConfig.h"
#include "RooFitResult.h"
#include "RooMsgService.h"
#include "AsymBifurcatedCB.cpp"
#include "RooFormulaVar.h"
#include "TH2.h"
#include "RooCrystalBall.h"
void fit_model(Int_t seed=1234, int ntoy=10, TString file_out="something") {
    using namespace RooFit;

    // Define the observable
    RooRealVar mass("mass", "Mass", 4700, 6400);

    // Parameters for the Gaussian signal
    RooRealVar mean("mean", "Mean of Gaussian", 5283.4, 5240, 5320);
    RooRealVar sigma("sigma","",84.8,70,120);
    RooRealVar deltaS("deltaS","",7.14,0.,20.);
    RooFormulaVar sigmaL("sigmaL","Left sigma", "@0-@1",RooArgList(sigma,deltaS));
    RooFormulaVar sigmaR("sigmaR", "Right sigma", "@0+@1",RooArgList(sigma,deltaS));    
    RooRealVar nL_("nL_","",0.8);
    RooRealVar aL_("aL_","",2.2);
    RooRealVar nR_("nR_","",4.2);
    RooRealVar aR_("aR_","",-1.585);//, -5.0 , -5.00e-01);
    AsymBifurcatedCB signal("signal","", mass, mean, sigmaL,sigmaR, aL_, nL_, aR_, nR_);
    //RooCrystalBall signal("signal","", mass, mean, sigma, aL_, nL_, aR_, nR_);
    //RooGaussian signal("signal", "Gaussian Signal", mass, mean, sigma);
    
    // Parameters for the combinatorial background
    RooRealVar a1("a1", "a1", -1.23e-4, -1./6400, 0.);
    RooPolynomial bkg1("bkg1", "Combinatorial Background", mass, RooArgList(a1));
    
    RooRealVar sigN("sigN", "Signal Ntion", 1.6735e+05, 0.0, 3.0e+5);
    RooRealVar bkg1N("bkg1N", "Background 1 Ntion", 3.55e+04, 0.0, 1e+5);    
        
    RooAddPdf model("model", "Signal + Background", RooArgList(signal, bkg1), RooArgList(sigN, bkg1N));

    Int_t total_evt = 1.6735e+05 + 3.55e+04 ;

    /// Fill the out put in tree
    
    TFile f(file_out,"RECREATE");
    TTree tree("fit_params", "Fit Parameters");
    Double_t mean_Val, mean_Err, sigma_Val, sigma_Err;
    Double_t a1_Val, a1_Err, sigN_Val,sigN_Err, bkg1N_Val, bkg1N_Err;
    Double_t fitEdm, aR_Val, aR_Err, deltaS_Val, deltaS_Err;
    Int_t seed_, fitQual,fitStatus;

    tree.Branch("mean_", &mean_Val, "mean_/D");
    tree.Branch("mean_Err", &mean_Err, "mean_Err/D");
    tree.Branch("sigma_", &sigma_Val, "sigma_/D");
    tree.Branch("sigma_Err", &sigma_Err, "sigma_Err/D");
    tree.Branch("a1_", &a1_Val, "a1_/D");
    tree.Branch("a1_Err", &a1_Err, "a1_Err/D");
    tree.Branch("aR_", &aR_Val, "aR_/D");
    tree.Branch("aR_Err", &aR_Err, "aR_Err/D");
    tree.Branch("deltaS_", &deltaS_Val, "deltaS_/D");
    tree.Branch("deltaS_Err", &deltaS_Err, "deltaS_Err/D");
    tree.Branch("sigN_", &sigN_Val, "sigN_/D");
    tree.Branch("sigN_Err", &sigN_Err, "sigN_Err/D");
    tree.Branch("bkg1N_", &bkg1N_Val, "bkg1N_/D");
    tree.Branch("bkg1N_Err", &bkg1N_Err, "bkg1N_Err/D");
    tree.Branch("seed", &seed_,"seed/I");
    tree.Branch("fitQual", &fitQual,"fitQual/I");
    tree.Branch("fitStatus", &fitStatus,"fitStatus/I");
    tree.Branch("fitEdm", &fitEdm,"fitEdm/D");

    Int_t genEvt;
    tree.Branch("genEvt", &genEvt,"genEvt/I");
    
    RooRandom::randomGenerator()->SetSeed(seed);

    for(int itoy=0; itoy<ntoy;itoy++){
      mean.setVal(5283.4);
      sigma.setVal(84.8);
      deltaS.setVal(7.14);      
      a1.setVal(-1.23e-4);
      sigN.setVal(1.6735e+05);
      bkg1N.setVal(3.55e+04);
      aR_.setVal(-1.585);
      
      cout<<"total "<<sigN.getVal()+bkg1N.getVal()<<endl;
      
      RooDataSet* data = model.generate(RooArgSet(mass), Extended());
      genEvt = data->sumEntries();
      cout<<" number of generated entry: "<<data->sumEntries()<<endl;
      
      RooFitResult* fit_result = model.fitTo(*data, Save(true),Extended(true), NumCPU(8));//PrefitDataFraction(0.7),Minimizer("Minuit", "hesse"), NumCPU(8));
      fit_result->Print("v");
      fit_result->Print("v");

      cout<<" fitQual: "<<fit_result->covQual()<<endl;

      mean_Val = mean.getVal();
      mean_Err = mean.getError();
      sigma_Val = sigma.getVal();
      sigma_Err = sigma.getError();
      a1_Val = a1.getVal();
      a1_Err = a1.getError();
      deltaS_Val = deltaS.getVal();
      deltaS_Err = deltaS.getError();      
      sigN_Val = sigN.getVal();
      sigN_Err = sigN.getError();
      bkg1N_Val = bkg1N.getVal();
      bkg1N_Err = bkg1N.getError();      
      aR_Val=aR_.getVal();
      aR_Err=aR_.getError();      
      seed_ = seed;
      fitQual = fit_result->covQual();
      fitStatus = fit_result->status();
      fitEdm = fit_result->edm();
      tree.Fill();
      // TCanvas* c = new TCanvas("c", "Fit Result", 800, 600);
      // RooPlot* frame = mass.frame();
      // data->plotOn(frame);
      // model.plotOn(frame);
      // model.plotOn(frame, Components(signal), LineStyle(kDashed),LineColor(kGreen));
      // model.plotOn(frame, Components(bkg1), LineStyle(kDashed), LineColor(kRed));
      // model.plotOn(frame, Components(argus1conv), LineStyle(kDashed), LineColor(kOrange));
      // model.plotOn(frame, Components(argus2conv), LineStyle(kDashed), LineColor(kOrange+6));
      // frame->Draw();
      // gPad->SetLogy();      
      // c->SaveAs("fit_result.png");
      // gPad->SetLogy(false);
      // //gStyle->SetOptStat(0);
      // TH2 *hcorr = fit_result->correlationHist();
      // c->Clear();
      // hcorr->Draw("colz");
      // c->SaveAs("fit_result_correlation.png");
      
    }

    tree.Write();
    f.Close();

}
void MassFit_study_more_onlyTwoComp(Int_t seed=1234, int ntoy=10, TString file_out="something")
{  
   RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
   RooMsgService::instance().setSilentMode(kTRUE);
   RooMsgService::instance().setStreamStatus(1,false);
   RooMsgService::instance().getStream(1).removeTopic(RooFit::Integration) ;
   RooMsgService::instance().getStream(1).removeTopic(RooFit::Minimization) ;
   RooMsgService::instance().getStream(1).removeTopic(RooFit::Fitting) ;
   RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration) ;
   RooMsgService::instance().getStream(1).removeTopic(RooFit::Optimization) ;
   RooMsgService::instance().getStream(1).removeTopic(RooFit::ObjectHandling) ;
   RooMsgService::instance().getStream(1).removeTopic(RooFit::Eval) ;
   RooAbsReal::defaultIntegratorConfig()->setEpsAbs(1e-6);
   RooAbsReal::defaultIntegratorConfig()->setEpsRel(1e-6);
   RooAbsReal::defaultIntegratorConfig()->getConfigSection("RooIntegrator1D").setRealValue("maxSteps", 50);
  
  TString file_prefit = Form("%s_fit_parameters_aRfixtwoComp_%i_toy%i.root", file_out.Data(), seed, ntoy);  
  fit_model(seed, ntoy, file_prefit);
}