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
0072   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0073     const double residual = (*resiter)[MuonResiduals1DOFFitter::kResid];
0074     const double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
0075     double weight = 1. / redchi2;
0076     if (!m_weightAlignment)
0077       weight = 1.;
0078
0079     if (!m_weightAlignment || TMath::Prob(redchi2 * 8, 8) < 0.99) {  // no spikes allowed
0080       if (fabs(residual) < 10.) {                                    // 10 cm
0081         resid_sum += weight * residual;
0082         resid_sum2 += weight * residual * residual;
0083         resid_N += weight;
0084       }
0085     }
0086   }
0087
0088   double resid_mean = resid_sum / resid_N;
0089   double resid_stdev = sqrt(resid_sum2 / resid_N - pow(resid_sum / resid_N, 2));
0090
0091   resid_sum = 0.;
0092   resid_sum2 = 0.;
0093   resid_N = 0.;
0094
0095   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0096     const double residual = (*resiter)[MuonResiduals1DOFFitter::kResid];
0097     const double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
0098     double weight = 1. / redchi2;
0099     if (!m_weightAlignment)
0100       weight = 1.;
0101
0102     if (!m_weightAlignment || TMath::Prob(redchi2 * 8, 8) < 0.99) {  // no spikes allowed
0103       if (fabs(residual - resid_mean) < 2.5 * resid_stdev) {
0104         resid_sum += weight * residual;
0105         resid_sum2 += weight * residual * residual;
0106         resid_N += weight;
0107       }
0108     }
0109   }
0110
0111   resid_mean = resid_sum / resid_N;
0112   resid_stdev = sqrt(resid_sum2 / resid_N - pow(resid_sum / resid_N, 2));
0113
0114   std::vector<int> num;
0115   std::vector<std::string> name;
0116   std::vector<double> start;
0117   std::vector<double> step;
0118   std::vector<double> low;
0119   std::vector<double> high;
0120
0121   if (fixed(kAlign)) {
0122     num.push_back(kAlign);
0123     name.push_back(std::string("Align"));
0124     start.push_back(0.);
0125     step.push_back(0.01 * resid_stdev);
0126     low.push_back(0.);
0127     high.push_back(0.);
0128   } else {
0129     num.push_back(kAlign);
0130     name.push_back(std::string("Align"));
0131     start.push_back(resid_mean);
0132     step.push_back(0.01 * resid_stdev);
0133     low.push_back(0.);
0134     high.push_back(0.);
0135   }
0136   num.push_back(kSigma);
0137   name.push_back(std::string("Sigma"));
0138   start.push_back(resid_stdev);
0139   step.push_back(0.01 * resid_stdev);
0140   low.push_back(0.);
0141   high.push_back(0.);
0142   if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
0143     num.push_back(kGamma);
0144     name.push_back(std::string("Gamma"));
0145     start.push_back(0.1 * resid_stdev);
0146     step.push_back(0.01 * resid_stdev);
0147     low.push_back(0.);
0148     high.push_back(0.);
0149   }
0150
0151   return dofit(&MuonResiduals1DOFFitter_FCN, num, name, start, step, low, high);
0152 }
0153
0154 double MuonResiduals1DOFFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
0155   sumofweights();
0156
0157   std::stringstream name_residual, name_residual_raw;
0158   name_residual << name << "_residual";
0159   name_residual_raw << name << "_residual_raw";
0160
0161   double min_residual = -100.;
0162   double max_residual = 100.;
0163   TH1F *hist_residual = dir->make<TH1F>(name_residual.str().c_str(), "", 100, min_residual, max_residual);
0164   TH1F *hist_residual_raw = dir->make<TH1F>(name_residual_raw.str().c_str(), "", 100, min_residual, max_residual);
0165
0166   name_residual << "_fit";
0167   TF1 *fit_residual = nullptr;
0168   if (residualsModel() == kPureGaussian) {
0169     fit_residual =
0170         new TF1(name_residual.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, min_residual, max_residual, 3);
0171     fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual) / 100.,
0172                                 10. * value(kAlign),
0173                                 10. * value(kSigma));
0174   } else if (residualsModel() == kPowerLawTails) {
0175     fit_residual =
0176         new TF1(name_residual.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, min_residual, max_residual, 4);
0177     fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual) / 100.,
0178                                 10. * value(kAlign),
0179                                 10. * value(kSigma),
0180                                 10. * value(kGamma));
0181   } else if (residualsModel() == kROOTVoigt) {
0182     fit_residual =
0183         new TF1(name_residual.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_residual, max_residual, 4);
0184     fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual) / 100.,
0185                                 10. * value(kAlign),
0186                                 10. * value(kSigma),
0187                                 10. * value(kGamma));
0188   } else if (residualsModel() == kGaussPowerTails) {
0189     fit_residual =
0190         new TF1(name_residual.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_residual, max_residual, 3);
0191     fit_residual->SetParameters(MuonResiduals1DOFFitter_sum_of_weights * (max_residual - min_residual) / 100.,
0192                                 10. * value(kAlign),
0193                                 10. * value(kSigma));
0194   } else {
0195     assert(false);
0196   }
0197
0198   fit_residual->SetLineColor(2);
0199   fit_residual->SetLineWidth(2);
0200   fit_residual->Write();
0201
0202   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0203     const double resid = (*resiter)[MuonResiduals1DOFFitter::kResid];
0204     const double redchi2 = (*resiter)[MuonResiduals1DOFFitter::kRedChi2];
0205     double weight = 1. / redchi2;
0206     if (!m_weightAlignment)
0207       weight = 1.;
0208
0209     if (!m_weightAlignment || TMath::Prob(redchi2 * 8, 8) < 0.99) {  // no spikes allowed
0210       hist_residual->Fill(10. * (resid + value(kAlign)), weight);
0211     }
0212
0213     hist_residual_raw->Fill(10. * resid);
0214   }
0215
0216   double chi2 = 0.;
0217   double ndof = 0.;
0218   for (int i = 1; i <= hist_residual->GetNbinsX(); i++) {
0219     double xi = hist_residual->GetBinCenter(i);
0220     double yi = hist_residual->GetBinContent(i);
0221     double yerri = hist_residual->GetBinError(i);
0222     double yth = fit_residual->Eval(xi);
0223     if (yerri > 0.) {
0224       chi2 += pow((yth - yi) / yerri, 2);
0225       ndof += 1.;
0226     }
0227   }
0228   ndof -= npar();
0229
0230   return (ndof > 0. ? chi2 / ndof : -1.);
0231 }