joint_fit_2018_tools.h File Reference

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "3FlavorAna/Prediction/PredictionAddRock.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "CAFAna/Analysis/CountingExperimentStatistics.h"
#include "3FlavorAna/Systs/JointAna2018Systs.h"
#include "CAFAna/Systs/XSecSysts.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Analysis/Style.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "Utilities/func/MathUtil.h"
#include "TCanvas.h"
#include "TBox.h"
#include "TLatex.h"
#include "TColor.h"
#include "TGraph.h"
#include "TVectorD.h"
#include "TF1.h"
#include "TLegend.h"
#include "TLine.h"
#include "TMarker.h"
#include "TStyle.h"
#include "TSystem.h"
#include "TGaxis.h"
#include "OscLib/IOscCalc.h"
#include "TFile.h"

Go to the source code of this file.

Classes
class	FitSinSqTheta23LowerOctant
class	FitSinSqTheta23UpperOctant

Functions
TFile *	GetGaussianSurface (TString par, TString hier)
void	CylindricalSmoothing (TH2 *h)
std::string	GetNuePredPath (std::string beam, bool isFakeData, bool NERSC)
const IPrediction *	GetNuePrediction2018 (std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
std::pair< TH1D *, double >	GetNueCosmics2018 (std::string beam, bool GetFromUPS=false, bool NERSC=false)
Spectrum *	GetNueData2018 (std::string beam)
std::vector< const IPrediction * >	GetNumuPredictions2018 (const int nq=4, bool useSysts=true, std::string beam="fhc", bool GetFromUPS=false, ENu2018ExtrapType numuExtrap=kNuMu, bool minimizeMemory=false, bool NERSC=false)
std::vector< std::pair< TH1D *, double > >	GetNumuCosmics2018 (const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
std::vector< Spectrum * >	GetNumuData2018 (const int nq=4, std::string beam="fhc")
Spectrum *	GetFakeData (const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
Spectrum *	GetMockData (const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
void	SetFakeCalc (osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
std::vector< const ISyst * >	GetJointFitSystematicList (bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool isFHC=true, bool isRHC=true, bool intersection=true)
std::vector< std::pair< const ISyst *, const ISyst * > >	GetCorrelations (bool isNue, bool isFHC)
void	DrawSliceCanvas (TString slicename, const double ymax, const double xmin=0, const double xmax=2., bool overlay=false)
TLegend *	SliceLegend (std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta, double low=0.7, bool isDmsq=false, bool overlay=false)
void	DrawContourCanvas (TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1)
TLegend *	ContourLegend (int hie, bool fillContour, bool fccorr, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor, bool bestFit)
TLegend *	FCOverlayContourLegend (int hie, bool fillContour, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor)

Variables
const IConstrainedFitVar *	kFitSinSqTheta23LowerOctant = new FitSinSqTheta23LowerOctant()
const IConstrainedFitVar *	kFitSinSqTheta23UpperOctant = new FitSinSqTheta23UpperOctant()

Function Documentation

TLegend* ContourLegend	(	int	hie,
		bool	fillContour,
		bool	fccorr,
		Int_t	kFillColor1,
		Int_t	kFillColor2,
		Int_t	kFillColor3,
		Int_t	kDarkColor,
		bool	bestFit
	)

Definition at line 643 of file joint_fit_2018_tools.h.

References APDHVSetting::dummy, ana::k2SigmaConfidenceStyle, ana::kBlessedLabelFontSize, and MECModelEnuComparisons::leg.

                                     {
  TLegend * leg = new TLegend(0.16,fccorr?0.19:0.17,0.65,0.28);
  leg->SetTextSize(kBlessedLabelFontSize);
  leg->SetFillStyle(0);
  leg->SetMargin(1.3*leg->GetMargin());
  leg->SetNColumns(3);

  TGraph * dummy = new TGraph;
  dummy->SetLineWidth(3);

  if(fillContour){
    dummy->SetLineColor(kDarkColor);
    dummy->SetFillColor(kFillColor1);
    leg->AddEntry(dummy->Clone(),"1#sigma","f");
    dummy->SetFillColor(kFillColor2);
    leg->AddEntry(dummy->Clone(),"2#sigma","f");
    dummy->SetFillColor(kFillColor3);
    leg->AddEntry(dummy->Clone(),"3#sigma","f");
  }
  else{
    leg->SetMargin(1.4*leg->GetMargin());
    dummy->SetLineColor(kFillColor1);
    dummy->SetLineStyle(kSolid);
    leg->AddEntry(dummy->Clone(),"1#sigma","l");
    dummy->SetLineStyle(k2SigmaConfidenceStyle);
    dummy->SetLineColor(kFillColor2);
    leg->AddEntry(dummy->Clone(),"2#sigma","l");
    dummy->SetLineStyle(kSolid);
    dummy->SetLineColor(kFillColor3);
    leg->AddEntry(dummy->Clone(),"3#sigma","l");
  }
  if(bestFit){
    leg->SetNColumns(4);
    //leg->SetTextSize(0.8*kBlessedLabelFontSize);
//  leg->SetColumnSeparation(0.08);
    leg->SetX2(0.75);
    dummy->SetMarkerStyle(43);
    dummy->SetMarkerSize(2);
    dummy->SetMarkerColor(kDarkColor);
    dummy->SetLineColor(kDarkColor);
    dummy->SetLineStyle(7);
    if(hie>0) leg->AddEntry(dummy->Clone(),"Best Fit","p");
    if(hie<0) leg->AddEntry((TObject*)0, "#color[0]{Best Fit}", "");
    dummy->SetLineStyle(1);
  }
  if(!fccorr) leg->SetHeader("No Feldman-Cousins");
  else leg->SetHeader("Feldman-Cousins");
  return leg;
}

void CylindricalSmoothing

(

TH2 *

h

)

Definition at line 63 of file joint_fit_2018_tools.h.

References MECModelEnuComparisons::i, calib::j, maxY, minY, APDHVSetting::temp, x1, submit_syst::x2, submit_syst::y, y1, and submit_syst::y2.

Referenced by joint_fit_2018_contours(), joint_fit_2019_contours(), nova_official_data_release(), and plot_joint_fit_2020_contours().

{
  TH2* temp = (TH2*)h->Clone();

  int NX = h->GetXaxis()->GetNbins();
  int NY = h->GetYaxis()->GetNbins();
  double minX = h->GetXaxis()->GetBinLowEdge(1);
  double maxX = h->GetXaxis()->GetBinUpEdge(NX);
  double widthX = h->GetXaxis()->GetBinWidth(1);
  double minY = h->GetYaxis()->GetBinLowEdge(1);
  double maxY = h->GetYaxis()->GetBinUpEdge(NY);
  double widthY = h->GetYaxis()->GetBinWidth(1);

  //  assert(NX % 2 == 0 && 
  //     "Histogram must have an even number of X bins for this to work");

  // Save the location we expect the end column of the original
  // TH2 to be after we do the folding
  int foldedBin = 1 + NX/2;
 
  // Fold the initial histogram to join the ends of the X axis
  for(int i = 1; i <= NX/2; i++) {
    for(int j = 1; j <= NY/2; j++) {
      // x1, y1 denotes a position in the left-hand side half plane
      // x2, y2 denotes a position in the right-hand side half plane
      int x1 = i+NX/2;
      int x2 = i%(NX/2+1)+2;
      int y1 = j;
      int y2 = j;
      double z1 = h->GetBinContent(x1, y1);
      double z2 = h->GetBinContent(x2, y2);
      temp->SetBinContent(x1, y1, z2);
      temp->SetBinContent(x2, y2, z1);
    }
  }
  // smooth both histograms
  temp->Smooth();
  h->Smooth();
  
  // Replace end columns
  for(int y = 1; y <= NY; y++) {
    double sZ1 = temp->GetBinContent(foldedBin, y);
    double sZ2 = temp->GetBinContent(foldedBin, y);
    h->SetBinContent(1, y, sZ1);
    h->SetBinContent(NX, y, sZ2);
  }
}

void DrawContourCanvas	(	TString	surfName,
		double	minx = 0,
		double	maxx = 2,
		double	miny = 0,
		double	maxy = 1
	)

Definition at line 606 of file joint_fit_2018_tools.h.

References file_size_ana::axes, demo5::c1, ana::kBlessedLabelFontSize, ana::kBlessedTitleFontSize, maxy, plotROC::title, ana::UniqueName(), and ana::XAxisDeltaCPLabels().

                                                                           {
  auto c1 = new TCanvas(ana::UniqueName().c_str());
  c1->SetFillStyle(0);
  c1->SetLeftMargin(0.14);
  c1->SetBottomMargin(0.15);
  TH2* axes = new TH2F();
  TString title;
  if(surfName.Contains("delta_th23")) title=";#delta_{CP};sin^{2}#theta_{23}";
  if(surfName.Contains("th13_delta")) title=";sin^{2}2#theta_{13};#delta_{CP}/#pi";
  if(surfName.Contains("th23_dm32"))  title=";sin^{2}#theta_{23};#Deltam^{2}_{32} (10^{-3}eV^{2})";
  axes = new TH2F("",title,100,minx,maxx,100,miny,maxy);
//  CenterTitles(axes);
  axes->GetXaxis()->CenterTitle();
  axes->GetYaxis()->CenterTitle(); 
  axes->Draw();
  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
  axes->GetXaxis()->SetTitleOffset(0.8);
  axes->GetYaxis()->SetTitleOffset(0.8);
  TGaxis::SetMaxDigits(3);
  if(surfName.Contains("th23_dm32")) {
    axes->GetYaxis()->SetDecimals() ;
    axes->GetYaxis()->SetTitleOffset(0.85);
    axes->GetYaxis()->SetLabelOffset(0.002);
    axes->GetYaxis()->SetLabelSize(0.058);
    axes->GetYaxis()->SetTitleSize(0.078);
    axes->GetXaxis()->SetTitleSize(.95*kBlessedTitleFontSize);
    axes->GetXaxis()->SetTitleOffset(0.86);
  }
  if(surfName.Contains("delta_th23"))  XAxisDeltaCPLabels(axes);
  c1->RedrawAxis();
}

void DrawSliceCanvas	(	TString	slicename,
		const double	ymax,
		const double	xmin = 0,
		const double	xmax = 2.,
		bool	overlay = false
	)

Definition at line 536 of file joint_fit_2018_tools.h.

References file_size_ana::axes, demo5::c1, ana::kBlessedLabelFontSize, ana::kBlessedTitleFontSize, allInOneTrainingPlots::overlay, plotROC::title, ana::UniqueName(), ana::XAxisDeltaCPLabels(), xmax, and make_mec_shifts_plots::xmin.

                                                                                          {
  auto c1 = new TCanvas(ana::UniqueName().c_str());
  c1->SetFillStyle(0);
  c1->SetLeftMargin(0.14);
  c1->SetBottomMargin(0.15);
  TString title;
  if(slicename.Contains("delta")) {
    title = ";delta;Significance (#sqrt{#Delta#chi^{2}})";
  }
  if(slicename.Contains("th23")) {
    title = ";sin^{2}#theta_{23};Significance (#sqrt{#Delta#chi^{2}})";
  }
  if(slicename.Contains("dmsq")) {
    title = ";|#Deltam^{2}_{32}| (10^{-3} eV^{2});Significance (#sqrt{#Delta#chi^{2}})";
  }
  if(slicename.Contains("th13")) {
    title = ";sin^{2}2#theta_{13};Significance (#sqrt{#Delta#chi^{2}})";
  }
  auto axes = new TH2F("",title, 100, xmin, xmax, 100, 0, ymax);
  if(slicename.Contains("fccorr") || overlay) axes->GetYaxis()->SetTitle("Significance (#sigma)");
  // CenterTitles(axes);
  axes->GetXaxis()->CenterTitle();
  axes->GetYaxis()->CenterTitle();
  axes->Draw();
  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
  axes->GetXaxis()->SetTitleOffset(0.8);
  axes->GetYaxis()->SetTitleOffset(0.75);
  if(slicename.Contains("delta")) XAxisDeltaCPLabels(axes);
  c1->RedrawAxis();
 }

TLegend* FCOverlayContourLegend	(	int	hie,
		bool	fillContour,
		Int_t	kFillColor1,
		Int_t	kFillColor2,
		Int_t	kFillColor3,
		Int_t	kDarkColor
	)

Definition at line 696 of file joint_fit_2018_tools.h.

References APDHVSetting::dummy, ana::k2SigmaConfidenceStyle, ana::kBlessedLabelFontSize, and MECModelEnuComparisons::leg.

{
  TLegend * leg = new TLegend(0.16,hie>0?.78:0.77,0.65,0.86);
  leg->SetTextSize(kBlessedLabelFontSize);
  leg->SetFillStyle(0);
  leg->SetMargin(1.3*leg->GetMargin());
  leg->SetNColumns(3);

  TGraph * dummy = new TGraph;
  dummy->SetLineWidth(3);

  if(fillContour){
    dummy->SetLineColor(kDarkColor);
    dummy->SetFillColor(kFillColor1);
    leg->AddEntry(dummy->Clone(),"1#sigma","f");
    dummy->SetFillColor(kFillColor2);
    leg->AddEntry(dummy->Clone(),"2#sigma","f");
    dummy->SetFillColor(kFillColor3);
    leg->AddEntry(dummy->Clone(),"3#sigma","f");
  }
  else{
    leg->SetMargin(1.4*leg->GetMargin());
    dummy->SetLineColor(kFillColor1);
    dummy->SetLineStyle(kSolid);
    leg->AddEntry(dummy->Clone(),"1#sigma","l");
    dummy->SetLineStyle(k2SigmaConfidenceStyle);
    dummy->SetLineColor(kFillColor2);
    leg->AddEntry(dummy->Clone(),"2#sigma","l");
    dummy->SetLineStyle(kSolid);
    dummy->SetLineColor(kFillColor3);
    leg->AddEntry(dummy->Clone(),"3#sigma","l");
  }

  leg->SetHeader("No Feldman-Cousins");
  return leg;
}

std::vector<std::pair<const ISyst*,const ISyst*> > GetCorrelations	(	bool	isNue,
		bool	isFHC
	)

Definition at line 489 of file joint_fit_2018_tools.h.

References om::cout, ana::GetJointFitSystematicList(), and runNovaSAM::ret.

                                                                                       {
  std::vector<std::pair<const ISyst*,const ISyst*> > ret ={};
  std::cout<<(isNue?"not nue ":"not numu ")<<(isFHC&&!isNue?" ":(isFHC?" FHC ":" RHC "));
  //to first approximation if they are not the same they are not correlated.
  //might have to be more complicated later
  auto badsyst = GetJointFitSystematicList(true, !isNue, isNue, isFHC, !isFHC, false);
  for (auto & s:badsyst) ret.push_back ({s,NULL});

  return ret;
}

Spectrum* GetFakeData	(	const IPrediction *	pred,
		osc::IOscCalc *	calc,
		const double	pot,
		TH1D *	cosmics = 0
	)

Definition at line 424 of file joint_fit_2018_tools.h.

References make_associated_cosmic_defs::cosmics, ana::Spectrum::FakeData(), ana::IPrediction::Predict(), and APDHVSetting::temp.

{
  auto temp = pred->Predict(calc).FakeData(pot);
  if(cosmics) temp += Spectrum(cosmics, pot,0);
  return new Spectrum(temp);
}

TFile* GetGaussianSurface	(	TString	par,
		TString	hier
	)

Definition at line 45 of file joint_fit_2018_tools.h.

References om::cout, allTimeWatchdog::endl, runNovaSAM::ret, and cvnie::subdir.

Referenced by nova_official_data_release().

{
  TString subdir = par.Contains("ssth23dmsq32")? "th23_dmsq" : "delta_th23";
  TString gaussDir = "/nova/ana/nu_e_ana/Ana2018/Results/RHC_and_FHC/contours/"+subdir+"/syst/";
  TString gaussName = "hist_contours_2018_joint_realData_both_only";
  gaussName += hier;
  if(par.Contains("ssth23dmsq32"))
    gaussName += "combo_dmsq_systs.root";
  else
    gaussName += "combo_systs.root";
  

  std::cout << "Opening " << gaussDir + gaussName << std::endl;
  TFile* ret = TFile::Open(gaussDir + gaussName);
  return ret;
}

std::vector<const ISyst*> GetJointFitSystematicList	(	bool	corrSysts,
		bool	nueExclusive = false,
		bool	numuExclusive = false,
		bool	isFHC = true,
		bool	isRHC = true,
		bool	intersection = true
	)

Definition at line 452 of file joint_fit_2018_tools.h.

References ana::AddJointAna2018NotCorrelSysts(), ana::AddNonLoadable2018Systs(), om::cout, allTimeWatchdog::endl, ana::getAllAna2018Systs(), isFHC, make_mec_shifts_plots::isRHC, ana::k2ndClassCurrs, ana::kAna2018NormFHC, ana::kAna2018NormRHC, ana::kBoth, ana::kFHC, ana::kJointAna2018, ana::kNueAna2018, ana::kNumuAna2018, ana::kRadCorrNue, ana::kRadCorrNuebar, ana::kRHC, and runNovaSAM::ret.

                                                                                                                                                                                       {
  std::vector<const ISyst*> ret ={} ;
  if(corrSysts){
    if(! intersection){
      if(nueExclusive){
        
         ret.push_back(&kRadCorrNue);
         ret.push_back(&kRadCorrNuebar);
         ret.push_back(&k2ndClassCurrs);
         AddJointAna2018NotCorrelSysts(ret, kNueAna2018, kBoth);
         AddNonLoadable2018Systs(ret, kNueAna2018);

      }
      if(numuExclusive){
        
         AddJointAna2018NotCorrelSysts(ret, kNumuAna2018, kBoth);
         if(isFHC) AddJointAna2018NotCorrelSysts(ret, kNueAna2018, kRHC);
         if(isRHC) AddJointAna2018NotCorrelSysts(ret, kNueAna2018, kFHC);
         if(isFHC) ret.push_back(&kAna2018NormRHC);
         if(isRHC) ret.push_back(&kAna2018NormFHC);
         AddNonLoadable2018Systs(ret, kNumuAna2018);

      }
    }
    else{
      if(nueExclusive) ret = getAllAna2018Systs(kNueAna2018);
      if(numuExclusive) ret = getAllAna2018Systs(kNumuAna2018);
      if(!nueExclusive && !numuExclusive) ret = getAllAna2018Systs(kJointAna2018, true, kBoth, true);
      //std::cout <<"I'm in else" <<std::endl;
      //if(!nueExclusive && !numuExclusive) ret = getAllAna2018Systs(kJointAna2018, true, true);
    }
  }
  std::cout << "Return list of  systematics: " ;
  for(auto & s:ret) std::cout << s->ShortName() << ", ";
  std::cout << std::endl;
  return ret;
}

Spectrum* GetMockData	(	const IPrediction *	pred,
		osc::IOscCalc *	calc,
		const double	pot,
		TH1D *	cosmics = 0
	)

Definition at line 432 of file joint_fit_2018_tools.h.

References make_associated_cosmic_defs::cosmics, ana::Spectrum::MockData(), ana::IPrediction::Predict(), and APDHVSetting::temp.

{
  auto temp = pred->Predict(calc).MockData(pot);
  if(cosmics) temp += Spectrum(cosmics, pot,0);
  return new Spectrum(temp);
}

std::pair<TH1D*, double > GetNueCosmics2018	(	std::string	beam,
		bool	GetFromUPS = false,
		bool	NERSC = false
	)

Definition at line 240 of file joint_fit_2018_tools.h.

References om::cout, allTimeWatchdog::endl, plot_validation_datamc::fname, cet::getenv(), make_syst_table_plots::h, ana::Spectrum::Integral(), ana::kAna2018FHCLivetime, ana::kAna2018RHCLivetime, ana::kLivetime, livetime, outtime(), runNovaSAM::release, std::sqrt(), string, and ana::Spectrum::ToTH1().

Referenced by bin_composition_pie_chart(), demoFitContours(), demoFitSlices(), fake_future_data(), getHists_FNEX(), joint_fit_2018_contours(), joint_fit_2018_slices(), jointsensitivity(), make_fc_mass_and_oct_nersc_2018(), make_fc_mh_nersc_2018(), make_fc_oct_nersc_2018(), make_fc_slices_nersc_2018(), make_fc_slices_nersc_2018_stats(), make_fc_surfaces_nersc_2018(), make_fc_surfaces_nersc_2018_stats(), MakeBkgEMPlot(), MakeCAFSensitivities_for_FNEX(), monoprob(), nuebar_signif(), plot_datamcpred(), sensitivity2018(), and syst_table_fit().

                                                                                                    {
  Spectrum* outtime;
  if (GetFromUPS){
    std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));            //update
    std::string fname = "/nue2017_cospred.root";                                      //update
    std::string name  = "cos";       
    outtime = LoadFromFile<Spectrum>(upsname+fname,name).release();
  }
  else{
    std::string name;
    std::string anaDir = "/nova/ana/nu_e_ana/Ana2018/";
    if(NERSC) anaDir = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
    if(beam=="rhc") name = "cosmic_spect_rhc";
    else name="cosmic_spect_fhc";
    outtime = LoadFromFile<Spectrum>
      (anaDir+"/Predictions/cosmic/cosmic_prediction_real_data.root", name ).release();
  }

  double livetime;
  if(beam=="rhc") livetime = kAna2018RHCLivetime ;
  else livetime = kAna2018FHCLivetime ;

  //double outN = outtime->Integral(outtime->Livetime(), 0, kLivetime);
  double outN = outtime->Integral(livetime, 0, kLivetime);

  std::cout << "\nAdding " << outtime->Integral(livetime, 0, kLivetime)
            << " cosmics from out-of-time." << std::endl;
  std::cout << "Which scale to " << outtime->Integral(livetime, 0,
                                                      kLivetime)
            << " in the spill window." << std::endl;

  std::cout << "uncertainty " << 1/sqrt(outN) << std::endl;

  TH1 *h = outtime->ToTH1(livetime, kBlack, kSolid, kLivetime);
//  assert( h->Integral() > 4 && "Nue cosmic prediction < 4. Something is wrong");    //RHC has ~0.6 cosmic events

  return {outtime->ToTH1(livetime, kBlack, kSolid, kLivetime), 1/sqrt(outN)};
}

Spectrum* GetNueData2018

(

std::string

beam

)

Definition at line 281 of file joint_fit_2018_tools.h.

References om::cout, allTimeWatchdog::endl, exit(), ana::Spectrum::Integral(), intime(), ana::kAna2018FHCLivetime, ana::kAna2018FHCPOT, ana::kAna2018RHCLivetime, ana::kAna2018RHCPOT, livetime, ana::Spectrum::Livetime(), POT, ana::Spectrum::POT(), and runNovaSAM::release.

Referenced by joint_fit_2018_contours(), joint_fit_2018_slices(), MakeCAFSensitivities_for_FNEX(), and nuebar_signif().

                                         {
  Spectrum* intime;
  
  if(beam=="fhc"){
    intime = LoadFromFile<Spectrum>
      ("/nova/ana/nu_e_ana/Ana2018/Data2018/spectra.root", "data_spect_fhc").release();
  }
  else if(beam=="rhc"){
    intime = LoadFromFile<Spectrum>
      ("/nova/ana/nu_e_ana/Ana2018/Data2018/spectra.root", "data_spect_rhc").release();
  }
  else{ std::cout << "Beam mode != fhc or rhc" << std::endl; exit(1);}
  double POT, livetime;
  if(beam == "rhc") {livetime = kAna2018RHCLivetime; POT = kAna2018RHCPOT;}
  if(beam == "fhc") {livetime = kAna2018FHCLivetime; POT = kAna2018FHCPOT;}

  std::cout << "Loaded data spectrum:\n"
            << "POT " << intime->POT()
            << " (expected " << POT  << ")\n"
            << "Livetime " << intime->Livetime()
            << " (expected " << livetime << ")\n"
            << "Integral " << intime->Integral(intime->POT())
            << std::endl;

  return intime;
}

const IPrediction* GetNuePrediction2018	(	std::string	decomp,
		osc::IOscCalc *	calc,
		bool	corrSysts,
		std::string	beam,
		bool	isFakeData,
		bool	GetFromUPS = false,
		bool	minimizeMemory = false,
		bool	NERSC = false
	)

Definition at line 136 of file joint_fit_2018_tools.h.

References POTSpillRate::beam, om::cout, dirName, allTimeWatchdog::endl, exit(), ana::getAllAna2018Systs(), cet::getenv(), GetNuePredPath(), ana::Spectrum::Integral(), ana::kAna2018FHCPOT, ana::kAna2018RHCPOT, ana::kCombo, ana::kFake, ana::kFHC, ana::kNoExtrap, ana::kNueAna2018, ana::kProportional, ana::kRHC, pot, ana::PredictionInterp::Predict(), runNovaSAM::release, rock(), and string.

Referenced by bicount_2018(), bin_composition_pie_chart(), demoFitContours(), demoFitSlices(), fake_future_data(), GetAllPredictions(), getHists_FNEX(), joint_fit_2018_contours(), joint_fit_2018_slices(), jointsensitivity(), make_fc_mass_and_oct_nersc_2018(), make_fc_mh_nersc_2018(), make_fc_oct_nersc_2018(), make_fc_slices_nersc_2018(), make_fc_slices_nersc_2018_stats(), make_fc_surfaces_nersc_2018(), make_fc_surfaces_nersc_2018_stats(), MakeCAFSensitivities_for_FNEX(), monoprob(), nuebar_signif(), plot_datamcpred(), sensitivity2018(), syst_table_fit(), systematics_extrap_comp_from_pred_interp(), and systematics_summary_from_pred_interp().

                                                            {

  ENu2018ExtrapType extrap = kProportional;
  if(decomp == "noextrap") extrap = kNoExtrap;
  if(decomp == "prop") extrap = kProportional;
  if(decomp == "combo") extrap = kCombo;
  if(decomp == "fakexp") extrap = kFake;//Hack 

  PredictionSystJoint2018 * pred_fid;
  std::string nue_pred_path = GetNuePredPath(beam, isFakeData, NERSC);
  // Interactive running
  if (!GetFromUPS){
    if(corrSysts){
      std::cout << "Getting Nue Prediction from " + nue_pred_path << endl;
      if(beam=="rhc") pred_fid = new PredictionSystJoint2018 (extrap, calc, 
                                                              getAllAna2018Systs( kNueAna2018, true, kRHC), 
                                                              beam, isFakeData, true, 
                                                              nue_pred_path, minimizeMemory);
      else if(beam=="fhc") pred_fid = new PredictionSystJoint2018 (extrap, calc,
                                                              getAllAna2018Systs( kNueAna2018, true, kFHC),
                                                              beam, isFakeData, true,
                                                              nue_pred_path, minimizeMemory);
      else{ std::cout << "Beam mode needs to be either fhc or rhc" << endl; exit(1);}
    }
    else pred_fid = new PredictionSystJoint2018 (extrap, calc,{}, beam, isFakeData,
                                                 true, nue_pred_path, minimizeMemory);
  }
  // Grid running for FC
  else{
    std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));                       // need an update
    std::string name = "/pred_nue_reduced_v9.root";                                              //need an update
    if(corrSysts){
      if(beam=="rhc") pred_fid = new PredictionSystJoint2018(extrap, calc,  
                                                             getAllAna2018Systs( kNueAna2018, true, kRHC),
                                                             beam, isFakeData, true,
                                                             upsname+name, minimizeMemory);
      else if(beam=="fhc") pred_fid = new PredictionSystJoint2018(extrap, calc,
                                                                  getAllAna2018Systs( kNueAna2018, true, kFHC),
                                                                  beam, isFakeData, true,
                                                                  upsname+name, minimizeMemory);
      else{ std::cout << "Beam mode should be either fhc or rhc" << endl; exit(1);}
    }
    else
      pred_fid = new PredictionSystJoint2018(extrap, calc, {},beam, isFakeData,
                                             true, upsname+name, minimizeMemory);
  }
  


  double pot;
  if(beam=="rhc") pot = kAna2018RHCPOT;
  else pot = kAna2018FHCPOT;
  std::cout << "Scale to "<<pot<<std::endl;

  std::cout << "Nue prediction (no rock)"
            << pred_fid->Predict(calc).Integral(pot) <<std::endl;
  // Load rock

  double nonsScale = 1./10.45;
  double swapScale = 1./12.91;
  //fix to the normal path later
  std::string dirName = "/nova/ana/nu_e_ana/Ana2018/Predictions/FDRock/";
  if(NERSC) dirName = std::string(getenv("FCHELPERANA2018_LIB_PATH"))+"/Predictions/FDRock/";
  std::string fName2 = "fdrock_nue2018.root";
  std::string rName;

  if(beam=="rhc") rName = "pred_FDRock_merg_rhc";
  else rName = "pred_FDRock_merg_fhc";
  if( decomp == "noextrap"){
    if(beam=="rhc") rName = "pred_FDRock_rhc";
    else rName = "pred_FDRock_fhc";
  }
  // Point to the rock file in the ups product for FC'ing
  if (GetFromUPS){
    dirName = std::string(getenv("FCHELPERANA2017_LIB_PATH"));                                  //need an update
    fName2  = "/fdrock_nue2017.root";                                                           //need an update
    rName = "pred_FDRock";
  }

  auto rock = LoadFromFile<PredictionNoExtrap>(dirName + fName2, rName).release();

  std::cout << "Nue prediction (rock only - not scaled) "
            << rock->Predict(calc).Integral(pot) <<std::endl;
  std::cout << "Adding Rock" << endl;
  auto nue_pred = new PredictionAddRock (pred_fid, rock, nonsScale, swapScale);

  std::cout << "Nue prediction (rock added) "
            <<  nue_pred->Predict(calc).Integral(pot) << std::endl;
  std::cout << "Pure rock events:  "
            <<  nue_pred->Predict(calc).Integral(pot) - pred_fid->Predict(calc).Integral(pot) << std::endl;
  
  
  return nue_pred;
}

std::string GetNuePredPath	(	std::string	beam,
		bool	isFakeData,
		bool	NERSC
	)

Definition at line 111 of file joint_fit_2018_tools.h.

References dir, cet::getenv(), isFHC, runNovaSAM::ret, and string.

Referenced by GetNuePrediction2018().

                                                                      {
  std::string ret;
  bool isFHC = true;
  if(beam == "rhc") isFHC = false;

  // default path name for Nue Preds constructor
  // to point to official files
  if(!NERSC) return "";

  // NERSC paths
  else {
    std::string dir = std::string(getenv("FCHELPERANA2018_LIB_PATH"))
      + "/Predictions/provisional_singles";
    if(isFHC && !isFakeData)
      ret = dir + "/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_FHCAllSyst_nueconcat_realNDData_2018-05-22.root";
    else if(!isFHC && !isFakeData)
      ret = dir + "/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_RHCAllSyst_nueconcat_realNDData_2018-05-22.root";
    else if(isFHC && isFakeData)
      ret = dir + "/NueProd4SystematicsOnFakeNDData_2018-05-03/pred_allxp_Nue2018Axis_full_FHCAllSysts_nueconcat_fakeNDData_05_03_2018_merged.root";
    else if(!isFHC & isFakeData)
      ret = dir + "/NueProd4SystematicsOnFakeNDData_2018-05-03/pred_allxp_Nue2018Axis_full_RHCAllSysts_nueconcat_fakeNDData_05_03_2018_merged.root";
  }
  return ret;
}

std::vector<std::pair <TH1D*, double> > GetNumuCosmics2018	(	const int	nq = 4,
		std::string	beam = "fhc",
		bool	GetFromUPS = false,
		bool	NERSC = false
	)

Definition at line 364 of file joint_fit_2018_tools.h.

References POTSpillRate::beam, om::cerr, om::cout, dir, allTimeWatchdog::endl, exit(), shutoffs::filename, cet::getenv(), make_syst_table_plots::h, MECModelEnuComparisons::i, std::sqrt(), string, art::to_string(), and ana::UniqueName().

Referenced by demoFitContours(), demoFitSlices(), fake_future_data(), joint_fit_2018_contours(), joint_fit_2018_slices(), jointsensitivity(), make_fc_mass_and_oct_nersc_2018(), make_fc_mh_nersc_2018(), make_fc_oct_nersc_2018(), make_fc_slices_nersc_2018(), make_fc_slices_nersc_2018_stats(), make_fc_surfaces_nersc_2018(), make_fc_surfaces_nersc_2018_stats(), MakeCAFSensitivities_for_FNEX(), sensitivity2018(), and syst_table_fit().

{
   cout<<"Beam is "<<beam<<endl;
  //cosmics
  //the versionn in cafana has float not double and it's a pain
  TFile *fcosm;
  if (GetFromUPS){
    std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));          //update
    fcosm = TFile::Open((upsname+"/cos_bkg_hists_v2.root").c_str());                //update
  }
  else{
    std::string dir = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
    if(NERSC) dir = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
    std::string filename;
    if (beam=="rhc") filename="cosmics_rhc__numu2018.root";
    else filename="cosmics_fhc__numu2018.root";
    fcosm = TFile::Open((dir+filename).c_str());
  }

  if(fcosm->IsZombie()) {std::cerr<< "bad cosmics\n"; exit(1);}
  std::vector <std::pair <TH1D*, double > > numu_cosmics;

  for (int i = 1; i <=nq; ++i) {
    auto h = (TH1D*) fcosm->Get((TString)"cosmics_q"+std::to_string(i))
      ->Clone(UniqueName().c_str());
    //auto hs = (TH1D*) fcosm->Get((TString)"q"+std::to_string(i)+"allsyst");        //update. no such stuff in numu files
    //numu_cosmics.push_back({h,hs->GetBinContent(1)});
    numu_cosmics.push_back({h,1/sqrt(h->Integral())}); 
    std::cout << "Cosmics all periods quantile " << i << " "
              << h->Integral()
              << " scale error "
              //<< hs->GetBinContent(1)
              << 1/sqrt(h->Integral())
              << "\n";
  }//end quantiles
  return numu_cosmics;
}

