{
  ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
  ROOT::Math::MinimizerOptions::SetDefaultStrategy(1);
  ROOT::Math::MinimizerOptions::SetDefaultPrintLevel(-1);
  ROOT::Math::MinimizerOptions::SetDefaultTolerance(0.001);

  RooMsgService::instance().getStream(1).removeTopic(RooFit::NumIntegration);
  RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);

  bool debug = false;
  bool useResDS = false;

  TFile *fIn = TFile::Open("aDS.root");
  RooWorkspace *w = (RooWorkspace *)fIn->Get("w");

  // The data
  RooDataSet *theDS = (RooDataSet*)w->data("AllRunsHM2");
  theDS->Print();
  double nExp = theDS->sumEntries();
  double sumw = nExp, sumw2 = theDS->sumEntriesW2();
  double nEq = sumw*sumw/sumw2;
  double R   = std::sqrt(sumw2/sumw);
  std::cout << " sumw = " << sumw << " +/- " << std::sqrt(sumw2)
	    << " nEq = " << nEq
	    << " naive rescaling factor for the error to go from as data to as MC " << R
	    << std::endl;

  // relative uncertainty for constrain
  double vsys = 0.1;
  // forced slope
  double fslo = 1.5;
  double dslo = -fslo/55.;
  // absolute uncertainty
  double esys = std::abs(dslo)*vsys;
  // double min/max slope
  double a0min = -10., a0max = 10.;

  // The observable
  RooRealVar myy = *((RooRealVar*)theDS->get(0)->find("myy"));

  // The weight
  RooRealVar *wei = theDS->weightVar();

  // A dataset where sumEntries is the effective number of events...
  RooDataSet *rDS;
  if (useResDS) {
    rDS = new RooDataSet("rDS","rDS",RooArgSet(myy,*wei),RooFit::WeightVar("wei"));
    for (int ie = 0; ie < theDS->numEntries(); ie++) {
      myy.setVal(theDS->get(ie)->getRealValue("myy"));
      rDS->add(myy,theDS->weight()*sumw/sumw2);
    }
    rDS->Print();
  }
  
  // Yields
  RooRealVar nb("nb","nb", nExp, 1, 10000);
  //nb.setConstant();

  // Adding a Gaussian to constrain the slope
  double nuimax = (-55*a0min/fslo-1.)/vsys;
  double nuimin = (-55*a0max/fslo-1.)/vsys;
  RooRealVar nui("nui","nui",0.,nuimin,nuimax);
  RooRealVar glo("glo","glo",0.,nuimin,nuimax);
  glo.setConstant();
  RooConstVar sn("sn","sn",1.);
  RooGaussian co("co","co",nui,glo,sn);

  // Parameterise the exponential slope with l = fsl * (1 + sys * n), n is the floating parameter
  RooRealVar sys("sys","sys",vsys,0.,1.);
  sys.setConstant();
  RooRealVar fsl("fsl","fsl",dslo,-10.,10.);
  fsl.setConstant();
  RooArgList aL(fsl);
  aL.add(sys);
  aL.add(nui);
  RooFormulaVar ca0("ca0","ca0","@0*(1+@1*@2)",aL);
  RooExponential fb("fb","fb",myy,ca0);
  RooAddPdf fsb("fsb","fsb",RooArgList(fb),RooArgList(nb)); // just to be able to quickly add a signal 

  // And the model with constraint
  RooProdPdf fsbc("fsbc","fsbc",fsb,co);

  for (bool b : { true, false }) {

    if (useResDS && b)
      continue;

    std::cout << "\n\n*********************************************" << std::endl;
    if (!useResDS) 
      std::cout << "Fit as if " << (b ? "MC power\n" : "Data stat\n");
    else
      std::cout << "Fit as if Data stat, but using the rescaled DS, so should correspond to MC power\n"; 
    
    // Unconstrained fit
     // start from initial value
    nb.setVal(nExp);
    nb.setError(0.);
    nui.setVal(0.);
    nui.setError(0.);
    std::cout << "\nUnconstrained fit, the fit parameter is n in slope = fixedSlope x (1 + s x n), s = " << vsys << " fixedSlope = " << dslo << std::endl;
    RooFitResult *fitR = fsb.fitTo(useResDS ? *rDS : *theDS, RooFit::Save(), RooFit::PrintLevel(-1), RooFit::Minimizer("Minuit2","Migrad"), RooFit::SumW2Error(b));
    double ucSl = dslo*(1+nui.getVal()*vsys);
    double uceSl = esys*nui.getError();
    if (debug)
      fitR->Print("v");
     std::cout << "  fitted nuisance = " << nui.getVal() << " +/- " << nui.getError() << std::endl;
    std::cout << "Corresponding slope = " << ucSl << " +/- " << uceSl << std::endl;
    
    // If I add a "measurement" -1.5/55 +/- 10%, I expect
    double s = 1./(1./esys/esys + 1./uceSl/uceSl);
    double sl = (dslo/esys/esys + ucSl/uceSl/uceSl)*s;
    s = TMath::Sqrt(s);
    std::cout << "\nAdding a measurement " << dslo << " +/- " << esys << ", expect a combined measurement of ";
    std::cout << " slope = " << sl << " +/- " << s << std::endl;

    // Now adding a constraint : if as if MC stat, how can this be properly handled ?
    std::cout << "\nConstrained fit\n";
    // start from initial value
    nb.setVal(nExp);
    nb.setError(0.);
    nui.setVal(0.);
    nui.setError(0.);
    RooFitResult *cfitR = fsbc.fitTo(useResDS ? *rDS : *theDS, RooFit::Save(), RooFit::PrintLevel(-1), RooFit::Minimizer("Minuit2","Migrad"), RooFit::SumW2Error(b));
    if (debug)
      cfitR->Print("v");
    std::cout << "  fitted nuisance = " << nui.getVal() << " +/- " << nui.getError() << std::endl;
    std::cout << "Corresponding slope = " << dslo*(1+nui.getVal()*vsys) << " +/- " << esys*nui.getError() << std::endl;
    
    delete fitR;
    delete cfitR;
  }

}