Classes | Functions
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().

433 {
434  osc::OscCalcPMNSOpt calc2;
435  ana::ResetOscCalcToDefault(&calc2);
436 
437  Spectrum fakeData(h1, kSecondAnaPOT, kSecondAnaLivetime);
438 
439  // Experiment
440  SingleSampleExperiment fakeExpt1(prediction, fakeData, hCosmic);
441 
442  std::vector<const IExperiment*> vExpt;
443  vExpt.push_back(&fakeExpt1);
444  vExpt.push_back(&kSolarConstraintsPDG2017);
445  vExpt.push_back(WorldReactorConstraint2017());
446  MultiExperiment fakeExpt(vExpt);
447 
448  std::vector<const IFitVar*> margPars_ssth23 = getNumuMarginalisedOscParam();
449  margPars_ssth23.push_back(&kFitDmSq32Scaled);
450 
451  MinuitFitter fitMaxMix(&fakeExpt, margPars_ssth23);
452  double MMLL = fitMaxMix.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();
453 
454  margPars_ssth23.push_back(&kFitSinSqTheta23);
455  MinuitFitter fit(&fakeExpt, margPars_ssth23);
456  double bestLL = fit.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();
457 
458  double bestX = pow(TMath::Sin(calc2.GetTh23()),2);
459  double bestY = 1e3 * calc2.GetDmsq32();
460 
461  double MMrej = 0;
462  if(MMLL > bestLL) MMrej = sqrt(MMLL-bestLL); // this should always happen, but just in case
463 
464  FitResult result = {bestX, bestY, bestLL, MMrej};
465 
466  return result;
467 }
ana::ResetOscCalcToDefault
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
cet::pow
constexpr T pow(T x)
Definition: pow.h:75
osc::_IOscCalcAdjustable::GetTh23
virtual T GetTh23() const
Definition: IOscCalc.h:89
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
osc::_OscCalcPMNSOpt
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
ana::WorldReactorConstraint2017
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
Definition: ReactorExperiment.cxx:112
osc::_IOscCalcAdjustable::GetDmsq32
virtual T GetDmsq32() const
Definition: IOscCalc.h:86
ana::kFitDmSq32Scaled
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
Definition: VarsWithPriors.cxx:13
h1
TH1F * h1
Definition: nslices_potspill.C:4
ana::MultiExperiment
Combine multiple component experiments.
Definition: MultiExperiment.h:11
ana::kSolarConstraintsPDG2017
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
Definition: SolarConstraints.h:56
ana::kFitSinSqTheta23
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
Definition: VarsWithPriors.cxx:9
ana::kSecondAnaPOT
const double kSecondAnaPOT
Definition: Exposures.h:87
ana::kSecondAnaLivetime
const double kSecondAnaLivetime
Definition: Exposures.h:113
ana::SingleSampleExperiment
Compare a single data spectrum to the MC + cosmics expectation.
Definition: SingleSampleExperiment.h:21
ana::getNumuMarginalisedOscParam
std::vector< const IFitVar * > getNumuMarginalisedOscParam()
Definition: FitVarsNumu.cxx:9
ana::MinuitFitter
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
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().

