ana::SBNOsc::Covariance Class Reference

#include <Covariance.h>

Inheritance diagram for ana::SBNOsc::Covariance:

Classes
class	EventSample

Public Member Functions
	Covariance ()
void	FileCleanup (TTree *eventTree)
void	Initialize (fhicl::ParameterSet *config)
void	ProcessEvent (const event::Event *event)
void	ProcessSubRun (const SubRun *subrun)
void	Finalize ()
void	GetCovs ()
void	Write ()
void	Scale ()
TMatrixDSym	CovarianceMatrix ()
Public Member Functions inherited from core::PostProcessorBase
	PostProcessorBase ()
virtual	~PostProcessorBase ()
void	Run (std::vector< std::string > filelist)
void	Initialize (char *config=NULL, const std::string &output_fname="", unsigned n_workers=1)

Public Attributes
std::vector< TH2D >	cov
	Covariance Matrix per variation. More...
std::vector< TH2D >	fcov
	Fractional Covariance Matrix per variation. More...
std::vector< TH2D >	corr
	Correlation Matrix per variation. More...

Private Member Functions
void	GetCovPerVariation (unsigned variation)

Private Attributes
std::vector< std::vector < std::string > >	fWeightKeys
std::vector< std::vector < std::string > >	fWeightKeysCC
std::vector< std::vector < std::string > >	fWeightKeysNC
int	fNumAltUnis
int	fAltUniOffset
std::string	fEnergyType
unsigned	fNVariations
double	fSelectionEfficiency
double	fBackgroundRejection
std::string	fOutputFile
bool	fSaveCentralValue
bool	fSaveUniverses
double	fWeightMax
unsigned	fSampleIndex
std::vector< EventSample >	fEventSamples

Additional Inherited Members
Protected Member Functions inherited from core::PostProcessorBase
virtual void	ProcessFileMeta (const FileMeta *filemeta)
virtual void	FileSetup (TFile *f, TTree *eventTree)
virtual void	InitializeThread ()
unsigned	NWorkers ()
unsigned	WorkerID ()
Protected Attributes inherited from core::PostProcessorBase
ProviderManager *	fProviderManager
	Interface for provider access. More...
int	fConfigExperimentID

Detailed Description

Definition at line 30 of file Covariance.h.

Constructor & Destructor Documentation

ana::SBNOsc::Covariance::Covariance ( )

inline

Definition at line 33 of file Covariance.h.

{}

Member Function Documentation

TMatrixDSym ana::SBNOsc::Covariance::CovarianceMatrix

(

)

Definition at line 407 of file Covariance.cxx.

                                         {
    TMatrixDSym E_mat(cov[0].GetNbinsX());
    for (int i = 0; i < cov[0].GetNbinsX(); i++) {
        for (int j = 0; j < cov[0].GetNbinsY(); j++) {
            E_mat[i][j] = 0.; // zero-initialize
            for (TH2D &cov_histo: cov)  E_mat[i][j] += cov_histo.GetBinContent(i+1, j+1); // add all systematics
            // add statistical uncertainty
            if (i == j) {
                unsigned stat_ind = i;
                float stat_uncertainty;
                // get index into CV bin to get stat uncertainty
                for (auto const &sample: fEventSamples) {
                    if (stat_ind < sample.fCentralValue->GetNbinsX()) {
                        stat_uncertainty = sample.fCentralValue->GetBinContent(1+stat_ind);
                        break;
                    }
                    else stat_ind -= sample.fCentralValue->GetNbinsX();
                }
                E_mat[i][i] += stat_uncertainty;
            }
        }
    }
    return E_mat;
}

void ana::SBNOsc::Covariance::FileCleanup ( TTree *  eventTree )

virtual

Any cleanup needed per file

Parameters

eventTree

the TTree associated with the sbncode event.

Reimplemented from core::PostProcessorBase.

Definition at line 257 of file Covariance.cxx.

                                             {
    // onto the next sample
    fSampleIndex ++;
}

void ana::SBNOsc::Covariance::Finalize ( )

inlinevirtual

Perform user-level finalization. Called after all events have been processed.

Reimplemented from core::PostProcessorBase.

Definition at line 40 of file Covariance.h.

{ Scale(); GetCovs(); Write(); }

