#include "TROOT.h"
#include "RooMsgService.h"
#include "RooRealVar.h"
#include "RooBDecay.h"
#include "RooAddPdf.h"
#include "RooCategory.h"
#include "RooTruthModel.h"
#include "RooArgSet.h"
#include "RooArgList.h"
#include "RooDataHist.h"
#include "RooGenericPdf.h"
#include "RooRealConstant.h"
#include "RooProdPdf.h"
#include "RooDataSet.h"
#include "TCanvas.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "RooRandom.h"
#include "RooNLLVar.h"
#include "RooMinuit.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "TString.h"
#include "RooHistPdf.h"
#include "RooErrorVar.h"
#include "RooEffProd.h"
#include "Roo1DTable.h"
#include "RooGaussModel.h"
#include "RooNumIntConfig.h"
#include "RooMCIntegrator.h"
#include "RooExtendPdf.h"
#include "RooExponential.h"
#include "RooDecay.h"
#include "RooGaussian.h"
#include "RooSimPdfBuilder.h"
#include "RooSimultaneous.h"
#include "RooGlobalFunc.h"
#include "RooEffProd.h"
#include "RooInt.h"
#include "RooAddModel.h"
#include "RooPolynomial.h"
#include "RooBVVSinhCoefVar.h"
#include "RooBVVSinCoefVar.h"
#include "RooBVVCoshCoefVar.h"
#include "RooBVVCosCoefVar.h"
#define mBs 5369.6
#define mW 50

using namespace RooFit;

RooDataSet *getExternal(RooRealVar& st,double mean=0.030, double width=0.005,int nEvents=-1)
{
  cout << "Generazione Pdf st con ";
  /*
  RooLandau stPdf("stPdf","distribution of per-event errors",st,RooRealConstant::value(mean),
		    RooRealConstant::value(width));
  cout << " Landau  "<<endl;
  */
  RooGaussian stPdf("stPdf","distribution of per-event errors",st,RooRealConstant::value(mean),
		    RooRealConstant::value(width));
  cout << " Gaussiana "<<endl;
   
  int NumEvents=10000;
  if(nEvents>0) NumEvents=nEvents;
  RooDataSet *data = stPdf.generate(st,NumEvents);
  return data;
}

RooProdPdf * makeBkgPromptPdf_byEv( const char *name,  const char *region,
                               const RooHistPdf& bkg_st_pdf,
                               RooRealVar  & parBkgPromptTimeMean,
                               RooRealVar  & parBkgPromptTimeSigma,
                               RooRealVar  & parBkgPromptMass,
                               RooRealVar  & mass,
                               RooRealVar  & time,  RooRealVar  & st,
                               RooRealVar  & cosTheta, 
                               RooRealVar  & cosPsi, 
                               RooRealVar  & phi,
                               RooCategory & tagFlav,
                               RooRealVar  & NbkgPrompt);

                               

RooProdPdf * makeBkgLongPdf_byEv( const char *name,  const char *region,
                             const RooHistPdf& bkg_st_pdf,
                             const RooResolutionModel& model,
                             RooRealVar  & parBkgLongTime1,
                             RooRealVar  & parBkgLongTime2,
                             RooRealVar  & fracBkgLongTime1,
                             RooRealVar  & parBkgLongMass,
                             RooRealVar  & mass,
                             RooRealVar  & time,   RooRealVar  & st,
                             RooRealVar  & cosTheta, 
                             RooRealVar  & cosPsi, 
                             RooRealVar  & phi,
                             RooCategory & tagFlav,
                             RooFormulaVar  & NbkgLong);



