Toy_analyses.C File Reference

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "3FlavorAna/Cuts/QuantileCuts.h"
#include "CAFAna/Core/EventList.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Experiment/AtmConstraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/ReactorExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/SolarConstraints.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Extrap/ModularExtrapComponent.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "3FlavorAna/Systs/EnergySysts.h"
#include "3FlavorAna/Systs/NumuSysts.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "CAFAna/Vars/FitVars.h"
#include "3FlavorAna/Vars/FitVarsNumu.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NumuVars.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "OscLib/OscCalcPMNSOpt.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TRandom3.h"
#include "TTree.h"

Go to the source code of this file.

Classes
struct	FitResult

Functions
TRandom3	r (0)
FitResult	getBestFit (IPrediction *prediction, TH1F *h1, TH1D *hCosmic)
FitResult	getBestFit (IPrediction *quant1, IPrediction *quant2, IPrediction *quant3, IPrediction *quant4, TH1F *h1, TH1F *h2, TH1F *h3, TH1F *h4, TH1D *hCosmic1, TH1D *hCosmic2, TH1D *hCosmic3, TH1D *hCosmic4)
std::vector< std::vector< double > >	QuantileVectorFromTH2 (TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles)
unsigned int	quantile (std::vector< std::vector< double >> v, float E, float hadEfrac)
void	getPoissonRandom (TH1D *s, TH1F *h)
void	getSArebin (TH1D *s, TH1D *h)
void	getSArebin (TH1F *s, TH1F *h)
void	Toy_analyses (unsigned int nExperiments, unsigned int volume=1)

Function Documentation

FitResult getBestFit	(	IPrediction *	prediction,
		TH1F *	h1,
		TH1D *	hCosmic
	)

Definition at line 432 of file Toy_analyses.C.

References ana::IFitter::Fit(), osc::_IOscCalcAdjustable< T >::GetDmsq32(), ana::getNumuMarginalisedOscParam(), osc::_IOscCalcAdjustable< T >::GetTh23(), ana::kFitDmSq32Scaled, ana::kFitSinSqTheta23, ana::IFitter::kQuiet, ana::kSecondAnaLivetime, ana::kSecondAnaPOT, ana::kSolarConstraintsPDG2017, cet::pow(), ana::ResetOscCalcToDefault(), fillBadChanDBTables::result, std::sqrt(), and ana::WorldReactorConstraint2017().

Referenced by Toy_analyses().

{  
  osc::OscCalcPMNSOpt calc2;
  ana::ResetOscCalcToDefault(&calc2);

  Spectrum fakeData(h1, kSecondAnaPOT, kSecondAnaLivetime);

  // Experiment                                                                                     
  SingleSampleExperiment fakeExpt1(prediction, fakeData, hCosmic);

  std::vector<const IExperiment*> vExpt;
  vExpt.push_back(&fakeExpt1);
  vExpt.push_back(&kSolarConstraintsPDG2017);
  vExpt.push_back(WorldReactorConstraint2017());
  MultiExperiment fakeExpt(vExpt);

  std::vector<const IFitVar*> margPars_ssth23 = getNumuMarginalisedOscParam();
  margPars_ssth23.push_back(&kFitDmSq32Scaled);

  MinuitFitter fitMaxMix(&fakeExpt, margPars_ssth23);
  double MMLL = fitMaxMix.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();

  margPars_ssth23.push_back(&kFitSinSqTheta23);
  MinuitFitter fit(&fakeExpt, margPars_ssth23);
  double bestLL = fit.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();

  double bestX = pow(TMath::Sin(calc2.GetTh23()),2);
  double bestY = 1e3 * calc2.GetDmsq32();

  double MMrej = 0;
  if(MMLL > bestLL) MMrej = sqrt(MMLL-bestLL); // this should always happen, but just in case

  FitResult result = {bestX, bestY, bestLL, MMrej};

  return result;
}

FitResult getBestFit	(	IPrediction *	quant1,
		IPrediction *	quant2,
		IPrediction *	quant3,
		IPrediction *	quant4,
		TH1F *	h1,
		TH1F *	h2,
		TH1F *	h3,
		TH1F *	h4,
		TH1D *	hCosmic1,
		TH1D *	hCosmic2,
		TH1D *	hCosmic3,
		TH1D *	hCosmic4
	)

Definition at line 384 of file Toy_analyses.C.

References ana::IFitter::Fit(), osc::_IOscCalcAdjustable< T >::GetDmsq32(), ana::getNumuMarginalisedOscParam(), osc::_IOscCalcAdjustable< T >::GetTh23(), ana::kFitDmSq32Scaled, ana::kFitSinSqTheta23, ana::IFitter::kQuiet, ana::kSecondAnaLivetime, ana::kSecondAnaPOT, ana::kSolarConstraintsPDG2017, cet::pow(), ana::ResetOscCalcToDefault(), fillBadChanDBTables::result, std::sqrt(), and ana::WorldReactorConstraint2017().

