Problem with splot when running rs301_splot.C

[jelov]

Hi

I’m trying to modify the tutorial code rs301_splot.C to run on my data.
But I always meet " *** Break *** segmentation violation "

I modify the rs301_splot.C code to simply use 1-d pdf to fit on input data, which generate by also by rs301_splot.C code.

So the addModel() become like this :

RooRealVar invMass(“invMass”, “M_{inv}”, lowRange, highRange,“GeV”);
RooRealVar isolation(“isolation”, “isolation”, 0., 20., “GeV”);

RooRealVar mZ(“mZ”, “Z Mass”, 91.2, lowRange, highRange);
RooRealVar sigmaZ(“sigmaZ”, “Width of Gaussian”, 2,0,10,“GeV”);
RooGaussian zModel(“zModel”, “Z+jets Model”, invMass, mZ, sigmaZ);
// we know Z mass
mZ.setConstant();

RooConstVar zIsolDecayConst(“zIsolDecayConst”,
“z isolation decay constant”, -1);
RooExponential zIsolationModel(“zIsolationModel”, “z isolation model”,
isolation, zIsolDecayConst);
// make the combined Z model
// RooProdPdf zModel(“zModel”, “4-d model for Z”,
// RooArgSet(mZModel, zIsolationModel));
// --------------------------------------
RooRealVar qcdMassDecayConst(“qcdMassDecayConst”,
“Decay const for QCD mass spectrum”,
-0.01, -100, 100,“1/GeV”);
RooExponential qcdModel(“qcdModel”, “qcd Mass Model”,
invMass, qcdMassDecayConst);

RooConstVar qcdIsolDecayConst(“qcdIsolDecayConst”,
“Et resolution constant”, -.1);
RooExponential qcdIsolationModel(“qcdIsolationModel”, “QCD isolation model”,
isolation, qcdIsolDecayConst);
// make the 2-d model
// RooProdPdf qcdModel(“qcdModel”, “2-d model for QCD”,
// RooArgSet(qcdMassModel, qcdIsolationModel));

// --------------------------------------
// combined model
RooRealVar zYield(“zYield”,“fitted yield for Z”,50 ,0.,1000) ;
RooRealVar qcdYield(“qcdYield”,“fitted yield for QCD”, 100 ,0.,1000) ;
RooAddPdf model(“model”,“z+qcd background models”,
RooArgList(zModel, qcdModel),
RooArgList(zYield,qcdYield));
model.graphVizTree(“fullModel.dot”);
ws->import(model);

i.e , I replace the 2-d model to 1-d model , then use this modified code to run on the data I generated in another code.
But it will terminate with " *** Break *** segmentation violation ".
this message won’t appear if I use same code but do not replace the 2-d model to 1-d model.

Does anyone know what cause and problem and how to solve it? Thanks.

[jblomer]

Hi,

I can reproduce the segfault but I’m not an expert on RooStats. @StephanH can you help?

Cheers,
Jakob

[StephanH]

Hello @jelov,

you make a 1-D model but ask for 2-D event generation. My interpreter seems to be a bit nicer than yours (may I ask which root version you are using?), because it shows the line where the error is:

Error in <HandleInterpreterException>: Trying to dereference null pointer or trying to call routine taking non-null arguments.
Execution of your code was aborted.
In file included from input_line_10:1:
/home/shageboe/root-dbg/rs301_splot.C:149:60: warning: null passed to a callee that requires a non-null argument [-Wnonnull]
   RooDataSet* data = model->generate(RooArgSet(*invMass, *isolation),nEvents);

Don’t ask it to generate 2-D data here by using:
RooDataSet* data = model->generate(RooArgSet(*invMass),nEvents);

Further down, you also have to remove the plotting in the isolation frame.

Cheers,
Stephan

[jelov]

Hi @StephanH

Thanks for your reply.
The root version I use is 6.02
I do not make 1-D model but ask for 2-D event generation.
What I did is :

1. use the original rs301_splot.C code (which is 2-D model) to generate 2-D data, and save this into file
2. use another code modified from rs301_splot.C , to read the data generated from step1, in this code, change the 2-d model to 1-d mode to fit the data (if I do not change to 1st order model, the original code works)

And why I need to remove the plotting in the isolation frame ? Only Data is ploton this frame , right?

[StephanH]