387 {
388  osc::OscCalcPMNSOpt calc2;
389  ana::ResetOscCalcToDefault(&calc2);
390 
391  Spectrum fakeData1(h1, kSecondAnaPOT, kSecondAnaLivetime);
392  Spectrum fakeData2(h2, kSecondAnaPOT, kSecondAnaLivetime);
393  Spectrum fakeData3(h3, kSecondAnaPOT, kSecondAnaLivetime);
394  Spectrum fakeData4(h4, kSecondAnaPOT, kSecondAnaLivetime);
395 
396  SingleSampleExperiment fakeExpt1 (quant1, fakeData1, hCosmic1);
397  SingleSampleExperiment fakeExpt2 (quant2, fakeData2, hCosmic2);
398  SingleSampleExperiment fakeExpt3 (quant3, fakeData3, hCosmic3);
399  SingleSampleExperiment fakeExpt4 (quant4, fakeData4, hCosmic4);
400 
401  std::vector<const IExperiment*> vExpt;
402  vExpt.push_back(&fakeExpt1);
403  vExpt.push_back(&fakeExpt2);
404  vExpt.push_back(&fakeExpt3);
405  vExpt.push_back(&fakeExpt4);
406  vExpt.push_back(&kSolarConstraintsPDG2017);
407  vExpt.push_back(WorldReactorConstraint2017());
408  MultiExperiment fakeExpt(vExpt);
409 
410  std::vector<const IFitVar*> margPars_ssth23 = getNumuMarginalisedOscParam();
411  margPars_ssth23.push_back(&kFitDmSq32Scaled);
412 
413  MinuitFitter fitMaxMix(&fakeExpt, margPars_ssth23);
414  double MMLL = fitMaxMix.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();
415 
416  margPars_ssth23.push_back(&kFitSinSqTheta23);
417  MinuitFitter fit(&fakeExpt, margPars_ssth23);
418  double bestLL = fit.Fit(&calc2, IFitter::kQuiet)->EvalMetricVal();
419 
420  double bestX = pow(TMath::Sin(calc2.GetTh23()),2);
421  double bestY = 1e3 * calc2.GetDmsq32();
422 
423  double MMrej = 0;
424  if(MMLL > bestLL) MMrej = sqrt(MMLL-bestLL); // this should always happen, but just in case
425 
426  FitResult result = {bestX, bestY, bestLL, MMrej};
427 
428  return result;
429 
430 }
ana::ResetOscCalcToDefault
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
h3
TH1F * h3
Definition: berger.C:36
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
cet::pow
constexpr T pow(T x)
Definition: pow.h:75
osc::_IOscCalcAdjustable::GetTh23
virtual T GetTh23() const
Definition: IOscCalc.h:89
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
osc::_OscCalcPMNSOpt
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
ana::WorldReactorConstraint2017
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
Definition: ReactorExperiment.cxx:112
osc::_IOscCalcAdjustable::GetDmsq32
virtual T GetDmsq32() const
Definition: IOscCalc.h:86
ana::kFitDmSq32Scaled
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
Definition: VarsWithPriors.cxx:13
h2
TH1F * h2
Definition: plot.C:45
h1
TH1F * h1
Definition: nslices_potspill.C:4
ana::MultiExperiment
Combine multiple component experiments.
Definition: MultiExperiment.h:11
ana::kSolarConstraintsPDG2017
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
Definition: SolarConstraints.h:56
ana::kFitSinSqTheta23
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
Definition: VarsWithPriors.cxx:9
ana::kSecondAnaPOT
const double kSecondAnaPOT
Definition: Exposures.h:87
ana::kSecondAnaLivetime
const double kSecondAnaLivetime
Definition: Exposures.h:113
ana::SingleSampleExperiment
Compare a single data spectrum to the MC + cosmics expectation.
Definition: SingleSampleExperiment.h:21
ana::getNumuMarginalisedOscParam
std::vector< const IFitVar * > getNumuMarginalisedOscParam()
Definition: FitVarsNumu.cxx:9
ana::MinuitFitter
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
void getPoissonRandom ( TH1D *  s,
TH1F *  h 
)

Definition at line 504 of file Toy_analyses.C.

References MECModelEnuComparisons::i, and r().

Referenced by Toy_analyses().

505 {
506  for(int i = 1; i <= s->GetNbinsX(); i++)
507  {
508  float fX = s->GetBinCenter(i);
509  float fY = r.Poisson(s->GetBinContent(i));
510 
511  h->Fill(fX, fY);
512  }
513 }
s
const XML_Char * s
Definition: expat.h:262
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
r
TRandom3 r(0)
void getSArebin ( TH1D *  s,
TH1D *  h 
)

Definition at line 515 of file Toy_analyses.C.

References MECModelEnuComparisons::i.

Referenced by Toy_analyses().

516 {
517  for(int i = 1; i <= s->GetNbinsX(); i++)
518  {
519  double fX = s->GetBinCenter(i);
520  double fY = s->GetBinContent(i);
521 
522  h->Fill(fX, fY);
523  }
524 }
s
const XML_Char * s
Definition: expat.h:262
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
void getSArebin ( TH1F *  s,
TH1F *  h 
)

Definition at line 526 of file Toy_analyses.C.

References MECModelEnuComparisons::i.

527 {
528  for(int i = 1; i <= s->GetNbinsX(); i++)
529  {
530  float fX = s->GetBinCenter(i);
531  float fY = s->GetBinContent(i);
532 
533  h->Fill(fX, fY);
534  }
535 }
s
const XML_Char * s
Definition: expat.h:262
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
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().

