16 const std::string MinimizerType =
"Minuit2";
23 const std::string MinimizerAlgo =
"Migrad";
25 MACH3LOG_INFO(
"Creating instance of Minimizer with {} and {}", MinimizerType, MinimizerAlgo);
27 minuit = std::unique_ptr<ROOT::Math::Minimizer>(
28 ROOT::Math::Factory::CreateMinimizer(MinimizerType.c_str(), MinimizerAlgo.c_str()));
48 const int NparsMinuitFull =
NPars;
53 minuit->SetFunction(fChi2);
57 minuit->SetTolerance(0.01);
58 minuit->SetMaxFunctionCalls(
fitMan->
raw()[
"General"][
"Minuit2"][
"NSteps"].as<
unsigned>());
59 minuit->SetMaxIterations(10000);
68 for(
int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
71 minuit->SetVariable(ParCounter, ((*it)->GetParName(i)), (*it)->GetParInit(i), (*it)->GetDiagonalError(i)/10);
72 minuit->SetVariableValue(ParCounter, (*it)->GetParInit(i));
74 if(!
fMirroring)
minuit->SetVariableLimits(ParCounter, (*it)->GetLowerBound(i), (*it)->GetUpperBound(i));
75 if((*it)->IsParameterFixed(i))
77 minuit->FixVariable(ParCounter);
83 for(
int i = 0; i < (*it)->GetNParameters(); ++i, ++ParCounter)
85 minuit->SetVariable(ParCounter, Form(
"%i_PCA", i), (*it)->GetPCAHandler()->GetParPropPCA(i), (*it)->GetPCAHandler()->GetEigenValuesMaster()[i]/10);
86 if((*it)->GetPCAHandler()->IsParameterFixedPCA(i))
88 minuit->FixVariable(ParCounter);
103 TVectorD* MinuitParValue =
new TVectorD(NparsMinuitFull);
104 TVectorD* MinuitParError =
new TVectorD(NparsMinuitFull);
105 TMatrixDSym* Postmatrix =
new TMatrixDSym(NparsMinuitFull);
107 for(
int i = 0; i < NparsMinuitFull; ++i)
109 (*MinuitParValue)(i) = 0;
110 (*MinuitParError)(i) = 0;
111 for(
int j = 0; j < NparsMinuitFull; ++j)
113 (*Postmatrix)(i,j) = 0;
114 (*Postmatrix)(i,j) =
minuit->CovMatrix(i,j);
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)
125 for(
int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
127 double ParVal = X[ParCounter];
131 if(ParVal < (*it)->GetLowerBound(i))
133 ParVal = (*it)->GetLowerBound(i) + ((*it)->GetLowerBound(i) - ParVal);
135 else if (ParVal > (*it)->GetUpperBound(i))
137 ParVal = (*it)->GetUpperBound(i) - ( ParVal - (*it)->GetUpperBound(i));
140 (*MinuitParValue)(ParCounter) = ParVal;
141 (*MinuitParError)(ParCounter) = err[ParCounter];
143 if((*it)->IsParameterFixed(i))
145 (*MinuitParError)(ParCounter) = (*it)->GetDiagonalError(i);
146 (*Postmatrix)(ParCounter,ParCounter) = (*MinuitParError)(ParCounter) * (*MinuitParError)(ParCounter);
153 TVectorD ParVals((*it)->GetNumParams());
154 TVectorD ParVals_PCA((*it)->GetNParameters());
156 TVectorD ErrorVals((*it)->GetNumParams());
157 TVectorD ErrorVals_PCA((*it)->GetNParameters());
159 TMatrixD MatrixVals((*it)->GetNumParams(), (*it)->GetNumParams());
160 TMatrixD MatrixVals_PCA((*it)->GetNParameters(), (*it)->GetNParameters());
164 const int StartVal = ParCounter;
165 for(
int i = 0; i < (*it)->GetNParameters(); ++i, ++ParCounter)
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)
172 MatrixVals_PCA(i,j) =
minuit->CovMatrix(ParCounter,ParCounterMatrix);
175 ParVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ParVals_PCA;
176 ErrorVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ErrorVals_PCA;
177 MatrixVals.Mult(((*it)->GetPCAHandler()->GetTransferMatrix()),MatrixVals_PCA);
179 ParCounter = StartVal;
181 for(
int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
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)
188 (*Postmatrix)(ParCounter,ParCounterMatrix) = MatrixVals(i,j);
191 if((*it)->GetPCAHandler()->IsParameterFixedPCA(i))
193 (*MinuitParError)(ParCounter) = (*it)->GetDiagonalError(i);
194 (*Postmatrix)(ParCounter,ParCounter) = (*MinuitParError)(ParCounter) * (*MinuitParError)(ParCounter);
200 MinuitParValue->Write(
"MinuitParValue");
201 MinuitParError->Write(
"MinuitParError");
202 Postmatrix->Write(
"Postmatrix");
203 delete MinuitParValue;
204 delete MinuitParError;
MaCh3Plotting::PlottingManager * man
void SaveOutput()
Save output and close files.
TFile * outputFile
Output.
manager * fitMan
The manager.
void SanitiseInputs()
Remove obsolete memory and make other checks before fit starts.
std::vector< ParameterHandlerBase * > systematics
Systematic holder.
Implementation of base Likelihood Fit class, it is mostly responsible for likelihood calculation whil...
int NParsPCA
Number of all parameters from all covariances in PCA base.
int NPars
Number of all parameters from all covariances.
void PrepareFit()
prepare output and perform sanity checks
bool fMirroring
Flag telling if mirroring is used or not.
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.
MinuitFit(manager *const fitMan)
Constructor.
void RunMCMC() override
Actual implementation of Minuit Fit algorithm.
virtual ~MinuitFit()
Destructor.
The manager class is responsible for managing configurations and settings.
YAML::Node const & raw()
Return config.