@jelov

What I did is :

1. use the original rs301_splot.C code (which is 2-D model) to generate 2-D data, and save this into file
2. use another code modified from rs301_splot.C , to read the data generated from step1, in this code, change the 2-d model to 1-d mode to fit the data (if I do not change to 1st order model, the original code works)

I tried to do the same following your description in the first post. I guess you understand that it’s hard to debug the macro if you describe something I cannot see and I have to re-engineer what you did. Following your description, I pasted the function you posted into the original script. Doing that, I found it to be crashing when asking for 2-D event generation of a 1-D model.

If you want me to give another try, you may post the complete macro and maybe also the data and tell me how you used it.

Cheers,
Stephan

[jelov]

Hi, code is in the following:

1, use original code to generate data and save it into file for later use

 #include "RooRealVar.h"
 #include "RooStats/SPlot.h"
 #include "RooDataSet.h"
 #include "RooRealVar.h"
 #include "RooGaussian.h"
 #include "RooExponential.h"
 #include "RooChebychev.h"
 #include "RooAddPdf.h"
 #include "RooProdPdf.h"
 #include "RooAddition.h"
 #include "RooProduct.h"
 #include "TCanvas.h"
 #include "RooAbsPdf.h"
 #include "RooFit.h"
 #include "RooFitResult.h"
 #include "RooWorkspace.h"
 #include "RooConstVar.h"
 // use this order for safety on library loading
 using namespace RooFit;
 using namespace RooStats;
 // see below for implementation
 void AddModel(RooWorkspace*);
 void AddData(RooWorkspace*);
 void DoSPlot(RooWorkspace*);
 void MakePlots(RooWorkspace*);
 void rs302_splot()
 {
    // Create a new workspace to manage the project.
    RooWorkspace* wspace = new RooWorkspace("myWS");
    // add the signal and background models to the workspace.
    // Inside this function you will find a description our model.
    AddModel(wspace);
    // add some toy data to the workspace
    AddData(wspace);
    // inspect the workspace if you wish
    //  wspace->Print();
    // do sPlot.
    //This wil make a new dataset with sWeights added for every event.
    DoSPlot(wspace);
    // Make some plots showing the discriminating variable and
    // the control variable after unfolding.
    MakePlots(wspace);
    // cleanup
    delete wspace;
 }
 //____________________________________
 void AddModel(RooWorkspace* ws){
    // Make models for signal (Higgs) and background (Z+jets and QCD)
    // In real life, this part requires an intelligent modeling
    // of signal and background -- this is only an example.
    // set range of observable
    Double_t lowRange = 00, highRange = 200;
    // make a RooRealVar for the observables
    RooRealVar invMass("invMass", "M_{inv}", lowRange, highRange,"GeV");
    RooRealVar isolation("isolation", "isolation", 0., 20., "GeV");
    // --------------------------------------
    // make 2-d model for Z including the invariant mass
    // distribution  and an isolation distribution which we want to
    // unfold from QCD.
    std::cout << "make z model" << std::endl;
    // mass model for Z
    RooRealVar mZ("mZ", "Z Mass", 91.2, lowRange, highRange);
    RooRealVar sigmaZ("sigmaZ", "Width of Gaussian", 2,0,10,"GeV");
    RooGaussian mZModel("mZModel", "Z+jets Model", invMass, mZ, sigmaZ);
    // we know Z mass
    mZ.setConstant();
    // we leave the width of the Z free during the fit in this example.
    // isolation model for Z.  Only used to generate toy MC.
    // the exponential is of the form exp(c*x).  If we want
    // the isolation to decay an e-fold every R GeV, we use
    // c = -1/R.
    RooConstVar zIsolDecayConst("zIsolDecayConst",
                               "z isolation decay  constant", -1);
    RooExponential zIsolationModel("zIsolationModel", "z isolation model",
                                  isolation, zIsolDecayConst);
    // make the combined Z model
    RooProdPdf zModel("zModel", "4-d model for Z",
                      RooArgSet(mZModel, zIsolationModel));
    // --------------------------------------
    // make QCD model
    std::cout << "make qcd model" << std::endl;
    // mass model for QCD.
    // the exponential is of the form exp(c*x).  If we want
    // the mass to decay an e-fold every R GeV, we use
    // c = -1/R.
    // We can leave this parameter free during the fit.
    RooRealVar qcdMassDecayConst("qcdMassDecayConst",
                               "Decay const for QCD mass spectrum",
                               -0.01, -100, 100,"1/GeV");
    RooExponential qcdMassModel("qcdMassModel", "qcd Mass Model",
                               invMass, qcdMassDecayConst);
    // isolation model for QCD.  Only used to generate toy MC
    // the exponential is of the form exp(c*x).  If we want
    // the isolation to decay an e-fold every R GeV, we use
    // c = -1/R.
    RooConstVar qcdIsolDecayConst("qcdIsolDecayConst",
                                  "Et resolution constant", -.1);
    RooExponential qcdIsolationModel("qcdIsolationModel", "QCD isolation model",
                                  isolation, qcdIsolDecayConst);
    // make the 2-d model
    RooProdPdf qcdModel("qcdModel", "2-d model for QCD",
       RooArgSet(qcdMassModel, qcdIsolationModel));
    // --------------------------------------
    // combined model
    // These variables represent the number of Z or QCD events
    // They will be fitted.
    RooRealVar zYield("zYield","fitted yield for Z",50 ,0.,1000) ;
    RooRealVar qcdYield("qcdYield","fitted yield for QCD", 100 ,0.,1000) ;
    // now make the combined model
    std::cout << "make full model" << std::endl;
    RooAddPdf model("model","z+qcd background models",
                   RooArgList(zModel, qcdModel),
                   RooArgList(zYield,qcdYield));
    // interesting for debugging and visualizing the model
    model.graphVizTree("fullModel.dot");
    std::cout << "import model" << std::endl;
    ws->import(model);
 }
 //____________________________________
 void AddData(RooWorkspace* ws){
   
     RooAbsData::setDefaultStorageType(RooAbsData::Tree);
     TFile *fout=new TFile("sPlotTest.root","recreate");
     fout->cd();
    // Add a toy dataset
    // how many events do we want?
    Int_t nEvents = 1000;
    // get what we need out of the workspace to make toy data
    RooAbsPdf* model = ws->pdf("model");
    RooRealVar* invMass = ws->var("invMass");
    RooRealVar* isolation = ws->var("isolation");
    // make the toy data
    std::cout << "make data set and import to workspace" << std::endl;
    RooDataSet* data = model->generate(RooArgSet(*invMass, *isolation),nEvents);
     const TTree *tree=data->tree();
     tree->Write();   
     fout->Close();
 
    // import data into workspace
    ws->import(*data, Rename("data"));
 }
 //____________________________________
 void DoSPlot(RooWorkspace* ws){
    std::cout << "Calculate sWeights" << std::endl;
    // get what we need out of the workspace to do the fit
    RooAbsPdf* model = ws->pdf("model");
    RooRealVar* zYield = ws->var("zYield");
    RooRealVar* qcdYield = ws->var("qcdYield");
    RooDataSet* data = (RooDataSet*) ws->data("data");
    // fit the model to the data.
    model->fitTo(*data, Extended() );
    // The sPlot technique requires that we fix the parameters
    // of the model that are not yields after doing the fit.
    //
    // This *could* be done with the lines below, however this is taken care of
    // by the RooStats::SPlot constructor (or more precisely the AddSWeight
    // method).
    //
    //RooRealVar* sigmaZ = ws->var("sigmaZ");
    //RooRealVar* qcdMassDecayConst = ws->var("qcdMassDecayConst");
    //sigmaZ->setConstant();
    //qcdMassDecayConst->setConstant();
    RooMsgService::instance().setSilentMode(true);
    // Now we use the SPlot class to add SWeights to our data set
    // based on our model and our yield variables
    RooStats::SPlot* sData = new RooStats::SPlot("sData","An SPlot",
                                              *data, model, RooArgList(*zYield,*qcdYield) );
    // Check that our weights have the desired properties
    std::cout << "Check SWeights:" << std::endl;
    std::cout << std::endl <<  "Yield of Z is "
                << zYield->getVal() << ".  From sWeights it is "
                << sData->GetYieldFromSWeight("zYield") << std::endl;
    std::cout << "Yield of QCD is "
                << qcdYield->getVal() << ".  From sWeights it is "
                << sData->GetYieldFromSWeight("qcdYield") << std::endl
                << std::endl;
    for(Int_t i=0; i < 10; i++)
       {
       std::cout << "z Weight   " << sData->GetSWeight(i,"zYield")
                   << "   qcd Weight   " << sData->GetSWeight(i,"qcdYield")
                   << "  Total Weight   " << sData->GetSumOfEventSWeight(i)
                   << std::endl;
       }
    std::cout << std::endl;
    // import this new dataset with sWeights
    std::cout << "import new dataset with sWeights" << std::endl;
    ws->import(*data, Rename("dataWithSWeights"));
 }
 void MakePlots(RooWorkspace* ws){
    // Here we make plots of the discriminating variable (invMass) after the fit
    // and of the control variable (isolation) after unfolding with sPlot.
    std::cout << "make plots" << std::endl;
    // make our canvas
    TCanvas* cdata = new TCanvas("sPlot","sPlot demo", 400, 600);
    cdata->Divide(1,3);
    // get what we need out of the workspace
    RooAbsPdf* model = ws->pdf("model");
    RooAbsPdf* zModel = ws->pdf("zModel");
    RooAbsPdf* qcdModel = ws->pdf("qcdModel");
    RooRealVar* isolation = ws->var("isolation");
    RooRealVar* invMass = ws->var("invMass");
    // note, we get the dataset with sWeights
    RooDataSet* data = (RooDataSet*) ws->data("dataWithSWeights");
    // this shouldn't be necessary, need to fix something with workspace
    // do this to set parameters back to their fitted values.
    model->fitTo(*data, Extended() );
    //plot invMass for data with full model and individual components overlaid
    //  TCanvas* cdata = new TCanvas();
    cdata->cd(1);
    RooPlot* frame = invMass->frame() ;
    data->plotOn(frame ) ;
    model->plotOn(frame) ;
    model->plotOn(frame,Components(*zModel),LineStyle(kDashed), LineColor(kRed)) ;
    model->plotOn(frame,Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;
    frame->SetTitle("Fit of model to discriminating variable");
    frame->Draw() ;
    // Now use the sWeights to show isolation distribution for Z and QCD.
    // The SPlot class can make this easier, but here we demonstrate in more
    // detail how the sWeights are used.  The SPlot class should make this
    // very easy and needs some more development.
    // Plot isolation for Z component.
    // Do this by plotting all events weighted by the sWeight for the Z component.
    // The SPlot class adds a new variable that has the name of the corresponding
    // yield + "_sw".
    cdata->cd(2);
    // create weighted data set
    RooDataSet * dataw_z = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"zYield_sw") ;
    RooPlot* frame2 = isolation->frame() ;
    dataw_z->plotOn(frame2, DataError(RooAbsData::SumW2) ) ;
    frame2->SetTitle("isolation distribution for Z");
    frame2->Draw() ;
    // Plot isolation for QCD component.
    // Eg. plot all events weighted by the sWeight for the QCD component.
    // The SPlot class adds a new variable that has the name of the corresponding
    // yield + "_sw".
    cdata->cd(3);
    RooDataSet * dataw_qcd = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"qcdYield_sw") ;
    RooPlot* frame3 = isolation->frame() ;
    dataw_qcd->plotOn(frame3,DataError(RooAbsData::SumW2) ) ;
    frame3->SetTitle("isolation distribution for QCD");
    frame3->Draw() ;
    //  cdata->SaveAs("SPlot.gif");
 }