481 {
482  float A17Ebins[20] = { 0, 0.75, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
483  1.9, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5 };
484  unsigned int quant = 0;
485 
486  unsigned int k = 0;
487 
488  while (A17Ebins[k] < E && k < 20)
489  // for(unsigned int k = 0; k < 20; k++)
490  {
491  unsigned int l = 0;
492  while (v.at(k)[l] < hadEfrac && l < v.at(k).size())
493  // for(unsigned int l = 0; l < v.at(k).size(); l++)
494  {
495  quant = l;
496  l++;
497  }
498  k++;
499  }
500 
501  return quant;
502 }
E
Float_t E
Definition: plot.C:20
submit_hadd.l
l
Definition: submit_hadd.py:89
registry_explorer.v
v
Definition: registry_explorer.py:22
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().

472  {
473 
474  std::vector<std::vector<double>> quantileBins = GetQuantileBins(quantileHist,
475  independentAxis, quantileAxis, numQuantiles);
476 
477  return quantileBins;
478 }
ana::GetQuantileBins
std::vector< std::vector< double > > GetQuantileBins(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
Returns a 2D vector First index is the independentAxis bin number Second index is the high bin edge f...
Definition: QuantileCuts.cxx:30
TRandom3 r ( 0  )

Referenced by util::HoughCalc::AddPoint(), genie::geometry::ROOTGeomAnalyzer::AdoptGeomVolSelector(), caf::SRProxySystController::AnyShifted(), art::AssnsNode< L, R, D >::AssnsNode(), osc::EarthModel::AveNe(), fnex::Spectrum::BaseVals(), 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(), ana::LogLikelihood(), 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(), ND_predictionbyMC_usingbincontent(), 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::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(), fnex::NuisanceParameter_Experiment::ToString(), fnex::NOvA_Experiment::ToString(), fnex::CorrSpec_SimpleExtrap::ToString(), fnex::CorrSpec_NoExtrapNuE::ToString(), fnex::CorrSpec_NoExtrapNuMu::ToString(), fnex::CorrSpec_MichelDecomp::ToString(), fnex::CorrSpec_PropDecomp::ToString(), fnex::CorrSpec_SimpleNuMuExtrap::ToString(), fnex::CorrSpec_BENDecomp::ToString(), fnex::Correlations::ToString(), Toy_analyses(), 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, plotROC::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(), and ana::Spectrum::ToTH1().

80 {
81  std::string inputDir = "/pnfs/nova/persistent/users/bzamoran/ToyPvalue/";
82 
83  TFile* fInputTree = TFile::Open(pnfs2xrootd(inputDir + "inputTree.root").c_str());
84  TTree* inputTree = (TTree*) fInputTree->Get("inputTree");
85  float TrueE, CCE, CCESA, HadEfrac, PPFX, CXW;
86  int isSA, isA17;
87  inputTree->SetBranchAddress("TrueE", &TrueE);
88  inputTree->SetBranchAddress("isSA", &isSA);
89  inputTree->SetBranchAddress("isA17", &isA17);
90  inputTree->SetBranchAddress("CCE", &CCE);
91  inputTree->SetBranchAddress("CCESA", &CCESA);
92  inputTree->SetBranchAddress("HadEfrac", &HadEfrac);
93  inputTree->SetBranchAddress("PPFX", &PPFX);
94  inputTree->SetBranchAddress("CXW", &CXW);
95 
96  std::cout << "Tree" << std::endl;
97 
98  TFile* fPred = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_AllQuants.root").c_str());
99  PredictionInterp* prediction_OR = LoadFrom<PredictionInterp>(fPred, "prediction_OR" ).release();
100  PredictionInterp* prediction_SA = LoadFrom<PredictionInterp>(fPred, "prediction_SA" ).release();
101  fPred->Close();
102 
103  std::cout << "Predictions 1" << std::endl;
104 
105  TFile* fPred_Q1 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant1.root").c_str());
106  TFile* fPred_Q2 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant2.root").c_str());
107  TFile* fPred_Q3 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant3.root").c_str());
108  TFile* fPred_Q4 = TFile::Open(pnfs2xrootd(inputDir + "Predictions/pred_3A_prod3_Quant4.root").c_str());
109 
110  PredictionInterp* prediction_Q1 = LoadFrom<PredictionInterp>(fPred_Q1, "prediction").release();
111  PredictionInterp* prediction_Q2 = LoadFrom<PredictionInterp>(fPred_Q2, "prediction").release();
112  PredictionInterp* prediction_Q3 = LoadFrom<PredictionInterp>(fPred_Q3, "prediction").release();
113  PredictionInterp* prediction_Q4 = LoadFrom<PredictionInterp>(fPred_Q4, "prediction").release();
114 
115  std::cout << "Predictions quant" << std::endl;
116 
117  TFile* fCosmics = TFile::Open(pnfs2xrootd(inputDir + "cos_bkg_hists_v2.root").c_str());
118  TH1D* hSpecCosmics_Q1 = (TH1D*) fCosmics->Get("q1all");
119  TH1D* hSpecCosmics_Q2 = (TH1D*) fCosmics->Get("q2all");
120  TH1D* hSpecCosmics_Q3 = (TH1D*) fCosmics->Get("q3all");
121  TH1D* hSpecCosmics_Q4 = (TH1D*) fCosmics->Get("q4all");
122 
123  hSpecCosmics_Q1->Scale(kSecondAnaLivetime / kAna2017Livetime ); // scale to 6e20
124  hSpecCosmics_Q2->Scale(kSecondAnaLivetime / kAna2017Livetime );
125  hSpecCosmics_Q3->Scale(kSecondAnaLivetime / kAna2017Livetime );
126  hSpecCosmics_Q4->Scale(kSecondAnaLivetime / kAna2017Livetime );
127 
128  TF1* fSA = new TF1("fSA","TMath::Exp(1.976*(x-2.26))/(1+TMath::Exp(1.976*(x-2.26)))",0,5);
129  TF1* fA17 = new TF1("fA17","TMath::Exp(-1.907*(x-1.967))/(1+TMath::Exp(-1.907*(x-1.967)))",0,5);
130  TF1* fBoth = new TF1("fBoth","1.-fA17-fSA",0,5);
131 
132  TF1* fAdjustSA = new TF1("fAdjustSA","TMath::Min(fSA/(1.-fSA),5.0)", 0.5, 5); // don't let it blow
133  TF1* fAdjustA17 = new TF1("fAdjustA17","fA17/(1.-fSA)", 0.5, 5);
134  TF1* fAdjustBoth = new TF1("fAdjustBoth","fBoth/(1.-fSA)", 0.5, 5);
135  TF1* fAdjustSAandBoth = new TF1("fAdjustSAandBoth","TMath::Min((1.-fA17)/(1.-fSA),5.0)", 0.5, 5);
136 
137  TH1D* hSpecCosmics_SA = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_SA");
138  hSpecCosmics_SA->Add(hSpecCosmics_Q2);
139  hSpecCosmics_SA->Add(hSpecCosmics_Q3);
140  hSpecCosmics_SA->Add(hSpecCosmics_Q4);
141  hSpecCosmics_SA->Multiply(fAdjustSA);
142  TH1D* hSpecCosmics_SA_rebin = new TH1D("hSpecCosmics_SA_rebin","",20,0,5);
143  getSArebin(hSpecCosmics_SA, hSpecCosmics_SA_rebin);
144 
145  TH1D* hSpecCosmics_SA_and_both = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_SA_and_both");
146  hSpecCosmics_SA_and_both->Add(hSpecCosmics_Q2);
147  hSpecCosmics_SA_and_both->Add(hSpecCosmics_Q3);
148  hSpecCosmics_SA_and_both->Add(hSpecCosmics_Q4);
149  hSpecCosmics_SA_and_both->Multiply(fAdjustSAandBoth);
150  TH1D* hSpecCosmics_SA_and_both_rebin = new TH1D("hSpecCosmics_SA_and_both_rebin","",20,0,5);
151  getSArebin(hSpecCosmics_SA_and_both, hSpecCosmics_SA_and_both_rebin);
152 
153  TH1D* hSpecCosmics_Both_Q1 = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_Both_Q1");
154  TH1D* hSpecCosmics_Both_Q2 = (TH1D*) hSpecCosmics_Q2->Clone("hSpecCosmics_Both_Q2");
155  TH1D* hSpecCosmics_Both_Q3 = (TH1D*) hSpecCosmics_Q3->Clone("hSpecCosmics_Both_Q3");
156  TH1D* hSpecCosmics_Both_Q4 = (TH1D*) hSpecCosmics_Q4->Clone("hSpecCosmics_Both_Q4");
157  hSpecCosmics_Both_Q1->Multiply(fAdjustBoth);
158  hSpecCosmics_Both_Q2->Multiply(fAdjustBoth);
159  hSpecCosmics_Both_Q3->Multiply(fAdjustBoth);
160  hSpecCosmics_Both_Q4->Multiply(fAdjustBoth);
161 
162  TH1D* hSpecCosmics_A17_Q1 = (TH1D*) hSpecCosmics_Q1->Clone("hSpecCosmics_A17_Q1");
163  TH1D* hSpecCosmics_A17_Q2 = (TH1D*) hSpecCosmics_Q2->Clone("hSpecCosmics_A17_Q2");
164  TH1D* hSpecCosmics_A17_Q3 = (TH1D*) hSpecCosmics_Q3->Clone("hSpecCosmics_A17_Q3");
165  TH1D* hSpecCosmics_A17_Q4 = (TH1D*) hSpecCosmics_Q4->Clone("hSpecCosmics_A17_Q4");
166  hSpecCosmics_A17_Q1->Multiply(fAdjustA17);
167  hSpecCosmics_A17_Q2->Multiply(fAdjustA17);
168  hSpecCosmics_A17_Q3->Multiply(fAdjustA17);
169  hSpecCosmics_A17_Q4->Multiply(fAdjustA17);
170 
171 
172  std::cout << "All's well with loading" << std::endl;
173 
174  // Get FD files for HadEFrac vs. neutrino energy distribution
175  std::string quantpath = "/pnfs/nova/persistent/analysis/numu/Ana2017/all_FD_histo_for_quantile_cuts.root";
176  TFile* quantfile = TFile::Open(pnfs2xrootd(quantpath).c_str());
177  TH2* quanthist = (TH2*)quantfile->Get("dir_FDSpec2D/FDSpec2D");
178  std::vector<std::vector<double>> QuantBins = QuantileVectorFromTH2(quanthist, kNumuCCOptimisedAxis, kHadEFracAxis, 4);
179 
180  // Now the fittting
181 
182  osc::OscCalcPMNSOpt calc;
183  ana::ResetOscCalcToDefault(&calc);
184  // calc.SetDmsq32(2.67e-3);
185  // calc.SetTh23(asin(sqrt(0.404))); // PRL-published Second Ana
186  calc.SetDmsq32(2.43071e-3);
187  calc.SetTh23(asin(sqrt(0.471974)));
188 
189  TFile* fOutput = new TFile(Form("Toy_Output_vol%i.root", volume),"RECREATE");
190  TTree* fTree = new TTree("T","Summary tree from p-value study");
191 
192  TH1F* hSA = new TH1F("hSA","",20,0,5);
193  TH1F* hSA_cosmic = new TH1F("hSA_cosmic","",20,0,5);
194 
195  float A17Ebins[20] = { 0, 0.75, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8,
196  1.9, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5 };
197 
198  TH1F* hA17 = new TH1F("hA17","",19, A17Ebins);
199  TH1F* hA17_Q1 = new TH1F("hA17_Q1","",19, A17Ebins);
200  TH1F* hA17_Q2 = new TH1F("hA17_Q2","",19, A17Ebins);
201  TH1F* hA17_Q3 = new TH1F("hA17_Q3","",19, A17Ebins);
202  TH1F* hA17_Q4 = new TH1F("hA17_Q4","",19, A17Ebins);
203 
204  TH1F* hA17_cosmic = new TH1F("hA17_cosmic","",19, A17Ebins);
205  TH1F* hA17_Q1_cosmic = new TH1F("hA17_Q1_cosmic","",19, A17Ebins);
206  TH1F* hA17_Q2_cosmic = new TH1F("hA17_Q2_cosmic","",19, A17Ebins);
207  TH1F* hA17_Q3_cosmic = new TH1F("hA17_Q3_cosmic","",19, A17Ebins);
208  TH1F* hA17_Q4_cosmic = new TH1F("hA17_Q4_cosmic","",19, A17Ebins);
209 
210  TH1F* hBoth = new TH1F("hBoth","",19, A17Ebins);
211 
212  TH1F* hBoth_cosmic_Q1 = new TH1F("hBoth_cosmic_Q1","",19, A17Ebins);
213  TH1F* hBoth_cosmic_Q2 = new TH1F("hBoth_cosmic_Q2","",19, A17Ebins);
214  TH1F* hBoth_cosmic_Q3 = new TH1F("hBoth_cosmic_Q3","",19, A17Ebins);
215  TH1F* hBoth_cosmic_Q4 = new TH1F("hBoth_cosmic_Q4","",19, A17Ebins);
216  TH1F* hBoth_cosmic = new TH1F("hBoth_cosmic","",19, A17Ebins);
217  TH1F* hBoth_cosmic_rebin = new TH1F("hA17_Both_cosmic_rebin","",20,0,5);
218 
219  double sigmaMaxSA, sigmaMaxA17, LLSA, LLA17;
220  double dmsq32SA, dmsq32A17, ssth23SA, ssth23A17;
221  unsigned int nPoints;
222 
223  fTree->Branch("sigmaMaxSA", &sigmaMaxSA);
224  fTree->Branch("sigmaMaxA17", &sigmaMaxA17);
225  fTree->Branch("LLSA", &LLSA);
226  fTree->Branch("LLA17", &LLA17);
227  fTree->Branch("dmsq32SA", &dmsq32SA);
228  fTree->Branch("dmsq32A17", &dmsq32A17);
229  fTree->Branch("ssth23SA", &ssth23SA);
230  fTree->Branch("ssth23A17", &ssth23A17);
231  fTree->Branch("Spectrum_SA", &hSA);
232  fTree->Branch("Cosmics_SA", &hSA_cosmic);
233  fTree->Branch("Spectrum_A17", &hA17);
234  fTree->Branch("Spectrum_A17_Q1", &hA17_Q1);
235  fTree->Branch("Spectrum_A17_Q2", &hA17_Q2);
236  fTree->Branch("Spectrum_A17_Q3", &hA17_Q3);
237  fTree->Branch("Spectrum_A17_Q4", &hA17_Q4);
238  fTree->Branch("Cosmics_A17", &hA17_cosmic);
239  fTree->Branch("Cosmics_A17_Q1", &hA17_Q1_cosmic);
240  fTree->Branch("Cosmics_A17_Q2", &hA17_Q2_cosmic);
241  fTree->Branch("Cosmics_A17_Q3", &hA17_Q3_cosmic);
242  fTree->Branch("Cosmics_A17_Q4", &hA17_Q4_cosmic);
243  fTree->Branch("Spectrum_Both", &hBoth);
244  fTree->Branch("Cosmics_Both", &hBoth_cosmic);
245 
246  unsigned int nEntries = inputTree->GetEntries();
247 
248  double MCpoints = prediction_OR->Predict(&calc).ToTH1(kSecondAnaPOT)->Integral();
249 
250  for(unsigned int i = 0; i < nExperiments; i++)
251  {
252  hSA->Reset();
253  hSA_cosmic->Reset();
254  hA17->Reset();
255  hA17_Q1->Reset();
256  hA17_Q2->Reset();
257  hA17_Q3->Reset();
258  hA17_Q4->Reset();
259  hA17_cosmic->Reset();
260  hA17_Q1_cosmic->Reset();
261  hA17_Q2_cosmic->Reset();
262  hA17_Q3_cosmic->Reset();
263  hA17_Q4_cosmic->Reset();
264  hBoth->Reset();
265  hBoth_cosmic_Q1->Reset();
266  hBoth_cosmic_Q2->Reset();
267  hBoth_cosmic_Q3->Reset();
268  hBoth_cosmic_Q4->Reset();
269  hBoth_cosmic->Reset();
270  hBoth_cosmic_rebin->Reset();
271 
272  nPoints = r.Poisson(MCpoints);
273 
274  unsigned int k = 0;
275  while(k <= nPoints)
276  {
277  unsigned int iEntry = r.Integer(nEntries);
278  inputTree->GetEntry(iEntry);
279 
280  double p = r.Uniform();
281  double prob = calc.P(14, 14, TrueE) * TMath::Min(CXW * PPFX * 0.5, 1.); // adjusting for weights
282 
283  if(p < prob && TrueE > 0.0 && TrueE < 5.0)
284  {
285  if(isSA)
286  {
287  hSA->Fill(CCE);
288  }
289  if(isA17)
290  {
291  int qi = quantile(QuantBins, CCE, HadEfrac);
292  if(qi == 0) hA17_Q1->Fill(CCE);
293  else if(qi == 1) hA17_Q2->Fill(CCE);
294  else if(qi == 2) hA17_Q3->Fill(CCE);
295  else hA17_Q4->Fill(CCE);
296  }
297  if(isA17 && isSA) hBoth->Fill(CCE);
298  if(isSA || isA17) k++;
299  }
300  } // loop over points
301 
302  // Protection against errors reading the tree
303  if(hSA->FindFirstBinAbove(1) == hSA->FindLastBinAbove(1)) continue;
304  if(hBoth->FindFirstBinAbove(1) == hBoth->FindLastBinAbove(1)) continue;
305  if(hA17_Q1->FindFirstBinAbove(1) == hA17_Q1->FindLastBinAbove(1)) continue;
306  if(hA17_Q2->FindFirstBinAbove(1) == hA17_Q2->FindLastBinAbove(1)) continue;
307  if(hA17_Q3->FindFirstBinAbove(1) == hA17_Q3->FindLastBinAbove(1)) continue;
308  if(hA17_Q4->FindFirstBinAbove(1) == hA17_Q4->FindLastBinAbove(1)) continue;
309 
310  // Now handle the cosmics
311  getPoissonRandom(hSpecCosmics_Both_Q1, hBoth_cosmic_Q1);
312  getPoissonRandom(hSpecCosmics_Both_Q2, hBoth_cosmic_Q2);
313  getPoissonRandom(hSpecCosmics_Both_Q3, hBoth_cosmic_Q3);
314  getPoissonRandom(hSpecCosmics_Both_Q4, hBoth_cosmic_Q4);
315 
316  hBoth_cosmic->Add(hBoth_cosmic_Q1);
317  hBoth_cosmic->Add(hBoth_cosmic_Q2);
318  hBoth_cosmic->Add(hBoth_cosmic_Q3);
319  hBoth_cosmic->Add(hBoth_cosmic_Q4);
320 
321  hBoth->Add(hBoth_cosmic);
322 
323  getSArebin(hBoth_cosmic, hBoth_cosmic_rebin);
324 
325  getPoissonRandom(hSpecCosmics_SA_rebin, hSA_cosmic);
326  hSA_cosmic->Add(hBoth_cosmic_rebin);
327  hSA->Add(hSA_cosmic);
328 
329  getPoissonRandom(hSpecCosmics_A17_Q1, hA17_Q1_cosmic);
330  getPoissonRandom(hSpecCosmics_A17_Q2, hA17_Q2_cosmic);
331  getPoissonRandom(hSpecCosmics_A17_Q3, hA17_Q3_cosmic);
332  getPoissonRandom(hSpecCosmics_A17_Q4, hA17_Q4_cosmic);
333 
334  hA17_Q1_cosmic->Add(hBoth_cosmic_Q1);
335  hA17_Q2_cosmic->Add(hBoth_cosmic_Q2);
336  hA17_Q3_cosmic->Add(hBoth_cosmic_Q3);
337  hA17_Q4_cosmic->Add(hBoth_cosmic_Q4);
338 
339  hA17_Q1->Add(hA17_Q1_cosmic);
340  hA17_Q2->Add(hA17_Q2_cosmic);
341  hA17_Q3->Add(hA17_Q3_cosmic);
342  hA17_Q4->Add(hA17_Q4_cosmic);
343 
344  hA17->Add(hA17_Q1);
345  hA17->Add(hA17_Q2);
346  hA17->Add(hA17_Q3);
347  hA17->Add(hA17_Q4);
348 
349  hA17_cosmic->Add(hA17_Q1_cosmic);
350  hA17_cosmic->Add(hA17_Q2_cosmic);
351  hA17_cosmic->Add(hA17_Q3_cosmic);
352  hA17_cosmic->Add(hA17_Q4_cosmic);
353 
354  // And finally make the fits
355  FitResult fitSA = getBestFit(prediction_SA, hSA, hSpecCosmics_SA_and_both_rebin);
356 
357  ssth23SA = fitSA.bestX;
358  dmsq32SA = fitSA.bestY;
359  LLSA = fitSA.LL;
360  sigmaMaxSA = fitSA.sigma;
361 
362  std::cout << "SA. rejection for expt " << i+1 << " " << sigmaMaxSA << std::endl;
363 
364  FitResult fitA17 = getBestFit(prediction_Q1, prediction_Q2, prediction_Q3, prediction_Q4,
365  hA17_Q1, hA17_Q2, hA17_Q3, hA17_Q4,
366  hSpecCosmics_Q1, hSpecCosmics_Q2, hSpecCosmics_Q3, hSpecCosmics_Q4);
367 
368  ssth23A17 = fitA17.bestX;
369  dmsq32A17 = fitA17.bestY;
370  LLA17 = fitA17.LL;
371  sigmaMaxA17 = fitA17.sigma;
372 
373  std::cout << "--- A17 rejection for expt " << i+1 << " " << sigmaMaxA17 << std::endl;
374 
375  fTree->Fill();
376 
377  } // loop over entries
378 
379  fTree->Write();
380  fOutput->Close();
381 
382 } // End of function
FitResult::LL
double LL
Definition: Toy_analyses.C:60
ana::PredictionInterp
Implements systematic errors by interpolation between shifted templates.
Definition: PredictionInterp.h:23
ana::ResetOscCalcToDefault
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
plotROC.nPoints
nPoints
Definition: plotROC.py:88
ana::Spectrum::ToTH1
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:209
getPoissonRandom
void getPoissonRandom(TH1D *s, TH1F *h)
Definition: Toy_analyses.C:504
p
const char * p
Definition: xmltok.h:285
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
ana::kHadEFracAxis
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
ana::pnfs2xrootd
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: Utilities.cxx:620
calc
osc::OscCalcDumb calc
Definition: make_nueFDprediction_kinematics.C:70
FitResult::sigma
double sigma
Definition: Toy_analyses.C:61
getBestFit
FitResult getBestFit(IPrediction *prediction, TH1F *h1, TH1D *hCosmic)
Definition: Toy_analyses.C:432
ana::kAna2017Livetime
const double kAna2017Livetime
Definition: Exposures.h:200
osc::_OscCalcPMNSOpt
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
osc::_OscCalcPMNSOpt::P
T P(int flavBefore, int flavAfter, double E) override
E in GeV; flavors as PDG codes (so, neg==>antinu)
Definition: OscCalcPMNSOpt.cxx:38
FitResult::bestX
double bestX
Definition: Toy_analyses.C:58
osc::_OscCalcPMNSOpt::SetTh23
void SetTh23(const T &th23) override
Definition: OscCalcPMNSOpt.h:42
ana::kNumuCCOptimisedAxis
const HistAxis kNumuCCOptimisedAxis("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kCCE)
HistAxis that implements optimised numuCCE from L. Vinton. See docdb 16332. This was close to &#39;custC&#39;...
Definition: HistAxes.h:24
quantile
unsigned int quantile(std::vector< std::vector< double >> v, float E, float hadEfrac)
Definition: Toy_analyses.C:480
om::cout
OStream cout
Definition: OStream.cxx:6
osc::_OscCalcPMNSOpt::SetDmsq32
void SetDmsq32(const T &dmsq32) override
Definition: OscCalcPMNSOpt.h:39
ana::PredictionInterp::Predict
Spectrum Predict(osc::IOscCalc *calc) const override
Definition: PredictionInterp.cxx:305
xml_schema::string
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
FitResult::bestY
double bestY
Definition: Toy_analyses.C:59
getSArebin
void getSArebin(TH1D *s, TH1D *h)
Definition: Toy_analyses.C:515
r
TRandom3 r(0)
QuantileVectorFromTH2
std::vector< std::vector< double > > QuantileVectorFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles)
Definition: Toy_analyses.C:469
e
Float_t e
Definition: plot.C:35
ana::kSecondAnaPOT
const double kSecondAnaPOT
Definition: Exposures.h:87
std::asin
T asin(T number)
Definition: d0nt_math.hpp:60
ana::kSecondAnaLivetime
const double kSecondAnaLivetime
Definition: Exposures.h:113
Generated on Wed Sep 23 2020 05:16:19 for NOvA by   doxygen 1.8.11