MuonAlignmentAlgorithms/src/MuonResidualsBfieldAngleFitter.cc

0001 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsBfieldAngleFitter.h"
0002
0003 static TMinuit *MuonResidualsBfieldAngleFitter_TMinuit;
0004
0005 void MuonResidualsBfieldAngleFitter::inform(TMinuit *tMinuit) { MuonResidualsBfieldAngleFitter_TMinuit = tMinuit; }
0006
0007 void MuonResidualsBfieldAngleFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
0008   MuonResidualsFitterFitInfo *fitinfo =
0009       (MuonResidualsFitterFitInfo *)(MuonResidualsBfieldAngleFitter_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)[MuonResidualsBfieldAngleFitter::kResidual];
0016     const double qoverpt = (*resiter)[MuonResidualsBfieldAngleFitter::kQoverPt];
0017     const double qoverpz = (*resiter)[MuonResidualsBfieldAngleFitter::kQoverPz];
0018
0019     double center = 0.;
0020     center += par[MuonResidualsBfieldAngleFitter::kAngle];
0021     center += par[MuonResidualsBfieldAngleFitter::kBfrompt] * qoverpt;
0022     center += par[MuonResidualsBfieldAngleFitter::kBfrompz] * qoverpz;
0023     center += par[MuonResidualsBfieldAngleFitter::kdEdx] * (1. / qoverpt / qoverpt + 1. / qoverpz / qoverpz) *
0024               (qoverpt > 0. ? 1. : -1.);
0025
0026     if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
0027       fval += -MuonResidualsFitter_logPureGaussian(residual, center, par[MuonResidualsBfieldAngleFitter::kSigma]);
0028     } else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
0029       fval += -MuonResidualsFitter_logPowerLawTails(
0030           residual, center, par[MuonResidualsBfieldAngleFitter::kSigma], par[MuonResidualsBfieldAngleFitter::kGamma]);
0031     } else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
0032       fval += -MuonResidualsFitter_logROOTVoigt(
0033           residual, center, par[MuonResidualsBfieldAngleFitter::kSigma], par[MuonResidualsBfieldAngleFitter::kGamma]);
0034     } else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
0035       fval += -MuonResidualsFitter_logGaussPowerTails(residual, center, par[MuonResidualsBfieldAngleFitter::kSigma]);
0036     } else {
0037       assert(false);
0038     }
0039   }
0040 }
0041
0042 bool MuonResidualsBfieldAngleFitter::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 qoverpt = (*resiter)[kQoverPt];
0052
0053     if (fabs(residual) < 0.1) {  // truncate at 100 mrad
0054       sum_x += residual;
0055       sum_xx += residual * residual;
0056       N++;
0057     }
0058   }
0059
0060   if (N < m_minHits)
0061     return false;
0062
0063   // truncated mean and stdev to seed the fit
0064   double mean = sum_x / double(N);
0065   double stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0066
0067   // refine the standard deviation calculation
0068   sum_x = 0.;
0069   sum_xx = 0.;
0070   N = 0;
0071   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0072     const double residual = (*resiter)[kResidual];
0073     if (mean - 1.5 * stdev < residual && residual < mean + 1.5 * stdev) {
0074       sum_x += residual;
0075       sum_xx += residual * residual;
0076       N++;
0077     }
0078   }
0079   mean = sum_x / double(N);
0080   stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0081
0082   sum_x = 0.;
0083   sum_xx = 0.;
0084   N = 0;
0085   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0086     const double residual = (*resiter)[kResidual];
0087     if (mean - 1.5 * stdev < residual && residual < mean + 1.5 * stdev) {
0088       sum_x += residual;
0089       sum_xx += residual * residual;
0090       N++;
0091     }
0092   }
0093   mean = sum_x / double(N);
0094   stdev = sqrt(sum_xx / double(N) - pow(sum_x / double(N), 2));
0095
0096   std::vector<int> parNum;
0097   std::vector<std::string> parName;
0098   std::vector<double> start;
0099   std::vector<double> step;
0100   std::vector<double> low;
0101   std::vector<double> high;
0102
0103   parNum.push_back(kAngle);
0104   parName.push_back(std::string("angle"));
0105   start.push_back(mean);
0106   step.push_back(0.1);
0107   low.push_back(0.);
0108   high.push_back(0.);
0109   parNum.push_back(kBfrompt);
0110   parName.push_back(std::string("bfrompt"));
0111   start.push_back(0.);
0112   step.push_back(0.1 * stdev / 0.05);
0113   low.push_back(0.);
0114   high.push_back(0.);
0115   parNum.push_back(kBfrompz);
0116   parName.push_back(std::string("bfrompz"));
0117   start.push_back(0.);
0118   step.push_back(0.1 * stdev / 0.05);
0119   low.push_back(0.);
0120   high.push_back(0.);
0121   parNum.push_back(kdEdx);
0122   parName.push_back(std::string("dEdx"));
0123   start.push_back(0.);
0124   step.push_back(0.1 * stdev / 0.05);
0125   low.push_back(0.);
0126   high.push_back(0.);
0127   parNum.push_back(kSigma);
0128   parName.push_back(std::string("sigma"));
0129   start.push_back(stdev);
0130   step.push_back(0.1 * stdev);
0131   low.push_back(0.);
0132   high.push_back(0.);
0133   if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
0134     parNum.push_back(kGamma);
0135     parName.push_back(std::string("gamma"));
0136     start.push_back(stdev);
0137     step.push_back(0.1 * stdev);
0138     low.push_back(0.);
0139     high.push_back(0.);
0140   }
0141
0142   return dofit(&MuonResidualsBfieldAngleFitter_FCN, parNum, parName, start, step, low, high);
0143 }
0144
0145 double MuonResidualsBfieldAngleFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
0146   std::stringstream raw_name, narrowed_name, qoverpt_name, qoverpz_name, psquared_name;
0147   raw_name << name << "_raw";
0148   narrowed_name << name << "_narrowed";
0149   qoverpt_name << name << "_qoverpt";
0150   qoverpz_name << name << "_qoverpz";
0151   psquared_name << name << "_psquared";
0152
0153   TH1F *raw_hist =
0154       dir->make<TH1F>(raw_name.str().c_str(), (raw_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
0155   TH1F *narrowed_hist = dir->make<TH1F>(
0156       narrowed_name.str().c_str(), (narrowed_name.str() + std::string(" (mrad)")).c_str(), 100, -100., 100.);
0157   TProfile *qoverpt_hist = dir->make<TProfile>(
0158       qoverpt_name.str().c_str(), (qoverpt_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
0159   TProfile *qoverpz_hist = dir->make<TProfile>(
0160       qoverpz_name.str().c_str(), (qoverpz_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
0161   TProfile *psquared_hist = dir->make<TProfile>(
0162       psquared_name.str().c_str(), (psquared_name.str() + std::string(" (mrad)")).c_str(), 100, -0.05, 0.05);
0163
0164   narrowed_name << "fit";
0165   qoverpt_name << "fit";
0166   qoverpz_name << "fit";
0167   psquared_name << "fit";
0168
0169   double scale_factor = double(numResiduals()) * (100. - -100.) / 100;  // (max - min)/nbins
0170
0171   TF1 *narrowed_fit = nullptr;
0172   if (residualsModel() == kPureGaussian) {
0173     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, -100., 100., 3);
0174     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
0175     narrowed_fit->Write();
0176   } else if (residualsModel() == kPowerLawTails) {
0177     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, -100., 100., 4);
0178     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
0179     narrowed_fit->Write();
0180   } else if (residualsModel() == kROOTVoigt) {
0181     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, -100., 100., 4);
0182     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000., value(kGamma) * 1000.);
0183     narrowed_fit->Write();
0184   } else if (residualsModel() == kGaussPowerTails) {
0185     narrowed_fit = new TF1(narrowed_name.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, -100., 100., 3);
0186     narrowed_fit->SetParameters(scale_factor, value(kAngle) * 1000., value(kSigma) * 1000.);
0187     narrowed_fit->Write();
0188   }
0189
0190   TF1 *qoverpt_fit = new TF1(qoverpt_name.str().c_str(), "[0]+x*[1]", -0.05, 0.05);
0191   qoverpt_fit->SetParameters(value(kAngle) * 1000., value(kBfrompt) * 1000.);
0192   qoverpt_fit->Write();
0193
0194   TF1 *qoverpz_fit = new TF1(qoverpz_name.str().c_str(), "[0]+x*[1]", -0.05, 0.05);
0195   qoverpz_fit->SetParameters(value(kAngle) * 1000., value(kBfrompz) * 1000.);
0196   qoverpz_fit->Write();
0197
0198   TF1 *psquared_fit = new TF1(psquared_name.str().c_str(), "[0]+[1]*x**2", -0.05, 0.05);
0199   psquared_fit->SetParameters(value(kAngle) * 1000., value(kdEdx) * 1000.);
0200   psquared_fit->Write();
0201
0202   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0203     const double raw_residual = (*resiter)[kResidual];
0204     const double qoverpt = (*resiter)[kQoverPt];
0205     const double qoverpz = (*resiter)[kQoverPz];
0206     const double psquared = (1. / qoverpt / qoverpt + 1. / qoverpz / qoverpz) * (qoverpt > 0. ? 1. : -1.);
0207
0208     double qoverpt_correction = value(kBfrompt) * qoverpt;
0209     double qoverpz_correction = value(kBfrompz) * qoverpz;
0210     double dEdx_correction = value(kdEdx) * psquared;
0211     double corrected_residual = raw_residual - qoverpt_correction - qoverpz_correction - dEdx_correction;
0212
0213     raw_hist->Fill(raw_residual * 1000.);
0214     narrowed_hist->Fill(corrected_residual * 1000.);
0215
0216     qoverpt_hist->Fill(qoverpt, (raw_residual - qoverpz_correction - dEdx_correction) * 1000.);
0217     qoverpz_hist->Fill(qoverpz, (raw_residual - qoverpt_correction - dEdx_correction) * 1000.);
0218     psquared_hist->Fill(psquared, (raw_residual - qoverpt_correction - qoverpz_correction) * 1000.);
0219   }
0220
0221   return 0.;
0222 }