2. read data generate from 1, use 1-d model to do the splot (not request it to generate 2-D event here, only read data):

 #include "RooRealVar.h"
 #include "RooStats/SPlot.h"
 #include "RooDataSet.h"
 #include "RooRealVar.h"
 #include "RooGaussian.h"
 #include "RooExponential.h"
 #include "RooChebychev.h"
 #include "RooAddPdf.h"
 #include "RooProdPdf.h"
 #include "RooAddition.h"
 #include "RooProduct.h"
 #include "TCanvas.h"
 #include "RooAbsPdf.h"
 #include "RooFit.h"
 #include "RooFitResult.h"
 #include "RooWorkspace.h"
 #include "RooConstVar.h"
 #include "TFile.h"
 #include "TTree.h"
 
 // use this order for safety on library loading
 using namespace RooFit;
 using namespace RooStats;
 // see below for implementation
 void AddModel(RooWorkspace*);
 void AddData(RooWorkspace*);
 void DoSPlot(RooWorkspace*);
 void MakePlots(RooWorkspace*);
 void rs305_splot()
 {
    // Create a new workspace to manage the project.
    RooWorkspace* wspace = new RooWorkspace("myWS");
    // add the signal and background models to the workspace.
    // Inside this function you will find a description our model.
    AddModel(wspace);
    // add some toy data to the workspace
    AddData(wspace);
    // inspect the workspace if you wish
    //  wspace->Print();
    // do sPlot.
    //This wil make a new dataset with sWeights added for every event.
    DoSPlot(wspace);
    // Make some plots showing the discriminating variable and
    // the control variable after unfolding.
    MakePlots(wspace);
    // cleanup
    delete wspace;
 }
 //____________________________________
 void AddModel(RooWorkspace* ws){
    // Make models for signal (Higgs) and background (Z+jets and QCD)
    // In real life, this part requires an intelligent modeling
    // of signal and background -- this is only an example.
    // set range of observable
    Double_t lowRange = 00, highRange = 200;
    // make a RooRealVar for the observables
    RooRealVar invMass("invMass", "M_{inv}", lowRange, highRange,"GeV");
    RooRealVar isolation("isolation", "isolation", 0., 20., "GeV");
    // --------------------------------------
    // make 2-d model for Z including the invariant mass
    // distribution  and an isolation distribution which we want to
    // unfold from QCD.
    std::cout << "make z model" << std::endl;
    // mass model for Z
    RooRealVar mZ("mZ", "Z Mass", 91.2, lowRange, highRange);
    RooRealVar sigmaZ("sigmaZ", "Width of Gaussian", 2,0,10,"GeV");
    RooGaussian zModel("zModel", "Z+jets Model", invMass, mZ, sigmaZ);
    // we know Z mass
    mZ.setConstant();
    // we leave the width of the Z free during the fit in this example.
    // isolation model for Z.  Only used to generate toy MC.
    // the exponential is of the form exp(c*x).  If we want
    // the isolation to decay an e-fold every R GeV, we use
    // c = -1/R.
    RooConstVar zIsolDecayConst("zIsolDecayConst",
                               "z isolation decay  constant", -1);
    RooExponential zIsolationModel("zIsolationModel", "z isolation model",
                                  isolation, zIsolDecayConst);
    // make the combined Z model
 //   RooProdPdf zModel("zModel", "4-d model for Z",
 //                     RooArgSet(mZModel, zIsolationModel));
    // --------------------------------------
    // make QCD model
    std::cout << "make qcd model" << std::endl;
    // mass model for QCD.
    // the exponential is of the form exp(c*x).  If we want
    // the mass to decay an e-fold every R GeV, we use
    // c = -1/R.
    // We can leave this parameter free during the fit.
    RooRealVar qcdMassDecayConst("qcdMassDecayConst",
                               "Decay const for QCD mass spectrum",
                               -0.01, -100, 100,"1/GeV");
    RooExponential qcdModel("qcdModel", "qcd Mass Model",
                               invMass, qcdMassDecayConst);
    // isolation model for QCD.  Only used to generate toy MC
    // the exponential is of the form exp(c*x).  If we want
    // the isolation to decay an e-fold every R GeV, we use
    // c = -1/R.
    RooConstVar qcdIsolDecayConst("qcdIsolDecayConst",
                                  "Et resolution constant", -.1);
    RooExponential qcdIsolationModel("qcdIsolationModel", "QCD isolation model",
                                  isolation, qcdIsolDecayConst);
    // make the 2-d model
 //   RooProdPdf qcdModel("qcdModel", "2-d model for QCD",
 //      RooArgSet(qcdMassModel, qcdIsolationModel));
 
    // --------------------------------------
    // combined model
    // These variables represent the number of Z or QCD events
    // They will be fitted.
    RooRealVar zYield("zYield","fitted yield for Z",50 ,0.,1000) ;
    RooRealVar qcdYield("qcdYield","fitted yield for QCD", 100 ,0.,1000) ;
    // now make the combined model
    std::cout << "make full model" << std::endl;
    RooAddPdf model("model","z+qcd background models",
                   RooArgList(zModel, qcdModel),
                   RooArgList(zYield,qcdYield));
    // interesting for debugging and visualizing the model
    model.graphVizTree("fullModel.dot");
    std::cout << "import model" << std::endl;
    ws->import(model);
 }
 //____________________________________
 void AddData(RooWorkspace* ws){
   
     RooAbsData::setDefaultStorageType(RooAbsData::Tree);
     TFile *fout=new TFile("sPlotTest.root","read");
     fout->cd();
     TTree *tree=(TTree*)fout->Get("modelData");
    // Add a toy dataset
    // how many events do we want?
    /* Int_t nEvents = 1000; */
    // get what we need out of the workspace to make toy data
    /* RooAbsPdf* model = ws->pdf("model"); */
    RooRealVar* invMass = ws->var("invMass");
    RooRealVar* isolation = ws->var("isolation");
    // make the toy data
    std::cout << "make data set and import to workspace" << std::endl;
    /* RooDataSet* data = model->generate(RooArgSet(*invMass, *isolation),nEvents); */
    /* RooDataSet data = new RooDataSet("data","data",RooArgSet(*invMass,*isolation),Import(*tree)); */
    RooDataSet datatemp("data","data",RooArgSet(*invMass,*isolation),Import(*tree));
     RooDataSet *data=&datatemp; 
 
     /* const TTree *tree=data->tree(); */
     /* tree->Write();    */
     // fout->Close();
 
    // import data into workspace
    ws->import(*data, Rename("data"));
 }
 //____________________________________
 void DoSPlot(RooWorkspace* ws){
    std::cout << "Calculate sWeights" << std::endl;
    // get what we need out of the workspace to do the fit
    RooAbsPdf* model = ws->pdf("model");
    RooRealVar* zYield = ws->var("zYield");
    RooRealVar* qcdYield = ws->var("qcdYield");
    RooDataSet* data = (RooDataSet*) ws->data("data");
    // fit the model to the data.
    model->fitTo(*data, Extended() );
    // The sPlot technique requires that we fix the parameters
    // of the model that are not yields after doing the fit.
    //
    // This *could* be done with the lines below, however this is taken care of
    // by the RooStats::SPlot constructor (or more precisely the AddSWeight
    // method).
    //
    //RooRealVar* sigmaZ = ws->var("sigmaZ");
    //RooRealVar* qcdMassDecayConst = ws->var("qcdMassDecayConst");
    //sigmaZ->setConstant();
    //qcdMassDecayConst->setConstant();
    RooMsgService::instance().setSilentMode(true);
    // Now we use the SPlot class to add SWeights to our data set
    // based on our model and our yield variables
    RooStats::SPlot* sData = new RooStats::SPlot("sData","An SPlot",
                                              *data, model, RooArgList(*zYield,*qcdYield) );
    // Check that our weights have the desired properties
    std::cout << "Check SWeights:" << std::endl;
    std::cout << std::endl <<  "Yield of Z is "
                << zYield->getVal() << ".  From sWeights it is "
                << sData->GetYieldFromSWeight("zYield") << std::endl;
    std::cout << "Yield of QCD is "
                << qcdYield->getVal() << ".  From sWeights it is "
                << sData->GetYieldFromSWeight("qcdYield") << std::endl
                << std::endl;
    for(Int_t i=0; i < 10; i++)
       {
       std::cout << "z Weight   " << sData->GetSWeight(i,"zYield")
                   << "   qcd Weight   " << sData->GetSWeight(i,"qcdYield")
                   << "  Total Weight   " << sData->GetSumOfEventSWeight(i)
                   << std::endl;
       }
    std::cout << std::endl;
    // import this new dataset with sWeights
    std::cout << "import new dataset with sWeights" << std::endl;
    ws->import(*data, Rename("dataWithSWeights"));
 }
 void MakePlots(RooWorkspace* ws){
    // Here we make plots of the discriminating variable (invMass) after the fit
    // and of the control variable (isolation) after unfolding with sPlot.
    std::cout << "make plots" << std::endl;
    // make our canvas
    TCanvas* cdata = new TCanvas("sPlot","sPlot demo", 400, 600);
    cdata->Divide(1,3);
    // get what we need out of the workspace
    RooAbsPdf* model = ws->pdf("model");
    RooAbsPdf* zModel = ws->pdf("zModel");
    RooAbsPdf* qcdModel = ws->pdf("qcdModel");
    RooRealVar* isolation = ws->var("isolation");
    RooRealVar* invMass = ws->var("invMass");
    // note, we get the dataset with sWeights
    RooDataSet* data = (RooDataSet*) ws->data("dataWithSWeights");
    // this shouldn't be necessary, need to fix something with workspace
    // do this to set parameters back to their fitted values.
    model->fitTo(*data, Extended() );
    //plot invMass for data with full model and individual components overlaid
    //  TCanvas* cdata = new TCanvas();
    cdata->cd(1);
    RooPlot* frame = invMass->frame() ;
    data->plotOn(frame ) ;
    model->plotOn(frame) ;
    model->plotOn(frame,Components(*zModel),LineStyle(kDashed), LineColor(kRed)) ;
    model->plotOn(frame,Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;
    frame->SetTitle("Fit of model to discriminating variable");
    frame->Draw() ;
    // Now use the sWeights to show isolation distribution for Z and QCD.
    // The SPlot class can make this easier, but here we demonstrate in more
    // detail how the sWeights are used.  The SPlot class should make this
    // very easy and needs some more development.
    // Plot isolation for Z component.
    // Do this by plotting all events weighted by the sWeight for the Z component.
    // The SPlot class adds a new variable that has the name of the corresponding
    // yield + "_sw".
    cdata->cd(2);
    // create weighted data set
    RooDataSet * dataw_z = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"zYield_sw") ;
    RooPlot* frame2 = isolation->frame() ;
    dataw_z->plotOn(frame2, DataError(RooAbsData::SumW2) ) ;
    frame2->SetTitle("isolation distribution for Z");
    frame2->Draw() ;
    // Plot isolation for QCD component.
    // Eg. plot all events weighted by the sWeight for the QCD component.
    // The SPlot class adds a new variable that has the name of the corresponding
    // yield + "_sw".
    cdata->cd(3);
    RooDataSet * dataw_qcd = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"qcdYield_sw") ;
    RooPlot* frame3 = isolation->frame() ;
    dataw_qcd->plotOn(frame3,DataError(RooAbsData::SumW2) ) ;
    frame3->SetTitle("isolation distribution for QCD");
    frame3->Draw() ;
    //  cdata->SaveAs("SPlot.gif");
 }