{
  osc::OscCalcPMNSOpt calc2;
  ana::ResetOscCalcToDefault(&calc2);

  Spectrum fakeData1(h1, kSecondAnaPOT, kSecondAnaLivetime);
  Spectrum fakeData2(h2, kSecondAnaPOT, kSecondAnaLivetime);
  Spectrum fakeData3(h3, kSecondAnaPOT, kSecondAnaLivetime);
  Spectrum fakeData4(h4, kSecondAnaPOT, kSecondAnaLivetime);

  SingleSampleExperiment fakeExpt1 (quant1, fakeData1, hCosmic1);
  SingleSampleExperiment fakeExpt2 (quant2, fakeData2, hCosmic2);
  SingleSampleExperiment fakeExpt3 (quant3, fakeData3, hCosmic3);
  SingleSampleExperiment fakeExpt4 (quant4, fakeData4, hCosmic4);

  std::vector<const IExperiment*> vExpt;
  vExpt.push_back(&fakeExpt1);
  vExpt.push_back(&fakeExpt2);
  vExpt.push_back(&fakeExpt3);
  vExpt.push_back(&fakeExpt4);
  vExpt.push_back(&kSolarConstraintsPDG2017);
  vExpt.push_back(WorldReactorConstraint2017());
  MultiExperiment fakeExpt(vExpt);

  std::vector<const IFitVar*> margPars_ssth23 = getNumuMarginalisedOscParam();
  margPars_ssth23.push_back(&kFitDmSq32Scaled);
 
  MinuitFitter fitMaxMix(&fakeExpt, margPars_ssth23);
  double MMLL = fitMaxMix.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();

  margPars_ssth23.push_back(&kFitSinSqTheta23);
  MinuitFitter fit(&fakeExpt, margPars_ssth23);
  double bestLL = fit.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();

  double bestX = pow(TMath::Sin(calc2.GetTh23()),2);
  double bestY = 1e3 * calc2.GetDmsq32();

  double MMrej = 0;
  if(MMLL > bestLL) MMrej = sqrt(MMLL-bestLL); // this should always happen, but just in case

  FitResult result = {bestX, bestY, bestLL, MMrej};

  return result;

}

void getPoissonRandom	(	TH1D *	s,
		TH1F *	h
	)

Definition at line 504 of file Toy_analyses.C.

References MECModelEnuComparisons::i, and r().

Referenced by Toy_analyses().

{
  for(int i = 1; i <= s->GetNbinsX(); i++)
  {
    float fX = s->GetBinCenter(i);
    float fY = r.Poisson(s->GetBinContent(i));

    h->Fill(fX, fY);
  }
}

void getSArebin	(	TH1D *	s,
		TH1D *	h
	)

Definition at line 515 of file Toy_analyses.C.

References MECModelEnuComparisons::i.

Referenced by Toy_analyses().

{
  for(int i = 1; i <= s->GetNbinsX(); i++)
  {
    double fX = s->GetBinCenter(i);
    double fY = s->GetBinContent(i);

    h->Fill(fX, fY);
  }
}

void getSArebin	(	TH1F *	s,
		TH1F *	h
	)

Definition at line 526 of file Toy_analyses.C.

References MECModelEnuComparisons::i.

{
  for(int i = 1; i <= s->GetNbinsX(); i++)
  {
    float fX = s->GetBinCenter(i);
    float fY = s->GetBinContent(i);

    h->Fill(fX, fY);
  }
}

unsigned int quantile	(	std::vector< std::vector< double >>	v,
		float	E,
		float	hadEfrac
	)

Definition at line 480 of file Toy_analyses.C.

References submit_hadd::l, and registry_explorer::v.

Referenced by MakeCovarSim(), PlotDataQA::MakePlots(), ana::MCMCSamples::NumSamples(), stan::mcmc::chains< RNG >::quantile(), ana::QuantileAndPIDCutsFromTH2(), ana::QuantileCuts(), ana::QuantileCutsFromTH2(), ana::MCMCSamples::QuantileLL(), stan::mcmc::chains< RNG >::quantiles(), ana::BayesianMarginal::QuantileThreshold(), ana::BayesianMarginal::ThresholdFromTH1(), Toy_analyses(), and ana::BayesianMarginal::Vars().

{
  float A17Ebins[20] = { 0, 0.75, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
                        1.9, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5 };
  unsigned int quant = 0;

  unsigned int k = 0;

  while (A17Ebins[k] < E && k < 20)
    //  for(unsigned int k = 0; k < 20; k++)
  {
    unsigned int l = 0;
    while (v.at(k)[l] < hadEfrac && l < v.at(k).size())
    //    for(unsigned int l = 0; l < v.at(k).size(); l++)
      {
        quant = l;
        l++;
      }
    k++;
  }

  return quant;
}

std::vector< std::vector< double > > QuantileVectorFromTH2	(	TH2 *	quantileHist,
		const HistAxis &	independentAxis,
		const HistAxis &	quantileAxis,
		const unsigned int &	numQuantiles
	)

Definition at line 469 of file Toy_analyses.C.

References ana::GetQuantileBins().

Referenced by Toy_analyses().

                                                                        {

  std::vector<std::vector<double>> quantileBins = GetQuantileBins(quantileHist, 
    independentAxis, quantileAxis, numQuantiles);

  return quantileBins;
}

TRandom3 r

(

0

)

