MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
LikelihoodFit.cpp
Go to the documentation of this file.
1#include "LikelihoodFit.h"
2
3// *******************
4// Run the Markov chain with all the systematic objects added
6// *******************
7 NPars = 0;
8 NParsPCA = 0;
9 fMirroring = GetFromManager<bool>(fitMan->raw()["General"]["Fitter"]["Mirroring"], false);
10 if(fMirroring) MACH3LOG_INFO("Mirroring enabled");
11}
12
13
14// *************************
15// Destructor: close the logger and output file
17// *************************
18
19}
20
21// *******************
23// *******************
24 // Save the settings into the output file
26
27 // Prepare the output branches
29
30 for (size_t s = 0; s < systematics.size(); ++s) {
31 NPars += systematics[s]->GetNumParams();
32 NParsPCA += systematics[s]->GetNParameters();
33 }
34
35 //KS: If PCA is note enabled NParsPCA == NPars
36 MACH3LOG_INFO("Total number of parameters {}", NParsPCA);
37}
38
39// *******************
40double LikelihoodFit::CalcChi2(const double* x) {
41// *******************
42 if (step % 10000 == 0) {
43 MACH3LOG_INFO("Iteration {}", step);
44 }
45
46 stepClock->Start();
47
48 int ParCounter = 0;
49 double llh = 0;
50 for (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
51 {
52 if(!(*it)->IsPCA())
53 {
54 std::vector<double> pars;
55 const int NumPar = (*it)->GetNumParams();
56 //KS: Avoid push back as they are slow
57 pars.resize(NumPar);
58 for(int i = 0; i < NumPar; ++i, ++ParCounter)
59 {
60 double ParVal = x[ParCounter];
61 //KS: Basically apply mirroring for parameters out of bounds
62 if(fMirroring)
63 {
64 if(ParVal < (*it)->GetLowerBound(i))
65 {
66 ParVal = (*it)->GetLowerBound(i) + ((*it)->GetLowerBound(i) - ParVal);
67 }
68 else if (ParVal > (*it)->GetUpperBound(i))
69 {
70 ParVal = (*it)->GetUpperBound(i) - ( ParVal - (*it)->GetUpperBound(i));
71 }
72 }
73 pars[i] = ParVal;
74 }
75 (*it)->SetParameters(pars);
76 }
77 else
78 {
79 std::vector<double> pars;
80 const int NumPar = (*it)->GetNParameters();
81 //KS: Avoid push back as they are slow
82 pars.resize(NumPar);
83 for(int i = 0; i < NumPar; ++i, ++ParCounter)
84 {
85 double ParVal = x[ParCounter];
86 //KS: Basically apply mirroring for parameters out of bounds
87 pars[i] = ParVal;
88 }
89 (*it)->GetPCAHandler()->SetParametersPCA(pars);
90 }
91 (*it)->AcceptStep();
92 }
93
94 // Loop over the systematics and propose the initial step
95 int stdIt = 0;
96 for (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it, ++stdIt)
97 {
98 //GetLikelihood will return LargeL if out of bounds, for minimizers this is not the problem, while calcLikelihood will return actual likelihood
99 syst_llh[stdIt] = (*it)->CalcLikelihood();
100 llh += syst_llh[stdIt];
101 #ifdef DEBUG
102 if (debug) debugFile << "LLH after " << systematics[stdIt]->GetName() << " " << llh << std::endl;
103 #endif
104 }
105 // Could multi-thread this
106 // But since sample reweight is multi-threaded it's probably better to do that
107 for (size_t i = 0; i < samples.size(); i++)
108 {
109 samples[i]->Reweight();
110 }
111
112 //DB for atmospheric event by event sample migration, need to fully reweight all samples to allow event passing prior to likelihood evaluation
113 for (size_t i = 0; i < samples.size(); i++) {
114 // Get the sample likelihoods and add them
115 sample_llh[i] = samples[i]->GetLikelihood();
116 llh += sample_llh[i];
117 #ifdef DEBUG
118 if (debug) debugFile << "LLH after sample " << i << " " << llh << std::endl;
119 #endif
120 }
121
122 // Save the proposed likelihood (class member)
123 logLProp = llh;
124 logLCurr = llh;
125 accProb = 1;
126
127 stepClock->Stop();
128 stepTime = stepClock->RealTime();
129
130 // Write step to output tree
131 outTree->Fill();
132
133 // Auto save the output
134 if (step % auto_save == 0) outTree->AutoSave();
135 step++;
136 accCount++;
137
138 llh = 2.0*llh;
139 return llh;
140}
141
MaCh3Plotting::PlottingManager * man
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
Base class for implementing fitting algorithms.
Definition: FitterBase.h:23
std::vector< SampleHandlerBase * > samples
Sample holder.
Definition: FitterBase.h:113
double logLProp
proposed likelihood
Definition: FitterBase.h:99
int accCount
counts accepted steps
Definition: FitterBase.h:103
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
double accProb
current acceptance prob
Definition: FitterBase.h:101
std::vector< double > sample_llh
store the llh breakdowns
Definition: FitterBase.h:108
double stepTime
Time of single step.
Definition: FitterBase.h:125
std::unique_ptr< TStopwatch > stepClock
tells how long single step/fit iteration took
Definition: FitterBase.h:123
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
TTree * outTree
Output tree with posteriors.
Definition: FitterBase.h:137
std::vector< ParameterHandlerBase * > systematics
Systematic holder.
Definition: FitterBase.h:118
int NParsPCA
Number of all parameters from all covariances in PCA base.
Definition: LikelihoodFit.h:30
LikelihoodFit(manager *const fitMan)
Constructor.
virtual ~LikelihoodFit()
Destructor.
int NPars
Number of all parameters from all covariances.
Definition: LikelihoodFit.h:28
void PrepareFit()
prepare output and perform sanity checks
bool fMirroring
Flag telling if mirroring is used or not.
Definition: LikelihoodFit.h:32
virtual double CalcChi2(const double *x)
Chi2 calculation over all included samples and syst objects.
The manager class is responsible for managing configurations and settings.
Definition: Manager.h:16
YAML::Node const & raw()
Return config.
Definition: Manager.h:41