//=============================================================================
/** 
 * @file CheckIfMinosChangesParrErr.C
 * @brief Check if the parameter errors in fits to weighted data depend
 * on the Minos option in fitTo()
 *
 * @author Tim-Philip Huecking, <huecking@hep.physik.uni-siegen.de>
 * 
 * Created:
 * @date 20-05-25 (TPH)
 *
 * Changed:
 *
 */
//=============================================================================

// system includes
#include<iomanip>

// ROOT includes 
#include <RooStats/SPlot.h>
#include <RooArgSet.h>

// local includes
#include "modelBuilder.C"
//#include "Config.C"

RooDataSet* getSWeightedData(ModelBuilder& model, RooDataSet* data);
void compareFitRes(RooFitResult* withMinos, RooFitResult* noMinos);

void CheckIfMinosChangesParrErr()
{
   // Model parameters are stored in Configurator
   Configurator cfg{};
   cfg.prntModelPars();
   // Build a model with 3 variables:
   //    - mass (bkg: Exponential, sig: Gaussian)
   //    - vexex (bkg: Exponential, sig: Exponential)
   //    - vgg (bkg: Gaussian, sig: Gaussian)
   ModelBuilder model{cfg};
   bool useRandomSeed = false;
   model.buildFullModel(useRandomSeed);
   
   // generate a data set 
   model.m_pdfFull->Print();
   model.m_nSigFit->Print();
   model.m_nBkgFit->Print();
   RooDataSet* data = model.m_pdfFull->generate(RooArgList(
            *model.m_mass.get(), *model.m_vgg.get()), Extended(kTRUE));

   // Get a data set weighted with the signal sWeights
   RooDataSet *sigSwData = getSWeightedData(model, data);

   // Now fit the signal Gaussian of vgg to it 
   model.m_pdfExtVggSig->Print();
   cout << "\n\n----- FIT WITH MINOS AND SUMW2ERROR ENABLED -----\n\n";
   RooFitResult* wFitWithMinos = model.m_pdfExtVggSig->fitTo( *sigSwData
         , Extended(kTRUE), SumW2Error(kTRUE), PrintLevel(0)
         , Verbose(kFALSE), Save(kTRUE), Minos(kTRUE));
   cout << "\n\n----- FIT WITH MINOS DISABLED AND SUMW2ERROR ENABLED -----\n\n";
   RooFitResult* wFitNoMinos = model.m_pdfExtVggSig->fitTo( *sigSwData
         , Extended(kTRUE), SumW2Error(kTRUE), PrintLevel(0)
         , Verbose(kFALSE), Save(kTRUE), Minos(kFALSE));


   compareFitRes(wFitWithMinos, wFitNoMinos);
   
   return;
}

void compareFitRes(RooFitResult* withMinos, RooFitResult* noMinos)
{
   cout << "\n\n----- COMPARE THE FIT RESULTS  -----\n\n";
   cout << "Have look at the error of vggSigStdDevFit in both cases" << endl;
   cout << "WITH MINOS AND SUMW2ERROR ENABLED:" <<endl;
   withMinos->Print(); 
   RooArgList pars = withMinos->floatParsFinal();
   for(int i=0; i<pars.getSize(); i++)
   {
      auto par = pars.at(i);
      double cov = withMinos->covarianceMatrix()[i][i];
      cout << setw(35) << 
         TString::Format("Covariance for %s: ", pars[i].GetName()).Data()
         << setw(13) << cov << setw(13) << "Sqrt(cov): " << setw(10) 
         << TMath::Sqrt(cov) << endl;
   }
   cout << "=> Parameter error and Sqrt of covariance matrix element do NOT agree" << endl;
   cout << "\nWITH MINOS DISABLED AND SUMW2ERROR ENABLED:" << endl;
   pars = noMinos->floatParsFinal();
   noMinos->Print();
   for(int i=0; i<pars.getSize(); i++)
   {
      auto par = pars.at(i);
      double cov = noMinos->covarianceMatrix()[i][i];
      cout << setw(35) << 
         TString::Format("Covariance for %s: ", pars[i].GetName()).Data()
         << setw(13) << cov << setw(13) << "Sqrt(cov): " << setw(10) 
         << TMath::Sqrt(cov) << endl;
   }
   cout << "=> Paramter error and Sqrt of covariance matrix element do agree" << endl;
   //cout << "pars" << pars[0].GetName() << endl;

   return;
}

RooDataSet* getSWeightedData(ModelBuilder& model, RooDataSet* data)
{
   // Do a fit in the discriminating variable, as it is needed for 
   // sPlot method
   RooFitResult* discrRes = model.m_pdfMassSum->fitTo(*data 
         , Extended(kTRUE), Save(kTRUE), Verbose(kFALSE), Minos(kFALSE));

   // add sWeights to data set
   auto splot = new RooStats::SPlot("splot", "splot");
   splot->SetSData(data);
   splot->AddSWeight( model.m_pdfMassSum.get() //model.m_pdfFull.get()
         , RooArgList(*model.m_nSigFit.get(), *model.m_nBkgFit.get())
         , RooArgSet{}, kFALSE);
   data->Print();
   delete splot;

   // now perform a fit to weighted events
   // Therefore create a data set weighted with signal sWeights
   auto sigSwData = new RooDataSet("sigSwData"
         , "Data weighted with signal sWeights", RooArgSet(*data->get())
         , Import(*data), WeightVar("nSigFit_sw"));
   return sigSwData;
}