void ana::SBNOsc::Covariance::GetCovPerVariation ( unsigned  variation )

private

Definition at line 266 of file Covariance.cxx.

                                                      {
    
    //// Get covariances, fractional covariances and correlation coefficients
    //// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    std::cout << std::endl << "Getting covs..." << std::endl;

    // get the total number of bins across all samples
    unsigned num_bins = 0;
    for (auto const &sample: fEventSamples) {
       num_bins += sample.fCentralValue->GetNbinsX();
    }
    

    std::string covariance_name = "Covariance "  + std::to_string(variation);
    std::string fractional_covariance_name = "Fractional Covariance "  + std::to_string(variation);
    // Covariance and fractional covariance
    cov.emplace_back(covariance_name.c_str(), covariance_name.c_str(), num_bins, 0, num_bins, num_bins, 0, num_bins);
    fcov.emplace_back(fractional_covariance_name.c_str(), fractional_covariance_name.c_str(), num_bins, 0, num_bins, num_bins, 0, num_bins);

    // indexes into covariance
    unsigned cov_i = 0;
    unsigned cov_j = 0;
    // iterate over samples and bins
    for (auto const &sample_i: fEventSamples) {
        for (unsigned sample_i_bin_index = 0; sample_i_bin_index < sample_i.fCentralValue->GetNbinsX(); sample_i_bin_index++) {
            // Get Central Value
            float i_central = sample_i.fCentralValue->GetBinContent(sample_i_bin_index+1);
            for (auto const &sample_j: fEventSamples) {
                for (unsigned sample_j_bin_index = 0; sample_j_bin_index < sample_j.fCentralValue->GetNbinsX(); sample_j_bin_index++) {
                    // Get Central Value
                    float j_central = sample_j.fCentralValue->GetBinContent(sample_j_bin_index+1);
                    // calculate covariance
                    float cov_value = 0.; 
                    for (int u = 0; u < fNumAltUnis; u++) {
                        // get variations
                        float i_variation = sample_i.fUniverses[variation][u]->GetBinContent(sample_i_bin_index+1);
                        float j_variation = sample_j.fUniverses[variation][u]->GetBinContent(sample_j_bin_index+1);
                        cov_value += (i_central - i_variation) * (j_central - j_variation);
                    }
                    // average
                    if (fNumAltUnis != 0) cov_value /= fNumAltUnis;

                    // calculate fractional covariance
                    float fcov_value;
                    if (i_central * j_central < 1e-6) {
                        fcov_value = 0.;
                    }
                    else {
                        fcov_value = cov_value / (i_central * j_central);
                    }

                    // save values
                    cov[variation].SetBinContent(cov_i + 1, cov_j + 1, cov_value);
                    fcov[variation].SetBinContent(cov_i + 1, cov_j + 1, fcov_value);
                    // update covariance position
                    cov_j ++;
                }
            }   
            // update covariance position
            cov_i ++;
            // reset j
            cov_j = 0;
        }
    } 

    std::string correlation_name = "Correlation " + std::to_string(variation);

    // Pearson Correlation Coefficients
    corr.emplace_back(correlation_name.c_str(), correlation_name.c_str(), num_bins, 0, num_bins, num_bins, 0, num_bins);
    
    for (int i = 0; i < num_bins; i++) {
        float covii = cov[variation].GetBinContent(i+1, i+1);
        for (int j = 0; j < num_bins; j++) {
            float covjj = cov[variation].GetBinContent(j+1, j+1);
            float covij = cov[variation].GetBinContent(i+1, j+1);

            float corrij;
            // handle case where covariance is 0
            if (covii*covjj < 1e-6) {
                corrij = 0.;
            }
            else {
                corrij = covij / TMath::Sqrt(covii*covjj);
            }
            corr[variation].SetBinContent(i+1, j+1, corrij);
        }
    }
    
    // Add bin labels
    std::vector <TH2D*> hists = {&cov[variation], &fcov[variation], &corr[variation]};
    unsigned bin = 0;
    for (auto const &sample: fEventSamples) {
        // Set label
        for (TH2D* hist : hists) {
            hist->GetXaxis()->SetBinLabel(1+bin, sample.fName.c_str());
            hist->GetYaxis()->SetBinLabel(1+bin, sample.fName.c_str());
        }
        // update bin positions
        bin += sample.fBins.size() - 1;

    }

    for (TH2D* hist : hists) { hist->GetXaxis()->LabelsOption("h"); hist->GetYaxis()->LabelsOption("v"); }
    
    std::cout << std::endl << "  Got covs." << std::endl;
    
}

