TwoSampleDecomp.cxx

Go to the documentation of this file.

#include "CAFAna/Decomp/TwoSampleDecomp.h"

#include <iostream>
#include <cmath>
#include <cassert>
#include <exception>
#include <iostream>
#include <string>
#include <sstream>

#include "CAFAna/Core/Spectrum.h"

#include "TH1D.h"
#include "TDirectory.h"
#include "TObjString.h"
#include "Utilities/func/MathUtil.h"
#include <iostream>

namespace ana
{

  TwoSampleDecomp::TwoSampleDecomp( 
    SpectrumLoaderBase& loaderMC,  ///< Spectrum loader for MC
    SpectrumLoaderBase& loaderData,///< Spectrum loader for data
    std::string label,             ///< Name of the thing
    const Binning& bins,           ///< Desired binning for spectra
    const Var& var,                ///< Variable to be plotted in spectra
    const Cut& firstSample,        ///< Cut to define samples, 
                                   ///< true: first, false: second 
    const Cut& isSignal,           ///< Signal or not? true:sig, false: bkg
    const Cut& anaCut,             ///< Cut to define selection (e.g. numu CC)
    const SystShifts& shift,
    const Var& weight              ///< Weight for events in spectrum
    ):TwoSampleDecomp(loaderMC, loaderData, HistAxis(label, bins, var), 
                      firstSample, isSignal, anaCut, shift, weight)
    {}



  TwoSampleDecomp::TwoSampleDecomp( 
    SpectrumLoaderBase& loaderMC,   ///< Spectrum loader for MC
    SpectrumLoaderBase& loaderData, ///< Spectrum loader for data
    const HistAxis& axis,           ///< Variable to be plotted in spectra
    const Cut& firstSample,         ///< Cut to define samples, 
                                    ///< true: first, false: second 
    const Cut& isSignal,            ///< Signal or not? true:sig, false: bkg
    const Cut& anaCut,              ///< Cut to define selection (e.g. numu CC)
    const SystShifts& shift,
    const Var& weight               ///< Weight for events in spectrum
    ):
  fLoaderMC  (loaderMC),
  fLoaderData(loaderData),
  fAxis(&axis),
  fFirstSig (loaderMC, axis,  firstSample &&  isSignal && anaCut, 
             shift, weight),
  fSecondSig(loaderMC, axis, !firstSample &&  isSignal && anaCut, 
             shift, weight),
  fFirstBkg (loaderMC, axis,  firstSample && !isSignal && anaCut, 
             shift, weight),
  fSecondBkg(loaderMC, axis, !firstSample && !isSignal && anaCut, 
             shift, weight),
  fTrueSig  (loaderMC, axis,                  isSignal && anaCut, 
             shift, weight),
  fTrueBkg  (loaderMC, axis,                 !isSignal && anaCut, 
             shift, weight),
  fFirst    (loaderData, axis,  firstSample            && anaCut, 
             shift, weight),
  fSecond   (loaderData, axis, !firstSample            && anaCut, 
             shift, weight),
  fSig      (Spectrum::Uninitialized()),
  fBkg      (Spectrum::Uninitialized()),
  fSigRatio (0),
  fBkgRatio (0),
  fDecomposed(false)
  {
    if(axis.NDimensions() != 1){
      std::cout << "Error: TwoSampleDecomp not equipped to handle multi-dimensional axis." << std::endl;
      abort();
    }
  }



  TwoSampleDecomp::TwoSampleDecomp():
  fLoaderMC  (kNullLoader),
  fLoaderData(kNullLoader),
  fFirstSig (Spectrum::Uninitialized()),
  fSecondSig(Spectrum::Uninitialized()),
  fFirstBkg (Spectrum::Uninitialized()),
  fSecondBkg(Spectrum::Uninitialized()),
  fTrueSig  (Spectrum::Uninitialized()),
  fTrueBkg  (Spectrum::Uninitialized()),
  fFirst    (Spectrum::Uninitialized()),
  fSecond   (Spectrum::Uninitialized()),
  fSig      (Spectrum::Uninitialized()),
  fBkg      (Spectrum::Uninitialized()),
  fSigRatio (0),
  fBkgRatio (0),
  fDecomposed(false)
  {  }


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

    dir = dir->mkdir(name.c_str()); // switch to subdir
    dir->cd();

    TObjString("TwoSampleDecomp").Write("type");
                