Referenced by util::HoughCalc::AddPoint(), genie::geometry::ROOTGeomAnalyzer::AdoptGeomVolSelector(), caf::SRProxySystController::AnyShifted(), art::AssnsNode< L, R, D >::AssnsNode(), osc::EarthModel::AveNe(), hough::Hough2P::BuildMap(), hough::MultiHough2P::BuildMap(), BuildSpectrum(), nova::database::Table::BulkInsertInDB(), genie::BardinIMDRadCorPXSec::C(), slid::ParticleIDAlg::CalcAsymIneria(), CalcFiberLoopCorr(), CalcRWithSysts(), CalcRWithSystsNus17(), CalcSignificance(), genie::QPMDISStrucFuncBase::Calculate(), CheckUnitarityLimit(), CheckVertexDistribution(), genie::NievesQELCCPXSec::CompareNievesTensors(), genie::HadronTransporter::Configure(), genie::EmpiricalMECPXSec2015::Configure(), genie::QPMDISPXSec::Configure(), genie::QPMDMDISPXSec::Configure(), genie::SmithMonizQELCCXSec::Configure(), genie::QELEventGeneratorSM::Configure(), genie::NuclearModelMap::Configure(), genie::SmithMonizQELCCPXSec::Configure(), nova::database::Table::CreateFromXML(), nova::dbi::Table::CreateFromXML(), nova::dbi::dBTable_(), nova::database::dBTable_(), genie::alvarezruso::AREikonalSolution::Deltamed(), genie::alvarezruso::AlvarezRusoCOHPiPDXSec::DeltaWidthPauliBlocked(), genie::alvarezruso::ARSampledNucleus::Density0FunctionFermiLiquid(), osc::EarthModel::DensityStacey(), util::HoughCalc::DistHit(), CLHEP::RandLandau::distributionName(), CLHEP::RandGaussQ::distributionName(), art::ResultsProducer::doBeginRun(), art::ResultsProducer::doEndRun(), evd::PlotView::Draw(), ana::CovMxSurface::DrawSurfacePoint(), DrawSurfacePoint(), genie::alvarezruso::AREikonalSolution::Element(), comi::CosmicMetrics::endJob(), comi::DataCheck::endJob(), boost::python::detail::proxy_links< Proxy, Container >::erase(), rsim::ExcessNoiseMaker::ExcessNoisePDF(), genie::BardinIMDRadCorPXSec::Fa(), om::Icons::FastForward(), genie::GiBUURESFormFactor::FormFactors::FFRes(), art::OutputModule::fileGranularity(), art::Assns< L, R, void >::fill_transients(), genie::alvarezruso::ARSampledNucleus::FillDensities(), calib::DCMTimingOffset::filter(), boost::python::detail::proxy_links< Proxy, Container >::find(), bpfit::Path::FindLayers(), hough::Hough2P::FindPeaks(), ana::PredictionInterp::FitComponent(), cet::sqlite::Connection::flush_no_throw(), om::Icons::FolderExplore(), fragmentAngularDistribution(), fragmentAngularDistributionGM(), fragmentEnergyDistributionDifferentAngles(), fragmentYields(), fragmentYieldsPlot(), stan::math::gaussian_dlm_obs_lpdf(), genie::CharmHadronization::GenerateCharmHadron(), genie::NucleonDecayPrimaryVtxGenerator::GenerateDecayedNucleonPosition(), GenerateEventsAtFixedInitState(), genie::flux::GJPARCNuFlux::GenerateNext(), genie::flux::GDk2NuFlux::GenerateNext(), genie::flux::GSimpleNtpFlux::GenerateNext(), genie::flux::GNuMIFlux::GenerateNext(), genie::NNBarOscPrimaryVtxGenerator::GenerateOscillatingNeutronPosition(), genie::VertexGenerator::GenerateVertex(), dbutils::hw::FEB::GetAllInstalled(), dbutils::hw::APD::GetAllInstalled(), novaddt::MoonShadow::GetAngularPrescaleDecision(), moonshadowana::MoonShadowAna::GetAngularPrescaleDecision(), moonshadowana::MoonShadowAnaHough::GetAngularPrescaleDecision(), GetCommandLineArgs(), genie::Algorithm::GetConfig(), getPoissonRandom(), GetRandomRay(), novaddt::MoonShadow::GetSmartPrescaleDecision(), novaddt::MichelETrigger::GetSmartPrescaleDecision(), moonshadowana::MoonShadowAna::GetSmartPrescaleDecision(), moonshadowana::MoonShadowAnaHough::GetSmartPrescaleDecision(), evgen::NuMISpillTimeStructure::GetTime(), genie::KPhaseSpace::GetTMaxDFR(), genie::HAIntranuke::HadronFateHA(), genie::HAIntranuke2018::HadronFateHA(), genie::HNIntranuke2018::HadronFateHN(), om::HeaderBar::HandleTimer(), om::Icons::Histo(), art::RunID::inRangeOrInvalid(), dbutils::hw::FEB::Installed(), dbutils::hw::APD::Installed(), inTorus(), stan::math::inv_Phi(), rawdata::RawDigit::IsMC(), novaddt::NNFastMMTrigger::LinFit(), zcl::FmmTrackerValidation::LinFit(), novaddt::UpMuTrigger::LinFit(), novaddt::UpMuAna::LinFit(), zcl::FastMMStudy::LinFit(), upmuana::UpMuAnalysis::LinFit(), upmuana::UpMuRecoAna::LinFit(), genie::NievesQELCCPXSec::LmunuAnumu(), keras::LayerConv2D::load_weights(), jmshower::RecoJMShower::LoadDedxHistograms(), nova::database::Table::LoadFromCSV(), nutools::dbi::Table::LoadFromCSV(), nova::database::Table::LoadFromSSV(), genie::GiBUURESFormFactor::FormFactors::LoadTables(), lu_backsubst(), lu_factorize(), lu_forwsubst(), LZ4_NbCommonBytes(), main(), make_extrap_figure(), cmf::CovarianceFitHelper::MakeCovarianceMatrixHistogram(), evdb::ColorScale::MakeHSVScale(), om::PlotOptions::MakeLabelText(), osc::EarthModel::MakeLayers(), trk::WindowTrackingAlg::MakeTrack(), cet::sqlite::max(), evd::SimulationDrawer::MCTruthVectors2D(), evd::SimulationDrawer::MCTruthVectors3D(), cet::sqlite::min(), genie::LocalFGM::ModelType(), hough::MultiHough2P::MultiMap(), util::HoughCalc::newNode(), nova::database::Table::NewRow(), nova::dbi::Table::NewRow(), nutools::dbi::Table::NewRow(), genie::alvarezruso::AlvarezRusoCOHPiPDXSec::NuclearCurrent(), CLHEP::operator!=(), fhicl::detail::it_value_get< std::complex< U >, typename tt::disable_if< std::is_reference< std::complex< U > >::value|| std::is_pointer< std::complex< U > >::value >::type >::operator()(), art::const_AssnsIter< L, R, D, Dir >::operator<(), CLHEP::operator<(), CLHEP::operator<<(), art::const_AssnsIter< L, R, D, Dir >::operator<=(), CLHEP::operator<=(), CLHEP::operator==(), art::const_AssnsIter< L, R, D, Dir >::operator>(), CLHEP::operator>(), art::const_AssnsIter< L, R, D, Dir >::operator>=(), CLHEP::operator>=(), INukeNucleonCorr::OutputFiles(), genie::NucDeExcitationSim::OxygenTargetSim(), nova::database::dBColumn_t::parse(), nova::dbi::dBRow_t::parse(), nova::database::dBRow_t::parse(), nova::dbi::dBTable::parse(), nova::database::dBTable::parse(), om::Icons::Pause(), genie::HAIntranuke::PiBounce(), genie::HAIntranuke2018::PiBounce(), om::Icons::Play(), TMVA::TMVAGlob::plot_logo(), art::DataViewImpl::PMValue::PMValue(), genie::HAIntranuke::PnBounce(), genie::HAIntranuke2018::PnBounce(), art::TimeTracker::postEndJob(), Print(), genie::Algorithm::Print(), genie::LocalFGM::ProbDistro(), geo::LiveGeometry::ProjectedDistance(), nova::dbi::RunHistory::PullAPDInfoFromDB(), nova::dbi::RunHistory::PullFEBInfoFromDB(), quantiles(), lem::dec::Forest::RandomOrthoMatrix(), keras::DataChunk2D::read_from_file(), art::EmptyEvent::readRun_(), jmshower::RecoJMShower::RecoShowers(), ana::RefineSeeds(), boost::python::detail::proxy_links< Proxy, Container >::remove(), util::HoughCalc::RemovePoint(), om::Icons::Repeat(), boost::python::detail::proxy_links< Proxy, Container >::replace(), om::Icons::Rewind(), stan::model::rvalue(), genie::alvarezruso::ARSampledNucleus::SamplePoint2(), art::SamplingInput::SamplingInput(), ana::FluxDecomp::SaveTo(), SaveToRootFile(), SelectInitState(), art::RangeSet::set_run(), evdb::ColorScale::SetBounds(), nova::dbi::RunHistory::SetFDRHCRunStart(), nova::dbi::RunHistory::SetFDRHCRunStop(), genie::Target::SetHitNucPosition(), nova::dbi::RunHistory::SetNDRHCRunStart(), nova::dbi::RunHistory::SetNDRHCRunStop(), novaddt::SetRho(), earms::GridSearch::SetStandardDirections(), rb::RecoHit::SetT(), nova::database::Table::SetValidityRange(), beamlinereco::CFDHitFinder< T >::SetWaveform(), gibuu::SimpleRecordList::SimpleRecordList(), genie::DISHadronicSystemGenerator::SimulateFormationZone(), art::EventRange::size(), osc::analytic::SolveCubic(), genie::utils::intranuke::StepParticle(), genie::utils::intranuke2018::StepParticle(), stub_reader(), stan::math::student_t_cdf(), stan::math::student_t_lccdf(), stan::math::student_t_lcdf(), TEST(), util::HoughCalc::TestHit(), testReconfigInCommonPool(), testReconfigInOwnedModules(), art::TFileService::TFileService(), Toy_analyses(), red::OverburdenCalculator::traceBack(), evdb::Colors::UnpackColorScale(), nuonecvntf::NuonECVNTF::vector_to_tensor(), regcvntf::RegCVNTF::vector_to_tensor_nue(), stan::math::von_mises_rng(), cmf::CMFDecorrelator::writeResults(), nova::database::Table::WriteToDB(), nova::dbi::Table::WriteToDB(), nutools::dbi::Table::WriteToDB(), genie::BardinIMDRadCorPXSec::XSec(), genie::NievesQELCCPXSec::XSec(), zheevq3(), testbeam::BeamlineEventDisplay::~BeamlineEventDisplay(), genie::GiBUURESFormFactor::FormFactors::~FormFactors(), and genie::flux::GDk2NuFluxXMLHelper::~GDk2NuFluxXMLHelper().