void ana::SBNOsc::Covariance::GetCovs

(

)

Definition at line 262 of file Covariance.cxx.

                         {
    for (unsigned i = 0; i < fNVariations; i++) GetCovPerVariation(i);
}

void ana::SBNOsc::Covariance::Initialize ( fhicl::ParameterSet *  config )

virtual

Perform user-level initialization.

Parameters

config

A configuration, as a JSON object.

Implements core::PostProcessorBase.

Definition at line 76 of file Covariance.cxx.

                                                     {
    // must have config
    assert(config != NULL);
    // Weight and universe stuff
    // WeightKey can either be the name of a single weight, or a name signifying a list of weights,
    // or a list of weights to be oscillated
    
    fhicl::ParameterSet pconfig = config->get<fhicl::ParameterSet>("Covariance");

    fWeightKeys = pconfig.get<std::vector<std::vector<std::string>>>("WeightKey");
    fWeightKeysCC = pconfig.get<std::vector<std::vector<std::string>>>("WeightKeyCC", {});
    fWeightKeysNC = pconfig.get<std::vector<std::vector<std::string>>>("WeightKeyNC", {});
    fNVariations = fWeightKeys.size();

    // all should have the same size
    assert(fWeightKeysCC.size() == 0 || fWeightKeysCC.size() == fWeightKeys.size());
    assert(fWeightKeysNC.size() == 0 || fWeightKeysNC.size() == fWeightKeys.size());

    // merge weight keys into CC/NC keys
    if (fWeightKeysCC.size() == 0) {
      fWeightKeysCC = fWeightKeys;
    }
    else {
      for (unsigned i = 0; i < fNVariations; i++) {
        fWeightKeysCC[i].insert(fWeightKeysCC[i].end(), fWeightKeys[i].begin(), fWeightKeys[i].end());
      }
    }
    if (fWeightKeysNC.size() == 0) {
      fWeightKeysNC = fWeightKeys;
    }
    else {
      for (unsigned i = 0; i < fNVariations; i++) {
        fWeightKeysNC[i].insert(fWeightKeysNC[i].end(), fWeightKeys[i].begin(), fWeightKeys[i].end());
      }
    }

    // number of universes to be used
    fNumAltUnis = pconfig.get<int>("NumAltUnis", 0);
    // and offset into the eventweight vector
    fAltUniOffset = pconfig.get<int>("AltUniOffset", 0);
    
    // Type of energy
    fEnergyType = pconfig.get<std::string>("EnergyType", "Reco");
    
    // Further selection and rejection 'efficiencies'
    fSelectionEfficiency = pconfig.get<double>("SelectionEfficiency", 1.0);
    fBackgroundRejection = pconfig.get<double>("BackgroundRejection", 0.0);

    // maximum weight for covariance calculation
    // If an event has any weight larger than this value, it will be 
    // thrown away for the purposes of covariance construction
    fWeightMax = pconfig.get<double>("WeightMax", -1.);
     
    // get event samples
    std::vector<fhicl::ParameterSet> configEventSamples = \
      config->get<std::vector<fhicl::ParameterSet> >("EventSamples");
    for (const fhicl::ParameterSet& sample : configEventSamples) {
      fEventSamples.emplace_back(sample, fNumAltUnis, fNVariations);
    }
    
    // set output directory
    fOutputFile = pconfig.get<std::string>("OutputFile", "");

    // whether to save things
    fSaveCentralValue = pconfig.get<bool>("SaveCentralValue", false);
    fSaveUniverses = pconfig.get<bool>("SaveUniverses", false);

    // start out at the zeroth sample
    fSampleIndex = 0;
}

void ana::SBNOsc::Covariance::ProcessEvent ( const event::Event *  event )

virtual

Process one event.

Parameters

event

The sbncode event for the current event

Implements core::PostProcessorBase.