std::vector<Spectrum * > GetNumuData2018	(	const int	nq = 4,
		std::string	beam = "fhc"
	)

Definition at line 404 of file joint_fit_2018_tools.h.

References POTSpillRate::beam, om::cerr, exit(), MakeMiniprodValidationCuts::f, file, MECModelEnuComparisons::i, ana::Spectrum::LoadFrom(), string, and art::to_string().

Referenced by joint_fit_2018_contours(), joint_fit_2018_slices(), and MakeCAFSensitivities_for_FNEX().

                                                                                 {
 
  std::string file;
  if(beam == "fhc") {file="/nova/ana/nu_mu_ana/Ana2018/Data/fd_dataspectrum_fhc_full__numu2018.root";}
  if(beam == "rhc") {file="/nova/ana/nu_mu_ana/Ana2018/Data/fd_dataspectrum_rhc_full__numu2018.root";}

  std::vector <Spectrum *> numu_data;

  for (int i = 1; i <=nq; ++i) {

    auto f= new TFile(file.c_str());
    if(f->IsZombie()) {std::cerr<< "bad data\n"; exit(1);}
    numu_data.push_back(Spectrum::LoadFrom(f, ("fd_data_q"+std::to_string(i)).c_str())).release();

  }//end quantiles
  return numu_data;
}

