Skip to content

Commit

Permalink
[fix] bootstrap with combined fit
Browse files Browse the repository at this point in the history
  • Loading branch information
@Marcellocosti
Marcellocosti committed Aug 20, 2024
1 parent 578e4c6 commit 12fca06
Show file tree
Hide file tree
Showing 2 changed files with 45 additions and 63 deletions.
80 changes: 40 additions & 40 deletions fempy/CombinedFitter.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,10 +40,13 @@ class GlobalChi2 {
fPars = pars;
fModelsSharedParsIdxs = sharedidxs;
DEBUG("fPars size: " << fPars.size());
for(int iPar=0; iPar<fPars.size(); iPar++) {
fChi2Functions.push_back(nullptr);
}
}

void AddChi2Function(ROOT::Math::IMultiGenFunction &func) {
fChi2Functions.push_back(&func);
void SetChi2Function(ROOT::Math::IMultiGenFunction &func, int idx) {
fChi2Functions[idx] = &func;
}

double operator() (const double *par) const {
Expand Down Expand Up @@ -87,7 +90,6 @@ class CombinedFitter {
// directly from the CorrelationFitter objects, but I was not able to do so
for(int iModel=0; iModel<this->fModels.size(); iModel++) {
fModelFuncts.push_back(fModels[iModel].GetFitFunction());
cout << "POINTER ADDRESS FMODELFUNCTS CONSTRUCTOR: " << &fModelFuncts[iModel] << endl;
}

SetTotalParameters();
Expand Down Expand Up @@ -137,7 +139,7 @@ class CombinedFitter {
// Data-fit difference extracting genuine correlation
if(difference) {
for(int iModel=0; iModel<this->fModels.size(); iModel++) {
hGenDiffOriginal.push_back(this->fModels[iModel].GetGenuineDifference());
hGenDiffOriginal.push_back(this->fModels[iModel].GetGenuineDifference(static_cast<TH1D *>(this->fModels[iModel].GetFitHisto())));
std::string titleWithModel = hGenDiffOriginal.back()->GetTitle();
hGenDiffOriginal.back()->SetTitle( (titleWithModel + Form("_Model%i", iModel)).c_str());
hGenDiffOriginal.back()->SetName( (titleWithModel + Form("_Model%i", iModel)).c_str());
Expand All @@ -148,47 +150,47 @@ class CombinedFitter {
double binGenuine = hGenDiffOriginal.back()->GetBinContent(iBin+1);
hBTHistos.push_back(new TH1D(Form("Model%i_Subtraction_%iMeV", iModel, binCenter),
Form("Model%i_Subtraction_%iMeV", iModel, binCenter),
10000, binGenuine - 0.1*binGenuine, binGenuine + 0.1*binGenuine));
500, binGenuine - 0.1*binGenuine, binGenuine + 0.1*binGenuine));
}
}
}

// Perform fits with sampled histos
DEBUG("Filling models");
std::vector<ROOT::Math::WrappedMultiTF1> wrappedFitFuncts;
ROOT::Fit::DataOptions opt;
std::vector<ROOT::Fit::DataRange> fitRanges;
std::vector<ROOT::Fit::BinData> fitBinData;

// Set fit range
DEBUG("Fit range");
for(int iModel=0; iModel<this->fModels.size(); iModel++) {
wrappedFitFuncts.push_back(ROOT::Math::WrappedMultiTF1(*fModelFuncts[iModel], 1));
fitRanges.push_back(ROOT::Fit::DataRange());
fitRanges.back().SetRange(fModels[iModel].GetLowFitRange(), fModels[iModel].GetUppFitRange());
fitBinData.push_back(ROOT::Fit::BinData(opt, fitRanges.back()));
}

for(int iTry=0; iTry<ntries; iTry++) {
std::vector<TH1D *> sampledHistos;
TH1D *sampledHistos[this->fModels.size()];
std::vector<ROOT::Fit::Chi2Function> chi2Functions;
GlobalChi2 global_chi2(fSingleModelsUniquePars, fSingleModelsSharedPars);

// Perform fits with sampled histos
DEBUG("Filling models");
std::vector<ROOT::Math::WrappedMultiTF1> wrappedFitFuncts;
ROOT::Fit::DataOptions opt;
std::vector<ROOT::Fit::DataRange> fitRanges;
std::vector<ROOT::Fit::BinData> fitBinData;

// Set fit range
DEBUG("Fit range");
for(int iModel=0; iModel<this->fModels.size(); iModel++) {
wrappedFitFuncts.push_back(ROOT::Math::WrappedMultiTF1(*fModelFuncts[iModel], 1));
fitRanges.push_back(ROOT::Fit::DataRange());
fitRanges.back().SetRange(fModels[iModel].GetLowFitRange(), fModels[iModel].GetUppFitRange());
fitBinData.push_back(ROOT::Fit::BinData(opt, fitRanges.back()));
}
int nPar = 0;
for(int iModel=0; iModel<this->fModels.size(); iModel++) {

sampledHistos.push_back(SampledHisto(static_cast<TH1D *>(fModels[iModel].GetFitHisto()),
fModels[iModel].GetUppFitRange(), iTry));
DEBUG("GetNbins histo: " << sampledHistos.back()->GetNbinsX());
ROOT::Fit::FillData(fitBinData[iModel], sampledHistos.back());
sampledHistos[iModel] = SampledHisto(static_cast<TH1D *>(fModels[iModel].GetFitHisto()),
fModels[iModel].GetUppFitRange(), iTry);
DEBUG("GetNbins histo: " << sampledHistos[iModel]->GetNbinsX());
ROOT::Fit::FillData(fitBinData[iModel], sampledHistos[iModel]);
nPar += fitBinData[iModel].Size();

DEBUG("Global chi2");
chi2Functions.push_back(ROOT::Fit::Chi2Function(fitBinData[iModel], wrappedFitFuncts[iModel]));

}

for(int iModel=0; iModel<this->fModels.size(); iModel++) {
global_chi2.AddChi2Function(chi2Functions[iModel]);
global_chi2.SetChi2Function(chi2Functions[iModel], iModel);
}

DEBUG("Fitter");
Expand All @@ -215,15 +217,15 @@ class CombinedFitter {
// Data-fit difference extracting genuine correlation
if(difference) {
for(int iModel=0; iModel<this->fModels.size(); iModel++) {
TH1D *hGenDiff = this->fModels[iModel].GetGenuineDifference();
TH1D *hGenDiff = this->fModels[iModel].GetGenuineDifference(sampledHistos[iModel]);
int nBinsFitCF = hGenDiff->GetNbinsX();
for(int iBin=0; iBin<nBinsFitCF; iBin++) {
hBTHistos[this->fInitPars.size() + iBin + iModel*nBinsFitCF]->Fill(hGenDiff->GetBinContent(iBin+1));
hBTHistos[this->fInitPars.size() + iModel * nBinsFitCF + iBin]->Fill(hGenDiff->GetBinContent(iBin+1));
}
delete sampledHistos[iModel];
delete hGenDiff;
}
}

sampledHistos.clear();
chi2Functions.clear();
}

Expand All @@ -232,9 +234,6 @@ class CombinedFitter {
double binWidth = this->fModels[iModel].GetFitHisto()->GetBinWidth(1);
int nBinsFitCF = fModels[iModel].GetUppFitRange()/binWidth;

// // Restore original fit results to evaluate the subtraction
// fModelFuncts[iModel]->SetFitResult(originalFitResult, fSingleModelsUniquePars[iModel].data());

TH1D *hYields = new TH1D(Form("hYields_Model%i_stat", iModel), Form("hYields_Model%i_stat", iModel), nBinsFitCF, 0, this->fModels[iModel].GetUppFitRange());
TH1D *hMeans = new TH1D(Form("hMeans_Model%i_stat", iModel), Form("hMeans_Model%i_stat", iModel), nBinsFitCF, 0, this->fModels[iModel].GetUppFitRange());
TH1D *hStdDevs = new TH1D(Form("hStdDevs_Model%i_stat", iModel), Form("hStdDevs_Model%i_stat", iModel), nBinsFitCF, 0, this->fModels[iModel].GetUppFitRange());
Expand All @@ -244,8 +243,9 @@ class CombinedFitter {
int iHistoIdx = this->fInitPars.size() + iModel * nBinsFitCF + iBin;
TF1 *gaus = new TF1("gaus", "gaus", hBTHistos[iHistoIdx]->GetBinLowEdge(1),
hBTHistos[iHistoIdx]->GetBinLowEdge(hBTHistos[iHistoIdx]->GetNbinsX()) + hBTHistos[iHistoIdx]->GetBinWidth(1));
cout << "Fitting iBin " << iBin+1 << " of model " << iModel << endl;
// cout << "HISTO TO BE FITTED LIMITS: " << hBTHistos[iHistoIdx]->GetBinLowEdge(1) << " " << hBTHistos[iHistoIdx]->GetBinLowEdge(hBTHistos[iHistoIdx]->GetNbinsX()) + hBTHistos[iHistoIdx]->GetBinWidth(1) << endl;
if((iBin+1)%10 == 0) {
cout << "Fitting iBin " << iBin+1 << " of model " << iModel << endl;
}
gaus->SetParameter(1, hBTHistos[iHistoIdx]->GetBinContent(hBTHistos[iHistoIdx]->GetMaximumBin()) );
gaus->SetParameter(2, hGenDiffOriginal[iModel]->GetBinContent(iBin+1) );
hBTHistos[iHistoIdx]->Fit(gaus, "SMRL+q", "");
Expand All @@ -254,7 +254,6 @@ class CombinedFitter {
hStdDevs->SetBinContent(iBin+1, gaus->GetParameter(2));
hRelUnc->SetBinContent(iBin+1, hStdDevs->GetBinContent(iBin+1) / hGenDiffOriginal[iModel]->GetBinContent(iBin+1));

// cout << "iBin" << iBin+1 << ", genuine " << hGenDiffOriginal[iModel]->GetBinContent(iBin+1) << endl;
hGenDiffOriginal[iModel]->SetBinError(iBin+1, hStdDevs->GetBinContent(iBin+1));

}
Expand Down Expand Up @@ -300,7 +299,7 @@ class CombinedFitter {
}
GlobalChi2 global_chi2(fSingleModelsUniquePars, fSingleModelsSharedPars);
for(int iChi2Fcn=0; iChi2Fcn<chi2Functions.size(); iChi2Fcn++) {
global_chi2.AddChi2Function(chi2Functions[iChi2Fcn]);
global_chi2.SetChi2Function(chi2Functions[iChi2Fcn], iChi2Fcn);
}

DEBUG("Fitter");
Expand Down Expand Up @@ -361,7 +360,7 @@ class CombinedFitter {
}
GlobalChi2 global_chi2(fSingleModelsUniquePars, fSingleModelsSharedPars);
for(int iChi2Fcn=0; iChi2Fcn<chi2Functions.size(); iChi2Fcn++) {
global_chi2.AddChi2Function(chi2Functions[iChi2Fcn]);
global_chi2.SetChi2Function(chi2Functions[iChi2Fcn], iChi2Fcn);
}

DEBUG("Fitter");
Expand All @@ -385,7 +384,8 @@ class CombinedFitter {
for(int iModel=0; iModel<this->fModels.size(); iModel++) {
fModelFuncts[iModel]->SetFitResult(result, fSingleModelsUniquePars[iModel].data());
fModels[iModel].SetFitFunction(fModelFuncts[iModel]);
hDifferences.push_back(this->fModels[iModel].GetGenuineDifference());
hDifferences.push_back(this->fModels[iModel].GetGenuineDifference(static_cast<TH1D *>(
this->fModels[iModel].GetFitHisto())));
}

return hDifferences;
Expand Down
28 changes: 5 additions & 23 deletions fempy/CorrelationFitter.hxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,12 +30,11 @@ class CorrelationFitter {
this->fFit = new TF1("fitFunction", "gaus", fFitRangeMin, fFitRangeMax);
this->fNPars = {0}; // The first parameter has index zero
this->fLambdaGen = 0;
cout << "POINTER ADDRESS FFIT CONSTRUCTOR: " << &this->fFit << endl;

}

void SetLambdaGen(double lambdagen) {
cout << "SETTING LAMBDA GENUINE TO: " << lambdagen << endl;
DEBUG("SETTING LAMBDA GENUINE TO: " << lambdagen);
this->fLambdaGen = lambdagen;
}

Expand Down Expand Up @@ -652,7 +651,7 @@ class CorrelationFitter {
double binWidth = this->fFitHist->GetBinWidth(1);
int nBinsFitCF = static_cast<int>(this->fFitRangeMax/binWidth);
TH1D *hSubtraction = new TH1D("hSubtraction", "hSubtraction", nBinsFitCF, 0, static_cast<int>(nBinsFitCF * fFitHist->GetBinWidth(1)));
TH1D *hDiffGenuineOriginal = GetGenuineDifference();
TH1D *hDiffGenuineOriginal = GetGenuineDifference(static_cast<TH1D *>(this->fFitHist));
for(int iBin=0; iBin<nBinsFitCF; iBin++) {
hSubtraction->SetBinContent(iBin+1, hDiffGenuineOriginal->GetBinContent(iBin+1) );
}
Expand Down Expand Up @@ -710,7 +709,7 @@ class CorrelationFitter {
hFitParsBT[iPar]->Fill(fitResults->Parameter(iPar));
}
hFitParsBT[this->fFit->GetNpar()]->Fill(fFit->GetChisquare() / fFit->GetNDF());
TH1D *hDiffGenuineTry = GetGenuineDifference();
TH1D *hDiffGenuineTry = GetGenuineDifference(sampledCF);
for(int iBin=0; iBin<nBinsFitCF; iBin++) {
hBinsCFDifferences[iBin]->Fill(hDiffGenuineTry->GetBinContent(iBin+1) );
}
Expand Down Expand Up @@ -755,8 +754,6 @@ class CorrelationFitter {
}

TH1D *GetfFitHisto() {
cout << "GET fFit HISTO" << endl;
cout << "POINTER ADDRESS fFit " << this->fFit << endl;
// TH1D *hfFit = new TH1D(Form("hCheck_fFit_%i", imodel), Form("hCheck_fFit_%i", imodel), 125, 0, 500);
TH1D *hfFit = new TH1D("hCheck_fFit", "hCheck_fFit", 125, 0, 500);
for(int iBin=0; iBin<hfFit->GetNbinsX(); iBin++) {
Expand All @@ -766,8 +763,6 @@ class CorrelationFitter {
}

TH1D *GetfFitHisto(TH1D *histo) {
cout << "GET fFit HISTO" << endl;
cout << "POINTER ADDRESS fFit " << this->fFit << endl;
TH1D *hfFit = static_cast<TH1D *>(histo->Clone("hCheck_fFit"));
hfFit->Reset("ICESM");
for(int iBin=0; iBin<hfFit->GetNbinsX(); iBin++) {
Expand All @@ -776,7 +771,7 @@ class CorrelationFitter {
return hfFit;
}

TH1D *GetGenuineDifference() {
TH1D *GetGenuineDifference(TH1D *histo) {

// Implemented formula for multiplicative minijet
// C_gen = (1 / (GlobNorm * lambda_gen * C_MJ) ) * (DATA - GlobNorm * (lambda_flat + other comps) )
Expand All @@ -787,25 +782,13 @@ class CorrelationFitter {
static_cast<int>(nBinsFitCF * fFitHist->GetBinWidth(1)));

TF1 *fMJtimesGlobNorm = this->GetBaseline(this->fGlobNorm);
TF1 *fMJtimesGlobNormtimesLambdaGen = this->GetBaseline(this->fGlobNorm, 1);

for(int iBin=0; iBin<hGenuineDiff->GetNbinsX(); iBin++) {
double binCenter = this->fFitHist->GetBinCenter(iBin+1);
// cout << "***************" << endl;
// cout << "Bin center: " << this->fFitHist->GetBinCenter(iBin+1) << endl;
// cout << "Setting bin content of iBin" << iBin+1 << " to " << 1 + (this->fFitHist->GetBinContent(iBin+1) - this->fFit->Eval(binCenter)) /
// (this->fLambdaGen * fMJtimesGlobNorm->Eval(binCenter)) << endl;
// cout << "Data value: " << this->fFitHist->GetBinContent(iBin+1) << endl;
// cout << "fFit value: " << this->fFit->Eval(binCenter) << endl;
// cout << "MJ value: " << fMJtimesGlobNorm->Eval(binCenter) << endl;
// cout << "***************" << endl;
// hGenuineDiff->SetBinContent(iBin+1, this->fFit->Eval(binCenter));
hGenuineDiff->SetBinContent(iBin+1,
1 + (this->fFitHist->GetBinContent(iBin+1) - this->fFit->Eval(binCenter)) /
1 + (histo->GetBinContent(iBin+1) - this->fFit->Eval(binCenter)) /
(this->fLambdaGen * fMJtimesGlobNorm->Eval(binCenter)));
}
cout << "POINTER ADDRESS FFIT: " << &this->fFit << endl;
cout << "CIAO" << endl;
return hGenuineDiff;
}

Expand Down Expand Up @@ -846,7 +829,6 @@ class CorrelationFitter {
// compute CF with gaussian sampling
TH1D *sampledCF = new TH1D(Form("hTry_%i", iTry), "", nBinsFitCF,
this->fFitHist->GetBinLowEdge(1), this->fFitRangeMax);
// cout << "CIAO4" << endl;
for(int iSampledBin=0; iSampledBin<sampledCF->GetNbinsX(); iSampledBin++) {
double CFvalue = this->fFitHist->GetBinContent(iSampledBin+1);
double CFerror = this->fFitHist->GetBinError(iSampledBin+1);
Expand Down

0 comments on commit 12fca06

Please sign in to comment.