RooProdPdf * makeSignalPdf_byEv( const char *name, const char *region,
                            const RooHistPdf& sig_st_pdf,
                            const RooResolutionModel& model,
                            RooAbsReal & tau, 
                            RooRealVar & deltaGs,
                            RooRealVar & R0, 
                            RooRealVar & RT, 
                            RooRealVar & phiS, 
                            RooRealVar & dpar, 
                            //                            RooFormulaVar & dpar, 
                            RooRealVar & dperp,
                            RooRealVar & mistag,
                            RooRealVar & deltaMs,
                            RooRealVar & meanMass, 
                            RooRealVar & sigmaMassTail,
                            RooRealVar & sigmaMassCore,
                            RooRealVar & fracMassCore,
                            RooRealVar & mass,
                            RooRealVar & time,   RooRealVar & st,
                            RooRealVar & cosTheta, 
                            RooRealVar & cosPsi, 
                            RooRealVar & phi,
                            RooCategory& tagFlav,
                            RooRealVar & Ns);

RooProdPdf * makeBkgPromptPdf_fix( const char *name, const char *region,
                               RooRealVar  & parBkgPromptTimeMean,
                               RooRealVar  & parBkgPromptTimeSigma,
                               RooRealVar  & parBkgPromptMass,
                               RooRealVar  & mass,
                               RooRealVar  & time,
                               RooRealVar  & cosTheta, 
                               RooRealVar  & cosPsi, 
                               RooRealVar  & phi,
                               RooCategory & tagFlav,
                               RooRealVar  & NbkgPrompt);

                               

RooProdPdf * makeBkgLongPdf_fix( const char *name,  const char *region,
                             const RooResolutionModel& model,
                             RooRealVar  & parBkgLongTime1,
                             RooRealVar  & parBkgLongTime2,
                             RooRealVar  & fracBkgLongTime1,
                             RooRealVar  & parBkgLongMass,
                             RooRealVar  & mass,
                             RooRealVar  & time,
                             RooRealVar  & cosTheta, 
                             RooRealVar  & cosPsi, 
                             RooRealVar  & phi,
                             RooCategory & tagFlav,
                             RooFormulaVar  & NbkgLong);



RooProdPdf * makeSignalPdf_fix( const char *name, const char *region,
                            const RooResolutionModel& model,
                            RooAbsReal & tau, 
                            RooRealVar & deltaGs,
                            RooRealVar & R0, 
                            RooRealVar & RT, 
                            RooRealVar & phiS, 
                            RooRealVar & dpar, 
                            RooRealVar & dperp,
                            RooRealVar & mistag,
                            RooRealVar & deltaMs,
                            RooRealVar & meanMass, 
                            RooRealVar & sigmaMassTail,
                            RooRealVar & sigmaMassCore,
                            RooRealVar & fracMassCore,
                            RooRealVar & mass,
                            RooRealVar & time,
                            RooRealVar & cosTheta, 
                            RooRealVar & cosPsi, 
                            RooRealVar & phi,
                            RooCategory& tagFlav,
                            RooRealVar & Ns);

////////////////////////////////++++++++++++++++++++++++++++++++++++////////////////////////////////++++++++++++++++++++++++++++++++++++

