MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
MinuitFit.cpp
Go to the documentation of this file.
1#include "MinuitFit.h"
2
3// *******************
4// Run the Minuit Fit with all the systematic objects added
6// *******************
8 // Minimizer type: determines the underlying implementation.
9 // Available types include:
10 // - "Minuit2" (recommended modern option)
11 // - "Minuit" (legacy)
12 // - "Fumili"
13 // - "GSLMultiMin" (for gradient-free minimization)
14 // - "GSLMultiFit"
15 // - "GSLSimAn" (Simulated Annealing)
16 const std::string MinimizerType = "Minuit2";
17 // Minimizer algorithm (specific to the selected type).
18 // For Minuit2, the following algorithms are available:
19 // - "Migrad" : gradient-based minimization (default)
20 // - "Simplex" : Nelder-Mead simplex method (derivative-free)
21 // - "Combined" : combination of Simplex and Migrad
22 // - "Scan" : parameter grid scan
23 const std::string MinimizerAlgo = "Migrad";
24
25 MACH3LOG_INFO("Creating instance of Minimizer with {} and {}", MinimizerType, MinimizerAlgo);
26
27 minuit = std::unique_ptr<ROOT::Math::Minimizer>(
28 ROOT::Math::Factory::CreateMinimizer(MinimizerType.c_str(), MinimizerAlgo.c_str()));
29}
30
31// *************************
32// Destructor: close the logger and output file
34// *************************
35}
36
37
38// *******************
39// Run the Minuit with all the systematic objects added
41// *******************
42 PrepareFit();
43
44 // Remove obsolete memory and make other checks before fit starts
46
47 //KS: For none PCA this will be equal to normal parameters
48 const int NparsMinuitFull = NPars;
49 const int NparsMinuit = NParsPCA;
50
51 //KS: Set SetFunction we will Minimize
52 ROOT::Math::Functor fChi2(this, &MinuitFit::CalcChi2, NparsMinuit);
53 minuit->SetFunction(fChi2);
54
55 //KS: add config or something
56 minuit->SetPrintLevel(2);
57 minuit->SetTolerance(0.01);
58 minuit->SetMaxFunctionCalls(fitMan->raw()["General"]["Minuit2"]["NSteps"].as<unsigned>());
59 minuit->SetMaxIterations(10000);
60
61 MACH3LOG_INFO("Preparing Minuit");
62 int ParCounter = 0;
63
64 for (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
65 {
66 if(!(*it)->IsPCA())
67 {
68 for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
69 {
70 //KS: Index, name, prior, step scale [different to MCMC],
71 minuit->SetVariable(ParCounter, ((*it)->GetParName(i)), (*it)->GetParInit(i), (*it)->GetDiagonalError(i)/10);
72 minuit->SetVariableValue(ParCounter, (*it)->GetParInit(i));
73 //KS: lower bound, upper bound, if Mirroring enabled then ignore
74 if(!fMirroring) minuit->SetVariableLimits(ParCounter, (*it)->GetLowerBound(i), (*it)->GetUpperBound(i));
75 if((*it)->IsParameterFixed(i))
76 {
77 minuit->FixVariable(ParCounter);
78 }
79 }
80 }
81 else
82 {
83 for(int i = 0; i < (*it)->GetNParameters(); ++i, ++ParCounter)
84 {
85 minuit->SetVariable(ParCounter, Form("%i_PCA", i), (*it)->GetPCAHandler()->GetParPropPCA(i), (*it)->GetPCAHandler()->GetEigenValuesMaster()[i]/10);
86 if((*it)->GetPCAHandler()->IsParameterFixedPCA(i))
87 {
88 minuit->FixVariable(ParCounter);
89 }
90 }
91 }
92 }
93
94 minuit->SetPrintLevel(2);
95
96 MACH3LOG_INFO("Starting MIGRAD");
97 minuit->Minimize();
98
99 MACH3LOG_INFO("Starting HESSE");
100 minuit->Hesse();
101 outputFile->cd();
102
103 TVectorD* MinuitParValue = new TVectorD(NparsMinuitFull);
104 TVectorD* MinuitParError = new TVectorD(NparsMinuitFull);
105 TMatrixDSym* Postmatrix = new TMatrixDSym(NparsMinuitFull);
106
107 for(int i = 0; i < NparsMinuitFull; ++i)
108 {
109 (*MinuitParValue)(i) = 0;
110 (*MinuitParError)(i) = 0;
111 for(int j = 0; j < NparsMinuitFull; ++j)
112 {
113 (*Postmatrix)(i,j) = 0;
114 (*Postmatrix)(i,j) = minuit->CovMatrix(i,j);
115 }
116 }
117
118 ParCounter = 0;
119 const double *X = minuit->X();
120 const double *err = minuit->Errors();
121 for (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
122 {
123 if(!(*it)->IsPCA())
124 {
125 for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
126 {
127 double ParVal = X[ParCounter];
128 //KS: Basically apply mirroring for parameters out of bounds
129 if(fMirroring)
130 {
131 if(ParVal < (*it)->GetLowerBound(i))
132 {
133 ParVal = (*it)->GetLowerBound(i) + ((*it)->GetLowerBound(i) - ParVal);
134 }
135 else if (ParVal > (*it)->GetUpperBound(i))
136 {
137 ParVal = (*it)->GetUpperBound(i) - ( ParVal - (*it)->GetUpperBound(i));
138 }
139 }
140 (*MinuitParValue)(ParCounter) = ParVal;
141 (*MinuitParError)(ParCounter) = err[ParCounter];
142 //KS: For fixed params HESS will not calculate error so we need to pass prior error
143 if((*it)->IsParameterFixed(i))
144 {
145 (*MinuitParError)(ParCounter) = (*it)->GetDiagonalError(i);
146 (*Postmatrix)(ParCounter,ParCounter) = (*MinuitParError)(ParCounter) * (*MinuitParError)(ParCounter);
147 }
148 }
149 }
150 else
151 {
152 //KS: We need to convert parameters from PCA to normal base
153 TVectorD ParVals((*it)->GetNumParams());
154 TVectorD ParVals_PCA((*it)->GetNParameters());
155
156 TVectorD ErrorVals((*it)->GetNumParams());
157 TVectorD ErrorVals_PCA((*it)->GetNParameters());
158
159 TMatrixD MatrixVals((*it)->GetNumParams(), (*it)->GetNumParams());
160 TMatrixD MatrixVals_PCA((*it)->GetNParameters(), (*it)->GetNParameters());
161
162 //First save them
163 //KS: This code is super convoluted as MaCh3 can store separate matrices while Minuit has one matrix. In future this will be simplified, keep it like this for now.
164 const int StartVal = ParCounter;
165 for(int i = 0; i < (*it)->GetNParameters(); ++i, ++ParCounter)
166 {
167 ParVals_PCA(i) = X[ParCounter];
168 ErrorVals_PCA(i) = err[ParCounter];
169 int ParCounterMatrix = StartVal;
170 for(int j = 0; j < (*it)->GetNParameters(); ++j, ++ParCounterMatrix)
171 {
172 MatrixVals_PCA(i,j) = minuit->CovMatrix(ParCounter,ParCounterMatrix);
173 }
174 }
175 ParVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ParVals_PCA;
176 ErrorVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ErrorVals_PCA;
177 MatrixVals.Mult(((*it)->GetPCAHandler()->GetTransferMatrix()),MatrixVals_PCA);
178
179 ParCounter = StartVal;
180 //KS: Now after going from PCA to normal let';s save it
181 for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
182 {
183 (*MinuitParValue)(ParCounter) = ParVals(i);
184 (*MinuitParError)(ParCounter) = std::fabs(ErrorVals(i));
185 int ParCounterMatrix = StartVal;
186 for(int j = 0; j < (*it)->GetNumParams(); ++j, ++ParCounterMatrix)
187 {
188 (*Postmatrix)(ParCounter,ParCounterMatrix) = MatrixVals(i,j);
189 }
190 //If fixed take prior
191 if((*it)->GetPCAHandler()->IsParameterFixedPCA(i))
192 {
193 (*MinuitParError)(ParCounter) = (*it)->GetDiagonalError(i);
194 (*Postmatrix)(ParCounter,ParCounter) = (*MinuitParError)(ParCounter) * (*MinuitParError)(ParCounter);
195 }
196 }
197 }
198 }
199
200 MinuitParValue->Write("MinuitParValue");
201 MinuitParError->Write("MinuitParError");
202 Postmatrix->Write("Postmatrix");
203 delete MinuitParValue;
204 delete MinuitParError;
205 delete Postmatrix;
206 // Save all the output
207 SaveOutput();
208}
209
MaCh3Plotting::PlottingManager * man
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
void SaveOutput()
Save output and close files.
Definition: FitterBase.cpp:221
TFile * outputFile
Output.
Definition: FitterBase.h:131
manager * fitMan
The manager.
Definition: FitterBase.h:92
void SanitiseInputs()
Remove obsolete memory and make other checks before fit starts.
Definition: FitterBase.cpp:213
std::vector< ParameterHandlerBase * > systematics
Systematic holder.
Definition: FitterBase.h:118
Implementation of base Likelihood Fit class, it is mostly responsible for likelihood calculation whil...
Definition: LikelihoodFit.h:6
int NParsPCA
Number of all parameters from all covariances in PCA base.
Definition: LikelihoodFit.h:30
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.
std::unique_ptr< ROOT::Math::Minimizer > minuit
Pointer to minimizer, which most often is Minuit.
Definition: MinuitFit.h:29
MinuitFit(manager *const fitMan)
Constructor.
Definition: MinuitFit.cpp:5
void RunMCMC() override
Actual implementation of Minuit Fit algorithm.
Definition: MinuitFit.cpp:40
virtual ~MinuitFit()
Destructor.
Definition: MinuitFit.cpp:33
The manager class is responsible for managing configurations and settings.
Definition: Manager.h:16
YAML::Node const & raw()
Return config.
Definition: Manager.h:41