Skip to content

Commit

Permalink
first draft at implementing weighted data
Browse files Browse the repository at this point in the history
  • Loading branch information
@cburgard
cburgard committed Jun 13, 2022
1 parent fd14c24 commit af44b7b
Show file tree
Hide file tree
Showing 2 changed files with 72 additions and 36 deletions.
2 changes: 2 additions & 0 deletions roofit/histfactory/inc/RooStats/HistFactory/HistoToWorkspaceFactoryFast.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ namespace RooStats{

struct Configuration {
bool binnedFitOptimization = true;
bool storeDataError = false;
};

HistoToWorkspaceFactoryFast() {}
Expand All @@ -60,6 +61,7 @@ namespace RooStats{
RooWorkspace* MakeSingleChannelModel( Measurement& measurement, Channel& channel );
RooWorkspace* MakeCombinedModel(std::vector<std::string>, std::vector<std::unique_ptr<RooWorkspace>>&);

static RooWorkspace* MakeCombinedModel( Measurement& measurement, const Configuration& config);
static RooWorkspace* MakeCombinedModel( Measurement& measurement );
static void PrintCovarianceMatrix(RooFitResult* result, RooArgSet* params,
std::string filename);
Expand Down
106 changes: 70 additions & 36 deletions roofit/histfactory/src/HistoToWorkspaceFactoryFast.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -240,6 +240,23 @@ namespace HistFactory{
}

RooWorkspace* HistoToWorkspaceFactoryFast::MakeCombinedModel( Measurement& measurement ) {

// This function takes a fully configured measurement
// which may contain several channels and returns
// a workspace holding the combined model
//
// This can be used, for example, within a script to produce
// a combined workspace on-the-fly
//
// This is a static function (for now) to make
// it a one-liner


Configuration config;
return MakeCombinedModel(measurement,config);
}

RooWorkspace* HistoToWorkspaceFactoryFast::MakeCombinedModel( Measurement& measurement , const Configuration& config) {

// This function takes a fully configured measurement
// which may contain several channels and returns
Expand All @@ -251,10 +268,10 @@ namespace HistFactory{
// This is a static function (for now) to make
// it a one-liner

RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::INFO, 0u, RooFit::ObjectHandling, false);
RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::INFO, 0u, RooFit::ObjectHandling, false);

// First, we create an instance of a HistFactory
HistoToWorkspaceFactoryFast factory( measurement );
HistoToWorkspaceFactoryFast factory( measurement , config);

// Loop over the channels and create the individual workspaces
vector<std::unique_ptr<RooWorkspace>> channel_workspaces;
Expand Down Expand Up @@ -1680,6 +1697,9 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
// THis works and is natural, but the memory size of the simultaneous dataset grows exponentially with channels
const char* weightName="weightVar";
proto->factory(TString::Format("%s[0,-1e10,1e10]",weightName));
const char* weightErrName="weightErr";
proto->factory(TString::Format("%s[0,-1e10,1e10]",weightErrName));

proto->defineSet("obsAndWeight",TString::Format("%s,%s",weightName,observablesStr.c_str()));

// New Asimov Generation: Use the code in the Asymptotic calculator
Expand All @@ -1695,9 +1715,15 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
if(TH1 const* mnominal = channel.GetData().GetHisto()) {
// This works and is natural, but the memory size of the simultaneous
// dataset grows exponentially with channels.
RooDataSet dataset{"obsData","",*proto->set("obsAndWeight"),weightName};
ConfigureHistFactoryDataset( dataset, *mnominal, *proto, fObsNameVec );
proto->import(dataset);
if(!fCfg.storeDataError){
RooDataSet dataset{"obsData","",*proto->set("obsAndWeight"),weightName};
ConfigureHistFactoryDataset( dataset, *mnominal, *proto, fObsNameVec );
proto->import(dataset);
} else {
RooDataSet dataset{"obsData","",*proto->set("obsAndWeight"), RooFit::WeightVar(weightName), RooFit::StoreError(*proto->var(weightErrName))};
ConfigureHistFactoryDataset( dataset, *mnominal, *proto, fObsNameVec );
proto->import(dataset);
}
} // End: Has non-null 'data' entry


Expand Down Expand Up @@ -1744,39 +1770,47 @@ RooArgList HistoToWorkspaceFactoryFast::createObservables(const TH1 *hist, RooWo
return;
}

TAxis const* ax = mnominal.GetXaxis();
TAxis const* ay = mnominal.GetYaxis();
TAxis const* az = mnominal.GetZaxis();

for (int i=1; i<=ax->GetNbins(); ++i) { // 1 or more dimension

double xval = ax->GetBinCenter(i);
proto.var( obsNameVec[0] )->setVal( xval );

if(obsNameVec.size()==1) {
double fval = mnominal.GetBinContent(i);
obsDataUnbinned.add( *proto.set("obsAndWeight"), fval );
} else { // 2 or more dimensions

for(int j=1; j<=ay->GetNbins(); ++j) {
double yval = ay->GetBinCenter(j);
proto.var( obsNameVec[1] )->setVal( yval );

if(obsNameVec.size()==2) {
double fval = mnominal.GetBinContent(i,j);
obsDataUnbinned.add( *proto.set("obsAndWeight"), fval );
} else { // 3 dimensions

for(int k=1; k<=az->GetNbins(); ++k) {
double zval = az->GetBinCenter(k);
proto.var( obsNameVec[2] )->setVal( zval );
double fval = mnominal.GetBinContent(i,j,k);
obsDataUnbinned.add( *proto.set("obsAndWeight"), fval );
TAxis const* ax = mnominal.GetXaxis();
TAxis const* ay = mnominal.GetYaxis();
TAxis const* az = mnominal.GetZaxis();

RooRealVar* err = proto.var("weightErr");

for (int i=1; i<=ax->GetNbins(); ++i) { // 1 or more dimension

double xval = ax->GetBinCenter(i);
proto.var( obsNameVec[0] )->setVal( xval );

if(obsNameVec.size()==1) {
double fval = mnominal.GetBinContent(i);
double ferr = mnominal.GetBinError(i);
err->setVal(ferr);
obsDataUnbinned.add( *proto.set("obsAndWeight"), fval );
} else { // 2 or more dimensions

for(int j=1; j<=ay->GetNbins(); ++j) {
double yval = ay->GetBinCenter(j);
proto.var( obsNameVec[1] )->setVal( yval );

if(obsNameVec.size()==2) {
double fval = mnominal.GetBinContent(i,j);
double ferr = mnominal.GetBinError(i,j);
err->setVal(ferr);
obsDataUnbinned.add( *proto.set("obsAndWeight"), fval );
} else { // 3 dimensions

for(int k=1; k<=az->GetNbins(); ++k) {
double zval = az->GetBinCenter(k);
proto.var( obsNameVec[2] )->setVal( zval );
double fval = mnominal.GetBinContent(i,j,k);
double ferr = mnominal.GetBinError(i,j,k);
err->setVal(ferr);
obsDataUnbinned.add( *proto.set("obsAndWeight"), fval );
}
}
}
}
}
}
}
}
}

void HistoToWorkspaceFactoryFast::GuessObsNameVec(const TH1* hist)
Expand Down

0 comments on commit af44b7b

Please sign in to comment.