FluxMultiverseSyst.cxx

Go to the documentation of this file.

#include "CAFAna/XSec/FluxMultiverseSyst.h"

#include <iostream>
#include <string>

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Systs/Systs.h"

#include "CAFAna/Core/LoadFromFile.h"

#include "TH1.h"
#include "TObjString.h"

using namespace std;

namespace ana
{
  //////
  FluxMultiverseSyst::FluxMultiverseSyst(SpectrumLoaderBase& loader,
                                         const HistAxis&     axis,
                                         const Cut&          cut,
                                         const SystShifts&   shift,
                                         const Var&          cv_wei,
                                         const std::vector<Var>& vuniv_wei):
    fNUniverses(0),fUpperSigma(NULL),fLowerSigma(NULL),fBoundariesFound(false),
    fNormUpperSigma(NULL), fNormLowerSigma(NULL){
    
    for(unsigned int i = 0; i < vuniv_wei.size(); i++){
      fSpectra.push_back( new Spectrum(loader, axis, cut, shift, vuniv_wei[i]));
    }
    fNUniverses = vuniv_wei.size();
    // last position is the cv
    fSpectra.push_back( new Spectrum(loader, axis, cut, shift, cv_wei) );
  }
    
  //////
  FluxMultiverseSyst::FluxMultiverseSyst(unsigned int nuniverses, 
                                       std::vector<std::unique_ptr<ana::Spectrum>>& spectra, bool fromfile): fLoadFromFile(fromfile)
  {
    fNUniverses = nuniverses;
    for(auto& it: spectra){
      fSpectra.push_back(it.release());
    }
    
  }

  FluxMultiverseSyst::FluxMultiverseSyst(std::vector<std::unique_ptr<ana::Spectrum>>& spectra){
    fNUniverses = spectra.size()-1;
     for(auto& it: spectra){
       fSpectra.push_back(it.release());
     }
  }
  
  FluxMultiverseSyst::~FluxMultiverseSyst(){
    if(!fLoadFromFile){
      for (auto it = fSpectra.begin(); it != fSpectra.end(); ++it) delete *it;
    }
    delete fUpperSigma;
    delete fLowerSigma;
    delete fNormUpperSigma;
    delete fNormLowerSigma;
  }
  
  Spectrum* FluxMultiverseSyst::CV(){
    return fSpectra[fNUniverses];
  }
  
  Spectrum* FluxMultiverseSyst::NormUpperSigma()
  {
    if(!fBoundariesFound) FindSigmaBoundaries();
    return fNormUpperSigma;
  }
  
  Spectrum* FluxMultiverseSyst::NormLowerSigma()
  {
    if(!fBoundariesFound) FindSigmaBoundaries();
    return fNormLowerSigma;
  }
  
  TGraphAsymmErrors* FluxMultiverseSyst::ToAreaNormalizedTH1(double pot, int col, int errCol,
                                                             float headroom, bool newaxis){
    if(col == -1){
      col = kRed;
      errCol = kRed-10;
    }
    else if(errCol == -1) errCol = col-7; // hopefully a lighter version
    
    TH1* cv =  CV()->ToTH1(pot);
    
    std::vector<TH1*> ups, dns;
    ups.push_back(fUpperSigma->ToTH1(pot));
    dns.push_back(fLowerSigma->ToTH1(pot));
    
    cv->SetLineColor(col);
    cv->GetXaxis()->CenterTitle();
    cv->GetYaxis()->CenterTitle();
    
    double yMax = cv->GetBinContent(cv->GetMaximumBin());
    
    TGraphAsymmErrors* g = new TGraphAsymmErrors;

    for(int binIdx = 0; binIdx < cv->GetNbinsX()+2; ++binIdx){
      const double y = cv->GetBinContent(binIdx);
      g->SetPoint(binIdx, cv->GetXaxis()->GetBinCenter(binIdx), y);

      const double w = cv->GetXaxis()->GetBinWidth(binIdx);

      double errUp = 0;
      double errDn = 0;

      for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
        double hi = ups[systIdx]->GetBinContent(binIdx)-y;
        double lo = y-dns[systIdx]->GetBinContent(binIdx);

        if(lo <= 0 && hi <= 0) std::swap(lo, hi);

        errUp += hi*hi;
        errDn += lo*lo;

        // TODO: what happens if they're both high or both low?
      } // end for systIdx


      g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
    } // end for i

    g->SetFillColor(errCol);
    //~ g->Draw("e2 same");
    g->GetYaxis()->SetRangeUser(0, headroom*yMax);
    cv->GetYaxis()->SetRangeUser(0, headroom*yMax);

