MichelDecompTest.C

Go to the documentation of this file.

#include "OscLib/OscCalcPMNSOpt.h"

#include "CAFAna/Cuts/Cuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "3FlavorAna/Decomp/MichelDecomp.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "3FlavorAna/Extrap/NueExtrap.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Vars/Vars.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"

#include "Utilities/func/MathUtil.h"

#include "TCanvas.h"
#include "TH2.h"
#include "TFile.h"
#include "TLegend.h"

#include <iostream>
#include <cmath>

using namespace ana;

void resetCalc(osc::IOscCalcAdjustable* calc)
{
  calc->SetL(810);
  calc->SetRho(2.75);
  calc->SetDmsq21(7.59e-5);
  calc->SetDmsq32(2.40e-3);
  calc->SetTh12(asin(sqrt(.87))/2);
  calc->SetTh13(asin(sqrt(.095))/2);
  calc->SetTh23(asin(sqrt(.95))/2);
  calc->SetdCP(0);
}

void MichelDecompTest()
{

  const Cut cutNDContNue( kNDContainNue );
  const Cut cutNDPIDNue( kLEM > .64 );
  const Cut cutNDNue( cutNDContNue && cutNDPIDNue );

  const Binning MBinning = Binning::Simple(8,0,8);

  const std::string farNonS = "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/000148/*/*fhc_nonswap*.caf.root";
  const std::string farSwap = "$NOVA_PROD/mc/FA14-10-28/genie/fd/caf/000148/*/*fhc_swap*.caf.root";
  const std::string nearNonS   = "$NOVA_PROD/mc/FA14-11-25/genie/nd/caf/*103/10378/*s101*.caf.root";

  SpectrumLoader loaderNonS(farNonS);
  SpectrumLoader loaderSwap(farSwap);
  SpectrumLoader loaderNDMC(nearNonS);
  
  const Binning bins = Binning::Simple(6, 1.25, 2.75);
  const HistAxis EAxis("Energy (GeV)", bins, kCaloE);

  SystShifts shift;
  shift.SetShift(new GenieRwgtTableSyst(rwgt::fReweightRvnCC2pi), +2);
  shift.SetShift(new GenieRwgtTableSyst(rwgt::fReweightRvnNC2pi), +2);
  shift.SetShift(new GenieRwgtTableSyst(rwgt::fReweightRvbarnCC2pi), +2);
  shift.SetShift(new GenieRwgtTableSyst(rwgt::fReweightRvbarnNC2pi), +2);

  Spectrum CCMC  (loaderNDMC,   EAxis, cutNDNue && kIsNumuCC);
  Spectrum NCMC  (loaderNDMC,   EAxis, cutNDNue && kIsNC);
  Spectrum NEMC  (loaderNDMC,   EAxis, cutNDNue && kIsBeamNue);

  Spectrum CCData(loaderNDMC, EAxis, cutNDNue && kIsNumuCC,  shift);
  Spectrum NCData(loaderNDMC, EAxis, cutNDNue && kIsNC,      shift);
  Spectrum NEData(loaderNDMC, EAxis, cutNDNue && kIsBeamNue, shift);

  // N Michels kNoShifted / nominal spectra
  Spectrum CCMCME  (EAxis.label, loaderNDMC,
                    EAxis.bins, EAxis.var,
                    MBinning, kNMichels,
                    cutNDNue && kIsNumuCC);
  Spectrum NCMCME  (EAxis.label, loaderNDMC,
                    EAxis.bins, EAxis.var,
                    MBinning, kNMichels,
                    cutNDNue && kIsNC);
  Spectrum NEMCME  (EAxis.label, loaderNDMC,
                    EAxis.bins, EAxis.var,
                    MBinning, kNMichels,
                    cutNDNue && kIsBeamNue);
  Spectrum CCDataME(EAxis.label, loaderNDMC,
                    EAxis.bins, EAxis.var,
                    MBinning, kNMichels,
                    cutNDNue && kIsNumuCC, shift);
  Spectrum NCDataME(EAxis.label, loaderNDMC,
                    EAxis.bins, EAxis.var,
                    MBinning, kNMichels,
                    cutNDNue && kIsNC, shift);
  Spectrum NEDataME(EAxis.label, loaderNDMC,
                    EAxis.bins, EAxis.var,
                    MBinning, kNMichels,
                    cutNDNue && kIsBeamNue, shift);

  MichelDecomp mDCMP(loaderNDMC, loaderNDMC, 
                       EAxis, cutNDNue, kNoShift, shift);
  NumuDecomp numuDCMP(loaderNDMC, loaderNDMC, 
                      EAxis, kNDContainNue && kNumuCC, kNoShift, shift);
  NueExtrap extrapDCMP(loaderNDMC,
                                                loaderSwap,
                                                loaderNonS, 
                                                kNullLoader,
                                                mDCMP,
                                                numuDCMP,
                                                EAxis, 
                                                EAxis, 
                                                cutNDPIDNue, 
                                                cutNDNue,
                                                kNDContainNue && kNumuCC,
                                                kNoShift);
  PredictionNoExtrap predNominal(loaderNonS, loaderSwap,
                                 EAxis.label, EAxis.bins, EAxis.var, cutNDPIDNue, kNoShift);
  
  PredictionExtrap predDCMP(&extrapDCMP);
  
  PredictionNoExtrap predShifts(loaderNonS, loaderSwap,
                                EAxis.label, EAxis.bins, EAxis.var, cutNDPIDNue, shift);


  // Decomp Loaders
  loaderNDMC.Go();
  loaderSwap.Go();
  loaderNonS.Go();

  Spectrum CCDCMP = mDCMP.NumuComponent() + mDCMP.AntiNumuComponent();
  Spectrum NCDCMP = mDCMP.NCTotalComponent();
  Spectrum NEDCMP = mDCMP.NueComponent()  + mDCMP.AntiNueComponent();

  osc::OscCalcPMNSOpt calc;
  resetCalc(&calc);

  double ourPOT = CCDataME.POT();
    
  Spectrum obsNominal = predNominal.Predict(&calc).FakeData(ourPOT);
  SingleSampleExperiment exptNominal(&predNominal, obsNominal);

  Spectrum obsShifts = predShifts.Predict(&calc).FakeData(ourPOT);
  SingleSampleExperiment exptShifts(&predNominal, obsShifts);

  Spectrum obsDCMP = predShifts.Predict(&calc).FakeData(ourPOT);
  SingleSampleExperiment exptDCMP(&predDCMP, obsDCMP);

  Spectrum obsBlah = predDCMP.Predict(&calc).FakeData(ourPOT);

  Spectrum bkgdNominal = obsNominal - predNominal.PredictComponent(&calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth);
  Spectrum bkgdShifted = obsShifts - predShifts.PredictComponent(&calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth);
  Spectrum bkgdDCMP = obsBlah - predDCMP.PredictComponent(&calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kBoth);

  TFile* outFile = new TFile("Rvn2pi.root","RECREATE");
  outFile->cd();

  CCMCME.ToTH2(ourPOT)->Write("ME_Dist_CCMC");
  NCMCME.ToTH2(ourPOT)->Write("ME_Dist_NCMC");
  NEMCME.ToTH2(ourPOT)->Write("ME_Dist_NuEMC");
  CCDataME.ToTH2(ourPOT)->Write("ME_Dist_CCData");
  NCDataME.ToTH2(ourPOT)->Write("ME_Dist_NCData");
  NEDataME.ToTH2(ourPOT)->Write("ME_Dist_NuEData");

  obsBlah.ToTH1(ourPOT)->Write("obsDCMP");
  obsNominal.ToTH1(ourPOT)->Write("obsNominal");
  obsShifts.ToTH1(ourPOT)->Write("obsShifted");

  bkgdDCMP.ToTH1(ourPOT)->Write("bkgdDCMP");
  bkgdNominal.ToTH1(ourPOT)->Write("bkgdNominal");
  bkgdShifted.ToTH1(ourPOT)->Write("bkgdShifted");

  CCMC.ToTH1(ourPOT)->Write("ND_CCMC");
  CCData.ToTH1(ourPOT)->Write("ND_CCData");
  CCDCMP.ToTH1(ourPOT)->Write("ND_CCDCMP");

  NCMC.ToTH1(ourPOT)->Write("ND_NCMC");
  NCData.ToTH1(ourPOT)->Write("ND_NCData");
  NCDCMP.ToTH1(ourPOT)->Write("ND_NCDCMP");

  NEMC.ToTH1(ourPOT)->Write("ND_NEMC");
  NEData.ToTH1(ourPOT)->Write("ND_NEData");
  NEDCMP.ToTH1(ourPOT)->Write("ND_NEDCMP");

}