std::vector<const IPrediction*> GetNumuPredictions2018	(	const int	nq = 4,
		bool	useSysts = true,
		std::string	beam = "fhc",
		bool	GetFromUPS = false,
		ENu2018ExtrapType	numuExtrap = kNuMu,
		bool	minimizeMemory = false,
		bool	NERSC = false
	)

Definition at line 310 of file joint_fit_2018_tools.h.

References POTSpillRate::beam, calc, om::cout, ana::DefaultOscCalc(), shutoffs::filename, plot_validation_datamc::fname, ana::getAllAna2018Systs(), cet::getenv(), ana::kNumuAna2018, and string.

Referenced by demoFitContours(), demoFitSlices(), fake_future_data(), GetAllPredictions(), joint_fit_2018_contours(), joint_fit_2018_slices(), jointsensitivity(), make_fc_mass_and_oct_nersc_2018(), make_fc_mh_nersc_2018(), make_fc_oct_nersc_2018(), make_fc_slices_nersc_2018(), make_fc_slices_nersc_2018_stats(), make_fc_surfaces_nersc_2018(), make_fc_surfaces_nersc_2018_stats(), MakeCAFSensitivities_for_FNEX(), sensitivity2018(), and syst_table_fit().

{
  auto calc = DefaultOscCalc();

  std::vector <const IPrediction*> numu_preds;
 
  std::string filename;
  std::string anaDir = "/nova/ana/nu_mu_ana/Ana2018/";
  if(NERSC) anaDir = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
  if(beam=="rhc") filename=anaDir + "/Predictions/pred_numuconcat_rhc__numu2018.root";
  else filename= anaDir + "/Predictions/pred_numuconcat_fhc__numu2018.root";

  std::cout << "\nLoading Numu Predictions\n" ;
  for (int quant = 0; quant < nq; ++quant){
    if (GetFromUPS){
      std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));     //update
      std::string fname = "/pred_numu_reduced_v5.root";                          //update
      if(!useSysts) {      
        numu_preds.push_back(new PredictionSystJoint2018(numuExtrap,
                                                         calc, beam,
                                                         quant+1,{},
                                                         upsname+fname,
                                                         minimizeMemory));         //update
      }
      else{
        numu_preds.push_back(new PredictionSystJoint2018(numuExtrap, calc, beam,        
                                                         quant+1,
                                                         getAllAna2018Systs( kNumuAna2018, true),
                                                         upsname+fname,
                                                         minimizeMemory));          //update 
      }
    }
    else{
      if(!useSysts) {
        numu_preds.push_back(new PredictionSystJoint2018(numuExtrap, calc, beam, quant+1, 
                                                         {}, filename, minimizeMemory));
      }
      else {
        numu_preds.push_back(new PredictionSystJoint2018(numuExtrap, calc, beam, quant+1,
                                                         getAllAna2018Systs( kNumuAna2018, true), 
                                                         filename, minimizeMemory));
      }
    }
  }
  return numu_preds;
}

