MakeFakeDataLLTest.C
Go to the documentation of this file.
1 /// MakeFakeData.C by J. Hewes <jhewes15@fnal.gov>
2 /// Modified by David Duenas <dueasdf@mail.uc.edu>
3 //#include "CAFAna/Prediction/CovarianceMatrix.h"
7 #include "Utilities/rootlogon.C"
8 #include "CAFAna/Core/Progress.h"
10 
11 #include "TH1.h"
12 #include "TCanvas.h"
13 #include "TMath.h"
14 #include "TLegend.h"
15 #include "TError.h"
16 
17 using namespace ana;
18 
20 
21  SetNus20Params(calc);
22 
23  if (n == 1) {
26  kFitDmSq41Sterile.SetValue(calc, 0.1);
27  } else if (n == 2) {
30  kFitDmSq41Sterile.SetValue(calc, 1.0);
31  } else if (n == 3) {
34  kFitDmSq41Sterile.SetValue(calc, 10.0);
35  } else {
36  throw std::runtime_error("Fake signal number not valid!");
37  }
38 } // function FakeSignal
39 
41 {
42  //MiniBooNe Best Fit values
43  Double_t deltm2_41=0.041;
44  Double_t theta_14=TMath::ASin(TMath::Sqrt(0.5000));
45  Double_t theta_24=TMath::ASin(TMath::Sqrt(0.9580));
46  Double_t theta_34=TMath::ASin(TMath::Sqrt(0.9580));
47  calc->SetDm(4, deltm2_41);
48  calc->SetAngle(1, 4, theta_14);
49  calc->SetAngle(2, 4, theta_24);
50  calc->SetAngle(3, 4, theta_34);
51 
52 }
53 
54 void MakeFakeData(bool rhc = false, bool stats = true)
55 {
56  try { // i know wrapping the entire macro in a try/catch is bad but here we are
57 
58  DontAddDirectory guard;
59  rootlogon();
60 
61  //Number of Universes
62  size_t nUniverses = 1000;
63 
65 
66  // TODO: confirm arguments to SetPrediction
69  nd.SetLivetime(-1);
70  nd.SetBinning(kNCNDBinning);
71  SetPrediction(nd, false); // Get ND prediction without systs
72 
76  fd.SetBinning(kNCFDBinning);
77  SetPrediction(fd, false); // Get FD prediction without systs
78 
79  TFile* f = TFile::Open("fhc_fake_data_lltest.root", "recreate");
80 
81  Spectrum* cosmic = LoadCosmic(kNusFHCFarDet, kUseCVMFS).release();
82  TH1D* hCosmic = cosmic->ToTH1(fd.GetLivetime(), kLivetime);
83 
84  // Predict three-flavour spectra here
85  SetNus20Params(calc);
86 
87  TH1D* hND3Flav = nd.GetPrediction()->Predict(calc).ToTH1(nd.GetPOT());
88  hND3Flav->GetYaxis()->SetMaxDigits(4);
89  hND3Flav->SetLineWidth(2);
90  hND3Flav->SetLineColor(kAzure-2);
91  hND3Flav->GetXaxis()->CenterTitle();
92  hND3Flav->GetYaxis()->CenterTitle();
93 
94  TH1D* hFD3Flav = fd.GetPrediction()->Predict(calc).ToTH1(fd.GetPOT());
95  hFD3Flav->GetYaxis()->SetMaxDigits(4);
96  hFD3Flav->SetLineWidth(2);
97  hFD3Flav->SetLineColor(kGreen+2);
98  hFD3Flav->GetXaxis()->CenterTitle();
99  hFD3Flav->GetYaxis()->CenterTitle();
100  hFD3Flav->Add(hCosmic);
101 
102  // Get list of active systematics
103  std::vector<std::string> systs;
104  for (const ISyst* syst : getAllXsecSysts_2018_RPAFix()) {
105  systs.push_back(syst->ShortName());
106  }
107  for (const ISyst* syst : LoadSysts(nd, kUseCVMFS)) {
108  systs.push_back(dynamic_cast<const NuISyst*>(syst)->BaseName());
109  }
110 
111  // Load matrix from manager
112  std::string mxPath = kNus20Path + "/covmx/covmxmanager.root";
113  TFile* mxFile = TFile::Open(pnfs2xrootd(mxPath).c_str(), "read");
114  CovMxManager* mgr = CovMxManager::LoadFrom(mxFile).release();
116  mx->Predict({nd,fd}, calc);
117 
118  // Get binning information
119  unsigned int nBins = mx->GetBinning().NBins();
120 
121  // Loop over every set of fake oscillation parameters
122  for (size_t u = 1; u < 4; ++u) {
123 
124  TDirectory* uDir = f->mkdir(Form("fs%zu", u));
125 
126  FakeSignal(calc, u);
127  mx->Predict({nd,fd}, calc);
128 
129  TH1D* hND4Flav = nd.GetPrediction()->Predict(calc).ToTH1(nd.GetPOT());
130  hND4Flav->GetYaxis()->SetMaxDigits(4);
131  hND4Flav->SetLineWidth(2);
132  hND4Flav->SetLineColor(kRed-3);
133  hND4Flav->GetXaxis()->CenterTitle();
134  hND4Flav->GetYaxis()->CenterTitle();
135 
136  TH1D* hFD4Flav = fd.GetPrediction()->Predict(calc).ToTH1(fd.GetPOT());
137  hFD4Flav->GetYaxis()->SetMaxDigits(4);
138  hFD4Flav->SetLineWidth(2);
139  hFD4Flav->SetLineColor(kRed-3);
140  hFD4Flav->GetXaxis()->CenterTitle();
141  hFD4Flav->GetYaxis()->CenterTitle();
142  hFD4Flav->Add(hCosmic);
143 
144  // Combine nominal 4-flavour spectra into a single histogram
145 
146  Progress* p = new Progress("Generating fake universes");
147  for (size_t i = 1; i <= nUniverses; ++i) {
148 
149  TCanvas* c3vs4 = new TCanvas(Form("c%zu", i), "", 1600, 1200);
150  c3vs4->Divide(3, 2);
151 
152  //ND Plots (NO Fluctuations)
153  c3vs4->cd(1);
154  TLegend l1(0.4,0.5,0.8,0.79);
155  l1.AddEntry(hND3Flav,"#splitline{ND Fake Data}{#splitline{3 Flavor }{NO Fluctuations}}","l");
156  hND3Flav->Draw("hist");
157  l1.Draw();
158 
159  //FD Plots (No Fluctuations)
160  c3vs4->cd(4);
161  TLegend l2(0.4,0.13,0.8,0.42);
162  l2.AddEntry(hFD3Flav,"#splitline{FD Fake Data}{#splitline{3 Flavor }{NO Fluctuations}}","l");
163  hFD3Flav->Draw("hist");
164  l2.Draw();
165 
166  TDirectory* dir = uDir->mkdir(Form("universe_%zu", i));
167  std::vector<Spectrum> specs = mx->GetUniverse(calc, {nd,fd});
168 
169  TH1D* hsND = specs[0].ToTH1(nd.GetPOT());
170  TH1D* hsFD = specs[1].ToTH1(fd.GetPOT());
171  hsFD->Add(hCosmic);
172 
173  //ND Fluctuation Plots
174  c3vs4->cd(2);
175  TLegend l3(0.4,0.5,0.8,0.79);
176  l3.AddEntry(hsND,"#splitline{ND Fake Data}{#splitline{4 Flavor }{Fluctuations}}","l");
177  hsND->GetYaxis()->SetMaxDigits(4);
178  hsND->SetLineWidth(2);
179  hsND->SetLineColor(kAzure-2);
180  hsND->Draw("hist");
181  hsND->GetXaxis()->CenterTitle();
182  hsND->GetYaxis()->CenterTitle();
183  if (hND4Flav->GetMaximum() > hsND->GetMaximum())
184  hsND->SetMaximum(hND4Flav->GetMaximum());
185  l3.AddEntry(hND4Flav, "#splitline{4 Flavor}{No fluctuations}", "l");
186  hND4Flav->Draw("hist same");
187  l3.Draw();
188 
189  //FD Fluctuation Plots
190  c3vs4->cd(5);
191  TLegend l4(0.4,0.13,0.8,0.42);
192  l4.AddEntry(hsFD,"#splitline{FD Fake Data}{#splitline{4 Flavor }{Fluctuations}}","l");
193  hsFD->SetLineWidth(2);
194  hsFD->SetLineColor(kGreen+2);
195  hsFD->Draw("hist");
196  hsFD->GetXaxis()->CenterTitle();
197  hsFD->GetYaxis()->CenterTitle();
198  if (hFD4Flav->GetMaximum() > hsFD->GetMaximum())
199  hsFD->SetMaximum(hFD4Flav->GetMaximum());
200  l4.AddEntry(hFD4Flav, "#splitline{4 Flavor}{No fluctuations}", "l");
201  hFD4Flav->Draw("hist same");
202  l4.Draw();
203 
204  //ratio Plots
205  c3vs4->cd(3);
206  TH1D* hratioND = (TH1D*)hND3Flav->Clone("hratioND");
207  hratioND->Divide(hsND);
208  hratioND->GetYaxis()->SetRangeUser(0,2);
209  hratioND->GetYaxis()->SetTitle("ND 3 Flavor/4 Flavor");
210  hratioND->SetLineColor(kBlack);
211  hratioND->Draw("hist");
212 
213  c3vs4->cd(6);
214  TH1D* hratioFD = (TH1D*)hFD3Flav->Clone("hratioFD");
215  hratioFD->Divide(hsFD);
216  hratioFD->GetYaxis()->SetRangeUser(0,2.5);
217  hratioFD->GetYaxis()->SetTitle("FD 3 Flavor/4 Flavor");
218  hratioFD->SetLineColor(kBlack);
219  hratioFD->Draw("hist");
220 
221  // saving Fake data with the SetNewParameters() values
222  specs[0].SaveTo(dir->mkdir(kNusFHCNearDet.GetTag().c_str()));
223  specs[1].SaveTo(dir->mkdir(kNusFHCFarDet.GetTag().c_str()));
224 
225  //Directory to store FakeData plots
226  gErrorIgnoreLevel = kError;
227  gSystem->Exec(Form("mkdir -p plots/FakeDataPlots/fs%zu", u));
228  c3vs4->SaveAs(Form("plots/FakeDataPlots/fs%zu/universe_%zu.png", u, i));
229  gErrorIgnoreLevel = kWarning;
230 
231  p->SetProgress(double(i)/double(nUniverses));
232 
233  delete c3vs4;
234  delete hsND;
235  delete hsFD;
236  delete hratioND;
237  delete hratioFD;
238  }
239  p->Done();
240 
241  }
242 
243  delete f;
244 
245  } catch(const std::exception& e) {
246  std::cout << std::endl << "Exception occurred! " << e.what() << std::endl;
247  }
248 
249 }
250 
int nBins
Definition: plotROC.py:16
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
enum BeamMode kRed
std::unique_ptr< covmx::CovarianceMatrix > GetCovarianceMatrix(std::vector< covmx::Sample > samples, std::vector< const ISyst * > systs)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
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:148
void MiniBooNEParameters(osc::OscCalcSterile *calc)
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
const char * p
Definition: xmltok.h:285
void SetPrediction(covmx::Sample &sample, bool systs=true)
Definition: Utilities.h:564
void SetLivetime(double l)
Definition: Sample.h:56
Adapt the PMNS_Sterile calculator to standard interface.
int stats(TString inputFilePath, Int_t firstRun, Int_t lastRun, Float_t thresh, TString myDet)
Definition: stats.C:13
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
std::string GetTag() const
Definition: Sample.cxx:94
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
void FakeSignal(osc::OscCalcSterile *calc, size_t n)
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
void MakeFakeData(bool rhc=false, bool stats=true)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
Sum up livetimes from individual cosmic triggers.
void SetNus20Params(osc::OscCalcSterile *calc)
Definition: Utilities.h:682
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
const covmx::Sample kNusFHCFarDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kFarDet)
OStream cout
Definition: OStream.cxx:6
double GetPOT() const
Definition: Sample.cxx:136
void rootlogon()
Definition: rootlogon.C:7
void SetAngle(int i, int j, double th)
std::vector< const ISyst * > getAllXsecSysts_2018_RPAFix()
const Binning kNCNDBinning
NC custom binning.
Definition: Binning.cxx:104
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
void SetDm(int i, double dm)
TDirectory * dir
Definition: macro.C:5
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
static std::unique_ptr< CovMxManager > LoadFrom(TDirectory *dir)
Standard implementation.
int NBins() const
Definition: Binning.h:29
const Binning kNCFDBinning
Definition: Binning.cxx:105
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::unique_ptr< Spectrum > LoadCosmic(covmx::Sample sample, bool cvmfs=true)
Get cosmics for a given sample.
Definition: Utilities.h:145
A simple ascii-art progress bar.
Definition: Progress.h:9
std::vector< const ISyst * > LoadSysts(covmx::Sample sample)
Get systematics for a given sample.
Definition: Utilities.h:524
void SetPOT(double d)
Definition: Sample.h:55
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
const double kAna2018FHCPOT
Definition: Exposures.h:207
double GetLivetime() const
Definition: Sample.cxx:142
Eigen::MatrixXd Predict(std::vector< Sample > samples, osc::IOscCalcAdjustable *calc, const SystShifts &systs=SystShifts::Nominal())
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
enum BeamMode kGreen
IPrediction * GetPrediction() const
Definition: Sample.cxx:148
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
Definition: Exposures.h:213
std::vector< Spectrum > GetUniverse(osc::OscCalcSterile *calc, std::vector< covmx::Sample > samples)
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...
enum BeamMode string