Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[RF] Broken weights after reducing RooDataSet created with RooAbsPdf::generate() #6951

Closed
1 task done
huebner-m opened this issue Dec 14, 2020 · 1 comment · Fixed by #10397
Closed
1 task done

Comments

@huebner-m
Copy link

  • Checked for duplicates

Describe the bug

When creating a RooDataSet from a PDF using RooAbsPdf::generate() and afterwards plotting it in a specific region using RooDataSet::reduce() the internal weights of the reduced RooDataSet are broken, i.e. they are all set to 1.

Expected behavior

The weights of the reduced RooDataSet should not change to 1 when calling RooDataSet::reduce().

To Reproduce

Run included macro with root broken_weights.C.

#include "RooStats/HistFactory/Measurement.h"
#include "RooStats/HistFactory/MakeModelAndMeasurementsFast.h"
#include "TFile.h"
#include "TROOT.h"
 
using namespace RooStats;
using namespace HistFactory;

void create_test_ws(){
   // taken from https://root.cern.ch/doc/v608/histfactory_2example_8C_source.html
   std::string InputFile = "example.root";
      // Create the measurement
   Measurement meas("meas", "meas");

   meas.SetOutputFilePrefix( "./results/example_UsingC" );
   meas.SetPOI( "SigXsecOverSM" );
   meas.AddConstantParam("alpha_syst1");
   meas.AddConstantParam("Lumi");

   meas.SetLumi( 1.0 );
   meas.SetLumiRelErr( 0.10 );
   meas.SetExportOnly( true );
   meas.SetBinHigh( 2 );

   // Create a channel

   Channel chan( "channel1" );
   chan.SetData( "data", InputFile );
   chan.SetStatErrorConfig( 0.05, "Poisson" );

   // Create another channel
   Channel chan2( "channel2" );
   chan2.SetData( "data", InputFile );
   chan2.SetStatErrorConfig( 0.05, "Poisson" );

   // Now, create some samples

   // Create the signal sample
   Sample signal( "signal", "signal", InputFile );
   signal.AddOverallSys( "syst1",  0.95, 1.05 );
   signal.AddNormFactor( "SigXsecOverSM", 1, 0, 3 );
   chan.AddSample( signal );
   chan2.AddSample( signal );

   // Background 1
   Sample background1( "background1", "background1", InputFile );
   background1.ActivateStatError( "background1_statUncert", InputFile );
   background1.AddOverallSys( "syst2", 0.95, 1.05  );
   chan.AddSample( background1 );
   chan2.AddSample( background1 );

   // Background 2
   Sample background2( "background2", "background2", InputFile );
   background2.ActivateStatError();
   background2.AddOverallSys( "syst3", 0.95, 1.05  );
   chan.AddSample( background2 );
   chan2.AddSample( background2 );

   // Done with this channel
   // Add it to the measurement:
   meas.AddChannel( chan );
   meas.AddChannel( chan2 );

   // Collect the histograms from their files,
   // print some output,
   meas.CollectHistograms();
   meas.PrintTree();

   // One can print XML code to an
   // output directory:
   // meas.PrintXML( "xmlFromCCode", meas.GetOutputFilePrefix() );

   // Now, do the measurement
   MakeModelAndMeasurementFast( meas );
}

void create_test_input(){
   // taken from https://root.cern.ch/doc/v608/makeExample_8C_source.html
   TFile* example = new TFile("example.root","RECREATE");
   TH1F* data = new TH1F("data","data", 2,1,2);
   TH1F* signal = new TH1F("signal","signal histogram (pb)", 2,1,2);
   TH1F* background1 = new TH1F("background1","background 1 histogram (pb)", 2,1,2);
   TH1F* background2 = new TH1F("background2","background 2 histogram (pb)", 2,1,2);
   TH1F* statUncert = new TH1F("background1_statUncert", "statUncert", 2,1,2);

   // run with 1 pb
   data->SetBinContent(1,122);
   data->SetBinContent(2,112);

   signal->SetBinContent(1,20);
   signal->SetBinContent(2,10);

   background1->SetBinContent(1,100);
   background2->SetBinContent(2,100);

   // A small statistical uncertainty
   statUncert->SetBinContent(1, .05);  // 5% uncertainty
   statUncert->SetBinContent(2, .05);  // 5% uncertainty

   example->Write();
   example->Close();
}

void broken_weights(){
        create_test_input();
        create_test_ws();
        TFile f("results/example_UsingC_combined_meas_model.root","READ");
        RooWorkspace* w = (RooWorkspace*)f.Get("combined");
        RooStats::ModelConfig mc = *(static_cast<RooStats::ModelConfig*>(w -> obj("ModelConfig")));
        RooSimultaneous* simPdf = dynamic_cast<RooSimultaneous*>((mc.GetPdf()));
        RooAbsPdf *pdf = mc.GetPdf();
        RooArgSet obsSet = *(mc.GetObservables());
        RooDataSet* data = static_cast<RooDataSet*>(w->data("obsData"));
        RooDataSet* toyData = pdf->generate(obsSet, RooFit::Extended());
        RooCategory* channelCat = (RooCategory*) (&simPdf->indexCat());
        TIterator* iter = channelCat->typeIterator() ;
        RooCatType* tt = NULL;
        for(int i=0; i<1; ++i){
                data->get(i);
                toyData->get(i);
                // FIXME: weights are broken in this reduced dataset
                RooAbsData *datatmp = toyData->reduce(Form("channelCat==channelCat::channel1"));
                datatmp->get(i);
                std::cout << "data weight: " << data->weight() << " toyData weight: " << toyData->weight() << " reduced toyData weight: " << datatmp->weight() << std::endl;
        }
}

This results in

data weight: 122 toyData weight: 99 reduced toyData weight: 1

while weights not equal to 1 are expected.

Setup

ROOT Version: 6.20/06
Built for linuxx8664gcc on Jun 10 2020, 06:10:57
From tags/v6-20-06@v6-20-06

Additional context

@guitargeek
Copy link
Contributor

Hi @huebner-m, sorry for the late reply. I understood the problem now, and it will be fixed in the 6.28 development cycle by this commit:

#10397

guitargeek added a commit to guitargeek/root that referenced this issue Apr 14, 2022
@guitargeek guitargeek changed the title Broken weights after reducing RooDataSet created with RooAbsPdf::generate() [RF] Broken weights after reducing RooDataSet created with RooAbsPdf::generate() Jun 30, 2022
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
No open projects
Development

Successfully merging a pull request may close this issue.

4 participants