int Gen_byEv(int argc, char *argv[]) {

  UInt_t rndSeed(5684);
  TString seq = "1";
  TString dirConfig = "../config";
  TString dirData   = "../";
  bool grid=true;       if(grid){gROOT->SetBatch(kTRUE); dirConfig="";  dirData=""; }

  TString configFile = dirConfig+"par.txt";
  TString ROOTFile = dirConfig+"PTerror.root";
  TString DATFile = dirData+"data_"+seq+".dat";
  cout << "Reading from. config "<< configFile<<" "<<endl;
  cout << "Reading PTerror distribution from "<< ROOTFile<<endl;
  
  double sig_sigma_st = 0.010;
  double sig_mean_st = 0.033;
  double bkg_sigma_st = 0.0105;
  double bkg_mean_st  = 0.036;
  
  //Define the observables
  double factor=1.;//80./31.5; //80./31.5;
  cout <<" Generating signal events with mean error "<<31.5*factor<<" fs"<<endl;

  RooRealVar * time  = new RooRealVar("time","proper time",-1,15.0,"ps");
  RooRealVar * st  = new RooRealVar("st","proper time error",0.005,factor*0.145,"ps");//Aggiunto errore sul tempo per file PT_error<0.80
  RooRealVar * mass = new RooRealVar("mass","invariant mass",mBs-3.*mW,mBs+3.*mW,"MeV/c^{2}");
  RooRealVar * cosTheta = new RooRealVar("cosTheta","cosTheta",-1.0,1.0);
  RooRealVar * cosPsi = new RooRealVar("cosPsi","cosPsi",-1.0,1.0);
  RooRealVar * phi = new RooRealVar("phi","phi",0.0,2*acos(-1.0));
  
  mass->setRange("SignalRegion",mBs-mW,mBs+mW);
  
  //Define variable
  RooRealVar * DM_s = new RooRealVar("DM_s","#DM_s",17.77,"ps^{-1}");
  RooRealVar * twobeta_s = new RooRealVar("twobeta_s","2#beta_s",0.0368,"rad");
  RooRealVar * gamma = new RooRealVar("gamma","#gamma",0.696,"ps");
  RooFormulaVar * tau = new RooFormulaVar("tau","#tau","1/gamma",RooArgSet(*gamma));
  RooRealVar * DG_s = new RooRealVar("DG_s","#DG_s",0.084,"ps^{-1}"); 
  
  //Longitudinal and transversal fractions
  RooRealVar * R_perp = new RooRealVar("R_perp","R_perp",0.233);
  RooRealVar * R_0 = new RooRealVar("R_0","R0",0.556);

  RooRealVar * d_perp = new RooRealVar("d_perp","dperp",2.91,"rad"); // set range to avoid unbiguity
  RooRealVar * d_par = new RooRealVar("d_par","d_par",-2.93,"rad"); // set range to avoid unbiguity
  
  //Mass resolution model
  // Signal mass parameter
  RooRealVar * meanMass  = new RooRealVar("meanMass","mean mass pdf",mBs);
  RooRealVar * sigmaMassTail = new RooRealVar ("sigmaMassTail","sigma mass core pdf",13.0);
  RooRealVar * sigmaMassCore = new RooRealVar ("sigmaMassCore","sigma mass core pdf",13.0);
  RooRealVar * fracMassCore = new RooRealVar ("fracMassCore","sigma mass core pdf",0.6,0,1);
 
 // Wrong tagging fraction
  RooRealVar * w = new RooRealVar("w","w",0.349);

  // Efficiency tagging fraction
  RooRealVar * Nt = new RooRealVar("Nt","Nt",40000);
  RooRealVar * Ns = new RooRealVar("Ns","Ns",10000);
  RooRealVar * Ngen = new RooRealVar("Ngen","Ngen",40000);
  RooRealVar * NbkgPrompt = new RooRealVar("NbkgPrompt","NbkgPrompt",25000);
  RooFormulaVar * NbkgLong = new RooFormulaVar("NbkgLong","NbkgLong","@0-@1-@2",RooArgSet(*Nt,*Ns,*NbkgPrompt)); //v1
  
  RooRealVar * parBkgPromptTimeMean = new RooRealVar("parBkgPromptTimeMean","parBkgPromptTimeMean",0.0,"ps");
  RooRealVar * parBkgPromptTimeSigma = new RooRealVar("parBkgPromptTimeSigma","parBkgPromptTimeSigma",0.06,"ps");

  RooRealVar * parBkgPromptMass = new RooRealVar("parBkgPromptMass","parBkgPromptMass",-0.0004);

  RooRealVar * parBkgLongTimeMean = new RooRealVar("parBkgLongTimeMean","parBkgLongTimeMean",0.0,"sigma units");
  RooRealVar * parBkgLongTimeSigmaCore = new RooRealVar("parBkgLongTimeSigmaCore","parBkgLongTimeSigmaCore",1,"sigma units");
  RooRealVar * parBkgLongTimeSigmaTail = new RooRealVar("parBkgLongTimeSigmaTail","parBkgLongTimeSigmaTail",1.61,"sigma units");
  RooRealVar * fracBkgLongTimeCore = new RooRealVar("fracBkgLongTimeCore","fracBkgLongTimeCore",0.6,0,1);

  RooRealVar * parBkgLongMass = new RooRealVar("parBkgLongMass","parBkgLongMass",-0.0025);
  RooRealVar * parBkgLongTime1 = new RooRealVar("parBkgLongTime1","parBkgLongTime1",1.15,"ps"); //v2
  RooRealVar * parBkgLongTime2 = new RooRealVar("parBkgLongTime2","parBkgLongTime2",0.15,"ps"); //v2
  RooRealVar * fracBkgLongTime1 = new RooRealVar("fracBkgLongTime1","fracBkgLongTime1",0.22,0,1); //v2


  // Resolution model
  RooRealVar * meanTime      = new RooRealVar("meanTime","meanTime",0.0);
  RooRealVar * sigmaTimeCore = new RooRealVar("sigmaTimeCore","sigmaTimeCore",1.);
  RooRealVar * sigmaTimeTail = new RooRealVar("sigmaTimeTail","sigmaTimeTail",2.);
  RooRealVar * fracTimeCore      = new RooRealVar("fracTimeCore","core fraction time resolution",0.76,0,1);
  
  //Tag category
  RooCategory * tagFlav = new RooCategory("tagFlav","tagFlav");
  tagFlav->defineType("B0",+1);
  tagFlav->defineType("B0Bar",-1);

   // Define a category with labels only
  RooCategory tagCat("tagCat","Tagging category") ;
  tagCat.defineType("c1",1) ;
  tagCat.defineType("c2",2) ;
.........
omitted
definition of *pdf_sig,*pdf_bkgLong,*pdf_bkgPrompt ....with some names
.........  
  RooAddPdf *pdf = new RooAddPdf("pdf","pdf",
                                 RooArgSet(*pdf_sig,*pdf_bkgLong,*pdf_bkgPrompt));
  
  
  // Instantiate a builder
  RooSimPdfBuilder b(*pdf) ;
  
  // Create an empty configuration object suitable for configuring 'sum'
  // A configuration object is a RooArgSet with RooStringVar objects
  RooArgSet& config = *b.createProtoBuildConfig() ;
  
  // Enter configuration data
  config.setStringValue("physModels", "pdf") ;    
  config.setStringValue("splitCats",  "tagCat") ; 
  config.setStringValue("pdf",       "tagCat : w,Nt,Ns,NbkgPrompt"); 
  RooDataSet *data = new RooDataSet("data",
                                    "data",
                                    RooArgList(RooArgSet(*mass,*time,*tagFlav,*cosTheta,*cosPsi,*phi,tagCat, *st)));

  // Build the PDF
  RooSimultaneous* simPdf  = (RooSimultaneous*)b.buildPdf(config,data) ;

  RooRandom::randomGenerator()->SetSeed(rndSeed);
......
omitted
......    
    data = simPdf->generate(RooArgSet(*mass,*time,*tagFlav,*cosTheta,*cosPsi,*phi,tagCat,*st),(int)Ntot);
    data->write("data_"+seq+".dat");
    return 1;
}