void SetFakeCalc	(	osc::IOscCalcAdjustable *	calc,
		double	ssth23 = -5,
		double	dmsq32 = -5,
		double	dCP_Pi = -5
	)

Definition at line 442 of file joint_fit_2018_tools.h.

References dmsq32, ana::kFitDeltaInPiUnits, ana::kFitDmSq32, ana::kFitSinSqTheta23, ana::FitVarWithPrior::SetValue(), ana::ConstrainedFitVarWithPrior::SetValue(), ana::FitDmSq32::SetValue(), and ssth23.

{
  if(dCP_Pi > -1) kFitDeltaInPiUnits.SetValue(calc, dCP_Pi);
  if(ssth23 > -1) kFitSinSqTheta23.SetValue(calc, ssth23);
  if(dmsq32 > -1) kFitDmSq32.SetValue(calc, dmsq32);
}

TLegend* SliceLegend	(	std::vector< std::pair< TGraph *, TString > > &	graphs,
		bool	isDelta,
		double	low = 0.7,
		bool	isDmsq = false,
		bool	overlay = false
	)

Definition at line 571 of file joint_fit_2018_tools.h.

References MECModelEnuComparisons::graphs, ana::kBlessedLabelFontSize, MECModelEnuComparisons::leg, allInOneTrainingPlots::overlay, and submit_syst::str.

