Problems importing to workspace

Roofit and RooStats
1

Hi,

I have a PDF & associated vars that I import to a RooWorkspace through a ModelConfig object and save to be accessed by another function. Simplified version:

        (lines defining many RooRealVars)
        RooFormulaVar BRconv("BRconv","@0*@1",RooArgList(BR,constBRrange));
        RooAbsPdf* bkg = new RooExponential("bkg","bkg", Mass,slope);
        RooAbsPdf* signal = new RooGaussian("signal", "signal", Mass, mean, sigma);
        RooAbsPdf* signal2 = new RooGaussian("signal2","gauss(x,mean,signal2)",ScaledMass, mean,sigma2);
        RooAddPdf sigmodel("sigmodel","sigmodel",RooArgList(*signal,*signal2), f2);
        RooAbsPdf* BsMass = new RooAddPdf("BsMass","BsMass",
                                   RooArgList(sigmodel,*bkg),RooArgList(BRconv,Nbkg));
        (lines making workspace, importing data, setting up priors and things...)
         w->defineSet("poi",RooArgSet(BR));
         w->defineSet("obs",RooArgSet(ScaledMass));
        // Make model
        ModelConfig model("model",w);
        model.SetPdf(*w->pdf("BsMass"));

When I then save it and open the workspace somewhere else, a fit to the data with the line

has BR fixed. Note this is not true if instead of BR in the mass fit, I use the number of events in the signal model. Does anyone know why this is?

Cheers,
Sean

2

Hi,

I am not sure I have understood your problem. Could you post eventually the minimal running code showing it ?

Thank you

Lorenzo

3

Hi Lorenzo,

Thank you for your reply. If you run the code below in association with your macro StandardHypoTestInvDemo, the fit in the beginning has BR constant.

Cheers,
Sean

        // Control options*********************
        double Bdmax=1.6e-6;
        double Bdmin=1.0e-7;
        int countmax=3;
        int numtoys=1000;
        // ************************************

        RooRealVar Mass("Mass","Mass(B_{s}^{0}) (MeV/c^{2})",5100.0,5600.0,"");

        RooRealVar mean("mean","mean", 5365.8); //variable
        RooRealVar sigma("Width","sigma", 13.6); //variable
        RooRealVar sigma2("Width2","sigma2", 29.5);
        RooRealVar slope("slope","slope", -0.000696); //variable

        RooRealVar constant("constant","constant",1.354e7);

        RooRealVar BR("BR","BR",Bdmin,Bdmin,Bdmax);
        RooFormulaVar BRconv("BRconv","@0*@1",RooArgList(BR,constant));

        // Constrain signal and background yields to be within errors
        double bkgMeanD = 74.6;
        double bkgErrD = TMath::Sqrt(74.6);
        RooRealVar bkgMean("bkgMean","bkgMean",bkgMeanD);
        RooRealVar bkgErr("bkgErr","bkgErr",bkgErrD);
        RooRealVar Nbkg("Nbkg","Nbkg",74.6,bkgMeanD-3.*bkgErrD,bkgMeanD+3.*bkgErrD);
        RooAbsPdf* gaussErrBkg = new RooGaussian("gaussErrBkg","gaussErrBkg",Nbkg,bkgMean,bkgErr);
        // Define pdf with Bd **********************************************************************
        RooAbsPdf* bkg = new RooExponential("bkg","bkg", Mass,slope);
        RooAbsPdf* signal = new RooGaussian("signal", "signal", Mass, mean, sigma);
        RooAbsPdf* signal2 = new RooGaussian("signal2","gauss(x,mean,signal2)",Mass, mean,sigma2);
        
        RooRealVar f2("f_{2}","f_{2}", 0.785, 0.00, 1.0);  f2.setConstant(true);
        RooAddPdf sigmodel("sigmodel","sigmodel",RooArgList(*signal,*signal2), f2);
        
        RooAbsPdf* BsMass = new RooAddPdf("BsMass","BsMass",
                                   RooArgList(sigmodel,*bkg),RooArgList(BRconv,Nbkg));
        RooDataSet* data = BsMass->generate(Mass,801);
        // Import pdf and data to a RooWorkspace
        RooWorkspace* w = new RooWorkspace("w");
        w->import(BR);
        w->import(*BsMass);
        w->import(*data);       
        w->import(*gaussErrBkg);        

        // Define our observable and parameters
        w->defineSet("obs",RooArgSet(Mass));
        //w->defineSet("poi",RooArgSet(NsigBd));
        w->defineSet("poi",RooArgSet(BR));
        
        // Make model
        ModelConfig sb_model("SB_model",w);
        sb_model.SetPdf(*w->pdf("BsMass"));
        sb_model.SetPriorPdf(*w->pdf("gaussErrBkg"));
        sb_model.SetObservables(*w->set("obs"));
        sb_model.SetParametersOfInterest(*w->set("poi"));
        w->import(sb_model);

        // save in to root file
        TFile* ftemp = new TFile("./savedWorkspace.root","RECREATE");
        w->Write();
        ftemp->Close();

        // Run the tool
   StandardHypoTestInvDemo("savedWorkspace.root","w","SB_model","","BsMassData",1,2,true,countmax,Bdmin,Bdmax,numtoys);
4

Hi,

The only problem I see in your model is that the range on BR is set too tight. Its best fit value is ~ 1.E-5 and you have set a maximum value of 1.6E-6. Otherwise I see it as a floating parameter in the fit.

One more thing, when defining the model, you should include in the global pdf also the Gaussian “gaussErrBkg”. Make then the full model as the product PDF of “gaussErrBkg” and “BsMass”

Lorenzo

5

Hi Lorenzo,

Even when I loosen the range, I see BR as a constant. Also, I am sure the range in the unsimplified version is correct.

Could this be due to a ROOT version? I am using 5.30.06.

Cheers,
Sean

6

Hi,
I see now the same problem as you when using 5.30. I think the cause of the problem is that you are defining a range in BR = [1.E-7,1.6E-6].
The macro, StandardHypoTestInvDemo, defines an alternate model (the background model) to be used in the test as the model with BR=0, which is outside the given range. This causes a problem when fitting because you have a value outside the range.
If you do not want to set th eminimum of BR to 0, you can try to define yourself the background model by doing:

      bModel = (ModelConfig*) sbModel->Clone();
      bModel->SetName(TString(modelSBName)+TString("_with_poi_0"));      
      RooRealVar * var = dynamic_cast<RooRealVar*>(bModel->GetParametersOfInterest()->first());
      if (!var) return 0;
      double oldval = var->getVal();
      //var->setVal(0);
      var->setVal(var->getMin() );
      bModel->SetSnapshot( RooArgSet(*var)  );
      var->setVal(oldval);

But as I said before the Hypothesis test inversion macro will probably not work anyway because also your upper value for the range is not right. The fit of the model should work.

Best Regards

Lorenzo

7

Hi Lorenzo,

Thank you for the assistance and advice. It took both setting the minimum to be zero and starting the scan at zero to work.

Cheers,
Sean