MuonAlignmentAlgorithms/src/MuonResidualsAngleFitter.cc

0001 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsAngleFitter.h"
0002
0003 static TMinuit *MuonResidualsAngleFitter_TMinuit;
0004
0005 void MuonResidualsAngleFitter::inform(TMinuit *tMinuit) { MuonResidualsAngleFitter_TMinuit = tMinuit; }
0006
0007 void MuonResidualsAngleFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
0008   MuonResidualsFitterFitInfo *fitinfo =
0009       (MuonResidualsFitterFitInfo *)(MuonResidualsAngleFitter_TMinuit->GetObjectFit());
0010   MuonResidualsFitter *fitter = fitinfo->fitter();
0011
0012   fval = 0.;
0013   for (std::vector<double *>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end();
0014        ++resiter) {
0015     const double residual = (*resiter)[MuonResidualsAngleFitter::kResidual];
0016     const double xangle = (*resiter)[MuonResidualsAngleFitter::kXAngle];
0017     const double yangle = (*resiter)[MuonResidualsAngleFitter::kYAngle];
0018
0019     double center = 0.;
0020     center += par[MuonResidualsAngleFitter::kAngle];
0021     center += par[MuonResidualsAngleFitter::kXControl] * xangle;
0022     center += par[MuonResidualsAngleFitter::kYControl] * yangle;
0023
0024     if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
0025       fval += -MuonResidualsFitter_logPureGaussian(residual, center, par[MuonResidualsAngleFitter::kSigma]);
0026     } else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
0027       fval += -MuonResidualsFitter_logPowerLawTails(
0028           residual, center, par[MuonResidualsAngleFitter::kSigma], par[MuonResidualsAngleFitter::kGamma]);
0029     } else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
0030       fval += -MuonResidualsFitter_logROOTVoigt(
0031           residual, center, par[MuonResidualsAngleFitter::kSigma], par[MuonResidualsAngleFitter::kGamma]);
0032     } else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
0033       fval += -MuonResidualsFitter_logGaussPowerTails(residual, center, par[MuonResidualsAngleFitter::kSigma]);
0034     } else {
0035       assert(false);
0036     }
0037   }
0038 }
0039
0040 bool MuonResidualsAngleFitter::fit(Alignable *ali) {
0041   initialize_table();  // if not already initialized
0042
0043   double sum_x = 0.;
0044   double sum_xx = 0.;
0045   int N = 0;
0046
0047   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0048     const double residual = (*resiter)[kResidual];
0049     //     const double xangle = (*resiter)[kXAngle];
0050     //     const double yangle = (*resiter)[kYAngle];
0051
0052     if (fabs(residual) < 0.1) {  // truncate at 100 mrad
0053       sum_x += residual;
0054       sum_xx += residual * residual;
0055       N++;
0056     }
0057   }
0058
0059   if (N < m_minHits)
0060     return false;
0061
0062   // truncated mean and stdev to seed the fit
0063   double mean = sum_x / double(N);
0064   double stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0065
0066   // refine the standard deviation calculation
0067   sum_x = 0.;
0068   sum_xx = 0.;
0069   N = 0;
0070   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0071     const double residual = (*resiter)[kResidual];
0072     if (mean - 1.5 * stdev < residual && residual < mean + 1.5 * stdev) {
0073       sum_x += residual;
0074       sum_xx += residual * residual;
0075       N++;
0076     }
0077   }
0078   mean = sum_x / double(N);
0079   stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0080
0081   sum_x = 0.;
0082   sum_xx = 0.;
0083   N = 0;
0084   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0085     const double residual = (*resiter)[kResidual];
0086     if (mean - 1.5 * stdev < residual && residual < mean + 1.5 * stdev) {
0087       sum_x += residual;
0088       sum_xx += residual * residual;
0089       N++;
0090     }
0091   }
0092   mean = sum_x / double(N);
0093   stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0094
0095   std::vector<int> parNum;
0096   std::vector<std::string> parName;
0097   std::vector<double> start;
0098   std::vector<double> step;
0099   std::vector<double> low;
0100   std::vector<double> high;
0101
0102   parNum.push_back(kAngle);
0103   parName.push_back(std::string("angle"));
0104   start.push_back(mean);
0105   step.push_back(0.1);
0106   low.push_back(0.);
0107   high.push_back(0.);
0108   parNum.push_back(kXControl);
0109   parName.push_back(std::string("xcontrol"));
0110   start.push_back(0.);
0111   step.push_back(0.1);
0112   low.push_back(0.);
0113   high.push_back(0.);
0114   parNum.push_back(kYControl);
0115   parName.push_back(std::string("ycontrol"));
0116   start.push_back(0.);
0117   step.push_back(0.1);
0118   low.push_back(0.);
0119   high.push_back(0.);
0120   parNum.push_back(kSigma);
0121   parName.push_back(std::string("sigma"));
0122   start.push_back(stdev);
0123   step.push_back(0.1 * stdev);
0124   low.push_back(0.);
0125   high.push_back(0.);
0126   if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
0127     parNum.push_back(kGamma);
0128     parName.push_back(std::string("gamma"));
0129     start.push_back(stdev);
0130     step.push_back(0.1 * stdev);
0131     low.push_back(0.);
0132     high.push_back(0.);
0133   }
0134
0135   return dofit(&MuonResidualsAngleFitter_FCN, parNum, parName, start, step, low, high);
0136 }
0137
0138 double MuonResidualsAngleFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
0139   std::stringstream raw_name, narrowed_name, xcontrol_name, ycontrol_name;
0140   raw_name << name << "_raw";
0141   narrowed_name << name << "_narrowed";
0142   xcontrol_name << name << "_xcontrol";
0143   ycontrol_name << name << "_ycontrol";
0144
0145   TH1F *raw_hist =
0146       dir->make<TH1F>(raw_name.str().c_str(), (raw_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
0147   TH1F *narrowed_hist = dir->make<TH1F>(
0148       narrowed_name.str().c_str(), (narrowed_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
0149   TProfile *xcontrol_hist = dir->make<TProfile>(
0150       xcontrol_name.str().c_str(), (xcontrol_name.str() + std::string(" (mrad)")).c_str(), 100, -1., 1.);
0151   TProfile *ycontrol_hist = dir->make<TProfile>(
0152       ycontrol_name.str().c_str(), (ycontrol_name.str() + std::string(" (mrad)")).c_str(), 100, -1., 1.);
0153
0154   narrowed_name << "fit";
0155   xcontrol_name << "fit";
0156   ycontrol_name << "fit";
0157
0158   double scale_factor = double(numResiduals()) * (100. - -100.) / 100;  // (max - min)/nbins
0159
0160   TF1 *narrowed_fit = nullptr;
0161   if (residualsModel() == kPureGaussian) {
0162     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, -100., 100., 3);
0163     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
0164     narrowed_fit->Write();
0165   } else if (residualsModel() == kPowerLawTails) {
0166     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, -100., 100., 4);
0167     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
0168     narrowed_fit->Write();
0169   } else if (residualsModel() == kROOTVoigt) {
0170     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, -100., 100., 4);
0171     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
0172     narrowed_fit->Write();
0173   } else if (residualsModel() == kGaussPowerTails) {
0174     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, -100., 100., 3);
0175     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
0176     narrowed_fit->Write();
0177   }
0178
0179   TF1 *xcontrol_fit = new TF1(xcontrol_name.str().c_str(), "[0]+x*[1]", -1., 1.);
0180   xcontrol_fit->SetParameters(value(kAngle) * 1000., value(kXControl) * 1000.);
0181   xcontrol_fit->Write();
0182
0183   TF1 *ycontrol_fit = new TF1(ycontrol_name.str().c_str(), "[0]+x*[1]", -1., 1.);
0184   ycontrol_fit->SetParameters(value(kAngle) * 1000., value(kYControl) * 1000.);
0185   ycontrol_fit->Write();
0186
0187   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0188     const double raw_residual = (*resiter)[kResidual];
0189     const double xangle = (*resiter)[kXAngle];
0190     const double yangle = (*resiter)[kYAngle];
0191
0192     double xangle_correction = value(kXControl) * xangle;
0193     double yangle_correction = value(kYControl) * yangle;
0194     double corrected_residual = raw_residual - xangle_correction - yangle_correction;
0195
0196     raw_hist->Fill(raw_residual * 1000.);
0197     narrowed_hist->Fill(corrected_residual * 1000.);
0198
0199     xcontrol_hist->Fill(xangle, (raw_residual - yangle_correction) * 1000.);
0200     ycontrol_hist->Fill(yangle, (raw_residual - xangle_correction) * 1000.);
0201   }
0202
0203   return 0.;
0204 }