File indexing completed on 2021-02-14 12:45:33
0001 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals1DOFFitter.h"
0002 #include "TMath.h"
0003
0004 static TMinuit *MuonResiduals1DOFFitter_TMinuit;
0005 static double MuonResiduals1DOFFitter_sum_of_weights;
0006 static double MuonResiduals1DOFFitter_number_of_hits;
0007 static bool MuonResiduals1DOFFitter_weightAlignment;
0008
0009 void MuonResiduals1DOFFitter::inform(TMinuit *tMinuit) { MuonResiduals1DOFFitter_TMinuit = tMinuit; }
0010
0011 void MuonResiduals1DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
0012 MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo *)(MuonResiduals1DOFFitter_TMinuit->GetObjectFit());
0013 MuonResidualsFitter *fitter = fitinfo->fitter();
0014
0015 fval = 0.;
0016 for (std::vector<double *>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end();
0017 ++resiter) {
0018 const double residual = (*resiter)[MuonResiduals1DOFFitter::kResid];
0019 const double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
0020
0021 const double residpeak = par[MuonResiduals1DOFFitter::kAlign];
0022 const double residsigma = par[MuonResiduals1DOFFitter::kSigma];
0023 const double residgamma = par[MuonResiduals1DOFFitter::kGamma];
0024
0025 double weight = (1. / redchi2) * MuonResiduals1DOFFitter_number_of_hits / MuonResiduals1DOFFitter_sum_of_weights;
0026 if (!MuonResiduals1DOFFitter_weightAlignment)
0027 weight = 1.;
0028
0029 if (!MuonResiduals1DOFFitter_weightAlignment || TMath::Prob(redchi2 * 8, 8) < 0.99) {
0030 if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
0031 fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
0032 } else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
0033 fval += -weight * MuonResidualsFitter_logPowerLawTails(residual, residpeak, residsigma, residgamma);
0034 } else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
0035 fval += -weight * MuonResidualsFitter_logROOTVoigt(residual, residpeak, residsigma, residgamma);
0036 } else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
0037 fval += -weight * MuonResidualsFitter_logGaussPowerTails(residual, residpeak, residsigma);
0038 } else {
0039 assert(false);
0040 }
0041 }
0042 }
0043 }
0044
0045 double MuonResiduals1DOFFitter::sumofweights() {
0046 MuonResiduals1DOFFitter_sum_of_weights = 0.;
0047 MuonResiduals1DOFFitter_number_of_hits = 0.;
0048 MuonResiduals1DOFFitter_weightAlignment = m_weightAlignment;
0049 for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0050 if (m_weightAlignment) {
0051 double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
0052 if (TMath::Prob(redchi2 * 8, 8) < 0.99) {
0053 MuonResiduals1DOFFitter_sum_of_weights += 1. / redchi2;
0054 MuonResiduals1DOFFitter_number_of_hits += 1.;
0055 }
0056 } else {
0057 MuonResiduals1DOFFitter_sum_of_weights += 1.;
0058 MuonResiduals1DOFFitter_number_of_hits += 1.;
0059 }
0060 }
0061 return MuonResiduals1DOFFitter_sum_of_weights;
0062 }
0063
0064 bool MuonResiduals1DOFFitter::fit(Alignable *ali) {
0065 initialize_table();
0066 sumofweights();
0067
0068 double resid_sum = 0.;
0069 double resid_sum2 = 0.;
0070 double resid_N = 0.;
0071 int N = 0;
0072
0073 for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0074 const double residual = (*resiter)[MuonResiduals1DOFFitter::kResid];
0075 const double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
0076 double weight = 1. / redchi2;
0077 if (!m_weightAlignment)
0078 weight = 1.;
0079
0080 if (!m_weightAlignment || TMath::Prob(redchi2 * 8, 8) < 0.99) {
0081 if (fabs(residual) < 10.) {
0082 resid_sum += weight * residual;
0083 resid_sum2 += weight * residual * residual;
0084 resid_N += weight;
0085 N++;
0086 }
0087 }
0088 }
0089
0090 double resid_mean = resid_sum / resid_N;
0091 double resid_stdev = sqrt(resid_sum2 / resid_N - pow(resid_sum / resid_N, 2));
0092
0093 resid_sum = 0.;
0094 resid_sum2 = 0.;
0095 resid_N = 0.;
0096
0097 for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0098 const double residual = (*resiter)[MuonResiduals1DOFFitter::kResid];
0099 const double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
0100 double weight = 1. / redchi2;
0101 if (!m_weightAlignment)
0102 weight = 1.;
0103
0104 if (!m_weightAlignment || TMath::Prob(redchi2 * 8, 8) < 0.99) {
0105 if (fabs(residual - resid_mean) < 2.5 * resid_stdev) {
0106 resid_sum += weight * residual;
0107 resid_sum2 += weight * residual * residual;
0108 resid_N += weight;
0109 }
0110 }
0111 }
0112
0113 resid_mean = resid_sum / resid_N;
0114 resid_stdev = sqrt(resid_sum2 / resid_N - pow(resid_sum / resid_N, 2));
0115
0116 std::vector<int> num;
0117 std::vector<std::string> name;
0118 std::vector<double> start;
0119 std::vector<double> step;
0120 std::vector<double> low;
0121 std::vector<double> high;
0122
0123 if (fixed(kAlign)) {
0124 num.push_back(kAlign);
0125 name.push_back(std::string("Align"));
0126 start.push_back(0.);
0127 step.push_back(0.01 * resid_stdev);
0128 low.push_back(0.);
0129 high.push_back(0.);
0130 } else {
0131 num.push_back(kAlign);
0132 name.push_back(std::string("Align"));
0133 start.push_back(resid_mean);
0134 step.push_back(0.01 * resid_stdev);
0135 low.push_back(0.);
0136 high.push_back(0.);
0137 }
0138 num.push_back(kSigma);
0139 name.push_back(std::string("Sigma"));
0140 start.push_back(resid_stdev);
0141 step.push_back(0.01 * resid_stdev);
0142 low.push_back(0.);
0143 high.push_back(0.);
0144 if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
0145 num.push_back(kGamma);
0146 name.push_back(std::string("Gamma"));
0147 start.push_back(0.1 * resid_stdev);
0148 step.push_back(0.01 * resid_stdev);
0149 low.push_back(0.);
0150 high.push_back(0.);
0151 }
0152
0153 return dofit(&MuonResiduals1DOFFitter_FCN, num, name, start, step, low, high);
0154 }
0155
0156 double MuonResiduals1DOFFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
0157 sumofweights();
0158
0159 std::stringstream name_residual, name_residual_raw;
0160 name_residual << name << "_residual";
0161 name_residual_raw << name << "_residual_raw";
0162
0163 double min_residual = -100.;
0164 double max_residual = 100.;
0165 TH1F *hist_residual = dir->make<TH1F>(name_residual.str().c_str(), "", 100, min_residual, max_residual);
0166 TH1F *hist_residual_raw = dir->make<TH1F>(name_residual_raw.str().c_str(), "", 100, min_residual, max_residual);
0167
0168 name_residual << "_fit";
0169 TF1 *fit_residual = nullptr;
0170 if (residualsModel() == kPureGaussian) {
0171 fit_residual =
0172 new TF1(name_residual.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, min_residual, max_residual, 3);
0173 fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual) / 100.,
0174 10. * value(kAlign),
0175 10. * value(kSigma));
0176 } else if (residualsModel() == kPowerLawTails) {
0177 fit_residual =
0178 new TF1(name_residual.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, min_residual, max_residual, 4);
0179 fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual) / 100.,
0180 10. * value(kAlign),
0181 10. * value(kSigma),
0182 10. * value(kGamma));
0183 } else if (residualsModel() == kROOTVoigt) {
0184 fit_residual =
0185 new TF1(name_residual.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_residual, max_residual, 4);
0186 fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual) / 100.,
0187 10. * value(kAlign),
0188 10. * value(kSigma),
0189 10. * value(kGamma));
0190 } else if (residualsModel() == kGaussPowerTails) {
0191 fit_residual =
0192 new TF1(name_residual.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_residual, max_residual, 3);
0193 fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual) / 100.,
0194 10. * value(kAlign),
0195 10. * value(kSigma));
0196 } else {
0197 assert(false);
0198 }
0199
0200 fit_residual->SetLineColor(2);
0201 fit_residual->SetLineWidth(2);
0202 fit_residual->Write();
0203
0204 for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0205 const double resid = (*resiter)[MuonResiduals1DOFFitter::kResid];
0206 const double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
0207 double weight = 1. / redchi2;
0208 if (!m_weightAlignment)
0209 weight = 1.;
0210
0211 if (!m_weightAlignment || TMath::Prob(redchi2 * 8, 8) < 0.99) {
0212 hist_residual->Fill(10. * (resid + value(kAlign)), weight);
0213 }
0214
0215 hist_residual_raw->Fill(10. * resid);
0216 }
0217
0218 double chi2 = 0.;
0219 double ndof = 0.;
0220 for (int i = 1; i <= hist_residual->GetNbinsX(); i++) {
0221 double xi = hist_residual->GetBinCenter(i);
0222 double yi = hist_residual->GetBinContent(i);
0223 double yerri = hist_residual->GetBinError(i);
0224 double yth = fit_residual->Eval(xi);
0225 if (yerri > 0.) {
0226 chi2 += pow((yth - yi) / yerri, 2);
0227 ndof += 1.;
0228 }
0229 }
0230 ndof -= npar();
0231
0232 return (ndof > 0. ? chi2 / ndof : -1.);
0233 }