MuonAlignmentAlgorithms/src/MuonResidualsPositionFitter.cc

0001 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsPositionFitter.h"
0002
0003 static TMinuit *MuonResidualsPositionFitter_TMinuit;
0004
0005 void MuonResidualsPositionFitter::inform(TMinuit *tMinuit) { MuonResidualsPositionFitter_TMinuit = tMinuit; }
0006
0007 void MuonResidualsPositionFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
0008   MuonResidualsFitterFitInfo *fitinfo =
0009       (MuonResidualsFitterFitInfo *)(MuonResidualsPositionFitter_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)[MuonResidualsPositionFitter::kResidual];
0016     const double angleerror = (*resiter)[MuonResidualsPositionFitter::kAngleError];
0017     const double trackangle = (*resiter)[MuonResidualsPositionFitter::kTrackAngle];
0018     const double trackposition = (*resiter)[MuonResidualsPositionFitter::kTrackPosition];
0019
0020     double center = 0.;
0021     center += par[MuonResidualsPositionFitter::kPosition];
0022     center += par[MuonResidualsPositionFitter::kZpos] * trackangle;
0023     center += par[MuonResidualsPositionFitter::kPhiz] * trackposition;
0024     center += par[MuonResidualsPositionFitter::kScattering] * angleerror;
0025
0026     if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
0027       fval += -MuonResidualsFitter_logPureGaussian(residual, center, par[MuonResidualsPositionFitter::kSigma]);
0028     } else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
0029       fval += -MuonResidualsFitter_logPowerLawTails(
0030           residual, center, par[MuonResidualsPositionFitter::kSigma], par[MuonResidualsPositionFitter::kGamma]);
0031     } else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
0032       fval += -MuonResidualsFitter_logROOTVoigt(
0033           residual, center, par[MuonResidualsPositionFitter::kSigma], par[MuonResidualsPositionFitter::kGamma]);
0034     } else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
0035       fval += -MuonResidualsFitter_logGaussPowerTails(residual, center, par[MuonResidualsPositionFitter::kSigma]);
0036     } else {
0037       assert(false);
0038     }
0039   }
0040 }
0041
0042 bool MuonResidualsPositionFitter::fit(Alignable *ali) {
0043   initialize_table();  // if not already initialized
0044
0045   double sum_x = 0.;
0046   double sum_xx = 0.;
0047   int N = 0;
0048
0049   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0050     const double residual = (*resiter)[kResidual];
0051     //    const double angleerror = (*resiter)[kAngleError];
0052     //    const double trackangle = (*resiter)[kTrackAngle];
0053     //    const double trackposition = (*resiter)[kTrackPosition];
0054
0055     if (fabs(residual) < 10.) {  // truncate at 100 mm
0056       sum_x += residual;
0057       sum_xx += residual * residual;
0058       N++;
0059     }
0060   }
0061
0062   if (N < m_minHits)
0063     return false;
0064
0065   // truncated mean and stdev to seed the fit
0066   double mean = sum_x / double(N);
0067   double stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0068
0069   // refine the standard deviation calculation
0070   sum_x = 0.;
0071   sum_xx = 0.;
0072   N = 0;
0073   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0074     const double residual = (*resiter)[kResidual];
0075     if (mean - 1.5 * stdev < residual && residual < mean + 1.5 * stdev) {
0076       sum_x += residual;
0077       sum_xx += residual * residual;
0078       N++;
0079     }
0080   }
0081   mean = sum_x / double(N);
0082   stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0083
0084   sum_x = 0.;
0085   sum_xx = 0.;
0086   N = 0;
0087   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0088     const double residual = (*resiter)[kResidual];
0089     if (mean - 1.5 * stdev < residual && residual < mean + 1.5 * stdev) {
0090       sum_x += residual;
0091       sum_xx += residual * residual;
0092       N++;
0093     }
0094   }
0095   mean = sum_x / double(N);
0096   stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0097
0098   std::vector<int> parNum;
0099   std::vector<std::string> parName;
0100   std::vector<double> start;
0101   std::vector<double> step;
0102   std::vector<double> low;
0103   std::vector<double> high;
0104
0105   parNum.push_back(kPosition);
0106   parName.push_back(std::string("position"));
0107   start.push_back(mean);
0108   step.push_back(0.1);
0109   low.push_back(0.);
0110   high.push_back(0.);
0111   parNum.push_back(kZpos);
0112   parName.push_back(std::string("zpos"));
0113   start.push_back(0.);
0114   step.push_back(0.1);
0115   low.push_back(0.);
0116   high.push_back(0.);
0117   parNum.push_back(kPhiz);
0118   parName.push_back(std::string("phiz"));
0119   start.push_back(0.);
0120   step.push_back(0.1);
0121   low.push_back(0.);
0122   high.push_back(0.);
0123   parNum.push_back(kScattering);
0124   parName.push_back(std::string("scattering"));
0125   start.push_back(0.);
0126   step.push_back(0.1 * 1000.);
0127   low.push_back(0.);
0128   high.push_back(0.);
0129   parNum.push_back(kSigma);
0130   parName.push_back(std::string("sigma"));
0131   start.push_back(stdev);
0132   step.push_back(0.1 * stdev);
0133   low.push_back(0.);
0134   high.push_back(0.);
0135   if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
0136     parNum.push_back(kGamma);
0137     parName.push_back(std::string("gamma"));
0138     start.push_back(stdev);
0139     step.push_back(0.1 * stdev);
0140     low.push_back(0.);
0141     high.push_back(0.);
0142   }
0143
0144   return dofit(&MuonResidualsPositionFitter_FCN, parNum, parName, start, step, low, high);
0145 }
0146
0147 double MuonResidualsPositionFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
0148   std::stringstream raw_name, narrowed_name, angleerror_name, trackangle_name, trackposition_name;
0149   raw_name << name << "_raw";
0150   narrowed_name << name << "_narrowed";
0151   angleerror_name << name << "_angleerror";
0152   trackangle_name << name << "_trackangle";
0153   trackposition_name << name << "_trackposition";
0154
0155   TH1F *raw_hist =
0156       dir->make<TH1F>(raw_name.str().c_str(), (raw_name.str() + std::string(" (mm)")).c_str(), 100, -100., 100.);
0157   TH1F *narrowed_hist = dir->make<TH1F>(
0158       narrowed_name.str().c_str(), (narrowed_name.str() + std::string(" (mm)")).c_str(), 100, -100., 100.);
0159   TProfile *angleerror_hist = dir->make<TProfile>(
0160       angleerror_name.str().c_str(), (angleerror_name.str() + std::string(" (mm)")).c_str(), 100, -30., 30.);
0161   TProfile *trackangle_hist = dir->make<TProfile>(
0162       trackangle_name.str().c_str(), (trackangle_name.str() + std::string(" (mm)")).c_str(), 100, -0.5, 0.5);
0163   TProfile *trackposition_hist = dir->make<TProfile>(
0164       trackposition_name.str().c_str(), (trackposition_name.str() + std::string(" (mm)")).c_str(), 100, -300., 300.);
0165
0166   angleerror_hist->SetAxisRange(-100., 100., "Y");
0167   trackangle_hist->SetAxisRange(-10., 10., "Y");
0168   trackposition_hist->SetAxisRange(-10., 10., "Y");
0169
0170   narrowed_name << "fit";
0171   angleerror_name << "fit";
0172   trackangle_name << "fit";
0173   trackposition_name << "fit";
0174
0175   double scale_factor = double(numResiduals()) * (100. - -100.) / 100;  // (max - min)/nbins
0176
0177   TF1 *narrowed_fit = nullptr;
0178   if (residualsModel() == kPureGaussian) {
0179     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, -100., 100., 3);
0180     narrowed_fit->SetParameters(scale_factor, value(kPosition) * 10., value(kSigma) * 10.);
0181     narrowed_fit->Write();
0182   } else if (residualsModel() == kPowerLawTails) {
0183     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, -100., 100., 4);
0184     narrowed_fit->SetParameters(scale_factor, value(kPosition) * 10., value(kSigma) * 10., value(kGamma) * 10.);
0185     narrowed_fit->Write();
0186   } else if (residualsModel() == kROOTVoigt) {
0187     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, -100., 100., 4);
0188     narrowed_fit->SetParameters(scale_factor, value(kPosition) * 10., value(kSigma) * 10., value(kGamma) * 10.);
0189     narrowed_fit->Write();
0190   } else if (residualsModel() == kGaussPowerTails) {
0191     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, -100., 100., 3);
0192     narrowed_fit->SetParameters(scale_factor, value(kPosition) * 10., value(kSigma) * 10.);
0193     narrowed_fit->Write();
0194   }
0195
0196   TF1 *angleerror_fit = new TF1(angleerror_name.str().c_str(), "[0]+x*[1]", -30., 30.);
0197   angleerror_fit->SetParameters(value(kPosition) * 10., value(kScattering) * 10. / 1000.);
0198   angleerror_fit->Write();
0199
0200   TF1 *trackangle_fit = new TF1(trackangle_name.str().c_str(), "[0]+x*[1]", -0.5, 0.5);
0201   trackangle_fit->SetParameters(value(kPosition) * 10., value(kZpos) * 10.);
0202   trackangle_fit->Write();
0203
0204   TF1 *trackposition_fit = new TF1(trackposition_name.str().c_str(), "[0]+x*[1]", -300., 300.);
0205   trackposition_fit->SetParameters(value(kPosition) * 10., value(kPhiz) * 10.);
0206   trackposition_fit->Write();
0207
0208   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0209     const double raw_residual = (*resiter)[kResidual];
0210     const double angleerror = (*resiter)[kAngleError];
0211     const double trackangle = (*resiter)[kTrackAngle];
0212     const double trackposition = (*resiter)[kTrackPosition];
0213
0214     double angleerror_correction = value(kScattering) * angleerror;
0215     double trackangle_correction = value(kZpos) * trackangle;
0216     double trackposition_correction = value(kPhiz) * trackposition;
0217
0218     double corrected_residual = raw_residual - angleerror_correction - trackangle_correction - trackposition_correction;
0219
0220     raw_hist->Fill(raw_residual * 10.);
0221     narrowed_hist->Fill(corrected_residual * 10.);
0222
0223     angleerror_hist->Fill(angleerror * 1000., (raw_residual - trackangle_correction - trackposition_correction) * 10.);
0224     trackangle_hist->Fill(trackangle, (raw_residual - angleerror_correction - trackposition_correction) * 10.);
0225     trackposition_hist->Fill(trackposition, (raw_residual - angleerror_correction - trackangle_correction) * 10.);
0226   }
0227
0228   return 0.;
0229 }