MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
FitterBase.cpp
Go to the documentation of this file.
1#include "FitterBase.h"
2
4#include "TRandom.h"
5#include "TStopwatch.h"
6#include "TTree.h"
7#include "TGraphAsymmErrors.h"
9
10#pragma GCC diagnostic ignored "-Wuseless-cast"
11
12// *************************
13// Initialise the manager and make it an object of FitterBase class
14// Now we can dump manager settings to the output file
16// *************************
17 //Get mach3 modes from manager
18 random = std::make_unique<TRandom3>(Get<int>(fitMan->raw()["General"]["Seed"], __FILE__, __LINE__));
19
20 // Counter of the accepted # of steps
21 accCount = 0;
22 step = 0;
23 stepStart = 0;
24
25 clock = std::make_unique<TStopwatch>();
26 stepClock = std::make_unique<TStopwatch>();
27 #ifdef DEBUG
28 // Fit summary and debug info
29 debug = GetFromManager<bool>(fitMan->raw()["General"]["Debug"], false, __FILE__ , __LINE__);
30 #endif
31
32 auto outfile = Get<std::string>(fitMan->raw()["General"]["OutputFile"], __FILE__ , __LINE__);
33 // Save output every auto_save steps
34 //you don't want this too often https://root.cern/root/html606/TTree_8cxx_source.html#l01229
35 auto_save = Get<int>(fitMan->raw()["General"]["MCMC"]["AutoSave"], __FILE__ , __LINE__);
36
37 #ifdef MULTITHREAD
38 //KS: TODO This should help with performance when saving entries to ROOT file. I didn't have time to validate hence commented out
39 //Based on other tests it is really helpful
40 //ROOT::EnableImplicitMT();
41 #endif
42 // Set the output file
43 outputFile = M3::Open(outfile, "RECREATE", __FILE__, __LINE__);
44 outputFile->cd();
45 // Set output tree
46 outTree = new TTree("posteriors", "Posterior_Distributions");
47 // Auto-save every 200MB, the bigger the better https://root.cern/root/html606/TTree_8cxx_source.html#l01229
48 outTree->SetAutoSave(-200E6);
49
50 FileSaved = false;
51 SettingsSaved = false;
52 OutputPrepared = false;
53
54 //Create TDirectory
55 CovFolder = outputFile->mkdir("CovarianceFolder");
56 outputFile->cd();
57 SampleFolder = outputFile->mkdir("SampleFolder");
58 outputFile->cd();
59
60 #ifdef DEBUG
61 // Prepare the output log file
62 if (debug) debugFile.open((outfile+".log").c_str());
63 #endif
64
65 TotalNSamples = 0;
66 fTestLikelihood = GetFromManager<bool>(fitMan->raw()["General"]["Fitter"]["FitTestLikelihood"], false, __FILE__ , __LINE__);
67}
68
69// *************************
70// Destructor: close the logger and output file
72// *************************
73 SaveOutput();
74
75 if(outputFile != nullptr) delete outputFile;
76 outputFile = nullptr;
77 MACH3LOG_DEBUG("Closing MaCh3 Fitter Engine");
78}
79
80// *******************
81// Prepare the output tree
83// *******************
84 if(SettingsSaved) return;
85
86 outputFile->cd();
87
88 TDirectory* MaCh3Version = outputFile->mkdir("MaCh3Engine");
89 MaCh3Version->cd();
90
91 if (std::getenv("MaCh3_ROOT") == NULL) {
92 MACH3LOG_ERROR("Need MaCh3_ROOT environment variable");
93 MACH3LOG_ERROR("Please remember about source bin/setup.MaCh3.sh");
94 throw MaCh3Exception(__FILE__ , __LINE__ );
95 }
96
97 if (std::getenv("MACH3") == NULL) {
98 MACH3LOG_ERROR("Need MACH3 environment variable");
99 throw MaCh3Exception(__FILE__ , __LINE__ );
100 }
101
102 std::string header_path = std::string(std::getenv("MACH3"));
103 header_path += "/version.h";
104 FILE* file = fopen(header_path.c_str(), "r");
105 //KS: It is better to use experiment specific header file. If given experiment didn't provide it we gonna use one given by Core MaCh3.
106 if (!file) {
107 header_path = std::string(std::getenv("MaCh3_ROOT"));
108 header_path += "/version.h";
109 } else {
110 fclose(file);
111 }
112
113 // EM: embed the cmake generated version.h file
114 TMacro versionHeader("version_header", "version_header");
115 versionHeader.ReadFile(header_path.c_str());
116 versionHeader.Write();
117
118 TNamed Engine(GetName(), GetName());
119 Engine.Write(GetName().c_str());
120
121 MaCh3Version->Write();
122 delete MaCh3Version;
123
124 outputFile->cd();
125
127
128 MACH3LOG_WARN("\033[0;31mCurrent Total RAM usage is {:.2f} GB\033[0m", MaCh3Utils::getValue("VmRSS") / 1048576.0);
129 MACH3LOG_WARN("\033[0;31mOut of Total available RAM {:.2f} GB\033[0m", MaCh3Utils::getValue("MemTotal") / 1048576.0);
130
131 MACH3LOG_INFO("#####Current Setup#####");
132 MACH3LOG_INFO("Number of covariances: {}", systematics.size());
133 for(unsigned int i = 0; i < systematics.size(); ++i)
134 MACH3LOG_INFO("{}: Cov name: {}, it has {} params", i, systematics[i]->GetName(), systematics[i]->GetNumParams());
135 MACH3LOG_INFO("Number of SampleHandlers: {}", samples.size());
136 for(unsigned int i = 0; i < samples.size(); ++i)
137 MACH3LOG_INFO("{}: SampleHandler name: {}, it has {} samples, {} OscChannels",i , samples[i]->GetTitle(), samples[i]->GetNsamples(), samples[i]->GetNOscChannels());
138
139 //TN: Have to close the folder in order to write it to disk before SaveOutput is called in the destructor
140 CovFolder->Close();
141 SampleFolder->Close();
142
143 SettingsSaved = true;
144}
145
146// *******************
147// Prepare the output tree
149// *******************
150 if(OutputPrepared) return;
151 //MS: Check if we are fitting the test likelihood, rather than T2K likelihood, and only setup T2K output if not
152 if(!fTestLikelihood)
153 {
154 // Check that we have added samples
155 if (!samples.size()) {
156 MACH3LOG_CRITICAL("No samples Found! Is this really what you wanted?");
157 //throw MaCh3Exception(__FILE__ , __LINE__ );
158 }
159
160 // Do we want to save proposal? This will break plotting scripts and is heave for disk space and step time. Only use when debugging
161 bool SaveProposal = GetFromManager<bool>(fitMan->raw()["General"]["SaveProposal"], false, __FILE__ , __LINE__);
162
163 // Prepare the output trees
164 for (ParameterHandlerBase *cov : systematics) {
165 cov->SetBranches(*outTree, SaveProposal);
166 }
167
168 outTree->Branch("LogL", &logLCurr, "LogL/D");
169 outTree->Branch("accProb", &accProb, "accProb/D");
170 outTree->Branch("step", &step, "step/I");
171 outTree->Branch("stepTime", &stepTime, "stepTime/D");
172
173 // Store individual likelihood components
174 // Using double means double as large output file!
175 sample_llh.resize(samples.size());
176 syst_llh.resize(systematics.size());
177
178 for (size_t i = 0; i < samples.size(); ++i) {
179 std::stringstream oss, oss2;
180 oss << "LogL_sample_" << i;
181 oss2 << oss.str() << "/D";
182 outTree->Branch(oss.str().c_str(), &sample_llh[i], oss2.str().c_str());
183 }
184
185 for (size_t i = 0; i < systematics.size(); ++i) {
186 std::stringstream oss, oss2;
187 oss << "LogL_systematic_" << systematics[i]->GetName();
188 oss2 << oss.str() << "/D";
189 outTree->Branch(oss.str().c_str(), &syst_llh[i], oss2.str().c_str());
190 }
191 }
192 else
193 {
194 outTree->Branch("LogL", &logLCurr, "LogL/D");
195 outTree->Branch("accProb", &accProb, "accProb/D");
196 outTree->Branch("step", &step, "step/I");
197 outTree->Branch("stepTime", &stepTime, "stepTime/D");
198 }
199
200 MACH3LOG_INFO("-------------------- Starting MCMC --------------------");
201 #ifdef DEBUG
202 if (debug) {
203 debugFile << "----- Starting MCMC -----" << std::endl;
204 }
205 #endif
206 // Time the progress
207 clock->Start();
208
209 OutputPrepared = true;
210}
211
212// *******************
214// *******************
215 for (size_t i = 0; i < samples.size(); ++i) {
216 samples[i]->CleanMemoryBeforeFit();
217 }
218}
219
220// *******************
222// *******************
223 if(FileSaved) return;
224 //Stop Clock
225 clock->Stop();
226
227 outputFile->cd();
228 outTree->Write();
229
230 MACH3LOG_INFO("{} steps took {:.2f} seconds to complete. ({:.2f}s / step).", step - stepStart, clock->RealTime(), clock->RealTime() / static_cast<double>(step - stepStart));
231 MACH3LOG_INFO("{} steps were accepted.", accCount);
232 #ifdef DEBUG
233 if (debug)
234 {
235 debugFile << "\n\n" << step << " steps took " << clock->RealTime() << " seconds to complete. (" << clock->RealTime() / step << "s / step).\n" << accCount<< " steps were accepted." << std::endl;
236 debugFile.close();
237 }
238 #endif
239
240 outputFile->Close();
241 FileSaved = true;
242}
243
244// *************************
245// Add SampleHandler object to the Markov Chain
247// *************************
248 //Check if the sample has a unique name
249 for (const auto &s : samples) {
250 if (s->GetTitle() == sample->GetTitle()) {
251 MACH3LOG_ERROR("SampleHandler with name '{}' already exists!", sample->GetTitle());
252 throw MaCh3Exception(__FILE__ , __LINE__ );
253 }
254 }
255 // Save additional info from samples
256 SampleFolder->cd();
257
259 TotalNSamples += sample->GetNsamples();
260 MACH3LOG_INFO("Adding {} object, with {} samples", sample->GetTitle(), sample->GetNsamples());
261 samples.push_back(sample);
262 outputFile->cd();
263}
264
265// *************************
266// Add flux systematics, cross-section systematics, ND systematics to the chain
268// *************************
269 MACH3LOG_INFO("Adding systematic object {}, with {} params", cov->GetName(), cov->GetNumParams());
270 systematics.push_back(cov);
271
272 CovFolder->cd();
273 std::vector<double> n_vec(cov->GetNumParams());
274 for (int i = 0; i < cov->GetNumParams(); ++i)
275 n_vec[i] = cov->GetParInit(i);
276
277 cov->GetCovMatrix()->Write(cov->GetName().c_str());
278
279 TH2D* CorrMatrix = cov->GetCorrelationMatrix();
280 CorrMatrix->Write((cov->GetName() + std::string("_Corr")).c_str());
281 delete CorrMatrix;
282
283 // If we have yaml config file for covariance let's save it
284 YAML::Node Config = cov->GetConfig();
285 if(!Config.IsNull())
286 {
287 TMacro ConfigSave = YAMLtoTMacro(Config, (std::string("Config_") + cov->GetName()));
288 ConfigSave.Write();
289 }
290
291 outputFile->cd();
292}
293
294// *******************
295void FitterBase::StartFromPreviousFit(const std::string& FitName) {
296// *******************
297 MACH3LOG_INFO("Getting starting position from {}", FitName);
298
299 TFile *infile = new TFile(FitName.c_str(), "READ");
300 TTree *posts = infile->Get<TTree>("posteriors");
301 int step_val = 0;
302 double log_val = M3::_LARGE_LOGL_;
303 posts->SetBranchAddress("step",&step_val);
304 posts->SetBranchAddress("LogL",&log_val);
305
306 for (size_t s = 0; s < systematics.size(); ++s)
307 {
308 TDirectory* CovarianceFolder = infile->Get<TDirectory>("CovarianceFolder");
309
310 std::string ConfigName = "Config_" + systematics[s]->GetName();
311 TMacro *ConfigCov = CovarianceFolder->Get<TMacro>(ConfigName.c_str());
312 // KS: Not every covariance uses yaml, if it uses yaml make sure they are identical
313 if (ConfigCov != nullptr) {
314 // Config which was in MCMC from which we are starting
315 YAML::Node CovSettings = TMacroToYAML(*ConfigCov);
316 // Config from currently used cov object
317 YAML::Node ConfigCurrent = systematics[s]->GetConfig();
318
319 if (!compareYAMLNodes(CovSettings, ConfigCurrent))
320 {
321 MACH3LOG_ERROR("Yaml configs in previous chain and current one are different", FitName);
322 throw MaCh3Exception(__FILE__ , __LINE__ );
323 }
324
325 delete ConfigCov;
326 }
327
328 CovarianceFolder->Close();
329 delete CovarianceFolder;
330
331 std::vector<double> branch_vals(systematics[s]->GetNumParams(), M3::_BAD_DOUBLE_);
332 for (int i = 0; i < systematics[s]->GetNumParams(); ++i) {
333 posts->SetBranchAddress(systematics[s]->GetParName(i).c_str(), &branch_vals[i]);
334 }
335 posts->GetEntry(posts->GetEntries()-1);
336
337 for (int i = 0; i < systematics[s]->GetNumParams(); ++i) {
338 if(branch_vals[i] == M3::_BAD_DOUBLE_)
339 {
340 MACH3LOG_ERROR("Parameter {} is unvitalised with value {}", i, branch_vals[i]);
341 MACH3LOG_ERROR("Please check more precisely chain you passed {}", FitName);
342 throw MaCh3Exception(__FILE__ , __LINE__ );
343 }
344 }
345 systematics[s]->SetParameters(branch_vals);
346 systematics[s]->AcceptStep();
347
348 MACH3LOG_INFO("Printing new starting values for: {}", systematics[s]->GetName());
349 systematics[s]->PrintNominalCurrProp();
350
351 // Resetting branch adressed to nullptr as we don't want to write into a delected vector out of scope...
352 for (int i = 0; i < systematics[s]->GetNumParams(); ++i) {
353 posts->SetBranchAddress(systematics[s]->GetParName(i).c_str(), nullptr);
354 }
355 }
356 logLCurr = log_val;
357
358 delete posts;
359 infile->Close();
360 delete infile;
361
362 for (size_t s = 0; s < systematics.size(); ++s) {
363 if(systematics[s]->GetDoAdaption()){ //Use separate throw matrix for xsec
364 systematics[s]->SetNumberOfSteps(step_val);
365 }
366 }
367}
368
369// *******************
370// Process the MCMC output to get postfit etc
372// *******************
373 if (fitMan == nullptr) return;
374
375 // Process the MCMC
376 if (GetFromManager<bool>(fitMan->raw()["General"]["ProcessMCMC"], false, __FILE__ , __LINE__)){
377 // Make the processor
378 MCMCProcessor Processor(std::string(outputFile->GetName()));
379
380 Processor.Initialise();
381 // Make the TVectorD pointers which hold the processed output
382 TVectorD *Central = nullptr;
383 TVectorD *Errors = nullptr;
384 TVectorD *Central_Gauss = nullptr;
385 TVectorD *Errors_Gauss = nullptr;
386 TVectorD *Peaks = nullptr;
387
388 // Make the postfit
389 Processor.GetPostfit(Central, Errors, Central_Gauss, Errors_Gauss, Peaks);
390 Processor.DrawPostfit();
391
392 // Make the TMatrix pointers which hold the processed output
393 TMatrixDSym *Covariance = nullptr;
394 TMatrixDSym *Correlation = nullptr;
395
396 // Make the covariance matrix
397 Processor.GetCovariance(Covariance, Correlation);
398 Processor.DrawCovariance();
399
400 std::vector<TString> BranchNames = Processor.GetBranchNames();
401
402 // Re-open the TFile
403 if (!outputFile->IsOpen()) {
404 MACH3LOG_INFO("Opening output again to update with means..");
405 outputFile = new TFile(Get<std::string>(fitMan->raw()["General"]["Output"]["Filename"], __FILE__, __LINE__).c_str(), "UPDATE");
406 }
407 Central->Write("PDF_Means");
408 Errors->Write("PDF_Errors");
409 Central_Gauss->Write("Gauss_Means");
410 Errors_Gauss->Write("Errors_Gauss");
411 Covariance->Write("Covariance");
412 Correlation->Write("Correlation");
413 }
414}
415
416// *************************
417// Run Drag Race
418void FitterBase::DragRace(const int NLaps) {
419// *************************
420 MACH3LOG_INFO("Let the Race Begin!");
421 MACH3LOG_INFO("All tests will be performed with {} threads", M3::GetNThreads());
422
423 // Reweight the MC
424 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs)
425 {
426 TStopwatch clockRace;
427 clockRace.Start();
428 for(int Lap = 0; Lap < NLaps; ++Lap) {
429 samples[ivs]->Reweight();
430 }
431 clockRace.Stop();
432 MACH3LOG_INFO("It took {:.4f} s to reweights {} times sample: {}", clockRace.RealTime(), NLaps, samples[ivs]->GetTitle());
433 MACH3LOG_INFO("On average {:.6f}", clockRace.RealTime()/NLaps);
434 }
435
436 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs)
437 {
438 TStopwatch clockRace;
439 clockRace.Start();
440 for(int Lap = 0; Lap < NLaps; ++Lap) {
441 samples[ivs]->GetLikelihood();
442 }
443 clockRace.Stop();
444 MACH3LOG_INFO("It took {:.4f} s to calculate GetLikelihood {} times sample: {}", clockRace.RealTime(), NLaps, samples[ivs]->GetTitle());
445 MACH3LOG_INFO("On average {:.6f}", clockRace.RealTime()/NLaps);
446 }
447 // Get vector of proposed steps. If we want to run LLH scan or something else after we need to revert changes after proposing steps multiple times
448 std::vector<std::vector<double>> StepsValuesBefore(systematics.size());
449 for (size_t s = 0; s < systematics.size(); ++s) {
450 StepsValuesBefore[s] = systematics[s]->GetProposed();
451 }
452 for (size_t s = 0; s < systematics.size(); ++s) {
453 TStopwatch clockRace;
454 clockRace.Start();
455 for(int Lap = 0; Lap < NLaps; ++Lap) {
456 systematics[s]->ProposeStep();
457 }
458 clockRace.Stop();
459 MACH3LOG_INFO("It took {:.4f} s to propose step {} times cov: {}", clockRace.RealTime(), NLaps, systematics[s]->GetName());
460 MACH3LOG_INFO("On average {:.6f}", clockRace.RealTime()/NLaps);
461 }
462 for (size_t s = 0; s < systematics.size(); ++s) {
463 systematics[s]->SetParameters(StepsValuesBefore[s]);
464 }
465
466 for (size_t s = 0; s < systematics.size(); ++s) {
467 TStopwatch clockRace;
468 clockRace.Start();
469 for(int Lap = 0; Lap < NLaps; ++Lap) {
470 systematics[s]->GetLikelihood();
471 }
472 clockRace.Stop();
473 MACH3LOG_INFO("It took {:.4f} s to calculate get likelihood {} times cov: {}", clockRace.RealTime(), NLaps, systematics[s]->GetName());
474 MACH3LOG_INFO("On average {:.6f}", clockRace.RealTime()/NLaps);
475 }
476 MACH3LOG_INFO("End of race");
477}
478
479// *************************
480bool FitterBase::GetScaneRange(std::map<std::string, std::vector<double>>& scanRanges) {
481// *************************
482 bool isScanRanges = false;
483 // YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D. Barrow
484 if(fitMan->raw()["LLHScan"]["ScanRanges"]){
485 YAML::Node scanRangesList = fitMan->raw()["LLHScan"]["ScanRanges"];
486 for (auto it = scanRangesList.begin(); it != scanRangesList.end(); ++it) {
487 std::string itname = it->first.as<std::string>();
488 std::vector<double> itrange = it->second.as<std::vector<double>>();
489 // Set the mapping as param_name:param_range
490 scanRanges[itname] = itrange;
491 }
492 isScanRanges = true;
493 } else {
494 MACH3LOG_INFO("There are no user-defined parameter ranges, so I'll use default param bounds for LLH Scans");
495 }
496 return isScanRanges;
497}
498
499// *************************
500bool FitterBase::CheckSkipParameter(const std::vector<std::string>& SkipVector, const std::string& ParamName) const {
501// *************************
502 bool skip = false;
503 for(unsigned int is = 0; is < SkipVector.size(); ++is)
504 {
505 if(ParamName.substr(0, SkipVector[is].length()) == SkipVector[is])
506 {
507 skip = true;
508 break;
509 }
510 }
511 return skip;
512}
513
514// *************************
515// Run LLH scan
517// *************************
518 // Save the settings into the output file
519 SaveSettings();
520
521 MACH3LOG_INFO("Starting LLH Scan");
522
523 //KS: Turn it on if you want LLH scan for each ND sample separately, which increase time significantly but can be useful for validating new samples or dials.
524 bool PlotLLHScanBySample = GetFromManager<bool>(fitMan->raw()["LLHScan"]["LLHScanBySample"], false, __FILE__ , __LINE__);;
525 auto SkipVector = GetFromManager<std::vector<std::string>>(fitMan->raw()["LLHScan"]["LLHScanSkipVector"], {}, __FILE__ , __LINE__);;
526
527 // Now finally get onto the LLH scan stuff
528 // Very similar code to MCMC but never start MCMC; just scan over the parameter space
529 std::vector<TDirectory *> Cov_LLH(systematics.size());
530 for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc )
531 {
532 std::string NameTemp = systematics[ivc]->GetName();
533 NameTemp = NameTemp.substr(0, NameTemp.find("_cov")) + "_LLH";
534 Cov_LLH[ivc] = outputFile->mkdir(NameTemp.c_str());
535 }
536
537 std::vector<TDirectory *> SampleClass_LLH(samples.size());
538 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
539 {
540 std::string NameTemp = samples[ivs]->GetTitle();
541 SampleClass_LLH[ivs] = outputFile->mkdir(NameTemp.c_str());
542 }
543
544 TDirectory *Sample_LLH = outputFile->mkdir("Sample_LLH");
545 TDirectory *Total_LLH = outputFile->mkdir("Total_LLH");
546
547 std::vector<TDirectory *>SampleSplit_LLH;
548 if(PlotLLHScanBySample)
549 {
550 SampleSplit_LLH.resize(TotalNSamples);
551 int SampleIterator = 0;
552 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
553 {
554 for(int is = 0; is < samples[ivs]->GetNsamples(); ++is )
555 {
556 SampleSplit_LLH[SampleIterator] = outputFile->mkdir((samples[ivs]->GetSampleName(is)+ "_LLH").c_str());
557 SampleIterator++;
558 }
559 }
560 }
561 // Number of points we do for each LLH scan
562 const int n_points = GetFromManager<int>(fitMan->raw()["LLHScan"]["LLHScanPoints"], 100, __FILE__ , __LINE__);
563 double nSigma = GetFromManager<int>(fitMan->raw()["LLHScan"]["LLHScanSigma"], 1., __FILE__, __LINE__);
564
565 // We print 5 reweights
566 const int countwidth = int(double(n_points)/double(5));
567
568 // YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D. Barrow
569 std::map<std::string, std::vector<double>> scanRanges;
570 const bool isScanRanges = GetScaneRange(scanRanges);
571
572 // Loop over the covariance classes
574 {
575 // Scan over all the parameters
576 // Get the number of parameters
577 int npars = cov->GetNumParams();
578 bool IsPCA = cov->IsPCA();
579 if (IsPCA) npars = cov->GetNParameters();
580 for (int i = 0; i < npars; ++i)
581 {
582 // Get the parameter name
583 std::string name = cov->GetParFancyName(i);
584 if (IsPCA) name += "_PCA";
585 // KS: Check if we want to skip this parameter
586 if(CheckSkipParameter(SkipVector, name)) continue;
587
588 // Get the parameter priors and bounds
589 double prior = cov->GetParInit(i);
590 if (IsPCA) prior = cov->GetPCAHandler()->GetParCurrPCA(i);
591
592 // Get the covariance matrix and do the +/- nSigma
593 // Set lower and upper bounds relative the prior
594 // Set the parameter ranges between which LLH points are scanned
595 double lower = prior - nSigma*cov->GetDiagonalError(i);
596 double upper = prior + nSigma*cov->GetDiagonalError(i);
597 // If PCA, transform these parameter values to the PCA basis
598 if (IsPCA) {
599 lower = prior - nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
600 upper = prior + nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
601 MACH3LOG_INFO("eval {} = {:.2f}", i, cov->GetPCAHandler()->GetEigenValues()(i));
602 MACH3LOG_INFO("prior {} = {:.2f}", i, prior);
603 MACH3LOG_INFO("lower {} = {:.2f}", i, lower);
604 MACH3LOG_INFO("upper {} = {:.2f}", i, upper);
605 MACH3LOG_INFO("nSigma = {:.2f}", nSigma);
606 }
607 // Implementation suggested by D. Barrow
608 // If param ranges are specified in scanRanges node, extract it from there
609 if(isScanRanges){
610 // Find matching entries through std::maps
611 auto it = scanRanges.find(name);
612 if (it != scanRanges.end() && it->second.size() == 2) { //Making sure the range is has only two entries
613 lower = it->second[0];
614 upper = it->second[1];
615 MACH3LOG_INFO("Found matching param name for setting specified range for {}", name);
616 MACH3LOG_INFO("Range for {} = [{:.2f}, {:.2f}]", name, lower, upper);
617 }
618 }
619
620 // Cross-section and flux parameters have boundaries that we scan between, check that these are respected in setting lower and upper variables
621 // This also applies for other parameters like osc, etc.
622 lower = std::max(lower, cov->GetLowerBound(i));
623 upper = std::min(upper, cov->GetUpperBound(i));
624 MACH3LOG_INFO("Scanning {} with {} steps, from [{:.2f} , {:.2f}], prior = {:.2f}", name, n_points, lower, upper, prior);
625
626 // Make the TH1D
627 auto hScan = std::make_unique<TH1D>((name + "_full").c_str(), (name + "_full").c_str(), n_points, lower, upper);
628 hScan->SetTitle((std::string("2LLH_full, ") + name + ";" + name + "; -2(ln L_{sample} + ln L_{xsec+flux} + ln L_{det})").c_str());
629 hScan->SetDirectory(nullptr);
630
631 auto hScanSam = std::make_unique<TH1D>((name + "_sam").c_str(), (name + "_sam").c_str(), n_points, lower, upper);
632 hScanSam->SetTitle((std::string("2LLH_sam, ") + name + ";" + name + "; -2(ln L_{sample})").c_str());
633 hScanSam->SetDirectory(nullptr);
634
635 std::vector<std::unique_ptr<TH1D>> hScanSample(samples.size());
636 std::vector<double> nSamLLH(samples.size());
637 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
638 {
639 std::string NameTemp = samples[ivs]->GetTitle();
640 hScanSample[ivs] = std::make_unique<TH1D>((name+"_"+NameTemp).c_str(), (name+"_" + NameTemp).c_str(), n_points, lower, upper);
641 hScanSample[ivs]->SetDirectory(nullptr);
642 hScanSample[ivs]->SetTitle(("2LLH_" + NameTemp + ", " + name + ";" + name + "; -2(ln L_{" + NameTemp +"})").c_str());
643 nSamLLH[ivs] = 0.;
644 }
645
646 std::vector<std::unique_ptr<TH1D>> hScanCov(systematics.size());
647 std::vector<double> nCovLLH(systematics.size());
648 for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc )
649 {
650 std::string NameTemp = systematics[ivc]->GetName();
651 NameTemp = NameTemp.substr(0, NameTemp.find("_cov"));
652 hScanCov[ivc] = std::make_unique<TH1D>((name + "_" + NameTemp).c_str(), (name + "_" + NameTemp).c_str(), n_points, lower, upper);
653 hScanCov[ivc]->SetDirectory(nullptr);
654 hScanCov[ivc]->SetTitle(("2LLH_" + NameTemp + ", " + name + ";" + name + "; -2(ln L_{" + NameTemp +"})").c_str());
655 nCovLLH[ivc] = 0.;
656 }
657
658 std::vector<TH1D*> hScanSamSplit;
659 std::vector<double> sampleSplitllh;
660 if(PlotLLHScanBySample)
661 {
662 int SampleIterator = 0;
663 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
664 {
665 hScanSamSplit.resize(TotalNSamples);
666 sampleSplitllh.resize(TotalNSamples);
667 for(int is = 0; is < samples[ivs]->GetNsamples(); ++is )
668 {
669 hScanSamSplit[SampleIterator] = new TH1D((name+samples[ivs]->GetSampleName(is)).c_str(), (name+samples[ivs]->GetSampleName(is)).c_str(), n_points, lower, upper);
670 hScanSamSplit[SampleIterator]->SetTitle((std::string("2LLH_sam, ") + name + ";" + name + "; -2(ln L_{sample})").c_str());
671 SampleIterator++;
672 }
673 }
674 }
675
676 // Scan over the parameter space
677 for (int j = 0; j < n_points; ++j)
678 {
679 if (j % countwidth == 0)
680 MaCh3Utils::PrintProgressBar(j, n_points);
681
682 // For PCA we have to do it differently
683 if (IsPCA) {
684 cov->GetPCAHandler()->SetParPropPCA(i, hScan->GetBinCenter(j+1));
685 } else {
686 // Set the parameter
687 cov->SetParProp(i, hScan->GetBinCenter(j+1));
688 }
689
690 // Reweight the MC
691 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
692 {
693 samples[ivs]->Reweight();
694 }
695 //Total LLH
696 double totalllh = 0;
697
698 // Get the -log L likelihoods
699 double samplellh = 0;
700
701 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
702 {
703 nSamLLH[ivs] = samples[ivs]->GetLikelihood();
704 samplellh += nSamLLH[ivs];
705 }
706
707 for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc )
708 {
709 nCovLLH[ivc] = systematics[ivc]->GetLikelihood();
710 totalllh += nCovLLH[ivc];
711 }
712
713 totalllh += samplellh;
714
715 if(PlotLLHScanBySample)
716 {
717 int SampleIterator = 0;
718 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
719 {
720 for(int is = 0; is < samples[ivs]->GetNsamples(); ++is)
721 {
722 sampleSplitllh[SampleIterator] = samples[ivs]->GetSampleLikelihood(is);
723 SampleIterator++;
724 }
725 }
726 }
727
728 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs ) {
729 hScanSample[ivs]->SetBinContent(j+1, 2*nSamLLH[ivs]);
730 }
731 for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc ) {
732 hScanCov[ivc]->SetBinContent(j+1, 2*nCovLLH[ivc]);
733 }
734
735 hScanSam->SetBinContent(j+1, 2*samplellh);
736 hScan->SetBinContent(j+1, 2*totalllh);
737
738 if(PlotLLHScanBySample)
739 {
740 int SampleIterator = 0;
741 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
742 {
743 for(int is = 0; is < samples[ivs]->GetNsamples(); ++is)
744 {
745 hScanSamSplit[SampleIterator]->SetBinContent(j+1, 2*sampleSplitllh[SampleIterator]);
746 SampleIterator++;
747 }
748 }
749 }
750 }
751 for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc )
752 {
753 Cov_LLH[ivc]->cd();
754 hScanCov[ivc]->Write();
755 }
756
757 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
758 {
759 SampleClass_LLH[ivs]->cd();
760 hScanSample[ivs]->Write();
761 }
762 Sample_LLH->cd();
763 hScanSam->Write();
764 Total_LLH->cd();
765 hScan->Write();
766
767 if(PlotLLHScanBySample)
768 {
769 int SampleIterator = 0;
770 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
771 {
772 for(int is = 0; is < samples[ivs]->GetNsamples(); ++is)
773 {
774 SampleSplit_LLH[SampleIterator]->cd();
775 hScanSamSplit[SampleIterator]->Write();
776 delete hScanSamSplit[SampleIterator];
777 SampleIterator++;
778 }
779 }
780 }
781
782 // Reset the parameters to their prior central values
783 if (IsPCA) {
784 cov->GetPCAHandler()->SetParPropPCA(i, prior);
785 } else {
786 cov->SetParProp(i, prior);
787 }
788 }//end loop over systematics
789 }//end loop covariance classes
790
791 for(unsigned int ivc = 0; ivc < systematics.size(); ++ivc )
792 {
793 Cov_LLH[ivc]->Write();
794 delete Cov_LLH[ivc];
795 }
796
797 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
798 {
799 SampleClass_LLH[ivs]->Write();
800 delete SampleClass_LLH[ivs];
801 }
802
803 Sample_LLH->Write();
804 delete Sample_LLH;
805
806 Total_LLH->Write();
807 delete Total_LLH;
808
809 if(PlotLLHScanBySample)
810 {
811 int SampleIterator = 0;
812 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs )
813 {
814 for(int is = 0; is < samples[ivs]->GetNsamples(); ++is )
815 {
816 SampleSplit_LLH[SampleIterator]->Write();
817 delete SampleSplit_LLH[SampleIterator];
818 SampleIterator++;
819 }
820 }
821 }
822}
823
824// *************************
825//LLH scan is good first estimate of step scale
827// *************************
828 TDirectory *Sample_LLH = outputFile->Get<TDirectory>("Sample_LLH");
829 MACH3LOG_INFO("Starting Get Step Scale Based On LLHScan");
830
831 if(Sample_LLH->IsZombie())
832 {
833 MACH3LOG_WARN("Couldn't find Sample_LLH, it looks like LLH scan wasn't run, will do this now");
834 RunLLHScan();
835 }
836
838 {
839 const int npars = cov->GetNumParams();
840 std::vector<double> StepScale(npars);
841 for (int i = 0; i < npars; ++i)
842 {
843 std::string name = cov->GetParFancyName(i);
844
845 StepScale[i] = cov->GetIndivStepScale(i);
846 TH1D* LLHScan = Sample_LLH->Get<TH1D>((name+"_sam").c_str());
847 if(LLHScan == nullptr)
848 {
849 MACH3LOG_WARN("Couldn't find LLH scan, for {}, skipping", name);
850 continue;
851 }
852 const double LLH_val = std::max(LLHScan->GetBinContent(1), LLHScan->GetBinContent(LLHScan->GetNbinsX()));
853 //If there is no sensitivity leave it
854 if(LLH_val < 0.001) continue;
855
856 // EM: assuming that the likelihood is gaussian, approximate sigma value is given by variation/sqrt(-2LLH)
857 // can evaluate this at any point, simple to evaluate it in the first bin of the LLH scan
858 // KS: We assume variation is 1 sigma, each dial has different scale so it becomes faff...
859 const double Var = 1.;
860 const double approxSigma = TMath::Abs(Var)/std::sqrt(LLH_val);
861
862 // Based on Ewan comment I just took the 1sigma width from the LLH, assuming it was Gaussian, but then had to also scale by 2.38/sqrt(N_params)
863 const double NewStepScale = approxSigma * 2.38/std::sqrt(npars);
864 StepScale[i] = NewStepScale;
865 MACH3LOG_DEBUG("Sigma: {}", approxSigma);
866 MACH3LOG_DEBUG("optimal Step Size: {}", NewStepScale);
867 }
868 cov->SetIndivStepScale(StepScale);
869 cov->SaveUpdatedMatrixConfig();
870 }
871}
872
873// *************************
874// Run 2D LLH scan
876// *************************
877 // Save the settings into the output file
878 SaveSettings();
879
880 MACH3LOG_INFO("Starting 2D LLH Scan");
881
882 TDirectory *Sample_2DLLH = outputFile->mkdir("Sample_2DLLH");
883 auto SkipVector = GetFromManager<std::vector<std::string>>(fitMan->raw()["LLHScan"]["LLHScanSkipVector"], {}, __FILE__ , __LINE__);;
884
885 // Number of points we do for each LLH scan
886 const int n_points = GetFromManager<int>(fitMan->raw()["LLHScan"]["2DLLHScanPoints"], 20, __FILE__ , __LINE__);
887 // We print 5 reweights
888 const int countwidth = int(double(n_points)/double(5));
889
890 std::map<std::string, std::vector<double>> scanRanges;
891 const bool isScanRanges = GetScaneRange(scanRanges);
892
893 const double nSigma = GetFromManager<int>(fitMan->raw()["LLHScan"]["LLHScanSigma"], 1., __FILE__, __LINE__);
894
895 // Loop over the covariance classes
897 {
898 // Scan over all the parameters
899 // Get the number of parameters
900 int npars = cov->GetNumParams();
901 bool IsPCA = cov->IsPCA();
902 if (IsPCA) npars = cov->GetNParameters();
903
904 for (int i = 0; i < npars; ++i)
905 {
906 std::string name_x = cov->GetParFancyName(i);
907 if (IsPCA) name_x += "_PCA";
908
909 // Get the parameter priors and bounds
910 double prior_x = cov->GetParInit(i);
911 if (IsPCA) prior_x = cov->GetPCAHandler()->GetParCurrPCA(i);
912
913 // Get the covariance matrix and do the +/- nSigma
914 // Set lower and upper bounds relative the prior
915 double lower_x = prior_x - nSigma*cov->GetDiagonalError(i);
916 double upper_x = prior_x + nSigma*cov->GetDiagonalError(i);
917 // If PCA, transform these parameter values to the PCA basis
918 if (IsPCA) {
919 lower_x = prior_x - nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
920 upper_x = prior_x + nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(i));
921 MACH3LOG_INFO("eval {} = {:.2f}", i, cov->GetPCAHandler()->GetEigenValues()(i));
922 MACH3LOG_INFO("prior {} = {:.2f}", i, prior_x);
923 MACH3LOG_INFO("lower {} = {:.2f}", i, lower_x);
924 MACH3LOG_INFO("upper {} = {:.2f}", i, upper_x);
925 MACH3LOG_INFO("nSigma = {:.2f}", nSigma);
926 }
927 // If param ranges are specified in scanRanges node, extract it from there
928 if(isScanRanges){
929 // Find matching entries through std::maps
930 auto it = scanRanges.find(name_x);
931 if (it != scanRanges.end() && it->second.size() == 2) { //Making sure the range is has only two entries
932 lower_x = it->second[0];
933 upper_x = it->second[1];
934 MACH3LOG_INFO("Found matching param name for setting specified range for {}", name_x);
935 MACH3LOG_INFO("Range for {} = [{:.2f}, {:.2f}]", name_x, lower_x, upper_x);
936 }
937 }
938
939 // Cross-section and flux parameters have boundaries that we scan between, check that these are respected in setting lower and upper variables
940 lower_x = std::max(lower_x, cov->GetLowerBound(i));
941 upper_x = std::min(upper_x, cov->GetUpperBound(i));
942 // KS: Check if we want to skip this parameter
943 if(CheckSkipParameter(SkipVector, name_x)) continue;
944
945 for (int j = 0; j < i; ++j)
946 {
947 std::string name_y = cov->GetParFancyName(j);
948 if (IsPCA) name_y += "_PCA";
949 // KS: Check if we want to skip this parameter
950 if(CheckSkipParameter(SkipVector, name_y)) continue;
951
952 // Get the parameter priors and bounds
953 double prior_y = cov->GetParInit(j);
954 if (IsPCA) prior_y = cov->GetPCAHandler()->GetParCurrPCA(j);
955
956 // Set lower and upper bounds relative the prior
957 double lower_y = prior_y - nSigma*cov->GetDiagonalError(j);
958 double upper_y = prior_y + nSigma*cov->GetDiagonalError(j);
959 // If PCA, transform these parameter values to the PCA basis
960 if (IsPCA) {
961 lower_y = prior_y - nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(j));
962 upper_y = prior_y + nSigma*std::sqrt((cov->GetPCAHandler()->GetEigenValues())(j));
963 MACH3LOG_INFO("eval {} = {:.2f}", i, cov->GetPCAHandler()->GetEigenValues()(j));
964 MACH3LOG_INFO("prior {} = {:.2f}", i, prior_y);
965 MACH3LOG_INFO("lower {} = {:.2f}", i, lower_y);
966 MACH3LOG_INFO("upper {} = {:.2f}", i, upper_y);
967 MACH3LOG_INFO("nSigma = {:.2f}", nSigma);
968 }
969 // If param ranges are specified in scanRanges node, extract it from there
970 if(isScanRanges){
971 // Find matching entries through std::maps
972 auto it = scanRanges.find(name_y);
973 if (it != scanRanges.end() && it->second.size() == 2) { //Making sure the range is has only two entries
974 lower_y = it->second[0];
975 upper_y = it->second[1];
976 MACH3LOG_INFO("Found matching param name for setting specified range for {}", name_y);
977 MACH3LOG_INFO("Range for {} = [{:.2f}, {:.2f}]", name_y, lower_y, upper_y);
978 }
979 }
980
981 // Cross-section and flux parameters have boundaries that we scan between, check that these are respected in setting lower and upper variables
982 lower_y = std::max(lower_y, cov->GetLowerBound(j));
983 upper_y = std::min(upper_y, cov->GetUpperBound(j));
984 MACH3LOG_INFO("Scanning X {} with {} steps, from {:.2f} - {:.2f}, prior = {}", name_x, n_points, lower_x, upper_x, prior_x);
985 MACH3LOG_INFO("Scanning Y {} with {} steps, from {:.2f} - {:.2f}, prior = {}", name_y, n_points, lower_y, upper_y, prior_y);
986
987 auto hScanSam = std::make_unique<TH2D>((name_x + "_" + name_y + "_sam").c_str(), (name_x + "_" + name_y + "_sam").c_str(),
988 n_points, lower_x, upper_x, n_points, lower_y, upper_y);
989 hScanSam->SetDirectory(nullptr);
990 hScanSam->GetXaxis()->SetTitle(name_x.c_str());
991 hScanSam->GetYaxis()->SetTitle(name_y.c_str());
992 hScanSam->GetZaxis()->SetTitle("2LLH_sam");
993
994 // Scan over the parameter space
995 for (int x = 0; x < n_points; ++x)
996 {
997 if (x % countwidth == 0)
998 MaCh3Utils::PrintProgressBar(x, n_points);
999
1000 for (int y = 0; y < n_points; ++y)
1001 {
1002 // For PCA we have to do it differently
1003 if (IsPCA) {
1004 cov->GetPCAHandler()->SetParPropPCA(i, hScanSam->GetXaxis()->GetBinCenter(x+1));
1005 cov->GetPCAHandler()->SetParPropPCA(j, hScanSam->GetYaxis()->GetBinCenter(y+1));
1006 } else {
1007 // Set the parameter
1008 cov->SetParProp(i, hScanSam->GetXaxis()->GetBinCenter(x+1));
1009 cov->SetParProp(j, hScanSam->GetYaxis()->GetBinCenter(y+1));
1010 }
1011 // Reweight the MC
1012 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs) {
1013 samples[ivs]->Reweight();
1014 }
1015
1016 // Get the -log L likelihoods
1017 double samplellh = 0;
1018 for(unsigned int ivs = 0; ivs < samples.size(); ++ivs) {
1019 samplellh += samples[ivs]->GetLikelihood();
1020 }
1021 hScanSam->SetBinContent(x+1, y+1, 2*samplellh);
1022 }// end loop over y points
1023 } // end loop over x points
1024
1025 Sample_2DLLH->cd();
1026 hScanSam->Write();
1027 // Reset the parameters to their prior central values
1028 if (IsPCA) {
1029 cov->GetPCAHandler()->SetParPropPCA(i, prior_x);
1030 cov->GetPCAHandler()->SetParPropPCA(j, prior_y);
1031 } else {
1032 cov->SetParProp(i, prior_x);
1033 cov->SetParProp(j, prior_y);
1034 }
1035 } //end loop over systematics y
1036 }//end loop over systematics X
1037 }//end loop covariance classes
1038 Sample_2DLLH->Write();
1039 delete Sample_2DLLH;
1040}
1041
1042// *************************
1044// *************************
1045 // Save the settings into the output file
1046 SaveSettings();
1047
1048 MACH3LOG_INFO("Starting Sigma Variation");
1049
1050 // Number of variations we want
1051 constexpr int numVar = 5;
1052 //-3 -1 0 +1 +3 sigma variation
1053 constexpr int sigmaArray[numVar] = {-3, -1, 0, 1, 3};
1054
1055 outputFile->cd();
1056
1057 //KS: this is only relevant if PlotByMode is turned on
1058 //Checking each mode is time consuming so we only consider one which are relevant for particular analysis
1059 constexpr int nRelevantModes = 2;
1060
1061 bool DoByMode = GetFromManager<int>(fitMan->raw()["General"]["DoByMode"], false, __FILE__ , __LINE__);
1062
1063 //KS: If true it will make additional plots with LLH sample contribution in each bin, should make it via config file...
1064 bool PlotLLHperBin = false;
1065 auto SkipVector = GetFromManager<std::vector<std::string>>(fitMan->raw()["LLHScan"]["LLHScanSkipVector"], {}, __FILE__ , __LINE__);;
1066
1068 {
1069 TMatrixDSym *Cov = cov->GetCovMatrix();
1070
1071 if(cov->IsPCA())
1072 {
1073 MACH3LOG_ERROR("Using PCAed matrix not implemented within sigma var code, I am sorry :(");
1074 throw MaCh3Exception(__FILE__ , __LINE__ );
1075 }
1076
1077 // Loop over xsec parameters
1078 for (int i = 0; i < cov->GetNumParams(); ++i)
1079 {
1080 // Get the parameter name
1081 std::string name = cov->GetParFancyName(i);
1082 // KS: Check if we want to skip this parameter
1083 if(CheckSkipParameter(SkipVector, name)) continue;
1084
1085 outputFile->cd();
1086 TDirectory* dirArryDial = outputFile->mkdir(name.c_str());
1087 std::vector<TDirectory*> dirArrySample(TotalNSamples);
1088
1089 int SampleIterator = 0;
1090 // Get each sample and how it's responded to our reweighted parameter
1091 for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
1092 {
1093 for(int k = 0; k < samples[ivs]->GetNsamples(); k++ )
1094 {
1095 std::string title = std::string(samples[ivs]->GetPDF(k)->GetName());
1096 dirArryDial->cd();
1097 dirArrySample[SampleIterator] = dirArryDial->mkdir(title.c_str());
1098 SampleIterator++;
1099 }
1100 }
1101
1102 // Get the initial value of ith parameter
1103 double init = cov->GetParInit(i);
1104
1105 std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_x(numVar);
1106 std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_y(numVar);
1107 std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_x_norm(numVar);
1108 std::vector<std::vector<std::unique_ptr<TH1D>>> sigmaArray_y_norm(numVar);
1109
1110 // Set up for single mode
1111 TH1D ****sigmaArray_mode_x = nullptr;
1112 TH1D ****sigmaArray_mode_y = nullptr;
1113 if (DoByMode)
1114 {
1115 sigmaArray_mode_x = new TH1D***[numVar]();
1116 sigmaArray_mode_y = new TH1D***[numVar]();
1117 }
1118
1119 MACH3LOG_INFO("{:<20}{}", name, "");
1120
1121 // Make asymmetric errors just in case
1122 for (int j = 0; j < numVar; ++j)
1123 {
1124 // New value = prior + variation*1sigma uncertainty
1125 double paramVal = cov->GetParInit(i)+sigmaArray[j]*std::sqrt((*Cov)(i,i));
1126
1127 // Check the bounds on the parameter
1128 paramVal = std::max(cov->GetLowerBound(i), std::min(paramVal, cov->GetUpperBound(i)));
1129
1130 // Set the parameter
1131 cov->SetParProp(i, paramVal);
1132 // And reweight the sample
1133 for(unsigned int ivs = 0; ivs < samples.size(); ivs++) {
1134 samples[ivs]->Reweight();
1135 }
1136
1137 sigmaArray_x[j].resize(TotalNSamples);
1138 sigmaArray_y[j].resize(TotalNSamples);
1139 sigmaArray_x_norm[j].resize(TotalNSamples);
1140 sigmaArray_y_norm[j].resize(TotalNSamples);
1141
1142 if (DoByMode)
1143 {
1144 sigmaArray_mode_x[j] = new TH1D**[TotalNSamples]();
1145 sigmaArray_mode_y[j] = new TH1D**[TotalNSamples]();
1146 }
1147
1148 SampleIterator = 0;
1149 // Get each sample and how it's responded to our reweighted parameter
1150 for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
1151 {
1152 for (int k = 0; k < samples[ivs]->GetNsamples(); ++k)
1153 {
1154 // Make a string of the double
1155 std::ostringstream ss;
1156 ss << paramVal;
1157 std::string parVarTitle = name + "_" + ss.str();
1158
1159 auto currSamp = M3::Clone<TH2Poly>(static_cast<TH2Poly*>(samples[ivs]->GetPDF(k)));
1160 // Set a descriptiv-ish title
1161 std::string title_long = std::string(currSamp->GetName())+"_"+parVarTitle;
1162
1163 // Enable the mode histograms AFTER addSelection is called
1164 //Get the 1d binning we want. Let's just use SetupBinning to get this as it already exists
1165 std::vector<double> xbins;
1166 std::vector<double> ybins;
1167 samples[ivs]->SetupBinning(M3::int_t(k), xbins, ybins);
1168
1169 //KS:here we loop over all reaction modes defined in "RelevantModes[nRelevantModes]"
1170 if (DoByMode)
1171 {
1172 //KS: this is only relevant if PlotByMode is turned on
1173 //Checking each mode is time consuming so we only consider one which are relevant for particular analysis
1174 MaCh3Modes_t RelevantModes[nRelevantModes] = {samples[ivs]->GetMaCh3Modes()->GetMode("CCQE"), samples[ivs]->GetMaCh3Modes()->GetMode("2p2h")};
1175
1176 sigmaArray_mode_x[j][SampleIterator] = new TH1D*[nRelevantModes]();
1177 sigmaArray_mode_y[j][SampleIterator] = new TH1D*[nRelevantModes]();
1178 // Now get the TH2D mode variations
1179 std::string mode_title_long;
1180
1181 for(int ir = 0; ir < nRelevantModes; ir++)
1182 {
1183 auto currSampMode = M3::Clone<TH2Poly>(static_cast<TH2Poly*>(samples[ivs]->GetPDFMode(k, RelevantModes[ir])));
1184
1185 mode_title_long = title_long + "_" + samples[ivs]->GetMaCh3Modes()->GetMaCh3ModeName(RelevantModes[ir]);
1186 currSampMode->SetNameTitle(mode_title_long.c_str(), mode_title_long.c_str());
1187 dirArrySample[SampleIterator]->cd();
1188 currSampMode->Write();
1189
1190 sigmaArray_mode_x[j][SampleIterator][ir] = PolyProjectionX(currSampMode.get(), (mode_title_long+"_xProj").c_str(), xbins);
1191 sigmaArray_mode_x[j][SampleIterator][ir]->SetDirectory(nullptr);
1192 sigmaArray_mode_y[j][SampleIterator][ir] = PolyProjectionY(currSampMode.get(), (mode_title_long+"_yProj").c_str(), ybins);
1193 sigmaArray_mode_y[j][SampleIterator][ir]->SetDirectory(nullptr);
1194 }
1195 }
1196
1197 //KS: This will give different results depending if data or Asimov, both have their uses
1198 if (PlotLLHperBin)
1199 {
1200 auto currLLHSamp = M3::Clone<TH2Poly>(static_cast<TH2Poly*>(samples[ivs]->GetPDF(k)));
1201 currLLHSamp->Reset("");
1202 currLLHSamp->Fill(0.0, 0.0, 0.0);
1203
1204 TH2Poly* MCpdf = static_cast<TH2Poly*>(samples[ivs]->GetPDF(k));
1205 TH2Poly* Datapdf = static_cast<TH2Poly*>(samples[ivs]->GetData(k));
1206 TH2Poly* W2pdf = samples[ivs]->GetW2(k);
1207
1208 for(int bin = 1; bin < currLLHSamp->GetNumberOfBins()+1; bin++)
1209 {
1210 const double mc = MCpdf->GetBinContent(bin);
1211 const double dat = Datapdf->GetBinContent(bin);
1212 const double w2 = W2pdf->GetBinContent(bin);
1213 currLLHSamp->SetBinContent(bin, samples[ivs]->GetTestStatLLH(dat, mc, w2));
1214 }
1215 currLLHSamp->SetNameTitle((title_long+"_LLH").c_str() ,(title_long+"_LLH").c_str());
1216 dirArrySample[SampleIterator]->cd();
1217 currLLHSamp->Write();
1218 }
1219
1220 // Project down onto x axis
1221 sigmaArray_x[j][SampleIterator] = std::unique_ptr<TH1D>(PolyProjectionX(currSamp.get(), (title_long+"_xProj").c_str(), xbins));
1222 sigmaArray_x[j][SampleIterator]->SetDirectory(nullptr);
1223 sigmaArray_x[j][SampleIterator]->GetXaxis()->SetTitle(currSamp->GetXaxis()->GetTitle());
1224 sigmaArray_y[j][SampleIterator] = std::unique_ptr<TH1D>(PolyProjectionY(currSamp.get(), (title_long+"_yProj").c_str(), ybins));
1225 sigmaArray_y[j][SampleIterator]->SetDirectory(nullptr);
1226 sigmaArray_y[j][SampleIterator]->GetXaxis()->SetTitle(currSamp->GetYaxis()->GetTitle());
1227
1228 sigmaArray_x_norm[j][SampleIterator] = M3::Clone<TH1D>(sigmaArray_x[j][SampleIterator].get());
1229 sigmaArray_x_norm[j][SampleIterator]->Scale(1., "width");
1230 sigmaArray_y_norm[j][SampleIterator] = M3::Clone<TH1D>(sigmaArray_y[j][SampleIterator].get());
1231 sigmaArray_y_norm[j][SampleIterator]->SetDirectory(nullptr);
1232 sigmaArray_y_norm[j][SampleIterator]->Scale(1., "width");
1233
1234 currSamp->SetNameTitle(title_long.c_str(), title_long.c_str());
1235 dirArrySample[k]->cd();
1236 currSamp->Write();
1237
1238 sigmaArray_x[j][k]->Write();
1239 sigmaArray_y[j][k]->Write();
1240 SampleIterator++;
1241 }//End loop over samples
1242 }
1243 } // End looping over variation
1244
1245 // Restore the parameter to prior value
1246 cov->SetParProp(i, init);
1247
1248 SampleIterator = 0;
1249 // Get each sample and how it's responded to our reweighted parameter
1250 for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
1251 {
1252 for (int k = 0; k < samples[ivs]->GetNsamples(); ++k)
1253 {
1254 std::string title = std::string(samples[ivs]->GetPDF(k)->GetName()) + "_" + name;
1255 auto var_x = MakeAsymGraph(sigmaArray_x[1][SampleIterator].get(), sigmaArray_x[2][SampleIterator].get(), sigmaArray_x[3][SampleIterator].get(), (title+"_X").c_str());
1256 auto var_y = MakeAsymGraph(sigmaArray_y[1][SampleIterator].get(), sigmaArray_y[2][SampleIterator].get(), sigmaArray_y[3][SampleIterator].get(), (title+"_Y").c_str());
1257
1258 auto var_x_norm = MakeAsymGraph(sigmaArray_x_norm[1][SampleIterator].get(), sigmaArray_x_norm[2][SampleIterator].get(), sigmaArray_x_norm[3][SampleIterator].get(), (title+"_X_norm").c_str());
1259 auto var_y_norm = MakeAsymGraph(sigmaArray_y_norm[1][SampleIterator].get(), sigmaArray_y_norm[2][SampleIterator].get(), sigmaArray_y_norm[3][SampleIterator].get(), (title+"_Y_norm").c_str());
1260
1261 dirArrySample[SampleIterator]->cd();
1262 var_x->Write();
1263 var_y->Write();
1264 var_x_norm->Write();
1265 var_y_norm->Write();
1266
1267 //KS: here we loop over all reaction modes defined in "RelevantModes[nRelevantModes]"
1268 if (DoByMode)
1269 {
1270 //KS: this is only relevant if PlotByMode is turned on
1271 //Checking each mode is time consuming so we only consider one which are relevant for particular analysis
1272 MaCh3Modes_t RelevantModes[nRelevantModes] = {samples[ivs]->GetMaCh3Modes()->GetMode("CCQE"), samples[ivs]->GetMaCh3Modes()->GetMode("2p2h")};
1273
1274 for(int ir = 0; ir < nRelevantModes;ir++)
1275 {
1276 auto var_mode_x = MakeAsymGraph(sigmaArray_mode_x[1][SampleIterator][ir], sigmaArray_mode_x[2][SampleIterator][ir], sigmaArray_mode_x[3][SampleIterator][ir], (title+"_"+samples[ivs]->GetMaCh3Modes()->GetMaCh3ModeName(RelevantModes[ir])+"_X").c_str());
1277 auto var_mode_y = MakeAsymGraph(sigmaArray_mode_y[1][SampleIterator][ir], sigmaArray_mode_y[2][SampleIterator][ir], sigmaArray_mode_y[3][SampleIterator][ir], (title+"_"+samples[ivs]->GetMaCh3Modes()->GetMaCh3ModeName(RelevantModes[ir])+"_Y").c_str());
1278
1279 dirArrySample[SampleIterator]->cd();
1280 var_mode_x->Write();
1281 var_mode_y->Write();
1282 } // end for nRelevantModes
1283 } // end if mode
1284
1285 SampleIterator++;
1286 }//End loop over samples(k)
1287 }//end looping over sample object
1288 SampleIterator = 0;
1289 for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
1290 {
1291 for (int k = 0; k < samples[ivs]->GetNsamples(); ++k)
1292 {
1293 dirArrySample[SampleIterator]->Close();
1294 delete dirArrySample[SampleIterator];
1295 SampleIterator++;
1296 }
1297 }
1298
1299 dirArryDial->Close();
1300 delete dirArryDial;
1301
1302 if (DoByMode)
1303 {
1304 for (int j = 0; j < numVar; ++j)
1305 {
1306 SampleIterator = 0;
1307 for(unsigned int ivs = 0; ivs < samples.size(); ivs++ )
1308 {
1309 for (int k = 0; k < samples[ivs]->GetNsamples(); ++k)
1310 {
1311 for(int ir = 0; ir < nRelevantModes;ir++)
1312 {
1313 delete sigmaArray_mode_x[j][SampleIterator][ir];
1314 delete sigmaArray_mode_y[j][SampleIterator][ir];
1315 }// end for nRelevantModes
1316 delete[] sigmaArray_mode_x[j][SampleIterator];
1317 delete[] sigmaArray_mode_y[j][SampleIterator];
1318 SampleIterator++;
1319 }// end for samples
1320 }
1321 }
1322 delete[] sigmaArray_mode_x;
1323 delete[] sigmaArray_mode_y;
1324 }
1325 } // end looping over xsec parameters (i)
1326 } // end looping over covarianceBase objects
1327}
#define _MaCh3_Safe_Include_Start_
KS: Avoiding warning checking for headers.
Definition: Core.h:106
#define _MaCh3_Safe_Include_End_
KS: Restore warning checking after including external headers.
Definition: Core.h:117
std::vector< TString > BranchNames
MaCh3Plotting::PlottingManager * man
TH1D * PolyProjectionX(TObject *poly, std::string TempName, const std::vector< double > &xbins, const bool computeErrors)
WP: Poly Projectors.
std::unique_ptr< TGraphAsymmErrors > MakeAsymGraph(TH1 *sigmaArrayLeft, TH1 *sigmaArrayCentr, TH1 *sigmaArrayRight, const std::string &title)
Used by sigma variation, check how 1 sigma changes spectra.
TH1D * PolyProjectionY(TObject *poly, std::string TempName, const std::vector< double > &ybins, const bool computeErrors)
WP: Poly Projectors.
#define MACH3LOG_CRITICAL
Definition: MaCh3Logger.h:26
#define MACH3LOG_DEBUG
Definition: MaCh3Logger.h:22
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:24
int MaCh3Modes_t
Enumerator of MaCh3Mode.
Definition: MaCh3Modes.h:14
TMacro YAMLtoTMacro(const YAML::Node &yaml_node, const std::string &name)
Convert a YAML node to a ROOT TMacro object.
Definition: YamlHelper.h:161
YAML::Node TMacroToYAML(const TMacro &macro)
KS: Convert a ROOT TMacro object to a YAML node.
Definition: YamlHelper.h:146
bool compareYAMLNodes(const YAML::Node &node1, const YAML::Node &node2)
Compare if yaml nodes are identical.
Definition: YamlHelper.h:179
void RunLLHScan()
Perform a 1D likelihood scan.
Definition: FitterBase.cpp:516
void AddSystObj(ParameterHandlerBase *cov)
This function adds a Covariance object to the analysis framework. The Covariance object will be utili...
Definition: FitterBase.cpp:267
void AddSampleHandler(SampleHandlerBase *sample)
This function adds a sample PDF object to the analysis framework. The sample PDF object will be utili...
Definition: FitterBase.cpp:246
std::unique_ptr< TRandom3 > random
Random number.
Definition: FitterBase.h:128
std::vector< SampleHandlerBase * > samples
Sample holder.
Definition: FitterBase.h:113
int stepStart
step start if restarting
Definition: FitterBase.h:105
bool CheckSkipParameter(const std::vector< std::string > &SkipVector, const std::string &ParamName) const
KS: Check whether we want to skip parameter using skip vector.
Definition: FitterBase.cpp:500
void ProcessMCMC()
Process MCMC output.
Definition: FitterBase.cpp:371
int accCount
counts accepted steps
Definition: FitterBase.h:103
virtual std::string GetName() const
Get name of class.
Definition: FitterBase.h:66
bool OutputPrepared
Checks if output prepared not repeat some operations.
Definition: FitterBase.h:149
void SaveOutput()
Save output and close files.
Definition: FitterBase.cpp:221
TFile * outputFile
Output.
Definition: FitterBase.h:131
void SaveSettings()
Save the settings that the MCMC was run with.
Definition: FitterBase.cpp:82
manager * fitMan
The manager.
Definition: FitterBase.h:92
unsigned int step
current state
Definition: FitterBase.h:95
void PrepareOutput()
Prepare the output file.
Definition: FitterBase.cpp:148
bool SettingsSaved
Checks if setting saved not repeat some operations.
Definition: FitterBase.h:147
double accProb
current acceptance prob
Definition: FitterBase.h:101
void GetStepScaleBasedOnLLHScan()
LLH scan is good first estimate of step scale.
Definition: FitterBase.cpp:826
virtual void StartFromPreviousFit(const std::string &FitName)
Allow to start from previous fit/chain.
Definition: FitterBase.cpp:295
bool FileSaved
Checks if file saved not repeat some operations.
Definition: FitterBase.h:145
std::vector< double > sample_llh
store the llh breakdowns
Definition: FitterBase.h:108
void RunSigmaVar()
Perform a 2D and 1D sigma var for all samples.
double stepTime
Time of single step.
Definition: FitterBase.h:125
std::unique_ptr< TStopwatch > clock
tells global time how long fit took
Definition: FitterBase.h:121
std::unique_ptr< TStopwatch > stepClock
tells how long single step/fit iteration took
Definition: FitterBase.h:123
TDirectory * CovFolder
Output cov folder.
Definition: FitterBase.h:133
TDirectory * SampleFolder
Output sample folder.
Definition: FitterBase.h:135
void DragRace(const int NLaps=100)
Calculates the required time for each sample or covariance object in a drag race simulation....
Definition: FitterBase.cpp:418
void Run2DLLHScan()
Perform a 2D likelihood scan.
Definition: FitterBase.cpp:875
unsigned int TotalNSamples
Total number of samples used.
Definition: FitterBase.h:115
double logLCurr
current likelihood
Definition: FitterBase.h:97
std::vector< double > syst_llh
systematic llh breakdowns
Definition: FitterBase.h:110
int auto_save
auto save every N steps
Definition: FitterBase.h:139
FitterBase(manager *const fitMan)
Constructor.
Definition: FitterBase.cpp:15
bool fTestLikelihood
Necessary for some fitting algorithms like PSO.
Definition: FitterBase.h:142
virtual ~FitterBase()
Destructor for the FitterBase class.
Definition: FitterBase.cpp:71
TTree * outTree
Output tree with posteriors.
Definition: FitterBase.h:137
void SanitiseInputs()
Remove obsolete memory and make other checks before fit starts.
Definition: FitterBase.cpp:213
bool GetScaneRange(std::map< std::string, std::vector< double > > &scanRanges)
YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D....
Definition: FitterBase.cpp:480
std::vector< ParameterHandlerBase * > systematics
Systematic holder.
Definition: FitterBase.h:118
Class responsible for processing MCMC chains, performing diagnostics, generating plots,...
Definition: MCMCProcessor.h:65
void Initialise()
Scan chain, what parameters we have and load information from covariance matrices.
const std::vector< TString > & GetBranchNames() const
Get the vector of branch names from root file.
void GetPostfit(TVectorD *&Central, TVectorD *&Errors, TVectorD *&Central_Gauss, TVectorD *&Errors_Gauss, TVectorD *&Peaks)
Get the post-fit results (arithmetic and Gaussian)
void DrawCovariance()
Draw the post-fit covariances.
void DrawPostfit()
Draw the post-fit comparisons.
void GetCovariance(TMatrixDSym *&Cov, TMatrixDSym *&Corr)
Get the post-fit covariances and correlations.
Custom exception class for MaCh3 errors.
Base class responsible for handling of systematic error parameters. Capable of using PCA or using ada...
Class responsible for handling implementation of samples used in analysis, reweighting and returning ...
virtual void SaveAdditionalInfo(TDirectory *Dir)
Store additional info in a chan.
The manager class is responsible for managing configurations and settings.
Definition: Manager.h:16
void SaveSettings(TFile *const OutputFile)
Add manager useful information's to TFile, in most cases to Fitter.
Definition: Manager.cpp:49
YAML::Node const & raw()
Return config.
Definition: Manager.h:41
int GetNumParams() const
Get total number of parameters.
TH2D * GetCorrelationMatrix()
KS: Convert covariance matrix to correlation matrix and return TH2D which can be used for fancy plott...
TMatrixDSym * GetCovMatrix() const
Return covariance matrix.
std::string GetName() const
Get name of covariance.
double GetParInit(const int i) const
Get prior parameter value.
YAML::Node GetConfig() const
Getter to return a copy of the YAML node.
virtual std::string GetTitle() const
virtual M3::int_t GetNsamples()
static constexpr const double _BAD_DOUBLE_
Default value used for double initialisation.
Definition: Core.h:43
TFile * Open(const std::string &Name, const std::string &Type, const std::string &File, const int Line)
Opens a ROOT file with the given name and mode.
static constexpr const double _LARGE_LOGL_
Large Likelihood is used it parameter go out of physical boundary, this indicates in MCMC that such s...
Definition: Core.h:70
int GetNThreads()
number of threads which we need for example for TRandom3
Definition: Monitor.cpp:322
int int_t
Definition: Core.h:29
void PrintProgressBar(const Long64_t Done, const Long64_t All)
KS: Simply print progress bar.
Definition: Monitor.cpp:212
int getValue(const std::string &Type)
CW: Get info like RAM.
Definition: Monitor.cpp:235