///////////////////===========================///////////////////===========================///////////////////===========================

int Fit_fix(int argc, char *argv[]) {

  UInt_t rndSeed(5684);
  TString seq = "1";
  TString dirConfig = "../config";
  TString dirData   = "../";
  bool grid=true;       if(grid){gROOT->SetBatch(kTRUE); dirConfig="";  dirData=""; }

  TString configFile = dirConfig+"par.txt";
  TString ROOTFile = dirConfig+"PTerror.root";
  TString DATFile = dirData+"data_"+seq+".dat";
  TString OUTFile = dirData+"result_"+seq+"_fix.root";
  cout << "Reading from config "<< configFile<<" "<<endl;


  //Define the observables
  RooRealVar * time  = new RooRealVar("time","proper time",-1,15.0,"ps");
  RooRealVar * mass = new RooRealVar("mass","invariant mass",mBs-3*mW,mBs+3*mW,"MeV/c^{2}");
  RooRealVar * cosTheta = new RooRealVar("cosTheta","cosTheta",-1.0,1.0);
  RooRealVar * cosPsi = new RooRealVar("cosPsi","cosPsi",-1.0,1.0);
  RooRealVar * phi = new RooRealVar("phi","phi",0.0,2*acos(-1.0));
  
  mass->setRange("SignalRegion",mBs-mW,mBs+mW);
  
  //Define variable
  RooRealVar * DM_s = new RooRealVar("DM_s","#DM_s",17.77,"ps^{-1}");
  RooRealVar * twobeta_s = new RooRealVar("twobeta_s","2#beta_s",0.0368,"rad");
  RooRealVar * gamma = new RooRealVar("gamma","#gamma",0.696,"ps");
  RooFormulaVar * tau = new RooFormulaVar("tau","#tau","1/gamma",RooArgSet(*gamma));
  RooRealVar * DG_s = new RooRealVar("DG_s","#DG_s",0.084,"ps^{-1}"); 
  
  //Longitudinal and transversal fractions
  RooRealVar * R_perp = new RooRealVar("R_perp","R_perp",0.233);
  RooRealVar * R_0 = new RooRealVar("R_0","R0",0.556);

  RooRealVar * d_perp = new RooRealVar("d_perp","dperp",2.91,"rad"); // set range to avoid unbiguity
  RooRealVar * d_par = new RooRealVar("d_par","d_par",-2.93,"rad"); // set range to avoid unbiguity
  
  //Mass resolution model
  // Signal mass parameter
  RooRealVar * meanMass  = new RooRealVar("meanMass","mean mass pdf",mBs);
  RooRealVar * sigmaMassTail = new RooRealVar ("sigmaMassTail","sigma mass core pdf",13.0);
  RooRealVar * sigmaMassCore = new RooRealVar ("sigmaMassCore","sigma mass core pdf",13.0);
  RooRealVar * fracMassCore = new RooRealVar ("fracMassCore","sigma mass core pdf",0.6,0,1);
 
 // Wrong tagging fraction
  RooRealVar * w = new RooRealVar("w","w",0.349);

  // Efficiency tagging fraction
  RooRealVar * Nt = new RooRealVar("Nt","Nt",40000);
  RooRealVar * Ns = new RooRealVar("Ns","Ns",10000);
  RooRealVar * Ngen = new RooRealVar("Ngen","Ngen",40000);
  RooRealVar * NbkgPrompt = new RooRealVar("NbkgPrompt","NbkgPrompt",25000);
  RooFormulaVar * NbkgLong = new RooFormulaVar("NbkgLong","NbkgLong","@0-@1-@2",RooArgSet(*Nt,*Ns,*NbkgPrompt));
  

  RooRealVar * parBkgPromptTimeMean = new RooRealVar("parBkgPromptTimeMean","parBkgPromptTimeMean",0.0,"ps");
  RooRealVar * parBkgPromptTimeSigma = new RooRealVar("parBkgPromptTimeSigma","parBkgPromptTimeSigma",0.06,"ps");

  RooRealVar * parBkgPromptMass = new RooRealVar("parBkgPromptMass","parBkgPromptMass",-0.0004);

  RooRealVar * parBkgLongTimeMean = new RooRealVar("parBkgLongTimeMean","parBkgLongTimeMean",0.0,"ps");
  RooRealVar * parBkgLongTimeSigmaCore = new RooRealVar("parBkgLongTimeSigmaCore","parBkgLongTimeSigmaCore",0.036,"sigma units");
  RooRealVar * parBkgLongTimeSigmaTail = new RooRealVar("parBkgLongTimeSigmaTail","parBkgLongTimeSigmaTail",0.058,"sigma units");
  RooRealVar * fracBkgLongTimeCore = new RooRealVar("fracBkgLongTimeCore","fracBkgLongTimeCore",0.6,0,1);

  RooRealVar * parBkgLongMass = new RooRealVar("parBkgLongMass","parBkgLongMass",-0.0025);
  RooRealVar * parBkgLongTime1 = new RooRealVar("parBkgLongTime1","parBkgLongTime1",1.15,"ps"); //v2
  RooRealVar * parBkgLongTime2 = new RooRealVar("parBkgLongTime2","parBkgLongTime2",0.15,"ps"); //v2
  RooRealVar * fracBkgLongTime1 = new RooRealVar("fracBkgLongTime1","fracBkgLongTime1",0.22,0,1); //v2


  // Resolution model
  RooRealVar * meanTime      = new RooRealVar("meanTime","meanTime",0.0);
  RooRealVar * sigmaTimeCore = new RooRealVar("sigmaTimeCore","sigmaTimeCore",0.030);
  RooRealVar * sigmaTimeTail = new RooRealVar("sigmaTimeTail","sigmaTimeTail",0.060);
  RooRealVar * fracTimeCore      = new RooRealVar("fracTimeCore","core fraction time resolution",0.76,0,1);
  
  //Tag category
  RooCategory * tagFlav = new RooCategory("tagFlav","tagFlav");
  tagFlav->defineType("B0",+1);
  tagFlav->defineType("B0Bar",-1);

   // Define a category with labels only
  RooCategory tagCat("tagCat","Tagging category") ;
  tagCat.defineType("c1",1) ;
  tagCat.defineType("c2",2) ;
  
  RooGaussModel * gaussTimeCore  = new RooGaussModel("gaussTimeCore","core gauss",*time,*meanTime,*sigmaTimeCore) ;
  RooGaussModel * gaussTimeTail  = new RooGaussModel("gaussTimeTail","tail gauss",*time,*meanTime,*sigmaTimeTail) ;
  
  cout<< "  Signal time resolution Model =======> DOUBLE Gaussian Ev fix"<<endl; 
  RooAddModel *resModel = new RooAddModel("resModel","core+tail gauss",
                                          RooArgList(*gaussTimeCore,*gaussTimeTail),*fracTimeCore) ;
  cout<< "  Long bkg time resolution Model =======> ";
  RooGaussModel * gaussBkgLongTimeCore  = new RooGaussModel("gaussBkgLongTimeCore","Bkg long core gauss",
                                                            *time,*parBkgLongTimeMean,*parBkgLongTimeSigmaCore) ;
  RooGaussModel * gaussBkgLongTimeTail  = new RooGaussModel("gaussBkgLongTimeTail","Bkg long tail gauss",
                                                            *time,*parBkgLongTimeMean,*parBkgLongTimeSigmaTail) ;

  cout<< "  DOUBLE Gaussian fixed "<<endl;
  RooAddModel *pdfBkgLongTimeResModel = new RooAddModel("pdfBkgLongTimeResModel","core+tail gauss resolution bkg long time",
                                         RooArgList(*gaussBkgLongTimeCore,*gaussBkgLongTimeTail),*fracBkgLongTimeCore) ;  // v1
  
  cout<< "  Prompt bkg time resolution Model =======> SINGLE Gaussian fix (independent)";
  RooProdPdf * pdf_bkgPrompt = makeBkgPromptPdf_fix("pdf_bkgPrompt", "SignalRegion",  w->setConstant(kFALSE);
                                                    *parBkgPromptTimeMean,*parBkgPromptTimeSigma,
                                                    *parBkgPromptMass,
                                                    *mass,
                                                    *time,
                                                    *cosTheta,
                                                    *cosPsi,
                                                    *phi,
                                                    *tagFlav,
                                                    *NbkgPrompt);

  RooProdPdf * pdf_bkgLong = makeBkgLongPdf_fix("pdf_bkgLong", "SignalRegion",
                                                *pdfBkgLongTimeResModel,
                                                *parBkgLongTime1,
                                                *parBkgLongTime2,
                                                *fracBkgLongTime1,
                                                *parBkgLongMass,
                                                *mass,
                                                *time,
                                                *cosTheta,
                                                *cosPsi,
                                                *phi,
                                                *tagFlav,
                                                *NbkgLong);
  
  
  RooProdPdf  * pdf_sig = makeSignalPdf_fix("pdf_sig", "SignalRegion",
                                            *resModel, 
                                            *tau,
                                            *DG_s,
                                            *R_0,
                                            *R_perp,
                                            *twobeta_s, 
                                            *d_par,
                                            *d_perp,
                                            *w,
                                            *DM_s,
                                            *meanMass,
                                            *sigmaMassTail,
                                            *sigmaMassCore,
                                            *fracMassCore,
                                            *mass,
                                            *time,
                                            *cosTheta,
                                            *cosPsi,
                                            *phi,
                                            *tagFlav,
                                            *Ns);
.........
omitted
definition of *pdf_sig,*pdf_bkgLong,*pdf_bkgPrompt ....with the same names as before
.........  
  
  RooAddPdf *pdf = new RooAddPdf("pdf","pdf",
                                 RooArgSet(*pdf_sig,*pdf_bkgLong,*pdf_bkgPrompt));


  
  // Instantiate a builder
  RooSimPdfBuilder b(*pdf) ;
  
  // Create an empty configuration object suitable for configuring 'sum'
  // A configuration object is a RooArgSet with RooStringVar objects
  RooArgSet& config = *b.createProtoBuildConfig() ;
  
  // Enter configuration data
  config.setStringValue("physModels", "pdf") ;    
  config.setStringValue("splitCats",  "tagCat") ; 
  config.setStringValue("pdf",       "tagCat : w,Nt,Ns,NbkgPrompt"); 
  RooDataSet *data = new RooDataSet("data",
                                    "data",
                                    RooArgList(RooArgSet(*mass,*time,*tagFlav,*cosTheta,*cosPsi,*phi,tagCat)));

  // Build the PDF
  RooSimultaneous* simPdf  = (RooSimultaneous*)b.buildPdf(config,data) ;

  RooRandom::randomGenerator()->SetSeed(rndSeed);

  
  RooArgSet* par = simPdf->getParameters(data) ;
  par->add(*Ngen);  
  par->add(*Nt);  
  par->add(*Ns);  
  par->add(*NbkgPrompt);  
  par->readFromFile(configFile.Data(), "READ", "Parameter"); 
  cout <<" Read parameters ============== "<<endl;
  par->Print("v");

    TString DATFile = dirData+"data_"+seq+".dat";
    cout << "Reading TOY data from "<< DATFile<<endl;
    data = RooDataSet::read(DATFile.Data(),RooArgSet(*mass,*time,*tagFlav,*cosTheta,*cosPsi,*phi,tagCat),"Q");

  RooNLLVar nll("nll","-log(L)",*simPdf,*data,Extended(kTRUE),Verbose(kTRUE));
  nll.Print("v");
  RooRealVar * offset = new RooRealVar("offset","offset",0);
  *offset = nll.getVal()+0.1;
  RooFormulaVar * L = new RooFormulaVar("L","@0-@1",RooArgList(nll,*offset));

  bool fit=true;
  
  if(fit) {
    TFile resFile(OUTFile.Data(),"RECREATE");
    RooMinuit m1(*L);    
    m1.setVerbose(kTRUE) ;
    m1.setProfile(1);
    m1.setPrintLevel(3);
    m1.setEps(1e-12);
    m1.setStrategy(1);
    m1.save()->Print("v");
    m1.simplex();
    m1.migrad();
    RooFitResult* r1 = m1.save();
    m1.hesse();
    if(m1.save()->covQual()!=3) {
      m1.migrad();
      m1.hesse();
    }
    RooFitResult* r2 = m1.save();
    cout << "fit result on NLL" << endl ;
    r1->Print("v") ;
    r2->Print("v") ;
    resFile.cd();
    r1->Write("",TObject::kSingleKey);
    r2->Write("",TObject::kSingleKey);
    resFile.Close();
  }
  
  return 1;
}

...... omitted .....

///////////////////===========================///////////////////===========================///////////////////===========================
////////////////////////////////++++++++++++++++++++++++++++++++++++////////////////////////////////++++++++++++++++++++++++++++++++++++