fit_mcmc_ND_systs.C
Go to the documentation of this file.
1 /*
2  * fit_mcmc_ND_systs_AllNumu.C:
3  * Perform fit to ND fake data with random x-sec systematic
4  * pulls applied via ThrowNDFDFakeData.C from previously generated predictions
5  * of "AllNumu" events.
6  *
7  * M. Dolce, michael.dolce@tufts.edu
8  * October, 2020.
9  */
10 
11 
12 #include "TMath.h"
13 
14 #include "CAFAna/Core/SystShifts.h"
15 #include "CAFAna/Fit/StanConfig.h"
16 #include "CAFAna/Fit/StanFitter.h"
19 
21 
22 #include "OscLib/OscCalcAnalytic.h"
23 #include "OscLib/OscCalcDMP.h"
24 
27 
28 
29 namespace xsec_systs
30 {
31  const std::vector<std::string> FSI_SYSTS
32  {
33  "hNFSI_MFP_2020",
34  "hNFSI_FateFracEV1_2020",
35  "hNFSI_FateFracEV2_2020",
36  "hNFSI_FateFracEV3_2020",
37  };
38 
39  const std::vector<std::string> QE_SYSTS
40  {
41  "ZNormCCQE",
42  "ZExpAxialFFSyst2020_EV1",
43  "ZExpAxialFFSyst2020_EV2",
44  "ZExpAxialFFSyst2020_EV3",
45  "ZExpAxialFFSyst2020_EV4",
46  };
47 
48  const std::vector<std::string> RES_SYSTS
49  {
50  "LowQ2RESSupp2020",
51  "MaCCRES",
52  "MvNCRES",
53  "MvCCRES",
54  "RPAShapesupp2020",
55  "RPAShapeenh2020",
56  };
57 
58  const std::vector<std::string> MEC_SYSTS
59  {
60  "MECDoubleGaussEnhSystSigmaQ0_1",
61  "MECDoubleGaussEnhSystCorr_1",
62  "MECDoubleGaussEnhSystMeanQ0_1",
63  "MECDoubleGaussEnhSystMeanQ3_2",
64  "MECDoubleGaussEnhSystSigmaQ3_2",
65  "MECDoubleGaussEnhSystNorm_1",
66  "MECDoubleGaussEnhSystNorm_2",
67  "MECDoubleGaussEnhSystMeanQ3_1",
68  "MECDoubleGaussEnhSystCorr_2",
69  "MECDoubleGaussEnhSystMeanQ0_2",
70  "MECDoubleGaussEnhSystSigmaQ0_2",
71  "MECDoubleGaussEnhSystSigmaQ3_1",
72  "MECDoubleGaussEnhSystBaseline",
73  };
74 
75  const std::vector<std::string> DIS_SYSTS
76  {
77  "FormZone2020",
78  "DISvpCC0pi_2020",
79  "DISvpCC1pi_2020",
80  "DISvpCC2pi_2020",
81  "DISvpCC3pi_2020",
82  "DISvpNC0pi_2020",
83  "DISvpNC1pi_2020",
84  "DISvpNC2pi_2020",
85  "DISvpNC3pi_2020",
86  "DISvnCC0pi_2020",
87  "DISvnCC1pi_2020",
88  "DISvnCC2pi_2020",
89  "DISvnCC3pi_2020",
90  "DISvnNC0pi_2020",
91  "DISvnNC1pi_2020",
92  "DISvnNC2pi_2020",
93  "DISvnNC3pi_2020",
94  "DISvbarpCC0pi_2020",
95  "DISvbarpCC1pi_2020",
96  "DISvbarpCC2pi_2020",
97  "DISvbarpCC3pi_2020",
98  "DISvbarpNC0pi_2020",
99  "DISvbarpNC1pi_2020",
100  "DISvbarpNC2pi_2020",
101  "DISvbarpNC3pi_2020",
102  "DISvbarnCC0pi_2020",
103  "DISvbarnCC1pi_2020",
104  "DISvbarnCC2pi_2020",
105  "DISvbarnCC3pi_2020",
106  "DISvbarnNC0pi_2020",
107  "DISvbarnNC1pi_2020",
108  "DISvbarnNC2pi_2020",
109  "DISvbarnNC3pi_2020",
110  };
111 }
112 
113 namespace test
114 {
115  // some 2020 test systematics from getAllXsecSysts_2020()
116  const std::vector<std::string> SYSTS_TO_THROW
117  {
118  "MvCCRES"
119  // "LowQ2RESSupp2020",
120  // "FormZone2020", // not ideal for large pulls
121  // "hNFSI_MFP_2020",
122  // "ZExpAxialFFSyst2020_EV1", unknown issues with ZExp syst currently
123  };
124 
125  void YouAreHere(const int & location)
126  {
127  std::cout << "You are at step number: " << location << "." << std::endl;
128  }
129 }
130 
131 // =====================================================================================================
132 void fit_mcmc_ND_systs(const std::string& fakeDataFile, // "$ana/mcmc/fake-data/<file.root>"
133  const std::string& fakeDataSpectraDir, // often "spectra"
134  const std::string& intType, // Int Type of systs to fit: QE, MEC, RES, DIS, TEST, ALL
135  bool allTopologies=false, // which topological prediction? AllNumu (false)
136  bool fitSysts=true, // true, almost always
137  bool fullWarmups=false, // true --> use full 500 warmup
138  bool gridSubmission=false, // use grid? true.
139  const std::string& dirPrefix=".", // often "$ana/mcmc/samples", "." for grid
140  const std::string& samplesFilename="ND_samples_3flavor.root")
141 // =====================================================================================================
142 {
143 
144  if (intType == "QE" || intType == "MEC" || intType == "DIS" ||
145  intType == "FSI" || intType == "RES" || intType == "TEST" || intType == "ALL")
146  { std::cout<< "Throwing Interaction Type: " << intType << "........" << std::endl;}
147  else
148  { std::cerr << "Abort. No Interaction type specified: QE, MEC, RES, DIS, FSI, TEST, ALL" << std::endl;
149  exit(1); }
150 
151  std::cout << "Saving MCMC samples to file: " << mcmc::FullFilename(dirPrefix, samplesFilename) << std::endl;
152 
153  std::unique_ptr<osc::IOscCalcAdjustable> calc;
154  std::unique_ptr<Spectrum> fakeData;
155 
156 // build my fit variables...
157 // no oscillation fit variables...this is a ND fit, only
158 
159  calc = std::make_unique<osc::OscCalcDMP>();
160  mcmc::ResetCalculator(*calc);
161  // auto predict = IPrediction::Predict(calc);
162 
163  std::vector<std::unique_ptr<ana::ISyst>> kDummyMECSysts = mcmc::LoadSystsCreateDummyMECSysts();
164 
165 
166  // try the fake data --- LoadFakeData() appears it can work for ND data too...
167  auto spectra = mcmc::LoadFakeData(fakeDataFile, fakeDataSpectraDir);
168 
169  // load the predictions and build a MultiExperiment out of them
170  std::vector<std::string> predNames {"pred_interp_EvElse",
171  "pred_interp_MuPiEtc",
172  "pred_interp_MuPr",
173  "pred_interp_MuPrEtc",
174  "pred_interp_Muon"};
175  if (!allTopologies)
176  {predNames.clear(); predNames.emplace_back("pred_interp_AllNumu");}
177 
178  std::string localPathStub = "/nova/ana/users/mdolce/mcmc/predictions";
179  std::string gridPathStub = "${CONDOR_DIR_INPUT}";
180  std::string finalPathStub;
181  (gridSubmission) ? finalPathStub = std::move(gridPathStub) : finalPathStub = std::move(localPathStub);
182 
183  auto preds = mcmc::LoadNDPreds({ { "fhc_nd",
184  Form("%s/predictions_systs_fit_fhc_doublegaussB_2020.root", finalPathStub.c_str() )} },
185  predNames);
186 // For screen sessions (no grid): "/nova/ana/users/mdolce/mcmc/predictions/predictions_systs_fit_fhc_doublegaussB_2020.root"
187 
188  auto expts = mcmc::BuildExperiments(spectra, preds);
190 
191  std::set<const ana::ISyst*> systs;
192  std::vector<const ana::ISyst*> intTypeSystVector;
193 
194  // use the newer Exponential class.
195  auto shifts = std::make_unique<ExponentialPriorSystShifts>();
196 
197  //add the systematics to fit and initially set them to 0.
198  if (fitSysts)
199  {
200  if (intType == "QE") {
201  for (const auto &systName : xsec_systs::QE_SYSTS) {
202  intTypeSystVector.emplace_back(Registry<ISyst>::ShortNameToPtr(systName));
203  shifts->SetShift(Registry<ISyst>::ShortNameToPtr(systName), 0.);
204  }
205  }
206 
207  else if (intType == "MEC"){
208  for (const auto &systName : xsec_systs::MEC_SYSTS) {
209  intTypeSystVector.emplace_back(Registry<ISyst>::ShortNameToPtr(systName));
210  shifts->SetShift(Registry<ISyst>::ShortNameToPtr(systName), 0.);
211  }
212  }
213 
214  else if (intType == "RES") {
215  for (const auto &systName : xsec_systs::RES_SYSTS) {
216  intTypeSystVector.emplace_back(Registry<ISyst>::ShortNameToPtr(systName));
217  shifts->SetShift(Registry<ISyst>::ShortNameToPtr(systName), 0.);
218  }
219  }
220 
221  else if (intType == "DIS") {
222  for (const auto &systName : xsec_systs::DIS_SYSTS) {
223  intTypeSystVector.emplace_back(Registry<ISyst>::ShortNameToPtr(systName));
224  shifts->SetShift(Registry<ISyst>::ShortNameToPtr(systName), 0.);
225  }
226  }
227 
228  else if (intType == "FSI") {
229  for (const auto &systName : xsec_systs::FSI_SYSTS) {
230  intTypeSystVector.emplace_back(Registry<ISyst>::ShortNameToPtr(systName));
231  shifts->SetShift(Registry<ISyst>::ShortNameToPtr(systName), 0.);
232  }
233  }
234 
235  else if (intType == "TEST") {
236  for (const auto &systName : test::SYSTS_TO_THROW) {
237  intTypeSystVector.emplace_back(Registry<ISyst>::ShortNameToPtr(systName));
238  shifts->SetShift(Registry<ISyst>::ShortNameToPtr(systName), 0.);
239  }
240  }
241 
242  else if (intType == "ALL") {systs = mcmc::SupportedSysts(expt.get());}
243  }
244 
245  // if doing "ALL" -- must be set outside of the loop
246  std::vector<const ana::ISyst*> allSystsVector(systs.begin(), systs.end());
247  for (const auto &syst : allSystsVector) {
248  const std::string& sname = syst->ShortName();
249  shifts->SetShift(Registry<ISyst>::ShortNameToPtr(sname), 0.);
250  }
251 
252 
253  // be clear about which vector of systematics is about to be fit...
254  std::vector<const ana::ISyst*> fitSystVector;
255  if (intType == "ALL") { fitSystVector = allSystsVector;}
256  else { fitSystVector = intTypeSystVector;}
257 
258 
259  StanConfig cfg;
260  cfg.num_warmup = 500;
261  cfg.num_samples = 1000;
262 
263  // from default of 0.8. Try to reduce the step size to minimize the divergences at boundary
264  //cfg.delta = 0.9;
265 
266  if (!fullWarmups) { cfg.num_samples = 1; cfg.num_warmup = 1; }
267  cfg.max_depth = 15;
268  cfg.verbosity = StanConfig::Verbosity::kQuiet; // kEverything;
269  // ND fit, only. No oscillation parameters to be fit via MCMC.
270  StanFitter fitter(expt.get(), {}, fitSystVector);
271  fitter.SetStanConfig(cfg);
272  fitter.Fit(calc.get(), *shifts);
273 // Test the derivatives of the fitting to ensure they are sensible
274 // ( "model" column of output is 0 ---> gradients are broken
275 // fitter.TestGradients(calc.get(), *shifts);
276 
277  std::unique_ptr<ana::MCMCSamples> warmup = fitter.ExtractSamples(ana::MemoryTupleWriter::WhichSamples::kWarmup);
278  std::unique_ptr<ana::MCMCSamples> samples = fitter.ExtractSamples(ana::MemoryTupleWriter::WhichSamples::kPostWarmup);
279 
280  mcmc::SaveToFile(mcmc::FullFilename(dirPrefix, samplesFilename),
281  fakeData.get(),
282  warmup.get(),
283  samples.get());
284 
285 
286 }
const std::vector< std::string > SYSTS_TO_THROW
std::string FullFilename(const std::string &dir, std::string file)
Definition: MCMC3FShared.h:53
#define location
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
std::unique_ptr< ana::IExperiment > BuildMultiExperiment(const std::vector< const ana::IExperiment * > &expts)
const std::vector< std::string > RES_SYSTS
void SaveToFile(const std::string &fileName, const ana::Spectrum *fakeData, const ana::MCMCSamples *warmup, const ana::MCMCSamples *samples, const ana::SystShifts *trueSystPulls, const osc::IOscCalcAdjustable *trueOscParams)
Configuration parameters for the Stan MCMC fitter, bundled up together for easy storage and parameter...
Definition: StanConfig.h:6
OStream cerr
Definition: OStream.cxx:7
std::vector< const ana::IExperiment * > ExptPtrs(const std::vector< std::unique_ptr< ana::IExperiment >> &exptVector)
Definition: MCMC3FShared.h:79
Verbosity verbosity
How verbose do you want me to be?
Definition: StanConfig.h:89
std::vector< std::unique_ptr< ana::ISyst > > LoadSystsCreateDummyMECSysts()
osc::OscCalcDumb calc
const std::vector< std::string > MEC_SYSTS
const std::vector< std::string > DIS_SYSTS
void SetStanConfig(const StanConfig &cfg)
Change the config used for Stan. See the StanConfig struct documentation for ideas.
Definition: StanFitter.h:224
const std::vector< std::string > FSI_SYSTS
void YouAreHere(const int &location)
expt
Definition: demo5.py:34
std::set< const ana::ISyst * > SupportedSysts(const ana::IExperiment *expt)
Which systs does this experiment support?
std::vector< std::unique_ptr< ana::IExperiment > > BuildExperiments(const std::map< std::string, ana::Spectrum > &fakeData, const std::vector< ana::predictions > &preds)
const std::vector< std::string > QE_SYSTS
int num_warmup
Number of initial steps in the Markov chain (used to enter the typical set). These are usually discar...
Definition: StanConfig.h:63
OStream cout
Definition: OStream.cxx:6
int max_depth
Depth of the binary tree used by the HMC when taking its leapfrog steps. Note: the number of leapfrog...
Definition: StanConfig.h:70
exit(0)
std::map< std::string, ana::Spectrum > LoadFakeData(const std::string &fakeDataFilename, const std::string &path)
void ResetCalculator(osc::IOscCalcAdjustable &calc)
std::vector< ana::predictions > LoadNDPreds(const std::vector< std::pair< std::string, std::string >> &filenames, const std::vector< std::string > &predObjNames)
Fitter type that bolts the Stan fitting tools onto CAFAna.
Definition: StanFitter.h:126
int num_samples
Number of steps in the Markov chain retained for analysis (after warmup).
Definition: StanConfig.h:64
void fit_mcmc_ND_systs(const std::string &fakeDataFile, const std::string &fakeDataSpectraDir, const std::string &intType, bool allTopologies=false, bool fitSysts=true, bool fullWarmups=false, bool gridSubmission=false, const std::string &dirPrefix=".", const std::string &samplesFilename="ND_samples_3flavor.root")
enum BeamMode string