    fFirstSig.SaveTo(dir, "FirstSig");
    fSecondSig.SaveTo(dir, "SecondSig");
    fFirstBkg.SaveTo(dir, "FirstBkg");
    fSecondBkg.SaveTo(dir, "SecondBkg");
    fTrueSig.SaveTo(dir, "TrueSig");
    fTrueBkg.SaveTo(dir, "TrueBkg");
    fFirst.SaveTo(dir, "First");
    fSecond.SaveTo(dir, "Second");
    fSig.SaveTo(dir, "Sig");
    fBkg.SaveTo(dir, "Bkg");

    dir->Write();
    delete dir;

    tmp->cd();
  }

  //----------------------------------------------------------------------
  std::unique_ptr<TwoSampleDecomp> TwoSampleDecomp::LoadFrom(TDirectory* dir, const std::string& name)
  {
    dir = dir->GetDirectory(name.c_str()); // switch to subdir
    assert(dir);

    std::unique_ptr<TwoSampleDecomp> ret(new TwoSampleDecomp);

    ret->fFirstSig = *Spectrum::LoadFrom(dir, "FirstSig");
    ret->fSecondSig = *Spectrum::LoadFrom(dir, "SecondSig");
    ret->fFirstBkg = *Spectrum::LoadFrom(dir, "FirstBkg");
    ret->fSecondBkg = *Spectrum::LoadFrom(dir, "SecondBkg");
    ret->fTrueSig = *Spectrum::LoadFrom(dir, "TrueSig");
    ret->fTrueBkg = *Spectrum::LoadFrom(dir, "TrueBkg");
    ret->fFirst = *Spectrum::LoadFrom(dir, "First");
    ret->fSecond = *Spectrum::LoadFrom(dir, "Second");
    ret->fSig = *Spectrum::LoadFrom(dir, "Sig");
    ret->fBkg = *Spectrum::LoadFrom(dir, "Bkg");    

    delete dir;

    return ret;
  }

  Spectrum TwoSampleDecomp::NCComponent()      const {std::cout << "TwoSampleDecomp::NCComponent() not implemented" << std::endl; abort();}

  Spectrum TwoSampleDecomp::NCAntiComponent()  const {std::cout << "TwoSampleDecomp::NCAntiComponent() not implemented" << std::endl; abort();}

  void TwoSampleDecomp::Decomp() const{
    if(fDecomposed) return;

    assert(fLoaderMC.Gone() 
          && "Call Go() on MC SpectrumLoader before TwoSampleDecomp::Decomp");
    assert(fLoaderData.Gone() 
          && "Call Go() on Data SpectrumLoader before TwoSampleDecomp::Decomp");

    RatioCalc();

    TH1D* hFirst = fFirst.ToTH1(fFirst.POT());
    TH1D* hSecond = fSecond.ToTH1(fSecond.POT());
    
    // Get the number of bins from the axis 
    int nBins = hFirst->GetNbinsX();

    // Make sure the histogram binning is sane 
    assert(hFirst->GetNbinsX() == nBins && 
       "Inconsistent binning across spectra.");
    assert(hSecond->GetNbinsX() == nBins && 
       "Inconsistent binning across spectra.");
    assert(fSigRatio->GetNbinsX() == nBins && 
       "MC and data have different numbers of bins.");
    assert(fBkgRatio->GetNbinsX() == nBins && 
       "MC and data have different numbers of bins.");

    std::stringstream title;
    title  << std::string(";") << hFirst->GetXaxis()->GetTitle() 
           << std::string(";Events");


    float axisMin = hFirst->GetBinLowEdge(1);
    float axisMax = hFirst->GetBinLowEdge(nBins+1);

    TH1D* hSig = new TH1D("Sig", title.str().c_str(), nBins, axisMin, axisMax);
    TH1D* hBkg = new TH1D("Bkg", title.str().c_str(), nBins, axisMin, axisMax);

    for(int i = 1; i <= nBins; ++i)
    {
      //These guys are what we're after, signal and background. Initialize them.  
      double nSig    = 0;
      double nBkg    = 0;
      double errSig  = 0;
      double errBkg  = 0;


      // Grab all of the things out of the histograms and store them so that 
      // the math looks a bit nicer.
      double rSig    = fSigRatio->GetBinContent(i);
      double rSigErr = fSigRatio->GetBinError(i);
      double rBkg    = fBkgRatio->GetBinContent(i);
      double rBkgErr = fBkgRatio->GetBinError(i);
      double nFirst  = hFirst->GetBinContent(i);
      double nSecond = hSecond->GetBinContent(i);

      // d2 is the square of the determinant of the matrix we're inverting.  
      // It shows up in the math a lot, so we'll store it for neatness.  
      double d2  = util::sqr(rSig - rBkg);


      //We're about to divide by d2 until the cows come home, so make sure it's
      //not zero.  When it is zero, it's probably because those bins were empty.  
      if  (d2 != 0) 
      {

        //Do the math in a way that stores the common terms in the error formula
        double sSig  = (1+rSig); 
        double sBkg  = (1+rBkg);
        nSig      = sSig * (nSecond - rBkg*nFirst) / (rSig - rBkg);
        nBkg      = sBkg * (nFirst*rSig - nSecond) / (rSig - rBkg);

        // Here's the error formula 
        errSig = sqrt((nSecond * util::sqr(sSig))/d2 + 
                    (nFirst*util::sqr(sSig) *util::sqr(rBkg))/d2 + 
                    util::sqr(rSigErr) * (-((sSig*(nSecond - nFirst*rBkg))/d2) + 
                    util::sqr(nSecond - nFirst*rBkg)/(rSig - rBkg)) + 
                    util::sqr(-((nFirst*(sSig))/(rSig - rBkg)) + 
                    (sSig*(nSecond - nFirst*rBkg))/d2) * util::sqr(rBkgErr));
        errBkg = sqrt((nSecond * util::sqr(sBkg))/d2 + 
                (nFirst*util::sqr(rSig) * util::sqr(sBkg))/d2 + 
                util::sqr(rSigErr)*util::sqr(-(((-nSecond + nFirst*rSig)*sBkg)/
                (rSig - rBkg)*2) + 
                (nFirst*sBkg)/(rSig - rBkg)) + 
                util::sqr((-nSecond + nFirst*rSig)/(rSig - rBkg) + 
                ((-nSecond + nFirst*rSig) *sBkg)/d2)*util::sqr(rBkgErr));
      }
      
      // Set the values in the resultant histograms, signal and background.  
      hSig->SetBinContent(i, nSig);
      hBkg->SetBinContent(i, nBkg);
      hSig->SetBinError(i, errSig);
      hBkg->SetBinError(i, errBkg);
    }

    fSig = Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(hSig->GetArray(), hSig->GetNbinsX()+2)),
                    HistAxis(fFirst.GetLabels(), fFirst.GetBinnings()),
                    fFirst.POT(), fFirst.Livetime());
    fBkg = Spectrum(Eigen::ArrayXd(Eigen::Map<Eigen::ArrayXd>(hBkg->GetArray(), hBkg->GetNbinsX()+2)),
                    HistAxis(fFirst.GetLabels(), fFirst.GetBinnings()),
                    fFirst.POT(), fFirst.Livetime());

    fDecomposed = true;
  }



  void TwoSampleDecomp::RatioCalc() const
  {


    // Check to make sure we haven't already computed ratios 
    // histogram pointers should be null
    assert(!fSigRatio && "Ratio histograms already constructed.");
    assert(!fBkgRatio && "Ratio histograms already constructed.");

    // Pull histograms out of the required Spectrum object
    TH1D* hFirstSig  = fFirstSig .ToTH1(fFirstSig.POT());        
    TH1D* hSecondSig = fSecondSig.ToTH1(fSecondSig.POT());      
    TH1D* hFirstBkg  = fFirstBkg .ToTH1(fFirstBkg.POT());   
    TH1D* hSecondBkg = fSecondBkg.ToTH1(fSecondBkg.POT());

    // Make the axis labels for the ratio histograms out of the 
    std::stringstream title;
    title  << std::string(";") << hFirstSig->GetXaxis()->GetTitle() 
           << std::string(";Ratio");

    // Get the number of bins from the axis 
    int nBins = hFirstSig->GetNbinsX();

    // Make sure the histogram binning is sane 
    assert(hFirstSig->GetNbinsX() == nBins &&
                                       "Inconsistent binning across spectra.");
    assert(hSecondSig->GetNbinsX() == nBins &&
                                       "Inconsistent binning across spectra.");
    assert(hFirstBkg->GetNbinsX() == nBins &&
                                       "Inconsistent binning across spectra.");
    assert(hSecondBkg->GetNbinsX() == nBins &&
                                       "Inconsistent binning across spectra.");

    const Binning bins = fAxis->GetBinnings()[0];

    // Construct the ratio histogram objects 
    fSigRatio = new TH1D("SigRatio", title.str().c_str(), nBins, 
                         bins.Min(), bins.Max());
    fBkgRatio = new TH1D("BkgRatio", title.str().c_str(), nBins, 
                         bins.Min(), bins.Max());


    // Loop over the bins in the histograms to compute the ratios bin-by-bin
    for (int i = 1; i <= nBins; ++i)
    {
      // Ratios for signal 
      if(hFirstSig->GetBinContent(i) != 0 and hSecondSig->GetBinContent(i) != 0)
      {
        // if non zero, make an object that does the math
        RatioError ratioError(hSecondSig->GetBinContent(i),  
                              hFirstSig->GetBinContent(i));
        // Set the bins 
        fSigRatio->SetBinContent(i, ratioError.fRatio);
        fSigRatio->SetBinError  (i, ratioError.fError);
      }
      // Repeat for background 
      if(hFirstBkg->GetBinContent(i) != 0 and hSecondBkg->GetBinContent(i) != 0)
      {
        RatioError ratioError(hSecondBkg->GetBinContent(i),  
                              hFirstBkg->GetBinContent(i));
        fBkgRatio->SetBinContent(i, ratioError.fRatio);
        fBkgRatio->SetBinError  (i, ratioError.fError);
      }
    }
  }


  void TwoSampleDecomp::DrawSigBkgOverlay()
  {
    Decomp();

    Spectrum data = fFirst + fSecond;

    TH1D* hData = data.ToTH1(data.POT());
    hData->Draw();
  
    TH1D* hSig = fSig.ToTH1(data.POT());
    TH1D* hBkg = fBkg.ToTH1(data.POT());
  
    hSig->SetLineColor(kGreen+2);
    hBkg->SetLineColor(kBlue-4);
    hSig->SetMarkerColor(kGreen+2);
    hBkg->SetMarkerColor(kBlue-4);
  
    hSig->Draw("same");
    hBkg->Draw("same");
  
  
    TH1D* hTrueSig = fTrueSig.ToTH1(data.POT());
    TH1D* hTrueBkg = fTrueBkg.ToTH1(data.POT());
  
    hTrueSig->SetLineColor(kGreen+2);
    hTrueBkg->SetLineColor(kBlue-4);
  
    hTrueSig->Draw("hist same");
    hTrueBkg->Draw("hist same");
  
    TH1D* hTrueTotal = (fTrueSig + fTrueBkg).ToTH1(data.POT());
  
    hTrueTotal->SetLineColor(kRed);
    hTrueTotal->Draw("hist same");
    
      
  }

  
  void TwoSampleDecomp::DrawTwoSamplesWithRatios()
  {
    Decomp();

    TH1D* hFirstSig = fFirstSig.ToTH1(fFirst.POT());
    TH1D* hSecondSig = fSecondSig.ToTH1(fFirst.POT());
    TH1D* hFirstBkg = fFirstBkg.ToTH1(fFirst.POT());
    TH1D* hSecondBkg = fSecondBkg.ToTH1(fFirst.POT());

    
    std::vector<TH1D*> hists = {hFirstSig, hSecondSig, hFirstBkg, hSecondBkg};

    float max = 0;
    for(const auto& h:hists){
      float m = h->GetBinContent(h->GetMaximumBin());
      if(m > max) max = m;
    }

    for(const auto& h:hists) h->SetMaximum(1.2*max);

    hFirstSig->SetLineColor(kRed);
    hSecondSig->SetLineColor(kBlue-4);
    hFirstBkg->SetLineColor(kRed);
    hSecondBkg->SetLineColor(kBlue-4);

    hFirstBkg->SetLineStyle(2);
    hSecondBkg->SetLineStyle(2);

    hFirstSig->Draw("hist");
    hSecondSig->Draw("hist same");
    hFirstBkg->Draw("hist same");
    hSecondBkg->Draw("hist same");
        
    
  }

  void TwoSampleDecomp::DrawRatios()
  {
    fSigRatio->Draw("");
    fBkgRatio->SetLineStyle(2);
    fBkgRatio->Draw("same");
  }



  RatioError::RatioError(double numerator, double denominator):
  fRatio(numerator/denominator),  // Ratio is the ratio
  // Error for gaussian counting errors  
  fError(fRatio * sqrt(1/ numerator + 1/denominator)) 
  {  }

}