{
  TLegend * leg;
  if(isDelta) leg = new TLegend(0.6, 0.6, 0.9,0.89);
  else {
    leg = new TLegend(0.18, 0.2, 0.4, 0.49);
  }
  leg->SetTextSize(0.8*kBlessedLabelFontSize);
  leg->SetFillStyle(0);
  for (auto &gr:graphs){
       TString str="";
    if(gr.second.Contains("O")){
      if(isDelta) leg->SetX1(0.57);
      if(isDmsq) leg->SetX1(0.18);
      if(!isDelta && !isDmsq) {leg->SetX1(0.6); leg->SetX2(0.89);}
      leg->SetMargin(0.18);
      str += gr.second.Contains("NH") ? "NH":"IH";
      str +=  gr.second.Contains("LO") ? " Lower": " Upper";
      str += " octant";
    }
    else{
      str="#splitline{";
      str += gr.second.Contains("NH") ? "Normal":"Inverted";
      str += "}{hierarchy}";    
      if(overlay) {
        str = "Uncorr. ";
        str += gr.second.Contains("NH") ? "NH" : "IH";
      }
    }
    leg->AddEntry(gr.first, str, "l");
  }
  return leg;
}

Variable Documentation

const IConstrainedFitVar* kFitSinSqTheta23LowerOctant = new FitSinSqTheta23LowerOctant()

Definition at line 514 of file joint_fit_2018_tools.h.

const IConstrainedFitVar* kFitSinSqTheta23UpperOctant = new FitSinSqTheta23UpperOctant()

Definition at line 528 of file joint_fit_2018_tools.h.