    //~ nom->Draw("hist ][ same");
    return g;
  }
  
  TH1* FluxMultiverseSyst::ToTH1()
  {
    return ToTH1(CV()->POT());
  }
  
  // return a TH1* with Sigma of all universes as the error band
  TH1* FluxMultiverseSyst::ToTH1(double pot){
    TH1* hCV = CV()->ToTH1(pot);
    vector<TH1*> hUniverses;
    map< int,float > cv_counts_in_bins;
    map< int, vector<float> > univ_counts_in_bins;

    unsigned int count = 0;
    for(auto& it: fSpectra){
      
      if(count!=fNUniverses){
        TH1* curhist = it->ToTH1(pot);
        hUniverses.push_back(curhist);
        int nbinsx = curhist->GetNbinsX();
        for(int i = 1; i <= nbinsx; i++){
          univ_counts_in_bins[i].push_back(curhist->GetBinContent(i));
        }
      }
      else if(count==fNUniverses){
        int nbinsx = hCV->GetNbinsX();
        for(int i = 1; i <= nbinsx; i++){
          cv_counts_in_bins[i] = hCV->GetBinContent(i);
        }
      }
      count++;
    } //end loop on spectra.
    
    for(auto& it: univ_counts_in_bins){
      float cv_val = hCV->GetBinContent( it.first);
      if(cv_val <=0)continue;
      float err = 0.0;
      for(unsigned int iuniv=0;iuniv<fNUniverses;iuniv++){
        err += pow(it.second[iuniv] - cv_val,2);
      }
      err /= float(fNUniverses);
      err = sqrt(err);
      hCV->SetBinError(it.first, err);
    }    
    hCV->SetLineColor(kRed);
    hCV->SetMarkerColor(kRed);
    
    return hCV;
  }
  
  // return upper boundary
  Spectrum* FluxMultiverseSyst::UpperSigma()
  {
    if(!fBoundariesFound) FindSigmaBoundaries();
    return fUpperSigma;
  }
  
  // function to fill the boudary spactra encircling the error band
  void FluxMultiverseSyst::FindSigmaBoundaries(){
    Eigen::ArrayXd aCV = CV()->GetEigen(CV()->POT());
    map< int,float > cv_events_in_bins;
    map<int, vector<float> > univ_events_in_bins;
    unsigned int count = 0;
    for(auto& it: fSpectra){
      if(count!=fNUniverses){
        Eigen::ArrayXd curhist = it->GetEigen(it->POT());
        int nbinsx = curhist.size()-2;
        for(int i = 1; i <= nbinsx; i++){
          univ_events_in_bins[i].push_back(curhist[i]);
        }
      }
      else if(count==fNUniverses){
        int nbinsx = aCV.size()-2;
        for(int i = 1; i <= nbinsx; i++){
          cv_events_in_bins[i] = aCV[i];
        }
      }
      count++;
    }
    
    Eigen::ArrayXd aUp = aCV;
    Eigen::ArrayXd aLow = aCV;
    int binIt = 1;
    for(auto& it: univ_events_in_bins){
      double cv_val = aCV[it.first];
      double err = 0.0;
      for(unsigned int iuniv=0;iuniv<fNUniverses;iuniv++){
        err += pow(it.second[iuniv] - cv_val,2);
      }
      err /= double(fNUniverses);
      err = sqrt(err);
      aUp[binIt] = cv_val+err;
      aLow[binIt] = cv_val-err;
      binIt++;
    }
    
    // release memory if they are already allocated
    delete fUpperSigma;
    delete fLowerSigma;
    delete fNormUpperSigma;
    delete fNormLowerSigma;
    fUpperSigma = new Spectrum(std::move(aUp),
                               HistAxis(CV()->GetLabels(), CV()->GetBinnings()),
                               fSpectra[0]->POT(), fSpectra[0]->Livetime());
    fLowerSigma = new Spectrum(std::move(aLow),
                               HistAxis(CV()->GetLabels(), CV()->GetBinnings()),
                               fSpectra[0]->POT(), fSpectra[0]->Livetime());
    fBoundariesFound = true;
    Spectrum cv = *CV();
    Spectrum tempMC = cv;
    fNormUpperSigma = new Spectrum(*fUpperSigma);
    fNormUpperSigma->OverridePOT(fNormUpperSigma->POT() * fNormUpperSigma->Integral(1) / tempMC.Integral(1));
    fNormLowerSigma = new Spectrum(*fLowerSigma);
    fNormLowerSigma->OverridePOT(fNormLowerSigma->POT() * fNormLowerSigma->Integral(1) / tempMC.Integral(1));
  }
  
  // return lower boundary
  Spectrum* FluxMultiverseSyst::LowerSigma()
  {
    if(!fBoundariesFound) FindSigmaBoundaries();
    return fLowerSigma;
  }

  //----------------------------------------------------------------------
  void FluxMultiverseSyst::SaveTo(TDirectory* dir, const std::string& name) const {
    TDirectory *tmp = gDirectory;

    dir = dir->mkdir(name.c_str()); // switch to subdir
    dir->cd();
    
   TObjString universes(Form("%lu",fSpectra.size()));
    universes.Write("nUniverses");

    unsigned int count = 0;
    for(auto& it: fSpectra){
      if(count != fNUniverses){
        it->SaveTo(dir, "flux_universe_"+std::to_string(count));
      }
      if(count == fNUniverses){
        it->SaveTo(dir, "flux_cv");
      }
      count++;
    }// end loop over spectra

    dir->Write();
    delete dir;

    tmp->cd();
  }// end of SaveTo
  
  
  std::unique_ptr<FluxMultiverseSyst> FluxMultiverseSyst::LoadFrom(TDirectory* dir, const std::string& name) {
    dir = dir->GetDirectory(name.c_str()); // switch to subdir
    assert(dir);

    TObjString* universes = (TObjString*)dir->Get("nUniverses");
    unsigned int nUniverses = universes->GetString().Atoi();

    std::vector<std::unique_ptr<ana::Spectrum>> spectra;
    for(unsigned int i = 0; i < nUniverses; i++){
      spectra.push_back(ana::LoadFrom<Spectrum>(dir, "flux_universe_"+std::to_string(i)));
    }
    spectra.push_back(ana::LoadFrom<Spectrum>(dir, "flux_cv"));

    delete dir;

    return std::make_unique<FluxMultiverseSyst>(nUniverses, spectra, true);
   
  }// end of LoadFrom
}