Skip to content

Commit

Permalink
Add a "true2d" histogram to the Universe class. This allows for correct
Browse files Browse the repository at this point in the history 
calculations of MC statistical covariances between true bin contents.
  • Loading branch information
@sjgardiner
sjgardiner committed Mar 3, 2024
1 parent f8b90ae commit eb18261
Show file tree
Hide file tree
Showing 3 changed files with 54 additions and 9 deletions.
29 changes: 27 additions & 2 deletions SystematicsCalculator.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,9 @@ void set_stats_and_dir( Universe& univ ) {

univ.hist_reco2d_->SetStats( false );
univ.hist_reco2d_->SetDirectory( nullptr );

univ.hist_true2d_->SetStats( false );
univ.hist_true2d_->SetDirectory( nullptr );
}

// Tests whether a string ends with another string. Taken from
Expand Down Expand Up @@ -506,21 +509,25 @@ void SystematicsCalculator::load_universes( TDirectoryFile& total_subdir ) {
TH2D* hist_2d = nullptr;
TH2D* hist_categ = nullptr;
TH2D* hist_reco2d = nullptr;
TH2D* hist_true2d = nullptr;

total_subdir.GetObject( (key + "_true").c_str(), hist_true );
total_subdir.GetObject( (key + "_reco").c_str(), hist_reco );
total_subdir.GetObject( (key + "_2d").c_str(), hist_2d );
total_subdir.GetObject( (key + "_categ").c_str(), hist_categ );
total_subdir.GetObject( (key + "_reco2d").c_str(), hist_reco2d );
total_subdir.GetObject( (key + "_true2d").c_str(), hist_true2d );

if ( !hist_true || !hist_reco || !hist_2d || !hist_categ || !hist_reco2d ) {
if ( !hist_true || !hist_reco || !hist_2d || !hist_categ
|| !hist_reco2d || !hist_true2d )
{
throw std::runtime_error( "Failed to retrieve histograms for the "
+ key + " universe" );
}

// Reconstruct the Universe object from the retrieved histograms
auto temp_univ = std::make_unique< Universe >( univ_name, univ_index,
hist_true, hist_reco, hist_2d, hist_categ, hist_reco2d );
hist_true, hist_reco, hist_2d, hist_categ, hist_reco2d, hist_true2d );

// Determine whether the current universe represents a detector
// variation or a reweightable variation. We'll use this information to
Expand Down Expand Up @@ -883,13 +890,17 @@ void SystematicsCalculator::build_universes( TDirectoryFile& root_tdir ) {
auto h_reco2d = get_object_unique_ptr< TH2D >(
(hist_name_prefix + "_reco2d"), *subdir );

auto h_true2d = get_object_unique_ptr< TH2D >(
(hist_name_prefix + "_true2d"), *subdir );

// Add their contributions to the owned histograms for the
// current Universe object
fake_data_universe_->hist_reco_->Add( h_reco.get() );
fake_data_universe_->hist_true_->Add( h_true.get() );
fake_data_universe_->hist_2d_->Add( h_2d.get() );
fake_data_universe_->hist_categ_->Add( h_categ.get() );
fake_data_universe_->hist_reco2d_->Add( h_reco2d.get() );
fake_data_universe_->hist_true2d_->Add( h_true2d.get() );

} // fake data sample

Expand Down Expand Up @@ -947,6 +958,9 @@ void SystematicsCalculator::build_universes( TDirectoryFile& root_tdir ) {
auto hist_reco2d = get_object_unique_ptr< TH2D >(
"unweighted_0_reco2d", *subdir );

auto hist_true2d = get_object_unique_ptr< TH2D >(
"unweighted_0_true2d", *subdir );

double temp_scale_factor = 1.;
if ( is_altCV ) {
// AltCV ntuple files are available for all runs, so scale
Expand All @@ -970,6 +984,7 @@ void SystematicsCalculator::build_universes( TDirectoryFile& root_tdir ) {
hist_2d->Scale( temp_scale_factor );
hist_categ->Scale( temp_scale_factor );
hist_reco2d->Scale( temp_scale_factor );
hist_true2d->Scale( temp_scale_factor );

// Add the scaled contents of these histograms to the
// corresponding histograms in the new Universe object
Expand All @@ -978,6 +993,7 @@ void SystematicsCalculator::build_universes( TDirectoryFile& root_tdir ) {
temp_univ_ptr->hist_2d_->Add( hist_2d.get() );
temp_univ_ptr->hist_categ_->Add( hist_categ.get() );
temp_univ_ptr->hist_reco2d_->Add( hist_reco2d.get() );
temp_univ_ptr->hist_true2d_->Add( hist_true2d.get() );

// Adjust the owned histograms to avoid auto-deletion problems
set_stats_and_dir( *temp_univ_ptr );
Expand Down Expand Up @@ -1095,13 +1111,17 @@ void SystematicsCalculator::build_universes( TDirectoryFile& root_tdir ) {
auto h_reco2d = get_object_unique_ptr< TH2D >(
(hist_name_prefix + "_reco2d"), *subdir );

auto h_true2d = get_object_unique_ptr< TH2D >(
(hist_name_prefix + "_true2d"), *subdir );

// Scale these histograms to the appropriate BNB data POT for
// the current run
h_reco->Scale( rw_scale_factor );
h_true->Scale( rw_scale_factor );
h_2d->Scale( rw_scale_factor );
h_categ->Scale( rw_scale_factor );
h_reco2d->Scale( rw_scale_factor );
h_true2d->Scale( rw_scale_factor );

// Add their contributions to the owned histograms for the
// current Universe object
Expand All @@ -1110,6 +1130,7 @@ void SystematicsCalculator::build_universes( TDirectoryFile& root_tdir ) {
universe.hist_2d_->Add( h_2d.get() );
universe.hist_categ_->Add( h_categ.get() );
universe.hist_reco2d_->Add( h_reco2d.get() );
universe.hist_true2d_->Add( h_true2d.get() );

} // universes indices

Expand Down Expand Up @@ -1205,6 +1226,7 @@ void SystematicsCalculator::save_universes( TDirectoryFile& out_tdf ) {
universe->hist_2d_->Write();
universe->hist_categ_->Write();
universe->hist_reco2d_->Write();
universe->hist_true2d_->Write();
}

// Save the alternate CV MC histograms
Expand All @@ -1217,6 +1239,7 @@ void SystematicsCalculator::save_universes( TDirectoryFile& out_tdf ) {
universe->hist_2d_->Write();
universe->hist_categ_->Write();
universe->hist_reco2d_->Write();
universe->hist_true2d_->Write();
}

// Save the reweightable systematic histograms
Expand All @@ -1229,6 +1252,7 @@ void SystematicsCalculator::save_universes( TDirectoryFile& out_tdf ) {
universe->hist_2d_->Write();
universe->hist_categ_->Write();
universe->hist_reco2d_->Write();
universe->hist_true2d_->Write();
}

}
Expand All @@ -1240,6 +1264,7 @@ void SystematicsCalculator::save_universes( TDirectoryFile& out_tdf ) {
fake_data_universe_->hist_2d_->Write();
fake_data_universe_->hist_categ_->Write();
fake_data_universe_->hist_reco2d_->Write();
fake_data_universe_->hist_true2d_->Write();
}

// Save the total BNB data POT for easy retrieval later
Expand Down
7 changes: 2 additions & 5 deletions TruthSystematicsCalculator.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,15 +50,12 @@ double TruthSystematicsCalculator::evaluate_observable( const Universe& univ,
double TruthSystematicsCalculator::evaluate_mc_stat_covariance(
const Universe& univ, int true_bin_a, int true_bin_b ) const
{
// TODO: revisit to get off-diagonal MC stat covariances
if ( true_bin_a != true_bin_b ) return 0.;

// ROOT histograms use one-based bin indices, so I correct for that here.
// Note that using the bin error (rather than the bin contents) enables a
// correct treatment for weighted events provided TH1::Sumw2() was called
// before filling the histogram.
double err = univ.hist_true_->GetBinError( true_bin_a + 1 );
double err2 = err * err;
double err = univ.hist_true2d_->GetBinError( true_bin_a + 1, true_bin_b + 1 );
double err2 = err*err;
return err2;
}

Expand Down
27 changes: 25 additions & 2 deletions UniverseMaker.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -265,6 +265,11 @@ class Universe {
"; reco bin number; reco bin number; counts", num_reco_bins, 0.,
num_reco_bins, num_reco_bins, 0., num_reco_bins );

hist_true2d_ = std::make_unique< TH2D >(
(hist_name_prefix + "_true2d").c_str(),
"; true bin number; true bin number; counts", num_true_bins, 0.,
num_true_bins, num_true_bins, 0., num_true_bins );

// Get the number of defined EventCategory values by checking the number
// of elements in the "label map" managed by the EventCategoryInterpreter
// singleton class
Expand All @@ -283,22 +288,25 @@ class Universe {
hist_2d_->Sumw2();
hist_categ_->Sumw2();
hist_reco2d_->Sumw2();
hist_true2d_->Sumw2();
}

// Note: the new Universe object takes ownership of the histogram
// pointers passed to this constructor
Universe( const std::string& universe_name,
size_t universe_index, TH1D* hist_true, TH1D* hist_reco, TH2D* hist_2d,
TH2D* hist_categ, TH2D* hist_reco2d )
TH2D* hist_categ, TH2D* hist_reco2d, TH2D* hist_true2d )
: universe_name_( universe_name ), index_( universe_index ),
hist_true_( hist_true ), hist_reco_( hist_reco ), hist_2d_( hist_2d ),
hist_categ_( hist_categ ), hist_reco2d_( hist_reco2d )
hist_categ_( hist_categ ), hist_reco2d_( hist_reco2d ),
hist_true2d_( hist_true2d )
{
hist_true_->SetDirectory( nullptr );
hist_reco_->SetDirectory( nullptr );
hist_2d_->SetDirectory( nullptr );
hist_categ_->SetDirectory( nullptr );
hist_reco2d_->SetDirectory( nullptr );
hist_true2d_->SetDirectory( nullptr );
}

std::unique_ptr< Universe > clone() const {
Expand All @@ -312,6 +320,7 @@ class Universe {
result->hist_2d_->Add( this->hist_2d_.get() );
result->hist_categ_->Add( this->hist_categ_.get() );
result->hist_reco2d_->Add( this->hist_reco2d_.get() );
result->hist_true2d_->Add( this->hist_true2d_.get() );

return result;
}
Expand All @@ -323,6 +332,7 @@ class Universe {
std::unique_ptr< TH2D > hist_2d_;
std::unique_ptr< TH2D > hist_categ_;
std::unique_ptr< TH2D > hist_reco2d_;
std::unique_ptr< TH2D > hist_true2d_;
};

class UniverseMaker {
Expand Down Expand Up @@ -682,6 +692,12 @@ void UniverseMaker::build_universes(
universe.hist_2d_->Fill( tb.bin_index_, rb.bin_index_,
tb.weight_ * rb.weight_ * safe_wgt );
} // reco bins

for ( const auto& other_tb : matched_true_bins ) {
universe.hist_true2d_->Fill( tb.bin_index_, other_tb.bin_index_,
tb.weight_ * other_tb.weight_ * safe_wgt );
} // true bins

} // true bins

for ( const auto& rb : matched_reco_bins ) {
Expand Down Expand Up @@ -711,6 +727,12 @@ void UniverseMaker::build_universes(
univ.hist_2d_->Fill( tb.bin_index_, rb.bin_index_,
tb.weight_ * rb.weight_ );
} // reco bins

for ( const auto& other_tb : matched_true_bins ) {
univ.hist_true2d_->Fill( tb.bin_index_, other_tb.bin_index_,
tb.weight_ * other_tb.weight_ );
} // true bins

} // true bins

for ( const auto& rb : matched_reco_bins ) {
Expand Down Expand Up @@ -871,6 +893,7 @@ void UniverseMaker::save_histograms(
univ.hist_true_->Write();
univ.hist_2d_->Write();
univ.hist_categ_->Write();
univ.hist_true2d_->Write();
}
} // universes
} // weight names
Expand Down

0 comments on commit eb18261

Please sign in to comment.