void Toy_analyses	(	unsigned int	nExperiments,
		unsigned int	volume = 1
	)

Definition at line 79 of file Toy_analyses.C.

References std::asin(), FitResult::bestX, FitResult::bestY, calc, om::cout, e, allTimeWatchdog::endl, getBestFit(), getPoissonRandom(), getSArebin(), MECModelEnuComparisons::i, ana::kAna2017Livetime, ana::kHadEFracAxis, ana::kNumuCCOptimisedAxis, ana::kSecondAnaLivetime, ana::kSecondAnaPOT, FitResult::LL, nPoints, osc::_OscCalcPMNSOpt< T >::P(), ana::pnfs2xrootd(), ana::PredictionInterp::Predict(), quantile(), QuantileVectorFromTH2(), r(), runNovaSAM::release, ana::ResetOscCalcToDefault(), osc::_OscCalcPMNSOpt< T >::SetDmsq32(), osc::_OscCalcPMNSOpt< T >::SetTh23(), FitResult::sigma, std::sqrt(), string, and ana::Spectrum::ToTH1().

{
  std::string inputDir = "/pnfs/nova/persistent/users/bzamoran/ToyPvalue/";

  TFile* fInputTree = TFile::Open(pnfs2xrootd(inputDir + "inputTree.root").c_str());
  TTree* inputTree = (TTree*) fInputTree->Get("inputTree");
  float TrueE, CCE, CCESA, HadEfrac, PPFX, CXW;
  int isSA, isA17;
  inputTree->SetBranchAddress("TrueE", &TrueE);
  inputTree->SetBranchAddress("isSA", &isSA);
  inputTree->SetBranchAddress("isA17", &isA17);
  inputTree->SetBranchAddress("CCE", &CCE);
  inputTree->SetBranchAddress("CCESA", &CCESA);
  inputTree->SetBranchAddress("HadEfrac", &HadEfrac);
  inputTree->SetBranchAddress("PPFX", &PPFX);
  inputTree->SetBranchAddress("CXW", &CXW);

  std::cout << "Tree" << std::endl;

  TFile* fPred = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_AllQuants.root").c_str());
  PredictionInterp* prediction_OR = LoadFrom<PredictionInterp>(fPred,  "prediction_OR" ).release();
  PredictionInterp* prediction_SA = LoadFrom<PredictionInterp>(fPred,  "prediction_SA" ).release();
  fPred->Close();

  std::cout << "Predictions 1" << std::endl;

  TFile* fPred_Q1 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant1.root").c_str());
  TFile* fPred_Q2 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant2.root").c_str());
  TFile* fPred_Q3 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant3.root").c_str());
  TFile* fPred_Q4 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant4.root").c_str());

  PredictionInterp* prediction_Q1 = LoadFrom<PredictionInterp>(fPred_Q1, "prediction").release();
  PredictionInterp* prediction_Q2 = LoadFrom<PredictionInterp>(fPred_Q2, "prediction").release();
  PredictionInterp* prediction_Q3 = LoadFrom<PredictionInterp>(fPred_Q3, "prediction").release();
  PredictionInterp* prediction_Q4 = LoadFrom<PredictionInterp>(fPred_Q4, "prediction").release();

  std::cout << "Predictions quant" << std::endl;

  TFile* fCosmics = TFile::Open(pnfs2xrootd(inputDir + "cos_bkg_hists_v2.root").c_str());
  TH1D* hSpecCosmics_Q1 = (TH1D*) fCosmics->Get("q1all");
  TH1D* hSpecCosmics_Q2 = (TH1D*) fCosmics->Get("q2all");
  TH1D* hSpecCosmics_Q3 = (TH1D*) fCosmics->Get("q3all");
  TH1D* hSpecCosmics_Q4 = (TH1D*) fCosmics->Get("q4all");

  hSpecCosmics_Q1->Scale(kSecondAnaLivetime / kAna2017Livetime ); // scale to 6e20
  hSpecCosmics_Q2->Scale(kSecondAnaLivetime / kAna2017Livetime );
  hSpecCosmics_Q3->Scale(kSecondAnaLivetime / kAna2017Livetime );
  hSpecCosmics_Q4->Scale(kSecondAnaLivetime / kAna2017Livetime );

  TF1* fSA = new TF1("fSA","TMath::Exp(1.976*(x-2.26))/(1+TMath::Exp(1.976*(x-2.26)))",0,5);
  TF1* fA17 = new TF1("fA17","TMath::Exp(-1.907*(x-1.967))/(1+TMath::Exp(-1.907*(x-1.967)))",0,5);
  TF1* fBoth = new TF1("fBoth","1.-fA17-fSA",0,5);

  TF1* fAdjustSA = new TF1("fAdjustSA","TMath::Min(fSA/(1.-fSA),5.0)", 0.5, 5); // don't let it blow
  TF1* fAdjustA17 = new TF1("fAdjustA17","fA17/(1.-fSA)", 0.5, 5);
  TF1* fAdjustBoth = new TF1("fAdjustBoth","fBoth/(1.-fSA)", 0.5, 5);
  TF1* fAdjustSAandBoth = new TF1("fAdjustSAandBoth","TMath::Min((1.-fA17)/(1.-fSA),5.0)", 0.5, 5);

  TH1D* hSpecCosmics_SA = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_SA");
  hSpecCosmics_SA->Add(hSpecCosmics_Q2);
  hSpecCosmics_SA->Add(hSpecCosmics_Q3);
  hSpecCosmics_SA->Add(hSpecCosmics_Q4);
  hSpecCosmics_SA->Multiply(fAdjustSA);
  TH1D* hSpecCosmics_SA_rebin = new TH1D("hSpecCosmics_SA_rebin","",20,0,5);
  getSArebin(hSpecCosmics_SA, hSpecCosmics_SA_rebin);

  TH1D* hSpecCosmics_SA_and_both = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_SA_and_both");
  hSpecCosmics_SA_and_both->Add(hSpecCosmics_Q2);
  hSpecCosmics_SA_and_both->Add(hSpecCosmics_Q3);
  hSpecCosmics_SA_and_both->Add(hSpecCosmics_Q4);
  hSpecCosmics_SA_and_both->Multiply(fAdjustSAandBoth);
  TH1D* hSpecCosmics_SA_and_both_rebin = new TH1D("hSpecCosmics_SA_and_both_rebin","",20,0,5);
  getSArebin(hSpecCosmics_SA_and_both, hSpecCosmics_SA_and_both_rebin);
  
  TH1D* hSpecCosmics_Both_Q1 = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_Both_Q1");
  TH1D* hSpecCosmics_Both_Q2 = (TH1D*) hSpecCosmics_Q2->Clone("hSpecCosmics_Both_Q2");
  TH1D* hSpecCosmics_Both_Q3 = (TH1D*) hSpecCosmics_Q3->Clone("hSpecCosmics_Both_Q3");
  TH1D* hSpecCosmics_Both_Q4 = (TH1D*) hSpecCosmics_Q4->Clone("hSpecCosmics_Both_Q4");
  hSpecCosmics_Both_Q1->Multiply(fAdjustBoth);
  hSpecCosmics_Both_Q2->Multiply(fAdjustBoth);
  hSpecCosmics_Both_Q3->Multiply(fAdjustBoth);
  hSpecCosmics_Both_Q4->Multiply(fAdjustBoth);

  TH1D* hSpecCosmics_A17_Q1 = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_A17_Q1");
  TH1D* hSpecCosmics_A17_Q2 = (TH1D*) hSpecCosmics_Q2->Clone("hSpecCosmics_A17_Q2");
  TH1D* hSpecCosmics_A17_Q3 = (TH1D*) hSpecCosmics_Q3->Clone("hSpecCosmics_A17_Q3");
  TH1D* hSpecCosmics_A17_Q4 = (TH1D*) hSpecCosmics_Q4->Clone("hSpecCosmics_A17_Q4");
  hSpecCosmics_A17_Q1->Multiply(fAdjustA17);
  hSpecCosmics_A17_Q2->Multiply(fAdjustA17);
  hSpecCosmics_A17_Q3->Multiply(fAdjustA17);
  hSpecCosmics_A17_Q4->Multiply(fAdjustA17);


  std::cout << "All's well with loading" << std::endl;

  // Get FD files for HadEFrac vs. neutrino energy distribution
  std::string quantpath = "/pnfs/nova/persistent/analysis/numu/Ana2017/all_FD_histo_for_quantile_cuts.root";
  TFile* quantfile = TFile::Open(pnfs2xrootd(quantpath).c_str());
  TH2* quanthist = (TH2*)quantfile->Get("dir_FDSpec2D/FDSpec2D");
  std::vector<std::vector<double>> QuantBins = QuantileVectorFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);

  // Now the fittting

  osc::OscCalcPMNSOpt calc;
  ana::ResetOscCalcToDefault(&calc);
  // calc.SetDmsq32(2.67e-3);
  // calc.SetTh23(asin(sqrt(0.404))); // PRL-published Second Ana
  calc.SetDmsq32(2.43071e-3);
  calc.SetTh23(asin(sqrt(0.471974)));

  TFile* fOutput = new TFile(Form("Toy_Output_vol%i.root", volume),"RECREATE");
  TTree* fTree = new TTree("T","Summary tree from p-value study");

  TH1F* hSA = new TH1F("hSA","",20,0,5);
  TH1F* hSA_cosmic = new TH1F("hSA_cosmic","",20,0,5);

  float A17Ebins[20] = { 0, 0.75, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
                        1.9, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5 };

  TH1F* hA17 = new TH1F("hA17","",19, A17Ebins);
  TH1F* hA17_Q1 = new TH1F("hA17_Q1","",19, A17Ebins);
  TH1F* hA17_Q2 = new TH1F("hA17_Q2","",19, A17Ebins);
  TH1F* hA17_Q3 = new TH1F("hA17_Q3","",19, A17Ebins);
  TH1F* hA17_Q4 = new TH1F("hA17_Q4","",19, A17Ebins);

  TH1F* hA17_cosmic = new TH1F("hA17_cosmic","",19, A17Ebins);
  TH1F* hA17_Q1_cosmic = new TH1F("hA17_Q1_cosmic","",19, A17Ebins);
  TH1F* hA17_Q2_cosmic = new TH1F("hA17_Q2_cosmic","",19, A17Ebins);
  TH1F* hA17_Q3_cosmic = new TH1F("hA17_Q3_cosmic","",19, A17Ebins);
  TH1F* hA17_Q4_cosmic = new TH1F("hA17_Q4_cosmic","",19, A17Ebins);

  TH1F* hBoth = new TH1F("hBoth","",19, A17Ebins);

  TH1F* hBoth_cosmic_Q1 = new TH1F("hBoth_cosmic_Q1","",19, A17Ebins);
  TH1F* hBoth_cosmic_Q2 = new TH1F("hBoth_cosmic_Q2","",19, A17Ebins);
  TH1F* hBoth_cosmic_Q3 = new TH1F("hBoth_cosmic_Q3","",19, A17Ebins);
  TH1F* hBoth_cosmic_Q4 = new TH1F("hBoth_cosmic_Q4","",19, A17Ebins);
  TH1F* hBoth_cosmic = new TH1F("hBoth_cosmic","",19, A17Ebins);
  TH1F* hBoth_cosmic_rebin = new TH1F("hA17_Both_cosmic_rebin","",20,0,5);

  double sigmaMaxSA, sigmaMaxA17, LLSA, LLA17;
  double dmsq32SA, dmsq32A17, ssth23SA, ssth23A17;
  unsigned int nPoints;

  fTree->Branch("sigmaMaxSA", &sigmaMaxSA);
  fTree->Branch("sigmaMaxA17", &sigmaMaxA17);
  fTree->Branch("LLSA", &LLSA);
  fTree->Branch("LLA17", &LLA17);
  fTree->Branch("dmsq32SA", &dmsq32SA);
  fTree->Branch("dmsq32A17", &dmsq32A17);
  fTree->Branch("ssth23SA", &ssth23SA);
  fTree->Branch("ssth23A17", &ssth23A17);
  fTree->Branch("Spectrum_SA", &hSA);
  fTree->Branch("Cosmics_SA", &hSA_cosmic);
  fTree->Branch("Spectrum_A17", &hA17);
  fTree->Branch("Spectrum_A17_Q1", &hA17_Q1);
  fTree->Branch("Spectrum_A17_Q2", &hA17_Q2);
  fTree->Branch("Spectrum_A17_Q3", &hA17_Q3);
  fTree->Branch("Spectrum_A17_Q4", &hA17_Q4);
  fTree->Branch("Cosmics_A17", &hA17_cosmic);
  fTree->Branch("Cosmics_A17_Q1", &hA17_Q1_cosmic);
  fTree->Branch("Cosmics_A17_Q2", &hA17_Q2_cosmic);
  fTree->Branch("Cosmics_A17_Q3", &hA17_Q3_cosmic);
  fTree->Branch("Cosmics_A17_Q4", &hA17_Q4_cosmic);
  fTree->Branch("Spectrum_Both", &hBoth);
  fTree->Branch("Cosmics_Both", &hBoth_cosmic);

  unsigned int nEntries = inputTree->GetEntries();

  double MCpoints = prediction_OR->Predict(&calc).ToTH1(kSecondAnaPOT)->Integral();

  for(unsigned int i = 0; i < nExperiments; i++)
  {
    hSA->Reset();
    hSA_cosmic->Reset();
    hA17->Reset();
    hA17_Q1->Reset();
    hA17_Q2->Reset();
    hA17_Q3->Reset();
    hA17_Q4->Reset();
    hA17_cosmic->Reset();
    hA17_Q1_cosmic->Reset();
    hA17_Q2_cosmic->Reset();
    hA17_Q3_cosmic->Reset();
    hA17_Q4_cosmic->Reset();
    hBoth->Reset();
    hBoth_cosmic_Q1->Reset();
    hBoth_cosmic_Q2->Reset();
    hBoth_cosmic_Q3->Reset();
    hBoth_cosmic_Q4->Reset();
    hBoth_cosmic->Reset();
    hBoth_cosmic_rebin->Reset();

    nPoints = r.Poisson(MCpoints);

    unsigned int k = 0;
    while(k <= nPoints)
    {
      unsigned int iEntry = r.Integer(nEntries);
      inputTree->GetEntry(iEntry);

      double p = r.Uniform();
      double prob = calc.P(14, 14, TrueE) * TMath::Min(CXW * PPFX * 0.5, 1.); // adjusting for weights

      if(p < prob && TrueE > 0.0 && TrueE < 5.0)
      { 
        if(isSA)
        {
          hSA->Fill(CCE);
        }
        if(isA17)
          {
            int qi = quantile(QuantBins, CCE, HadEfrac);
            if(qi == 0) hA17_Q1->Fill(CCE);
            else if(qi == 1) hA17_Q2->Fill(CCE);
            else if(qi == 2) hA17_Q3->Fill(CCE);
            else hA17_Q4->Fill(CCE);
          }
        if(isA17 && isSA) hBoth->Fill(CCE);
        if(isSA || isA17) k++;
      }
    } // loop over points

    // Protection against errors reading the tree
    if(hSA->FindFirstBinAbove(1) == hSA->FindLastBinAbove(1)) continue;
    if(hBoth->FindFirstBinAbove(1) == hBoth->FindLastBinAbove(1)) continue;
    if(hA17_Q1->FindFirstBinAbove(1) == hA17_Q1->FindLastBinAbove(1)) continue;
    if(hA17_Q2->FindFirstBinAbove(1) == hA17_Q2->FindLastBinAbove(1)) continue;
    if(hA17_Q3->FindFirstBinAbove(1) == hA17_Q3->FindLastBinAbove(1)) continue;
    if(hA17_Q4->FindFirstBinAbove(1) == hA17_Q4->FindLastBinAbove(1)) continue;

    // Now handle the cosmics
    getPoissonRandom(hSpecCosmics_Both_Q1, hBoth_cosmic_Q1);
    getPoissonRandom(hSpecCosmics_Both_Q2, hBoth_cosmic_Q2);
    getPoissonRandom(hSpecCosmics_Both_Q3, hBoth_cosmic_Q3);
    getPoissonRandom(hSpecCosmics_Both_Q4, hBoth_cosmic_Q4);

    hBoth_cosmic->Add(hBoth_cosmic_Q1);
    hBoth_cosmic->Add(hBoth_cosmic_Q2);
    hBoth_cosmic->Add(hBoth_cosmic_Q3);
    hBoth_cosmic->Add(hBoth_cosmic_Q4);

    hBoth->Add(hBoth_cosmic);

    getSArebin(hBoth_cosmic, hBoth_cosmic_rebin);

    getPoissonRandom(hSpecCosmics_SA_rebin, hSA_cosmic);
    hSA_cosmic->Add(hBoth_cosmic_rebin);
    hSA->Add(hSA_cosmic);

    getPoissonRandom(hSpecCosmics_A17_Q1, hA17_Q1_cosmic);
    getPoissonRandom(hSpecCosmics_A17_Q2, hA17_Q2_cosmic);
    getPoissonRandom(hSpecCosmics_A17_Q3, hA17_Q3_cosmic);
    getPoissonRandom(hSpecCosmics_A17_Q4, hA17_Q4_cosmic);

    hA17_Q1_cosmic->Add(hBoth_cosmic_Q1);
    hA17_Q2_cosmic->Add(hBoth_cosmic_Q2);
    hA17_Q3_cosmic->Add(hBoth_cosmic_Q3);
    hA17_Q4_cosmic->Add(hBoth_cosmic_Q4);

    hA17_Q1->Add(hA17_Q1_cosmic);
    hA17_Q2->Add(hA17_Q2_cosmic);
    hA17_Q3->Add(hA17_Q3_cosmic);
    hA17_Q4->Add(hA17_Q4_cosmic);

    hA17->Add(hA17_Q1);
    hA17->Add(hA17_Q2);
    hA17->Add(hA17_Q3);
    hA17->Add(hA17_Q4);

    hA17_cosmic->Add(hA17_Q1_cosmic);
    hA17_cosmic->Add(hA17_Q2_cosmic);
    hA17_cosmic->Add(hA17_Q3_cosmic);
    hA17_cosmic->Add(hA17_Q4_cosmic);

    // And finally make the fits
    FitResult fitSA = getBestFit(prediction_SA, hSA, hSpecCosmics_SA_and_both_rebin);

    ssth23SA   = fitSA.bestX;
    dmsq32SA   = fitSA.bestY;
    LLSA       = fitSA.LL;
    sigmaMaxSA = fitSA.sigma;

    std::cout << "SA. rejection for expt " << i+1 << " " << sigmaMaxSA << std::endl;

    FitResult fitA17 = getBestFit(prediction_Q1, prediction_Q2, prediction_Q3, prediction_Q4,
                           hA17_Q1, hA17_Q2, hA17_Q3, hA17_Q4,
                           hSpecCosmics_Q1, hSpecCosmics_Q2, hSpecCosmics_Q3, hSpecCosmics_Q4);

    ssth23A17   = fitA17.bestX;
    dmsq32A17   = fitA17.bestY;
    LLA17       = fitA17.LL;
    sigmaMaxA17 = fitA17.sigma;

    std::cout << "--- A17 rejection for expt " << i+1 << " " << sigmaMaxA17 << std::endl;
    
    fTree->Fill();

  } // loop over entries

    fTree->Write();
    fOutput->Close();

} // End of function