Definition at line 151 of file Covariance.cxx.

                                                     {
    // iterate over each interaction in the event
    for (int n = 0; n < event->reco.size(); n++) {
        unsigned truth_ind = event->reco[n].truth_index;

        // Get energy
        int idx = event->reco[n].truth_index;
        if (idx < 0 || idx > event->truth.size()) {
          continue;
        }
        double true_nuE = event->truth[idx].neutrino.energy;
        double nuE; 
        if (fEnergyType == "CCQE") {
            nuE = event->truth[truth_ind].neutrino.eccqe;
        } else if (fEnergyType == "True") {
            nuE = true_nuE;
        } else if (fEnergyType == "Reco") {
            nuE = event->reco[n].reco_energy;
        }
        nuE *= fEventSamples[fSampleIndex].fScaleEnergy;
    
        // Apply selection (or rejection) efficiencies
        int isCC = event->truth[truth_ind].neutrino.iscc;
        // double wgt = isCC*(fSelectionEfficiency) + (1-isCC)*(1 - fBackgroundRejection);
        // apply scaling from fScaleFactor
        // wgt *= fEventSamples[fSampleIndex].fScaleFactor;
        // apply uniform weights
        // for (auto const &key: fUniformWeights) {
        //     wgt *= event->truth[truth_ind].weights.at(key)[0];
        // }
        double wgt = event->reco[n].weight;
        // apply POT scaling if configured
        if (fEventSamples[fSampleIndex].fScalePOT > 0) {
          wgt *= fEventSamples[fSampleIndex].fScalePOT / fEventSamples[fSampleIndex].fPOT;
        }

        if (event->truth[truth_ind].neutrino.iscc) {
          wgt *= fEventSamples[fSampleIndex].fScaleCC;
        }
        else {
          wgt *= fEventSamples[fSampleIndex].fScaleNC;
        }
    
        std::vector<std::vector <float>> uweights; 
        // Get weights for each alternative universe
        if (isCC) {
          for (std::vector<std::string> &weight_keys: fWeightKeysCC) {
            uweights.push_back(
              GetUniWeights(event->truth[truth_ind].weights, weight_keys, fNumAltUnis, fAltUniOffset)
            );
          }
        }
        else {
          for (std::vector<std::string> &weight_keys: fWeightKeysNC) {
            uweights.push_back(
              GetUniWeights(event->truth[truth_ind].weights, weight_keys, fNumAltUnis, fAltUniOffset)
            );
          }
        }

        // see if weight is too big
        if (fWeightMax > 0) {
          double max_weight = 0;
          for (std::vector<float> &weights: uweights) {
            float this_max_weight = *std::max_element(weights.begin(), weights.end());
            if (this_max_weight > max_weight) max_weight = this_max_weight;
          }
          if (max_weight > fWeightMax) {
            std::cout << "Weight (" << max_weight << ") over cap (" << fWeightMax << 
              ") for event with energy (" << nuE << ") in sample " << fEventSamples[fSampleIndex].fName << std::endl;
            continue;
          }
        }
    
        // Add to base and bkg count histogram
        fEventSamples[fSampleIndex].fCentralValue->Fill(nuE, wgt);
    
        // Fill alternative universe histograms
        for (unsigned variation = 0; variation < uweights.size(); variation++) {
          for (int u = 0; u < uweights[variation].size(); u++) {
              fEventSamples[fSampleIndex].fUniverses[variation][u]->Fill(nuE, wgt*uweights[variation][u]);
          }
        }
    }
}

void ana::SBNOsc::Covariance::ProcessSubRun ( const SubRun *  subrun )

virtual

Reimplemented from core::PostProcessorBase.

Definition at line 147 of file Covariance.cxx.

                                                   {
  fEventSamples[fSampleIndex].fPOT += subrun->totgoodpot;
}

void ana::SBNOsc::Covariance::Scale

(

)

