TrivialCrossSectionAnalysis.cxx

Go to the documentation of this file.

#include "CAFAna/XSec/TrivialCrossSectionAnalysis.h"
#include "CAFAna/XSec/Multiverse.h"
#include "CAFAna/XSec/TrivialSignalEstimator.h"
#include "TObjString.h"
#include "CAFAna/Core/Ratio.h"
#include "CAFAna/XSec/Flux.h"

#include "TH2.h"
#include "TH3.h"
#include "TVectorD.h"

namespace ana
{
  ///////////////////////////////////////////////////////////////////////
  void 
  TrivialCrossSectionAnalysis::InitSystematicSpectra(XSecSysts::Syst_t syst,
                                                     SpectrumLoaderBase * loader)
  {
    assert( CanFillSpectra() &&
            "TrivialCrossSectionAnalysis object not yet ready to setup systematic spectra." );

    fXSecSpectra.insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
                        (syst,
                         new CrossSectionSpectra(*loader,
                                                 *fRecoAxis,
                                                 *fTruthAxis,
                                                 *fSelectionCut,
                                                 *fSignalCut,
                                                 new TrivialSignalEstimator(loader,
                                                                            fRecoAxis,
                                                                            *fSelectionCut,
                                                                            fBkgdCuts,
                                                                            *fSystShifts[syst],
                                                                            *fWeights[syst]),
                                                 fFidMax,
                                                 fFidMin,
                                                 *fENuBinning,
                                                 fPDG,
                                                 *fWeights[syst],
                                                 *fNuTruthWeights[syst],
                                                 *fSystShifts[syst])));
  }
  
  ///////////////////////////////////////////////////////////////////////
  void 
  TrivialCrossSectionAnalysis::InitSystematicSpectra(XSecSysts::Syst_t syst,
                                                     SpectrumLoaderBase * loader,
                                                     const Var * weight,
                                                     const NuTruthVar * weight_nt,
                                                     const SystShifts * shift)
  {
    SetWeight(syst, weight);
    SetNuTruthWeight(syst, weight_nt);
    SetSystShifts(syst, shift);
    InitSystematicSpectra(syst, loader);
  }

  ///////////////////////////////////////////////////////////////////////
  void
  TrivialCrossSectionAnalysis::InitGenieUniverses(SpectrumLoaderBase * loader,
                                                   std::vector<SystShifts> genie_shifts)
  {
    assert( CanFillSpectra() &&
            "CrossSectionAnalysis object not yet ready to setup genie multiverses." );
    fGenieMVEffNum = new Multiverse(*loader,
                                    HistAxisFromNuTruthHistAxis(*fTruthAxis),
                                    *fSelectionCut && CutFromNuTruthCut(*fSignalCut),
                                    genie_shifts,
                                    *fWeights[XSecSysts::kNominal]);
    fGenieMVEffDenom = new Multiverse(*loader,
                                      *fTruthAxis,
                                      *fSignalCut,
                                      genie_shifts,
                                      *fNuTruthWeights[XSecSysts::kNominal]);
    fHasGenieMultiverses = true;
  }
                                                   
  ///////////////////////////////////////////////////////////////////////
  void
  TrivialCrossSectionAnalysis::InitPPFXUniverses(SpectrumLoaderBase * loader,
                                                  std::vector<Var> ppfx_weights,
                                                  std::vector<NuTruthVar> ppfx_weights_nt)
  {
    assert( CanFillSpectra() &&
            "CrossSectionAnalysis object not yet ready to setup PPFX multiverses.");
    std::vector<Spectrum *> mv_flux;
    for(int i = 0; i < (int) ppfx_weights.size(); i++) {
      mv_flux.push_back(DeriveFlux(*loader,
                                   *fENuBinning,
                                   fPDG,
                                   &fFidMin,
                                   &fFidMax,
                                   *fSystShifts[XSecSysts::kNominal],
                                   ppfx_weights_nt[i]));
    }
    fPPFXMVFlux = new Multiverse(mv_flux);
    fPPFXMVEffNum = new Multiverse(*loader,
                                   HistAxisFromNuTruthHistAxis(*fTruthAxis),
                                   *fSelectionCut && CutFromNuTruthCut(*fSignalCut),
                                   *fSystShifts[XSecSysts::kNominal],
                                   ppfx_weights,
                                   *fWeights[XSecSysts::kNominal]);
    fPPFXMVEffDenom = new Multiverse(*loader,
                                     *fTruthAxis,
                                     *fSignalCut,
                                     *fSystShifts[XSecSysts::kNominal],
                                     ppfx_weights_nt,
                                     *fNuTruthWeights[XSecSysts::kNominal]);

    // make sure we transferred ownership completely to the multiverse object
    for(int i = 0; i < (int) ppfx_weights.size(); i++) {
      mv_flux[i] = NULL;
    }
    fHasPPFXMultiverses = true;
  }

  ///////////////////////////////////////////////////////////////////////
  TH1 *
  TrivialCrossSectionAnalysis::Purity(XSecSysts::Syst_t syst) 
  {
    TH1 * num = ToHist(*fXSecSpectra[syst]->GetSelectedSignalRecoSpace(), fPOT);
    TH1 * denom = ToHist(*fXSecSpectra[syst]->GetSelectedAllRecoSpace(), fPOT);
    num->Divide(num, denom, 1, 1, "B"); // calculate binomial errors
    return num;
  }

  ///////////////////////////////////////////////////////////////////////
  // basic constructor
  TrivialCrossSectionAnalysis::TrivialCrossSectionAnalysis(const Cut & selection_cut,
                                                           const HistAxis * reco_axis,
                                                           const NuTruthCut & signal_cut,
                                                           const NuTruthHistAxis * truth_axis,
                                                           TVector3 fid_max,
                                                           TVector3 fid_min,
                                                           const Binning * enu_binning,
                                                           std::vector<Cut> bkgd_cuts,
                                                           int pdg,
                                                           bool is_differential)
    : fBkgdCuts(bkgd_cuts)
  {
    this->fPDG = pdg;
    this->fENuBinning = enu_binning;
    this->fFidMax = fid_max;
    this->fFidMin = fid_min;
    this->fTarget = TargetCount(fFidMin, fFidMax, 1e6);
    this->fSelectionCut = new Cut(selection_cut);
    this->fSignalCut = new NuTruthCut(signal_cut);
    this->fRecoAxis = reco_axis;
    this->fTruthAxis = truth_axis;
    this->fNDims = fRecoAxis->NDimensions();
    this->fIsDifferential = is_differential;
  }
  
  ///////////////////////////////////////////////////////////////////////
  std::pair<TH1 *, TH1 *>
  TrivialCrossSectionAnalysis::TotalErrors()
  {
    // load file systematic shifts
    std::vector<XSecSysts::Syst_t> systs = XSecSysts::GetAllSystematics();
    std::vector<TH1 *> shifts;
    for(XSecSysts::Syst_t syst : systs) {
      shifts.push_back(RelativeUncertainty(syst));
    }
    // genie and ppfx shifts
    std::pair<TH1 *, TH1 *> flux_error = FluxError();
    std::pair<TH1 *, TH1 *> genie_error = GenieError();
    TH1 * flux_up = flux_error.first;
    TH1 * flux_down = flux_error.second;
    TH1 * genie_up = genie_error.first;
    TH1 * genie_down = genie_error.second;
    
    TH1 * nominal = fXSec[XSecSysts::kNominal];
    TH1 * errup = (TH1*) nominal->Clone();
    TH1 * errdown = (TH1*) nominal->Clone();
    double up, down, delta;

    for(int i = 1; i <= nominal->GetNbinsX(); i++) {
      for(int j = 1; j <= nominal->GetNbinsY(); j++) {
        for(int k = 1; k <= nominal->GetNbinsZ(); k++) {
          double totalup = 0;
          double totaldown = 0;
          // process file systematics
          for(int ishift = 0; ishift < (int) shifts.size(); ishift++) {
            delta = 
              shifts[ishift]->GetBinContent(i, j, k) - nominal->GetBinContent(i, j, k);
            if(delta > 0) {
              up = delta;
              down = 0;
            }
            else {
              up = 0;
              down = delta;
            }
            totalup += up * up;
            totaldown += down * down;
          }
          // now genie and ppfx universes
          up = genie_up->GetBinContent(i, j, k) - nominal->GetBinContent(i, j, k);
          down = genie_down->GetBinContent(i, j, k) - nominal->GetBinContent(i, j, k);
          totalup += up * up;
          totaldown += down * down;

          up = flux_up->GetBinContent(i, j, k) - nominal->GetBinContent(i, j, k);
          down = flux_down->GetBinContent(i, j, k) - nominal->GetBinContent(i, j, k);
          totalup += up * up;
          totaldown += down * down;
          errup->SetBinContent(i, j, k, totalup);
          errdown->SetBinContent(i, j, k, totaldown);
        }
      }
    }
    return std::pair<TH1 *, TH1 *> (errup, errdown);
  }

  ///////////////////////////////////////////////////////////////////////
  TH1 *
  TrivialCrossSectionAnalysis::RelativeUncertainty(XSecSysts::Syst_t syst)
  {
    TH1 * rel = (TH1*) fXSec[XSecSysts::kNominal]->Clone();
    double shift, cv, delta;
    for(int i = 1; i <= rel->GetNbinsX(); i++) {
      for(int j = 1; j <= rel->GetNbinsY(); j++) {
        for(int k = 1; k <= rel->GetNbinsZ(); k++) {
          cv = fXSec[XSecSysts::kNominal]->GetBinContent(i, j, k);
          shift = fXSec[syst]->GetBinContent(i, j, k);    
          delta = fabs(shift - cv);
          rel->SetBinContent(i, j, k, delta / cv);
        }
      }
    }
    return rel;
  }

  ///////////////////////////////////////////////////////////////////////
  std::pair<TH1 *, TH1 *>
  TrivialCrossSectionAnalysis::GenieError()
  {
    // xsec is proportional to 1/efficiency so efficiency denom goes up in the 
    // numerator
    Multiverse * genie_xsec = new Multiverse(*fGenieMVEffDenom);
    genie_xsec->Divide(*fGenieMVEffNum, true);
    *genie_xsec /= *Flux(XSecSysts::kNominal);
    *genie_xsec *= *UnfoldedSignal(XSecSysts::kNominal);
    return std::pair<TH1 *, TH1 *>
      (genie_xsec->GetPlusOneSigmaShift(fXSec[XSecSysts::kNominal]),
       genie_xsec->GetMinusOneSigmaShift(fXSec[XSecSysts::kNominal]));
  }

  ///////////////////////////////////////////////////////////////////////
  std::pair<TH1 *, TH1 *>
  TrivialCrossSectionAnalysis::FluxError()
  {
    // xsec is proportional to 1/efficiency so efficiency denom goes up in the 
    // numerator
    Multiverse * ppfx_xsec = new Multiverse(*fPPFXMVEffDenom);
    ppfx_xsec->Divide(*fPPFXMVEffNum, true);
    ppfx_xsec->Divide(*fPPFXMVFlux);
    *ppfx_xsec *= *UnfoldedSignal(XSecSysts::kNominal);

    return std::pair<TH1 *, TH1 *>
      (ppfx_xsec->GetPlusOneSigmaShift(fXSec[XSecSysts::kNominal]),
       ppfx_xsec->GetMinusOneSigmaShift(fXSec[XSecSysts::kNominal]));
  }

  ///////////////////////////////////////////////////////////////////////
  TH1 *
  TrivialCrossSectionAnalysis::Flux(XSecSysts::Syst_t syst)
  {
    return ToHist(*fXSecSpectra[syst]->GetFlux(), fPOT);
  }

  ///////////////////////////////////////////////////////////////////////
  TH1 *
  TrivialCrossSectionAnalysis::Efficiency(XSecSysts::Syst_t syst)
  {
    TH1 * num = ToHist(*fXSecSpectra[syst]->GetSelectedSignalTrueSpace(), fPOT);
    TH1 * denom = ToHist(*fXSecSpectra[syst]->GetTrueSignalTrueSpace(), fPOT);
    num->Divide(num, denom, 1, 1, "B"); // calculate binomial errors
    return num;
  }

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

    TObjString * ptag = (TObjString*) dir->Get("type");
    assert(ptag);
    
    const TString tag = ptag->GetString();
    assert(tag == "TrivialCrossSectionAnalysis" && "Type does not match TrivialCrossSectionAnalysis");
    delete ptag;

    Spectrum * data = Spectrum::LoadFrom(dir, "fData").release();
    std::map<XSecSysts::Syst_t, CrossSectionSpectra *> xsec_spectra;
    std::vector<XSecSysts::Syst_t> systs 
      = XSecSysts::GetAllSystematics();
    
    for(XSecSysts::Syst_t syst : systs) {
      xsec_spectra.
        insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *>
               (syst, 
                CrossSectionSpectra::LoadFrom(dir, XSecSysts::LabelFromSystematic(syst).Data()).release()));
    }

    TVector3 * fidmax = (TVector3*) dir->Get("fFidMax");
    TVector3 * fidmin = (TVector3*) dir->Get("fFidMin");
    TVectorD * ndims  = (TVectorD*) dir->Get("ndims");
    Multiverse * genie_mv_eff_num = NULL;
    Multiverse * genie_mv_eff_denom = NULL;
    Multiverse * ppfx_mv_eff_num = NULL;
    Multiverse * ppfx_mv_eff_denom= NULL;
    Multiverse * ppfx_mv_flux = NULL;
    bool isdifferential = false;
    if(dir->GetDirectory("fIsDifferential"))
      isdifferential = true;
    if(dir->GetDirectory("fGenieMVEffNum")) {
      genie_mv_eff_num = 
        Multiverse::LoadFrom(dir, "fGenieMVEffNum").release();
    }
    if(dir->GetDirectory("fGenieMVEffDenom")) {
      genie_mv_eff_denom = 
        Multiverse::LoadFrom(dir, "fGenieMVEffDenom").release();
    }
    if(dir->GetDirectory("fPPFXMVFlux")) {
      ppfx_mv_flux = 
        Multiverse::LoadFrom(dir, "fPPFXMVFlux").release();
    }
    if(dir->GetDirectory("fPPFXMVEffNum")) {
      ppfx_mv_eff_num = 
        Multiverse::LoadFrom(dir, "fPPFXMVEffNum").release();
    }
    if(dir->GetDirectory("fPPFXMVEffDenom")) {
      ppfx_mv_eff_denom = 
        Multiverse::LoadFrom(dir, "fPPFXMVEffDenom").release();
    }

    std::unique_ptr<TrivialCrossSectionAnalysis> ret
      ( new TrivialCrossSectionAnalysis(data,
                                        xsec_spectra,
                                        ppfx_mv_flux,
                                        ppfx_mv_eff_num,
                                        ppfx_mv_eff_denom,
                                        genie_mv_eff_num,
                                        genie_mv_eff_denom,
                                        fidmax,
                                        fidmin,
                                        ndims,
                                        isdifferential));
    double POT = ret->fData->POT();
    ret->SetPOT(POT);
    for(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *> xsec : xsec_spectra) 
      if(xsec.second)
        xsec.second = NULL;
    data = NULL;
    genie_mv_eff_num = NULL;
    genie_mv_eff_denom = NULL;
    ppfx_mv_flux = NULL;
    ppfx_mv_eff_num = NULL;
    ppfx_mv_eff_denom = NULL;
    fidmax = NULL;
    fidmin = NULL;
    ndims = NULL;

    delete dir;

    return ret;
  }
      
  ///////////////////////////////////////////////////////////////////////
  void
  TrivialCrossSectionAnalysis::SaveTo(TDirectory * dir, const std::string& name) const
  {
    TDirectory * tmp = gDirectory;

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

    TObjString("TrivialCrossSectionAnalysis").Write("type");
    if(fIsDifferential)
      TObjString("true").Write("fIsDifferential");
    for(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *> xsec : fXSecSpectra) {
      if(xsec.second)
        xsec.second->SaveTo(dir, XSecSysts::LabelFromSystematic(xsec.first).Data());
    }
    fFidMax.Write("fFidMax");
    fFidMin.Write("fFidMin");
    TVectorD ndims(1);
    ndims[0] = fNDims;
    ndims.Write("ndims");
    if(fHasGenieMultiverses) {
      fGenieMVEffNum->SaveTo(dir, "fGenieMVEffNum");
      fGenieMVEffDenom->SaveTo(dir, "fGenieMVEffDenom");
    }
    if(fHasPPFXMultiverses) {
      fPPFXMVEffNum->SaveTo(dir, "fPPFXMVEffNum");
      fPPFXMVEffDenom->SaveTo(dir, "fPPFXMVEffDenom");
      fPPFXMVFlux->SaveTo(dir, "fPPFXMVFlux");
    }
    if(HasData())
      fData->SaveTo(dir, "fData");

    dir->Write();
    delete dir;

    tmp->cd();
  }

  ///////////////////////////////////////////////////////////////////////
  ICrossSectionAnalysis * 
  TrivialCrossSectionAnalysis::Clone() const
  {
    return new TrivialCrossSectionAnalysis(*this);
  }

  ///////////////////////////////////////////////////////////////////////
  // Make sure things get properly allocated during assignment
  TrivialCrossSectionAnalysis &
  TrivialCrossSectionAnalysis::operator=(const TrivialCrossSectionAnalysis & rhs)
  {
    if ( this == &rhs )
      return *this;
    std::for_each(rhs.fXSecSpectra.begin(), rhs.fXSecSpectra.end(),
                  [this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*> xsec) 
                  {
                    this->fXSecSpectra
                      .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
                              (xsec.first,
                               new CrossSectionSpectra(*xsec.second)));
                  });

    std::for_each(rhs.fWeights.begin(), rhs.fWeights.end(),
                  [this](std::pair<XSecSysts::Syst_t, const Var *> weight) 
                  {
                    this->fWeights
                      .insert(std::pair<XSecSysts::Syst_t, const Var *>
                              (weight.first,
                               Copy(weight.second)));
                  });
                  
    std::for_each(rhs.fSystShifts.begin(), rhs.fSystShifts.end(),
                  [this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift) 
                  {
                    this->fSystShifts
                      .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
                              (shift.first,
                               Copy(shift.second)));
                  });
                  
    this->fTarget = rhs.fTarget;
    this->fFidMax = rhs.fFidMax;
    this->fFidMin = rhs.fFidMin;
    this->fSelectionCut   = Copy(rhs.fSelectionCut);
    this->fRecoAxis       = Copy(rhs.fRecoAxis);
    this->fSignalCut      = Copy(rhs.fSignalCut);
    this->fTruthAxis      = Copy(rhs.fTruthAxis);
    this->fUnfoldMethod   = rhs.fUnfoldMethod;
    this->fUnfoldReg      = rhs.fUnfoldReg;
    this->fNDims          = rhs.fNDims;
    this->fIsDifferential = rhs.fIsDifferential;
    if(rhs.fPPFXXSec.size() != 0) {
      std::for_each(rhs.fPPFXXSec.begin(), rhs.fPPFXXSec.end(),
                    [this](TH1 * ppfx)
                    {
                      this->fPPFXXSec.push_back((TH1*) ppfx->Clone());
                    });
    }
    if(rhs.fGenieXSec.size() != 0) {
      std::for_each(rhs.fGenieXSec.begin(), rhs.fGenieXSec.end(),
                    [this](TH1 * genie)
                    {
                      this->fGenieXSec.push_back((TH1*) genie->Clone());
                    });
    }

    this->fGenieMVEffNum   = new Multiverse(*rhs.fGenieMVEffNum);
    this->fGenieMVEffDenom = new Multiverse(*rhs.fGenieMVEffDenom);
    this->fPPFXMVFlux      = new Multiverse(*rhs.fPPFXMVFlux);
    this->fPPFXMVEffNum    = new Multiverse(*rhs.fPPFXMVEffNum);
    this->fPPFXMVEffDenom  = new Multiverse(*rhs.fPPFXMVEffDenom);
    this->fData            = new Spectrum(*rhs.fData);
    return *this;
  }

  ///////////////////////////////////////////////////////////////////////
  // Move assignment
  TrivialCrossSectionAnalysis &
  TrivialCrossSectionAnalysis::operator=(TrivialCrossSectionAnalysis && rhs)
  {
    if ( this == &rhs )
      return *this;
    std::for_each(rhs.fXSecSpectra.begin(), rhs.fXSecSpectra.end(),
                  [this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*> xsec) 
                  {
                    this->fXSecSpectra
                      .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
                              (xsec.first,
                               xsec.second));
                    xsec.second = NULL;
                  });

    std::for_each(rhs.fWeights.begin(), rhs.fWeights.end(),
                  [this](std::pair<XSecSysts::Syst_t, const Var *> weight) 
                  {
                    this->fWeights
                      .insert(std::pair<XSecSysts::Syst_t, const Var *>
                              (weight.first,
                               weight.second));
                    weight.second = NULL;
                  });
                  
    std::for_each(rhs.fSystShifts.begin(), rhs.fSystShifts.end(),
                  [this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift) 
                  {
                    this->fSystShifts
                      .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
                              (shift.first,
                               shift.second));
                    shift.second = NULL;
                  });
                  
    this->fTarget = rhs.fTarget;
    this->fFidMax = rhs.fFidMax;
    this->fFidMin = rhs.fFidMin;
    this->fSelectionCut = rhs.fSelectionCut; rhs.fSelectionCut = NULL;
    this->fRecoAxis     = rhs.fRecoAxis;     rhs.fRecoAxis     = NULL;
    this->fSignalCut    = rhs.fSignalCut;    rhs.fSignalCut    = NULL;
    this->fTruthAxis    = rhs.fTruthAxis;    rhs.fTruthAxis    = NULL;
    this->fUnfoldMethod = rhs.fUnfoldMethod;
    this->fUnfoldReg    = rhs.fUnfoldReg;
    this->fNDims        = rhs.fNDims;
    this->fIsDifferential = rhs.fIsDifferential;
    if(rhs.fPPFXXSec.size() != 0) {
      std::for_each(rhs.fPPFXXSec.begin(), rhs.fPPFXXSec.end(),
                    [this](TH1 * ppfx)
                    {
                      this->fPPFXXSec.push_back(ppfx);
                      ppfx = NULL;
                    });
    }
    if(rhs.fGenieXSec.size() != 0) {
      std::for_each(rhs.fGenieXSec.begin(), rhs.fGenieXSec.end(),
                    [this](TH1 * genie)
                    {
                      this->fGenieXSec.push_back(genie);
                      genie = NULL;
                    });
    }

    this->fGenieMVEffNum   = rhs.fGenieMVEffNum;   rhs.fGenieMVEffNum   = NULL;
    this->fGenieMVEffDenom = rhs.fGenieMVEffDenom; rhs.fGenieMVEffDenom = NULL;
    this->fPPFXMVFlux      = rhs.fPPFXMVFlux;      rhs.fPPFXMVFlux      = NULL;
    this->fPPFXMVEffNum    = rhs.fPPFXMVEffNum;    rhs.fPPFXMVEffNum    = NULL;
    this->fPPFXMVEffDenom  = rhs.fPPFXMVEffDenom;  rhs.fPPFXMVEffDenom  = NULL;
    this->fData            = rhs.fData;            rhs.fData            = NULL;
    return *this;
  }

  ///////////////////////////////////////////////////////////////////////
  // Move constructor
  TrivialCrossSectionAnalysis::TrivialCrossSectionAnalysis(TrivialCrossSectionAnalysis && rhs)
  {
    std::for_each(rhs.fXSecSpectra.begin(), rhs.fXSecSpectra.end(),
                  [this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*> xsec) 
                  {
                    this->fXSecSpectra
                      .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
                              (xsec.first,
                               xsec.second));
                    xsec.second = NULL;
                  });

    std::for_each(rhs.fWeights.begin(), rhs.fWeights.end(),
                  [this](std::pair<XSecSysts::Syst_t, const Var *> weight) 
                  {
                    this->fWeights
                      .insert(std::pair<XSecSysts::Syst_t, const Var *>
                              (weight.first,
                               weight.second));
                    weight.second = NULL;
                  });
                  
    std::for_each(rhs.fSystShifts.begin(), rhs.fSystShifts.end(),
                  [this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift) 
                  {
                    this->fSystShifts
                      .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
                              (shift.first,
                               shift.second));
                    shift.second = NULL;
                  });
                  
    this->fTarget = rhs.fTarget;
    this->fFidMax = rhs.fFidMax;
    this->fFidMin = rhs.fFidMin;
    this->fSelectionCut = rhs.fSelectionCut; rhs.fSelectionCut = NULL;
    this->fRecoAxis     = rhs.fRecoAxis;     rhs.fRecoAxis     = NULL;
    this->fSignalCut    = rhs.fSignalCut;    rhs.fSignalCut    = NULL;
    this->fTruthAxis    = rhs.fTruthAxis;    rhs.fTruthAxis    = NULL;
    this->fUnfoldMethod = rhs.fUnfoldMethod;
    this->fUnfoldReg    = rhs.fUnfoldReg;
    this->fNDims        = rhs.fNDims;
    this->fIsDifferential = rhs.fIsDifferential;
    if(rhs.fPPFXXSec.size() != 0) {
      std::for_each(rhs.fPPFXXSec.begin(), rhs.fPPFXXSec.end(),
                    [this](TH1 * ppfx)
                    {
                      this->fPPFXXSec.push_back(ppfx);
                      ppfx = NULL;
                    });
    }
    if(rhs.fGenieXSec.size() != 0) {
      std::for_each(rhs.fGenieXSec.begin(), rhs.fGenieXSec.end(),
                    [this](TH1 * genie)
                    {
                      this->fGenieXSec.push_back(genie);
                      genie = NULL;
                    });
    }

    this->fGenieMVEffNum   = rhs.fGenieMVEffNum;   rhs.fGenieMVEffNum   = NULL;
    this->fGenieMVEffDenom = rhs.fGenieMVEffDenom; rhs.fGenieMVEffDenom = NULL;
    this->fPPFXMVFlux      = rhs.fPPFXMVFlux;      rhs.fPPFXMVFlux      = NULL;
    this->fPPFXMVEffNum    = rhs.fPPFXMVEffNum;    rhs.fPPFXMVEffNum    = NULL;
    this->fPPFXMVEffDenom  = rhs.fPPFXMVEffDenom;  rhs.fPPFXMVEffDenom  = NULL;
    this->fData            = rhs.fData;            rhs.fData            = NULL;
  }

  ///////////////////////////////////////////////////////////////////////
  // Copy constructor
  TrivialCrossSectionAnalysis::TrivialCrossSectionAnalysis(const TrivialCrossSectionAnalysis & rhs)
  {
    std::for_each(rhs.fXSecSpectra.begin(), rhs.fXSecSpectra.end(),
                  [this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*> xsec) 
                  {
                    this->fXSecSpectra
                      .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
                              (xsec.first,
                               new CrossSectionSpectra(*xsec.second)));
                  });

    std::for_each(rhs.fWeights.begin(), rhs.fWeights.end(),
                  [this](std::pair<XSecSysts::Syst_t, const Var *> weight) 
                  {
                    this->fWeights
                      .insert(std::pair<XSecSysts::Syst_t, const Var *>
                              (weight.first,
                               Copy(weight.second)));
                  });
                  
    std::for_each(rhs.fSystShifts.begin(), rhs.fSystShifts.end(),
                  [this](std::pair<XSecSysts::Syst_t, const SystShifts *> shift) 
                  {
                    this->fSystShifts
                      .insert(std::pair<XSecSysts::Syst_t, const SystShifts *>
                              (shift.first,
                               Copy(shift.second)));
                  });
                  
    this->fTarget = rhs.fTarget;
    this->fFidMax = rhs.fFidMax;
    this->fFidMin = rhs.fFidMin;
    this->fSelectionCut = Copy(rhs.fSelectionCut);
    this->fSignalCut    = Copy(rhs.fSignalCut);
    this->fRecoAxis     = Copy(rhs.fRecoAxis);
    this->fTruthAxis    = Copy(rhs.fTruthAxis);
    this->fUnfoldMethod = rhs.fUnfoldMethod;
    this->fUnfoldReg    = rhs.fUnfoldReg;
    this->fNDims        = rhs.fNDims;
    this->fIsDifferential = rhs.fIsDifferential;
    if(rhs.fPPFXXSec.size() != 0) {
      std::for_each(rhs.fPPFXXSec.begin(), rhs.fPPFXXSec.end(),
                    [this](TH1 * ppfx)
                    {
                      this->fPPFXXSec.push_back((TH1*) ppfx->Clone());
                    });
    }
    if(rhs.fGenieXSec.size() != 0) {
      std::for_each(rhs.fGenieXSec.begin(), rhs.fGenieXSec.end(),
                    [this](TH1 * genie)
                    {
                      this->fGenieXSec.push_back((TH1*) genie->Clone());
                    });
    }
    this->fGenieMVEffNum = new Multiverse(*rhs.fGenieMVEffNum);
    this->fGenieMVEffDenom = new Multiverse(*rhs.fGenieMVEffDenom);
    this->fPPFXMVFlux = new Multiverse(*rhs.fPPFXMVFlux);
    this->fPPFXMVEffNum = new Multiverse(*rhs.fPPFXMVEffNum);
    this->fPPFXMVEffDenom = new Multiverse(*rhs.fPPFXMVEffDenom);
    this->fData = new Spectrum(*rhs.fData);
    }

  
  ///////////////////////////////////////////////////////////////////////
  TrivialCrossSectionAnalysis::TrivialCrossSectionAnalysis(Spectrum * data,
                                                           std::map<XSecSysts::Syst_t, CrossSectionSpectra*> _fXSecSpectra,
                                                           Multiverse * ppfx_mv_flux,
                                                           Multiverse * ppfx_mv_eff_num,
                                                           Multiverse * ppfx_mv_eff_denom,
                                                           Multiverse * genie_mv_eff_num,
                                                           Multiverse * genie_mv_eff_denom,
                                                           TVector3 * _fFidMax,
                                                           TVector3 * _fFidMin,
                                                           TVectorD * ndims,
                                                           bool isdifferential)
  {
    this->fData = data;
    std::for_each(_fXSecSpectra.begin(), _fXSecSpectra.end(),
                  [this](std::pair<XSecSysts::Syst_t, CrossSectionSpectra*> xsec) 
                  {
                    this->fXSecSpectra
                      .insert(std::pair<XSecSysts::Syst_t, CrossSectionSpectra*>
                              (xsec.first,
                               xsec.second));
                  });
    this->fFidMax = *_fFidMax;
    this->fFidMin = *_fFidMin;
    this->fNDims  = (*ndims)[0];
    this->fIsDifferential = isdifferential;
    this->fGenieMVEffNum = genie_mv_eff_num;
    this->fGenieMVEffDenom = genie_mv_eff_denom;
    this->fPPFXMVFlux = ppfx_mv_flux;
    this->fPPFXMVEffNum = ppfx_mv_eff_num;
    this->fPPFXMVEffDenom = ppfx_mv_eff_denom;

    for(std::pair<XSecSysts::Syst_t, CrossSectionSpectra *> xsec : _fXSecSpectra) 
      if(xsec.second != NULL)
        xsec.second = NULL;
    genie_mv_eff_num = NULL;
    genie_mv_eff_denom = NULL;
    ppfx_mv_flux = NULL;
    ppfx_mv_eff_num = NULL;
    ppfx_mv_eff_denom = NULL;
  }

  ///////////////////////////////////////////////////////////////////////
  TrivialCrossSectionAnalysis::~TrivialCrossSectionAnalysis()
  {
    for(int ipxsec = 0; ipxsec < (int) fPPFXXSec.size(); ipxsec++) {
      delete fPPFXXSec[ipxsec];
    }
    for(int igxsec = 0; igxsec < (int) fGenieXSec.size(); igxsec++) {
      delete fGenieXSec[igxsec];
    }
    delete fData;
    delete fGenieMVEffNum;
    delete fGenieMVEffDenom;
    delete fPPFXMVEffNum;
    delete fPPFXMVEffDenom;
    delete fPPFXMVFlux;
  }
}