MakeSurfaceJoint.C
Go to the documentation of this file.
1 
2 #include "CAFAna/Core/Progress.h"
3 #include "CAFAna/Core/Sample.h"
10 #include "CAFAna/Analysis/Calcs.h"
12 #include "CAFAna/Analysis/CovMxSurface.h"
13 
14 #include "CAFAna/covmx/Utilities.h"
15 
16 using namespace ana;
17 
19 {
20  DontAddDirectory guard;
21 
22  // Define exposures
23  double potND = kAna2018SensitivityFHCNDPOT;
24  double potFD = kAna2018FHCPOT;
25  double livetimeFD = kAna2018FHCLivetime;
26 
27  /////////////////////////////////////////////////////////////////////
28  // //
29  // ARGUMENT PARSING //
30  // //
31  /////////////////////////////////////////////////////////////////////
32 
33  std::vector<covmx::Sample> samples;
34  std::vector<bool> unoscillated(samples.size(), false);
37  samples.push_back(covmx::Sample(covmx::kNC, pol, det));
38  unoscillated.push_back(det == covmx::Detector::kNearDet);
39  samples.push_back(covmx::Sample(covmx::Analysis::kCCNue, pol, det));
40  unoscillated.push_back(det == covmx::Detector::kNearDet);
42  samples.push_back(covmx::Sample(covmx::Analysis::kCCNumu, pol, det, q));
43  unoscillated.push_back(det == covmx::Detector::kNearDet);
44  }
45  }
46  }
47  for (size_t i = 0; i < samples.size(); ++i) {
48  std::cout << "Sample: " << samples[i].GetTag() << ", disabling oscillations? " << unoscillated[i] << std::endl;
49  }
50 
51  /////////////////////////////////////////////////////////////////////
52  // //
53  // PREPARE INPUTS //
54  // //
55  /////////////////////////////////////////////////////////////////////
56 
57  // Create sterile oscillation calculator
59  SetNus20Params(calc);
60 
61  // Load predictions
62  std::vector<IPrediction*> preds;
63  for (covmx::Sample s : samples) preds.push_back(LoadPrediction(s, "systs", kUseCVMFS));
64 
65  // Load cosmic spectra
66  std::vector<Spectrum*> cosmic;
67  for (auto sample : samples)
68  cosmic.push_back(LoadCosmic(sample, kUseCVMFS));
69 
70  // Load data
71  std::vector<Spectrum*> data;
72  // if (fakedata_dir) {
73  // for (auto sample : samples)
74  // data.push_back(Spectrum::LoadFrom(fakedata_dir->GetDirectory(sample.GetName().c_str())).release());
75  // }
76  // else {
77  for (size_t i = 0; i < preds.size(); ++i) {
78  double pot = (samples[i].detector == covmx::Detector::kNearDet) ? potND : potFD;
79  data.push_back(new Spectrum(preds[i]->Predict(calc).FakeData(pot)));
80  if (samples[i].detector == covmx::kFarDet)
81  data.back()->OverrideLivetime(livetimeFD);
82  if (cosmic[i]) *data[i] += *cosmic[i];
83  }
84  // }
85 
86  /////////////////////////////////////////////////////////////////////
87  // //
88  // MAKE SURFACE //
89  // //
90  /////////////////////////////////////////////////////////////////////
91 
92  // Create experiment
93  OscCovMxExperiment* expt = nullptr;
94  // if (mx) expt = new OscCovMxExperiment(preds, mx, data, cosmic);
95  // else expt = new OscCovMxExperiment(preds, data, cosmic);
96  expt = new OscCovMxExperiment(preds, data, cosmic);
97  auto noosccalc = DefaultSterileCalc(4);
98  for (int i = 2; i <= 4; ++i) calc->SetDm(i, 0);
99  expt->DisableOscillations(unoscillated, noosccalc);
100 
101  std::vector<const IChiSExperiment*> expts;
102  expts.push_back(new const GaussianConstraint(&kFitSinSq2Theta13Sterile, 0.082, 0.004));
103  expts.push_back(expt);
104  MultiExperiment multiExpt(expts);
105 
106  // Set seed values and fit vars
107  std::vector<const IFitVar*> fitVars;
108  std::map<const IFitVar*, std::vector<double> > seedValues;
109  std::vector<double> fit_seeds;
110 
111  fitVars.push_back(&kFitDelta13InPiUnitsSterile);
112  seedValues[&kFitDelta13InPiUnitsSterile] = { 1.5 };
113 
114  fitVars.push_back(&kFitSinSq2Theta13Sterile);
115  seedValues[&kFitSinSq2Theta13Sterile] = { 0.082 };
116 
117  // Get max length of fit var name, for cout padding
118  TH2D* varHists[fitVars.size()];
119  unsigned int padWidth = (*std::max_element(fitVars.begin(), fitVars.end(),
120  [] (const IFitVar* lhs, const IFitVar* rhs) {
121  return lhs->ShortName().length() < rhs->ShortName().length();
122  }))->ShortName().length();
123 
124  Binning xbins(Binning::Simple(50, 0.3, 0.7));
125  Binning ybins(Binning::Simple(50, 2e-3, 3.3e-3));
126 
127  //size_t part = 1275;
128  //size_t njobs = 2500;
129  size_t part = 10;
130  size_t njobs = 20;
131  if (RunningOnGrid()) {
132  part = JobNumber();
133  njobs = 25;
134  }
135 
136  CovMxSurface surf(xbins, ybins, &multiExpt, calc,
138  fitVars, seedValues, {}, SystShifts::Nominal(),
139  true, nullptr, preds, data, cosmic, part, njobs);
140 
141  // Open output file
142  TFile* f = TFile::Open("surface_joint.root", "recreate");
143  surf.SaveTo(f->mkdir("surface"));
144  delete f;
145 
146 }
size_t JobNumber()
What&#39;s the process number for a grid job?
Definition: Utilities.cxx:366
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void FakeData()
Definition: rootlogon.C:156
bool RunningOnGrid()
Is this a grid (condor) job?
Definition: Utilities.cxx:359
A simple Gaussian constraint on an arbitrary IFitVar.
Adapt the PMNS_Sterile calculator to standard interface.
const FitDmSq32Sterile kFitDmSq32Sterile
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
void MakeSurfaceJoint()
const XML_Char const XML_Char * data
Definition: expat.h:268
const XML_Char * s
Definition: expat.h:262
TString part[npart]
Definition: Style.C:32
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
expt
Definition: demo5.py:34
Sum up livetimes from individual cosmic triggers.
void SetNus20Params(osc::OscCalcSterile *calc)
Definition: Utilities.h:682
const int xbins
Definition: MakePlots.C:82
#define pot
const FitSinSq2Theta13Sterile kFitSinSq2Theta13Sterile
Compare a single data spectrum to the MC + cosmics expectation.
const FitDelta13InPiUnitsSterile kFitDelta13InPiUnitsSterile
std::vector< float > Spectrum
Definition: Constants.h:610
OStream cout
Definition: OStream.cxx:6
const int ybins
Combine multiple component experiments.
void SetDm(int i, double dm)
const std::string & ShortName() const
Definition: IFitVar.h:36
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
const FitSinSqTheta23Sterile kFitSinSqTheta23Sterile
std::unique_ptr< Spectrum > LoadCosmic(covmx::Sample sample, bool cvmfs=true)
Get cosmics for a given sample.
Definition: Utilities.h:145
Interface definition for fittable variables.
Definition: IFitVar.h:16
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const double kAna2018FHCLivetime
Definition: Exposures.h:213
surf
Definition: demo5.py:35