[StephanH]

Hi @jelov,

no, please not like this. Please attach the files because the forum changes all the quotation marks and seems to remove the pointer *. I really don’t want to fix all the errors that occur due to what the forum tries to do to your code.

[jelov]

Hi @StephanH

Ah, I do not know there is a funtion I could upload files, the code I use is as following,

1. use original code to generate data and save it into file for later use
   rs302_splot.C (10.4 KB)

2. read data generate from 1, use 1-d model to do the splot (not request it to generate 2-D event here, only read data):
   rs305_splot.C (10.7 KB)

Thanks.

[StephanH]

I immediately get in line #135
[#0] WARNING:InputArguments -- RooAbsData::tree(modelData) WARNING: is not of StorageType::Tree. Use export_tree() instead or convert to tree storage.

So it looks like your data is not in a tree, but in a vector. See e.g.
https://root.cern.ch/doc/master/classRooAbsData.html#af56b00a9f46570acc940af933a7d32b3

The solution in the warning message seen above is outdated. You could use GetClonedTree() to get a tree.

---

If you fix that, you will run into problems in rs305:
In AddModel(), you are not saving/importing the variable “isolation”. Add
ws->Print();
at the end of the function to see your workspace. For me, it prints:

RooWorkspace(myWS) myWS contents

variables
---------
(invMass,mZ,qcdMassDecayConst,qcdYield,sigmaZ,zYield)

p.d.f.s
-------
RooAddPdf::model[ zYield * zModel + qcdYield * qcdModel ] = 0.245274
RooExponential::qcdModel[ x=invMass c=qcdMassDecayConst ] = 0.367879
RooGaussian::zModel[ x=invMass mean=mZ sigma=sigmaZ ] = 6.25215e-05

None of the PDFs contains “isolation”. Therefore, when you ask for it, you just get a nullptr. Anyway, you should always check the pointer if you ask for a variable or PDF, e.g.:

RooRealVar* isolation = ws->var("isolation");
if (isolation == nullptr) {
PROBLEM!
}

I hope that now you can figure out the details of your model.

Cheers,
Stephan