Definition at line 237 of file Covariance.cxx.

                       {
  for (int sample_i = 0; sample_i < fEventSamples.size(); sample_i++) {
    for (int i = 1; i < fEventSamples[sample_i].fCentralValue->GetNbinsX()+1; i++) {
      //std::cout << "pre sample: " << sample_i << " bin content: " << fEventSamples[sample_i].fCentralValue->GetBinContent(i);
      fEventSamples[sample_i].fCentralValue->SetBinContent(i,
        fEventSamples[sample_i].fEnergyBinScale[i-1] *  fEventSamples[sample_i].fCentralValue->GetBinContent(i));
      //std::cout << "post sample: " << sample_i << " bin content: " << fEventSamples[sample_i].fCentralValue->GetBinContent(i);
    }
    for (int variation = 0; variation < fNVariations; variation++) {
      for (int uni = 0; uni < fNumAltUnis; uni++) {
        for (int i = 1; i < fEventSamples[sample_i].fUniverses[variation][uni]->GetNbinsX(); i++) {
          double content = fEventSamples[sample_i].fUniverses[variation][uni]->GetBinContent(i);
          fEventSamples[sample_i].fUniverses[variation][uni]->SetBinContent(i, content * fEventSamples[sample_i].fEnergyBinScale[i-1]);
        }
      }
    }
  }
}

void ana::SBNOsc::Covariance::Write

(

)

Definition at line 375 of file Covariance.cxx.

                       {
    // if output file not set, don't write
    if (fOutputFile == "") return;

    // Write Covariances
    TFile *output = TFile::Open(fOutputFile.c_str(), "recreate");
    assert(output && output->IsOpen());
    
    for (TH2D &h: cov) h.Write();
    for (TH2D &h: fcov) h.Write();
    for (TH2D &h: corr) h.Write();

    // write histos
    if (fSaveCentralValue) {
      for (auto const &sample: fEventSamples) {
        sample.fCentralValue->Write();
      }
    }
    if (fSaveUniverses) {
        for (auto const &sample: fEventSamples) {
            for (unsigned i =0; i < fNVariations; i++) {
              for (TH1D *h: sample.fUniverses[i]) {
                  h->Write();
              }
            }
        }
    }
    // close file
    output->Close(); 
}

Member Data Documentation

std::vector<TH2D> ana::SBNOsc::Covariance::corr

Correlation Matrix per variation.

Definition at line 53 of file Covariance.h.

std::vector<TH2D> ana::SBNOsc::Covariance::cov

Covariance Matrix per variation.

Definition at line 51 of file Covariance.h.

int ana::SBNOsc::Covariance::fAltUniOffset

private

Definition at line 81 of file Covariance.h.

double ana::SBNOsc::Covariance::fBackgroundRejection

private

Definition at line 88 of file Covariance.h.

std::vector<TH2D> ana::SBNOsc::Covariance::fcov

Fractional Covariance Matrix per variation.

Definition at line 52 of file Covariance.h.

std::string ana::SBNOsc::Covariance::fEnergyType

private

Definition at line 82 of file Covariance.h.

std::vector<EventSample> ana::SBNOsc::Covariance::fEventSamples

private

Definition at line 101 of file Covariance.h.

int ana::SBNOsc::Covariance::fNumAltUnis

private

Definition at line 80 of file Covariance.h.

unsigned ana::SBNOsc::Covariance::fNVariations

private

Definition at line 85 of file Covariance.h.

std::string ana::SBNOsc::Covariance::fOutputFile

private

Definition at line 89 of file Covariance.h.

unsigned ana::SBNOsc::Covariance::fSampleIndex

private

Definition at line 98 of file Covariance.h.

bool ana::SBNOsc::Covariance::fSaveCentralValue

private

Definition at line 91 of file Covariance.h.

bool ana::SBNOsc::Covariance::fSaveUniverses

private

Definition at line 92 of file Covariance.h.

double ana::SBNOsc::Covariance::fSelectionEfficiency

private

Definition at line 87 of file Covariance.h.

std::vector<std::vector<std::string> > ana::SBNOsc::Covariance::fWeightKeys

private

Definition at line 77 of file Covariance.h.

std::vector<std::vector<std::string> > ana::SBNOsc::Covariance::fWeightKeysCC

private

Definition at line 78 of file Covariance.h.

std::vector<std::vector<std::string> > ana::SBNOsc::Covariance::fWeightKeysNC

private

Definition at line 79 of file Covariance.h.

double ana::SBNOsc::Covariance::fWeightMax

private

Definition at line 95 of file Covariance.h.

---

The documentation for this class was generated from the following files:

• Covariance.h
• Covariance.cxx