Skip to content

Commit 78bf753

Browse files
nadrinonadrino
authored
Merge pull request #170 from nadrino/feature/banffFit
Feature/banff fit
2 parents 8dfc49b + e6346b8 commit 78bf753

File tree

15 files changed

+373
-64
lines changed

15 files changed

+373
-64
lines changed

‎src/FitParameters/include/DialSet.h

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,9 @@ class DialSet {
3535
void setOwner(const FitParameter* owner_);
3636
void setConfig(const nlohmann::json &config_);
3737

38+
void setMinDialResponse(double minDialResponse_){ _minDialResponse_ = minDialResponse_; }
39+
void setMaxDialResponse(double maxDialResponse_){ _maxDialResponse_ = maxDialResponse_; }
40+
3841
void initialize();
3942

4043
// Getters

‎src/FitParameters/include/FitParameter.h

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@ class FitParameter {
4040
void setParameterIndex(int parameterIndex);
4141
void setParameterValue(double parameterValue);
4242
void setPriorValue(double priorValue);
43+
void setThrowValue(double throwValue);
4344
void setStdDevValue(double stdDevValue);
4445
void setDialSetConfig(const nlohmann::json &jsonConfig_);
4546
void setParameterDefinitionConfig(const nlohmann::json &config_);
@@ -67,6 +68,7 @@ class FitParameter {
6768
double getParameterValue() const;
6869
double getStdDevValue() const;
6970
double getPriorValue() const;
71+
double getThrowValue() const;
7072
std::vector<DialSet> &getDialSetList();
7173
double getMinValue() const;
7274
double getMaxValue() const;
@@ -89,6 +91,7 @@ class FitParameter {
8991
int _parameterIndex_{-1}; // to get the right definition in the json config (in case "name" is not specified)
9092
double _parameterValue_{};
9193
double _priorValue_{};
94+
double _throwValue_{std::nan("UNSET")};
9295
double _stdDevValue_{};
9396
double _minValue_{std::nan("UNSET")};
9497
double _maxValue_{std::nan("UNSET")};

‎src/FitParameters/src/FitParameter.cpp

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,9 @@ void FitParameter::setParameterValue(double parameterValue) {
7171
void FitParameter::setPriorValue(double priorValue) {
7272
_priorValue_ = priorValue;
7373
}
74+
void FitParameter::setThrowValue(double throwValue){
75+
_throwValue_ = throwValue;
76+
}
7477
void FitParameter::setStdDevValue(double stdDevValue) {
7578
_stdDevValue_ = stdDevValue;
7679
}
@@ -188,6 +191,9 @@ double FitParameter::getStdDevValue() const {
188191
double FitParameter::getPriorValue() const {
189192
return _priorValue_;
190193
}
194+
double FitParameter::getThrowValue() const{
195+
return _throwValue_;
196+
}
191197
PriorType::PriorType FitParameter::getPriorType() const {
192198
return _priorType_;
193199
}

‎src/FitParameters/src/FitParameterSet.cpp

Lines changed: 27 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -170,6 +170,16 @@ void FitParameterSet::prepareFitParameters(){
170170
_originalParBuffer_ = std::make_shared<TVectorD>(_strippedCovarianceMatrix_->GetNrows() );
171171
_eigenParBuffer_ = std::make_shared<TVectorD>(_strippedCovarianceMatrix_->GetNrows() );
172172

173+
// LogAlert << "Disabling par/dial limits" << std::endl;
174+
// for( auto& par : _parameterList_ ){
175+
// par.setMinValue(std::nan(""));
176+
// par.setMaxValue(std::nan(""));
177+
// for( auto& dialSet : par.getDialSetList() ){
178+
// dialSet.setMinDialResponse(std::nan(""));
179+
// dialSet.setMaxDialResponse(std::nan(""));
180+
// }
181+
// }
182+
173183
// Put original parameters to the prior
174184
for( auto& par : _parameterList_ ){
175185
par.setValueAtPrior();
@@ -277,7 +287,7 @@ void FitParameterSet::throwFitParameters(double gain_){
277287

278288
LogThrowIf(_strippedCovarianceMatrix_==nullptr, "No covariance matrix provided")
279289

280-
if( not _useEigenDecompInFit_ ){
290+
// if( not _useEigenDecompInFit_ ){
281291
LogInfo << "Throwing parameters for " << _name_ << " using Cholesky matrix" << std::endl;
282292

283293
if( _choleskyMatrix_ == nullptr ){
@@ -294,23 +304,30 @@ void FitParameterSet::throwFitParameters(double gain_){
294304
if( par.isEnabled() and not par.isFixed() and not par.isFree() ){
295305
iFit++;
296306
LogInfo << "Throwing par " << par.getTitle() << ": " << par.getParameterValue();
297-
par.setParameterValue( par.getPriorValue() + gain_ * throws[iFit] );
307+
par.setThrowValue(par.getPriorValue() + gain_ * throws[iFit]);
308+
par.setParameterValue( par.getThrowValue() );
298309
LogInfo << " → " << par.getParameterValue() << std::endl;
299310
}
300311
else{
301312
LogWarning << "Skipping parameter: " << par.getTitle() << std::endl;
302313
}
303314
}
304-
}
305-
else{
306-
LogInfo << "Throwing eigen parameters for " << _name_ << std::endl;
315+
// }
316+
// else{
317+
// LogInfo << "Throwing eigen parameters for " << _name_ << std::endl;
318+
// for( auto& eigenPar : _eigenParameterList_ ){
319+
// if( eigenPar.isFixed() ){ LogWarning << "Eigen parameter #" << eigenPar.getParameterIndex() << " is fixed. Not throwing" << std::endl; continue; }
320+
// eigenPar.setThrowValue(eigenPar.getPriorValue() + gain_ * gRandom->Gaus(0, eigenPar.getStdDevValue()));
321+
// eigenPar.setParameterValue( eigenPar.getThrowValue() );
322+
// }
323+
// this->propagateEigenToOriginal();
324+
// }
325+
326+
if( _useEigenDecompInFit_ ){
327+
this->propagateOriginalToEigen();
307328
for( auto& eigenPar : _eigenParameterList_ ){
308-
if( eigenPar.isFixed() ){ LogWarning << "Eigen parameter #" << eigenPar.getParameterIndex() << " is fixed. Not throwing" << std::endl; continue; }
309-
eigenPar.setParameterValue(
310-
eigenPar.getPriorValue() + gain_ * gRandom->Gaus(0, eigenPar.getStdDevValue())
311-
);
329+
eigenPar.setThrowValue(eigenPar.getParameterValue());
312330
}
313-
this->propagateEigenToOriginal();
314331
}
315332

316333
}

‎src/FitSamples/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@ set( SRCFILES
66
src/FitSample.cpp
77
src/PhysicsEvent.cpp
88
src/PlotGenerator.cpp
9+
src/JointProbability.cpp
910
)
1011

1112
if( USE_STATIC_LINKS )

‎src/FitSamples/include/FitSampleSet.h

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -8,6 +8,7 @@
88
#include "FitSample.h"
99
#include "FitParameterSet.h"
1010
#include "Likelihoods.hh"
11+
#include "JointProbability.h"
1112

1213
#include "GenericToolbox.h"
1314
#include "nlohmann/json.hpp"
@@ -36,7 +37,7 @@ class FitSampleSet {
3637
const std::vector<FitSample> &getFitSampleList() const;
3738
std::vector<FitSample> &getFitSampleList();
3839
const nlohmann::json &getConfig() const;
39-
const std::shared_ptr<CalcLLHFunc> &getLikelihoodFunctionPtr() const;
40+
const std::shared_ptr<JointProbability::JointProbability> &getJointProbabilityFct() const;
4041

4142
//Core
4243
bool empty() const;
@@ -54,7 +55,7 @@ class FitSampleSet {
5455
nlohmann::json _config_;
5556

5657
std::vector<FitSample> _fitSampleList_;
57-
std::shared_ptr<CalcLLHFunc> _likelihoodFunctionPtr_{nullptr};
58+
std::shared_ptr<JointProbability::JointProbability> _jointProbabilityPtr_{nullptr};
5859
std::vector<std::string> _eventByEventDialLeafList_;
5960

6061
};

‎src/FitSamples/include/JointProbability.h

Lines changed: 54 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,54 @@
1+
//
2+
// Created by Adrien BLANCHET on 23/06/2022.
3+
//
4+
5+
#ifndef GUNDAM_JOINTPROBABILITY_H
6+
#define GUNDAM_JOINTPROBABILITY_H
7+
8+
#include "FitSample.h"
9+
10+
11+
12+
namespace JointProbability{
13+
14+
class JointProbability {
15+
16+
public:
17+
JointProbability() = default;
18+
virtual ~JointProbability() = default;
19+
20+
// two choices -> either override bin by bin llh or global eval function
21+
virtual double eval(const FitSample& sample_, int bin_){ return 0; }
22+
virtual double eval(const FitSample& sample_){
23+
double out{0};
24+
int nBins = int(sample_.getBinning().getBinsList().size());
25+
for( int iBin = 1 ; iBin <= nBins ; iBin++ ){ out += this->eval(sample_, iBin); }
26+
return out;
27+
}
28+
};
29+
30+
class PoissonLLH : public JointProbability{
31+
double eval(const FitSample& sample_, int bin_) override;
32+
};
33+
34+
class BarlowLLH : public JointProbability{
35+
double eval(const FitSample& sample_, int bin_) override;
36+
private:
37+
double rel_var, b, c, beta, mc_hat, chi2;
38+
};
39+
40+
class BarlowLLH_BANFF_OA2020 : public JointProbability{
41+
double eval(const FitSample& sample_, int bin_) override;
42+
};
43+
class BarlowLLH_BANFF_OA2021 : public JointProbability{
44+
double eval(const FitSample& sample_, int bin_) override;
45+
};
46+
47+
}
48+
49+
50+
51+
52+
53+
54+
#endif //GUNDAM_JOINTPROBABILITY_H

‎src/FitSamples/src/FitSampleSet.cpp

Lines changed: 9 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,6 @@ void FitSampleSet::reset() {
2626
_isInitialized_ = false;
2727
_config_.clear();
2828

29-
_likelihoodFunctionPtr_ = nullptr;
3029
_fitSampleList_.clear();
3130

3231
_eventByEventDialLeafList_.clear();
@@ -101,18 +100,11 @@ void FitSampleSet::initialize() {
101100

102101
std::string llhMethod = JsonUtils::fetchValue(_config_, "llhStatFunction", "PoissonLLH");
103102
LogInfo << "Using \"" << llhMethod << "\" LLH function." << std::endl;
104-
if( llhMethod == "PoissonLLH" ){
105-
_likelihoodFunctionPtr_ = std::make_shared<PoissonLLH>();
106-
}
107-
else if( llhMethod == "BarlowLLH" ){
108-
_likelihoodFunctionPtr_ = std::make_shared<BarlowLLH>();
109-
}
110-
else if( llhMethod == "BarlowLLH_OA2020_Bad" ){
111-
_likelihoodFunctionPtr_ = std::make_shared<BarlowOA2020BugLLH>();
112-
}
113-
else{
114-
LogThrow("Unknown LLH Method: " << llhMethod)
115-
}
103+
if ( llhMethod == "PoissonLLH" ){ _jointProbabilityPtr_ = std::make_shared<JointProbability::PoissonLLH>(); }
104+
else if( llhMethod == "BarlowLLH" ) { _jointProbabilityPtr_ = std::make_shared<JointProbability::BarlowLLH>(); }
105+
else if( llhMethod == "BarlowLLH_BANFF_OA2020" ) { _jointProbabilityPtr_ = std::make_shared<JointProbability::BarlowLLH_BANFF_OA2020>(); }
106+
else if( llhMethod == "BarlowLLH_BANFF_OA2021" ) { _jointProbabilityPtr_ = std::make_shared<JointProbability::BarlowLLH_BANFF_OA2021>(); }
107+
else{ LogThrow("Unknown LLH Method: " << llhMethod); }
116108

117109
_isInitialized_ = true;
118110
}
@@ -126,29 +118,20 @@ std::vector<FitSample> &FitSampleSet::getFitSampleList() {
126118
const nlohmann::json &FitSampleSet::getConfig() const {
127119
return _config_;
128120
}
129-
const std::shared_ptr<CalcLLHFunc> &FitSampleSet::getLikelihoodFunctionPtr() const {
130-
return _likelihoodFunctionPtr_;
121+
const std::shared_ptr<JointProbability::JointProbability> &FitSampleSet::getJointProbabilityFct() const{
122+
return _jointProbabilityPtr_;
131123
}
132124

133125
bool FitSampleSet::empty() const {
134126
return _fitSampleList_.empty();
135127
}
136128
double FitSampleSet::evalLikelihood() const{
137129
double llh = 0.;
138-
for( auto& sample : _fitSampleList_ ){
139-
llh += this->evalLikelihood(sample);
140-
}
130+
for( auto& sample : _fitSampleList_ ){ llh += this->evalLikelihood(sample); }
141131
return llh;
142132
}
143133
double FitSampleSet::evalLikelihood(const FitSample& sample_) const{
144-
double sampleLlh = 0;
145-
for( int iBin = 1 ; iBin <= sample_.getMcContainer().histogram->GetNbinsX() ; iBin++ ){
146-
sampleLlh += (*_likelihoodFunctionPtr_)(
147-
sample_.getMcContainer().histogram->GetBinContent(iBin),
148-
std::pow(sample_.getMcContainer().histogram->GetBinError(iBin), 2),
149-
sample_.getDataContainer().histogram->GetBinContent(iBin));
150-
}
151-
return sampleLlh;
134+
return _jointProbabilityPtr_->eval(sample_);
152135
}
153136

154137
void FitSampleSet::copyMcEventListToDataContainer(){

0 commit comments

Comments
 (0)