MuonAlignmentAlgorithms/src/MuonResiduals1DOFFitter.cc

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) {  // no spikes allowed
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) {  // no spikes allowed
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();  // if not already initialized
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) {  // no spikes allowed
0081       if (fabs(residual) < 10.) {                                    // 10 cm
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) {  // no spikes allowed
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) {  // no spikes allowed
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 }