Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:45:35

0001 // $Id: MuonResidualsFitter.cc,v 1.11 2011/04/15 21:48:21 khotilov Exp $
0002 
0003 #ifdef STANDALONE_FITTER
0004 #include "MuonResidualsFitter.h"
0005 #else
0006 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResidualsFitter.h"
0007 #endif
0008 
0009 #include <fstream>
0010 #include <set>
0011 #include "TMath.h"
0012 #include "TH1.h"
0013 #include "TF1.h"
0014 #include "TRobustEstimator.h"
0015 
0016 // all global variables begin with "MuonResidualsFitter_" to avoid
0017 // namespace clashes (that is, they do what would ordinarily be done
0018 // with a class structure, but Minuit requires them to be global)
0019 
0020 const double MuonResidualsFitter_gsbinsize = 0.01;
0021 const double MuonResidualsFitter_tsbinsize = 0.1;
0022 const int MuonResidualsFitter_numgsbins = 500;
0023 const int MuonResidualsFitter_numtsbins = 500;
0024 
0025 bool MuonResidualsFitter_table_initialized = false;
0026 double MuonResidualsFitter_lookup_table[MuonResidualsFitter_numgsbins][MuonResidualsFitter_numtsbins];
0027 
0028 static TMinuit *MuonResidualsFitter_TMinuit;
0029 
0030 // fit function
0031 double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma) {
0032   sigma = fabs(sigma);
0033   static const double cgaus = 0.5 * log(2. * M_PI);
0034   return (-pow(residual - center, 2) * 0.5 / sigma / sigma) - cgaus - log(sigma);
0035 }
0036 
0037 // TF1 interface version
0038 Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par) {
0039   return par[0] * exp(MuonResidualsFitter_logPureGaussian(xvec[0], par[1], par[2]));
0040 }
0041 
0042 // 2D Gauss fit function
0043 double MuonResidualsFitter_logPureGaussian2D(double x, double y, double x0, double y0, double sx, double sy, double r) {
0044   sx = fabs(sx);
0045   sy = fabs(sy);
0046   static const double c2gaus = log(2. * M_PI);
0047   double one_r2 = 1. - r * r;
0048   double dx = x - x0;
0049   double dy = y - y0;
0050   return -0.5 / one_r2 * (pow(dx / sx, 2) + pow(dy / sy, 2) - 2. * r * dx / sx * dy / sy) - c2gaus -
0051          log(sx * sy * sqrt(one_r2));
0052 }
0053 
0054 double MuonResidualsFitter_compute_log_convolution(
0055     double toversigma, double gammaoversigma, double max, double step, double power) {
0056   if (gammaoversigma == 0.)
0057     return (-toversigma * toversigma / 2.) - log(sqrt(2 * M_PI));
0058 
0059   double sum = 0.;
0060   double uplus = 0.;
0061   double integrandplus_last = 0.;
0062   double integrandminus_last = 0.;
0063   for (double inc = 0.; uplus < max; inc += step) {
0064     double uplus_last = uplus;
0065     uplus = pow(inc, power);
0066 
0067     double integrandplus = exp(-pow(uplus - toversigma, 2) / 2.) / (uplus * uplus / gammaoversigma + gammaoversigma) /
0068                            M_PI / sqrt(2. * M_PI);
0069     double integrandminus = exp(-pow(-uplus - toversigma, 2) / 2.) / (uplus * uplus / gammaoversigma + gammaoversigma) /
0070                             M_PI / sqrt(2. * M_PI);
0071 
0072     sum += integrandplus * (uplus - uplus_last);
0073     sum += 0.5 * fabs(integrandplus_last - integrandplus) * (uplus - uplus_last);
0074 
0075     sum += integrandminus * (uplus - uplus_last);
0076     sum += 0.5 * fabs(integrandminus_last - integrandminus) * (uplus - uplus_last);
0077 
0078     integrandplus_last = integrandplus;
0079     integrandminus_last = integrandminus;
0080   }
0081   return log(sum);
0082 }
0083 
0084 // fit function
0085 double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma) {
0086   sigma = fabs(sigma);
0087   double gammaoversigma = fabs(gamma / sigma);
0088   double toversigma = fabs((residual - center) / sigma);
0089 
0090   int gsbin0 = int(floor(gammaoversigma / MuonResidualsFitter_gsbinsize));
0091   int gsbin1 = int(ceil(gammaoversigma / MuonResidualsFitter_gsbinsize));
0092   int tsbin0 = int(floor(toversigma / MuonResidualsFitter_tsbinsize));
0093   int tsbin1 = int(ceil(toversigma / MuonResidualsFitter_tsbinsize));
0094 
0095   bool gsisbad = (gsbin0 >= MuonResidualsFitter_numgsbins || gsbin1 >= MuonResidualsFitter_numgsbins);
0096   bool tsisbad = (tsbin0 >= MuonResidualsFitter_numtsbins || tsbin1 >= MuonResidualsFitter_numtsbins);
0097 
0098   if (gsisbad || tsisbad) {
0099     return log(gammaoversigma / M_PI) - log(toversigma * toversigma + gammaoversigma * gammaoversigma) - log(sigma);
0100   } else {
0101     double val00 = MuonResidualsFitter_lookup_table[gsbin0][tsbin0];
0102     double val01 = MuonResidualsFitter_lookup_table[gsbin0][tsbin1];
0103     double val10 = MuonResidualsFitter_lookup_table[gsbin1][tsbin0];
0104     double val11 = MuonResidualsFitter_lookup_table[gsbin1][tsbin1];
0105 
0106     double val0 = val00 + ((toversigma / MuonResidualsFitter_tsbinsize) - tsbin0) * (val01 - val00);
0107     double val1 = val10 + ((toversigma / MuonResidualsFitter_tsbinsize) - tsbin0) * (val11 - val10);
0108 
0109     double val = val0 + ((gammaoversigma / MuonResidualsFitter_gsbinsize) - gsbin0) * (val1 - val0);
0110 
0111     return val - log(sigma);
0112   }
0113 }
0114 
0115 // TF1 interface version
0116 Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par) {
0117   return par[0] * exp(MuonResidualsFitter_logPowerLawTails(xvec[0], par[1], par[2], par[3]));
0118 }
0119 
0120 double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma) {
0121   return log(TMath::Voigt(residual - center, fabs(sigma), fabs(gamma) * 2.));
0122 }
0123 
0124 Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par) {
0125   return par[0] * TMath::Voigt(xvec[0] - par[1], fabs(par[2]), fabs(par[3]) * 2.);
0126 }
0127 
0128 double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma) {
0129   double x = residual - center;
0130   double s = fabs(sigma);
0131   double m = 2. * s;
0132   static const double e2 = exp(-2.), sqerf = sqrt(2. * M_PI) * erf(sqrt(2.));
0133   double a = pow(m, 4) * e2;
0134   double n = sqerf * s + 2. * a * pow(m, -3) / 3.;
0135 
0136   if (fabs(x) < m)
0137     return -x * x / (2. * s * s) - log(n);
0138   else
0139     return log(a) - 4. * log(fabs(x)) - log(n);
0140 }
0141 
0142 Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par) {
0143   return par[0] * exp(MuonResidualsFitter_logGaussPowerTails(xvec[0], par[1], par[2]));
0144 }
0145 
0146 double MuonResidualsFitter_integrate_pureGaussian(double low, double high, double center, double sigma) {
0147   static const double isqr2 = 1. / sqrt(2.);
0148   return (erf((high + center) * isqr2 / sigma) - erf((low + center) * isqr2 / sigma)) * exp(0.5 / sigma / sigma) * 0.5;
0149 }
0150 
0151 MuonResidualsFitter::MuonResidualsFitter(int residualsModel, int minHits, int useResiduals, bool weightAlignment)
0152     : m_residualsModel(residualsModel),
0153       m_minHits(minHits),
0154       m_useResiduals(useResiduals),
0155       m_weightAlignment(weightAlignment),
0156       m_printLevel(0),
0157       m_strategy(1),
0158       m_cov(1),
0159       m_loglikelihood(0.) {
0160   if (m_residualsModel != kPureGaussian && m_residualsModel != kPowerLawTails && m_residualsModel != kROOTVoigt &&
0161       m_residualsModel != kGaussPowerTails && m_residualsModel != kPureGaussian2D)
0162     throw cms::Exception("MuonResidualsFitter") << "unrecognized residualsModel";
0163 }
0164 
0165 MuonResidualsFitter::~MuonResidualsFitter() {
0166   for (std::vector<double *>::const_iterator residual = residuals_begin(); residual != residuals_end(); ++residual) {
0167     delete[](*residual);
0168   }
0169 }
0170 
0171 void MuonResidualsFitter::fix(int parNum, bool dofix) {
0172   assert(0 <= parNum && parNum < npar());
0173   if (m_fixed.empty())
0174     m_fixed.resize(npar(), false);
0175   m_fixed[parNum] = dofix;
0176 }
0177 
0178 bool MuonResidualsFitter::fixed(int parNum) {
0179   assert(0 <= parNum && parNum < npar());
0180   if (m_fixed.empty())
0181     return false;
0182   else
0183     return m_fixed[parNum];
0184 }
0185 
0186 void MuonResidualsFitter::fill(double *residual) {
0187   m_residuals.push_back(residual);
0188   m_residuals_ok.push_back(true);
0189 }
0190 
0191 double MuonResidualsFitter::covarianceElement(int parNum1, int parNum2) {
0192   assert(0 <= parNum1 && parNum1 < npar());
0193   assert(0 <= parNum2 && parNum2 < npar());
0194   assert(m_cov.GetNcols() == npar());  // m_cov might have not yet been resized to account for proper #parameters
0195   return m_cov(parNum2parIdx(parNum1), parNum2parIdx(parNum2));
0196 }
0197 
0198 long MuonResidualsFitter::numsegments() {
0199   long num = 0;
0200   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter)
0201     num++;
0202   return num;
0203 }
0204 
0205 void MuonResidualsFitter::initialize_table() {
0206   if (MuonResidualsFitter_table_initialized || residualsModel() != kPowerLawTails)
0207     return;
0208   MuonResidualsFitter_table_initialized = true;
0209 
0210   std::ifstream convolution_table("convolution_table.txt");
0211   if (convolution_table.is_open()) {
0212     int numgsbins = 0;
0213     int numtsbins = 0;
0214     double tsbinsize = 0.;
0215     double gsbinsize = 0.;
0216 
0217     convolution_table >> numgsbins >> numtsbins >> tsbinsize >> gsbinsize;
0218     if (numgsbins != MuonResidualsFitter_numgsbins || numtsbins != MuonResidualsFitter_numtsbins ||
0219         tsbinsize != MuonResidualsFitter_tsbinsize || gsbinsize != MuonResidualsFitter_gsbinsize) {
0220       throw cms::Exception("MuonResidualsFitter") << "convolution_table.txt has the wrong bin width/bin size.  Throw "
0221                                                      "it away and let the fitter re-create the file.\n";
0222     }
0223 
0224     for (int gsbin = 0; gsbin < MuonResidualsFitter_numgsbins; gsbin++) {
0225       for (int tsbin = 0; tsbin < MuonResidualsFitter_numtsbins; tsbin++) {
0226         int read_gsbin = 0;
0227         int read_tsbin = 0;
0228         double val = 0.;
0229 
0230         convolution_table >> read_gsbin >> read_tsbin >> val;
0231         if (read_gsbin != gsbin || read_tsbin != tsbin) {
0232           throw cms::Exception("MuonResidualsFitter")
0233               << "convolution_table.txt is out of order.  Throw it away and let the fitter re-create the file.\n";
0234         }
0235         MuonResidualsFitter_lookup_table[gsbin][tsbin] = val;
0236       }
0237     }
0238     convolution_table.close();
0239   } else {
0240     std::ofstream convolution_table2("convolution_table.txt");
0241 
0242     if (!convolution_table2.is_open())
0243       throw cms::Exception("MuonResidualsFitter") << "Couldn't write to file convolution_table.txt\n";
0244 
0245     convolution_table2 << MuonResidualsFitter_numgsbins << " " << MuonResidualsFitter_numtsbins << " "
0246                        << MuonResidualsFitter_tsbinsize << " " << MuonResidualsFitter_gsbinsize << std::endl;
0247 
0248     std::cout << "Initializing convolution look-up table (takes a few minutes)..." << std::endl;
0249 
0250     for (int gsbin = 0; gsbin < MuonResidualsFitter_numgsbins; gsbin++) {
0251       double gammaoversigma = double(gsbin) * MuonResidualsFitter_gsbinsize;
0252       std::cout << "    gsbin " << gsbin << "/" << MuonResidualsFitter_numgsbins << std::endl;
0253       for (int tsbin = 0; tsbin < MuonResidualsFitter_numtsbins; tsbin++) {
0254         double toversigma = double(tsbin) * MuonResidualsFitter_tsbinsize;
0255 
0256         // 1e-6 errors (out of a value of ~0.01) with max=100, step=0.001, power=4 (max=1000 does a little better with the tails)
0257         MuonResidualsFitter_lookup_table[gsbin][tsbin] =
0258             MuonResidualsFitter_compute_log_convolution(toversigma, gammaoversigma);
0259 
0260         // <10% errors with max=20, step=0.005, power=4 (faster computation for testing)
0261         // MuonResidualsFitter_lookup_table[gsbin][tsbin] = MuonResidualsFitter_compute_log_convolution(toversigma, gammaoversigma, 100., 0.005, 4.);
0262 
0263         convolution_table2 << gsbin << " " << tsbin << " " << MuonResidualsFitter_lookup_table[gsbin][tsbin]
0264                            << std::endl;
0265       }
0266     }
0267     convolution_table2.close();
0268     std::cout << "Initialization done!" << std::endl;
0269   }
0270 }
0271 
0272 bool MuonResidualsFitter::dofit(void (*fcn)(int &, double *, double &, double *, int),
0273                                 std::vector<int> &parNum,
0274                                 std::vector<std::string> &parName,
0275                                 std::vector<double> &start,
0276                                 std::vector<double> &step,
0277                                 std::vector<double> &low,
0278                                 std::vector<double> &high) {
0279   MuonResidualsFitterFitInfo *fitinfo = new MuonResidualsFitterFitInfo(this);
0280 
0281   MuonResidualsFitter_TMinuit = new TMinuit(npar());
0282   // MuonResidualsFitter_TMinuit->SetPrintLevel(m_printLevel);
0283   MuonResidualsFitter_TMinuit->SetPrintLevel();
0284   MuonResidualsFitter_TMinuit->SetObjectFit(fitinfo);
0285   MuonResidualsFitter_TMinuit->SetFCN(fcn);
0286   inform(MuonResidualsFitter_TMinuit);
0287 
0288   std::vector<int>::const_iterator iNum = parNum.begin();
0289   std::vector<std::string>::const_iterator iName = parName.begin();
0290   std::vector<double>::const_iterator istart = start.begin();
0291   std::vector<double>::const_iterator istep = step.begin();
0292   std::vector<double>::const_iterator ilow = low.begin();
0293   std::vector<double>::const_iterator ihigh = high.begin();
0294 
0295   //MuonResidualsFitter_TMinuit->SetPrintLevel(-1);
0296 
0297   for (; iNum != parNum.end(); ++iNum, ++iName, ++istart, ++istep, ++ilow, ++ihigh) {
0298     MuonResidualsFitter_TMinuit->DefineParameter(*iNum, iName->c_str(), *istart, *istep, *ilow, *ihigh);
0299     if (fixed(*iNum))
0300       MuonResidualsFitter_TMinuit->FixParameter(*iNum);
0301   }
0302 
0303   double arglist[10];
0304   int ierflg;
0305   int smierflg;  //second MIGRAD ierflg
0306 
0307   // chi^2 errors should be 1.0, log-likelihood should be 0.5
0308   for (int i = 0; i < 10; i++)
0309     arglist[i] = 0.;
0310   arglist[0] = 0.5;
0311   ierflg = 0;
0312   smierflg = 0;
0313   MuonResidualsFitter_TMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
0314   if (ierflg != 0) {
0315     delete MuonResidualsFitter_TMinuit;
0316     delete fitinfo;
0317     return false;
0318   }
0319 
0320   // set strategy = 2 (more refined fits)
0321   for (int i = 0; i < 10; i++)
0322     arglist[i] = 0.;
0323   arglist[0] = m_strategy;
0324   ierflg = 0;
0325   MuonResidualsFitter_TMinuit->mnexcm("SET STR", arglist, 1, ierflg);
0326   if (ierflg != 0) {
0327     delete MuonResidualsFitter_TMinuit;
0328     delete fitinfo;
0329     return false;
0330   }
0331 
0332   bool try_again = false;
0333 
0334   // minimize
0335   for (int i = 0; i < 10; i++)
0336     arglist[i] = 0.;
0337   arglist[0] = 50000;
0338   ierflg = 0;
0339   MuonResidualsFitter_TMinuit->mnexcm("MIGRAD", arglist, 1, ierflg);
0340   if (ierflg != 0)
0341     try_again = true;
0342 
0343   // just once more, if needed (using the final Minuit parameters from the failed fit; often works)
0344   if (try_again) {
0345     for (int i = 0; i < 10; i++)
0346       arglist[i] = 0.;
0347     arglist[0] = 50000;
0348     MuonResidualsFitter_TMinuit->mnexcm("MIGRAD", arglist, 1, smierflg);
0349   }
0350 
0351   Double_t fmin, fedm, errdef;
0352   Int_t npari, nparx, istat;
0353   MuonResidualsFitter_TMinuit->mnstat(fmin, fedm, errdef, npari, nparx, istat);
0354 
0355   if (istat != 3) {
0356     for (int i = 0; i < 10; i++)
0357       arglist[i] = 0.;
0358     ierflg = 0;
0359     MuonResidualsFitter_TMinuit->mnexcm("HESSE", arglist, 0, ierflg);
0360   }
0361 
0362   // read-out the results
0363   m_loglikelihood = -fmin;
0364 
0365   m_value.clear();
0366   m_error.clear();
0367   for (int i = 0; i < npar(); i++) {
0368     double v, e;
0369     MuonResidualsFitter_TMinuit->GetParameter(i, v, e);
0370     m_value.push_back(v);
0371     m_error.push_back(e);
0372   }
0373   m_cov.ResizeTo(npar(), npar());
0374   MuonResidualsFitter_TMinuit->mnemat(m_cov.GetMatrixArray(), npar());
0375 
0376   delete MuonResidualsFitter_TMinuit;
0377   delete fitinfo;
0378   if (smierflg != 0)
0379     return false;
0380   return true;
0381 }
0382 
0383 void MuonResidualsFitter::write(FILE *file, int which) {
0384   long rows = numResiduals();
0385   int cols = ndata();
0386   int whichcopy = which;
0387 
0388   fwrite(&rows, sizeof(long), 1, file);
0389   fwrite(&cols, sizeof(int), 1, file);
0390   fwrite(&whichcopy, sizeof(int), 1, file);
0391 
0392   double *likeAChecksum = new double[cols];
0393   double *likeAChecksum2 = new double[cols];
0394   for (int i = 0; i < cols; i++) {
0395     likeAChecksum[i] = 0.;
0396     likeAChecksum2[i] = 0.;
0397   }
0398 
0399   for (std::vector<double *>::const_iterator residual = residuals_begin(); residual != residuals_end(); ++residual) {
0400     fwrite((*residual), sizeof(double), cols, file);
0401     for (int i = 0; i < cols; i++) {
0402       if (fabs((*residual)[i]) > likeAChecksum[i])
0403         likeAChecksum[i] = fabs((*residual)[i]);
0404       if (fabs((*residual)[i]) < likeAChecksum2[i])
0405         likeAChecksum2[i] = fabs((*residual)[i]);
0406     }
0407   }  // end loop over residuals
0408 
0409   // the idea is that mal-formed doubles are likely to be huge values (or tiny values)
0410   // because the exponent gets screwed up; we want to check for that
0411   fwrite(likeAChecksum, sizeof(double), cols, file);
0412   fwrite(likeAChecksum2, sizeof(double), cols, file);
0413 
0414   delete[] likeAChecksum;
0415   delete[] likeAChecksum2;
0416 }
0417 
0418 void MuonResidualsFitter::read(FILE *file, int which) {
0419   long rows = -100;
0420   int cols = -100;
0421   int readwhich = -100;
0422 
0423   fread(&rows, sizeof(long), 1, file);
0424   fread(&cols, sizeof(int), 1, file);
0425   fread(&readwhich, sizeof(int), 1, file);
0426 
0427   if (cols != ndata() || rows < 0 || readwhich != which) {
0428     throw cms::Exception("MuonResidualsFitter")
0429         << "temporary file is corrupted (which = " << which << " readwhich = " << readwhich << " rows = " << rows
0430         << " cols = " << cols << ")\n";
0431   }
0432 
0433   double *likeAChecksum = new double[cols];
0434   double *likeAChecksum2 = new double[cols];
0435   for (int i = 0; i < cols; i++) {
0436     likeAChecksum[i] = 0.;
0437     likeAChecksum2[i] = 0.;
0438   }
0439 
0440   m_residuals.reserve(rows);
0441   for (long row = 0; row < rows; row++) {
0442     double *residual = new double[cols];
0443     fread(residual, sizeof(double), cols, file);
0444     fill(residual);
0445     for (int i = 0; i < cols; i++) {
0446       if (fabs(residual[i]) > likeAChecksum[i])
0447         likeAChecksum[i] = fabs(residual[i]);
0448       if (fabs(residual[i]) < likeAChecksum2[i])
0449         likeAChecksum2[i] = fabs(residual[i]);
0450     }
0451   }  // end loop over records in file
0452 
0453   double *readChecksum = new double[cols];
0454   double *readChecksum2 = new double[cols];
0455   fread(readChecksum, sizeof(double), cols, file);
0456   fread(readChecksum2, sizeof(double), cols, file);
0457 
0458   for (int i = 0; i < cols; i++) {
0459     if (fabs(likeAChecksum[i] - readChecksum[i]) > 1e-10 ||
0460         fabs(1. / likeAChecksum2[i] - 1. / readChecksum2[i]) > 1e10) {
0461       throw cms::Exception("MuonResidualsFitter")
0462           << "temporary file is corrupted (which = " << which << " rows = " << rows << " likeAChecksum "
0463           << likeAChecksum[i] << " != readChecksum " << readChecksum[i] << " "
0464           << " likeAChecksum2 " << likeAChecksum2[i] << " != readChecksum2 " << readChecksum2[i] << ")\n";
0465     }
0466   }
0467 
0468   m_residuals_ok.resize(numResiduals(), true);
0469 
0470   delete[] likeAChecksum;
0471   delete[] likeAChecksum2;
0472 }
0473 
0474 void MuonResidualsFitter::plotsimple(std::string name, TFileDirectory *dir, int which, double multiplier) {
0475   double window = 100.;
0476   if (which == 0)
0477     window = 2. * 30.;
0478   else if (which == 1)
0479     window = 2. * 30.;
0480   else if (which == 2)
0481     window = 2. * 20.;
0482   else if (which == 3)
0483     window = 2. * 50.;
0484 
0485   TH1F *hist = dir->make<TH1F>(name.c_str(), "", 200, -window, window);
0486   for (std::vector<double *>::const_iterator r = residuals_begin(); r != residuals_end(); ++r)
0487     hist->Fill(multiplier * (*r)[which]);
0488 }
0489 
0490 void MuonResidualsFitter::plotweighted(
0491     std::string name, TFileDirectory *dir, int which, int whichredchi2, double multiplier) {
0492   double window = 100.;
0493   if (which == 0)
0494     window = 2. * 30.;
0495   else if (which == 1)
0496     window = 2. * 30.;
0497   else if (which == 2)
0498     window = 2. * 20.;
0499   else if (which == 3)
0500     window = 2. * 50.;
0501 
0502   TH1F *hist = dir->make<TH1F>(name.c_str(), "", 200, -window, window);
0503   for (std::vector<double *>::const_iterator r = residuals_begin(); r != residuals_end(); ++r) {
0504     double weight = 1. / (*r)[whichredchi2];
0505     if (TMath::Prob(1. / weight * 12, 12) < 0.99)
0506       hist->Fill(multiplier * (*r)[which], weight);
0507   }
0508 }
0509 
0510 void MuonResidualsFitter::computeHistogramRangeAndBinning(int which, int &nbins, double &a, double &b) {
0511   // first, make a numeric array while discarding some crazy outliers
0512   double *data = new double[numResiduals()];
0513   int n = 0;
0514   for (std::vector<double *>::const_iterator r = m_residuals.begin(); r != m_residuals.end(); r++)
0515     if (fabs((*r)[which]) < 50.) {
0516       data[n] = (*r)[which];
0517       n++;
0518     }
0519 
0520   // compute "3 normal sigma" and regular interquantile ranges
0521   const int n_quantiles = 7;
0522   double probabilities[n_quantiles] = {0.00135, 0.02275, 0.25, 0.5, 0.75, 0.97725, 0.99865};  // "3 normal sigma"
0523   //double probabilities[n_quantiles] = {0.02275, 0.25, 0.75, 0.97725}; // "2 normal sigma"
0524   double quantiles[n_quantiles];
0525   std::sort(data, data + n);
0526   TMath::Quantiles(n, n_quantiles, data, quantiles, probabilities, true, nullptr, 7);
0527   delete[] data;
0528   double iqr = quantiles[4] - quantiles[2];
0529 
0530   // estimate optimal bin size according to Freedman-Diaconis rule
0531   double hbin = 2 * iqr / pow(n, 1. / 3);
0532 
0533   a = quantiles[1];
0534   b = quantiles[5];
0535   nbins = (int)((b - a) / hbin + 3.);  // add extra safety margin of 3
0536 
0537   std::cout << "   quantiles: ";
0538   for (int i = 0; i < n_quantiles; i++)
0539     std::cout << quantiles[i] << " ";
0540   std::cout << std::endl;
0541   //cout<<"n="<<select_count<<" quantiles ["<<quantiles[1]<<", "<<quantiles[2]<<"]  IQR="<<iqr
0542   //  <<"  full range=["<<minx<<","<<maxx<<"]"<<"  2 normal sigma quantile range = ["<<quantiles[0]<<", "<<quantiles[3]<<"]"<<endl;
0543   std::cout << "   optimal h=" << hbin << " nbins=" << nbins << std::endl;
0544 }
0545 
0546 void MuonResidualsFitter::histogramChi2GaussianFit(int which, double &fit_mean, double &fit_sigma) {
0547   int nbins;
0548   double a, b;
0549   computeHistogramRangeAndBinning(which, nbins, a, b);
0550   if (a == b || a > b) {
0551     fit_mean = a;
0552     fit_sigma = 0;
0553     return;
0554   }
0555 
0556   TH1D *hist = new TH1D("htmp", "", nbins, a, b);
0557   for (std::vector<double *>::const_iterator r = m_residuals.begin(); r != m_residuals.end(); ++r)
0558     hist->Fill((*r)[which]);
0559 
0560   // do simple chi2 gaussian fit
0561   TF1 *f1 = new TF1("f1", "gaus", a, b);
0562   f1->SetParameter(0, hist->GetEntries());
0563   f1->SetParameter(1, 0);
0564   f1->SetParameter(2, hist->GetRMS());
0565   hist->Fit("f1", "RQ");
0566   // hist->Fit(f1,"RQ");
0567 
0568   fit_mean = f1->GetParameter(1);
0569   fit_sigma = f1->GetParameter(2);
0570   std::cout << " h(" << nbins << "," << a << "," << b << ") mu=" << fit_mean << " sig=" << fit_sigma << std::endl;
0571 
0572   delete f1;
0573   delete hist;
0574 }
0575 
0576 // simple non-turned ellipsoid selection
0577 // THIS WAS selectPeakResiduals_simple, but I changed it to use this simple function as default
0578 // void MuonResidualsFitter::selectPeakResiduals_simple(double nsigma, int nvar, int *vars)
0579 void MuonResidualsFitter::selectPeakResiduals(double nsigma, int nvar, int *vars) {
0580   // does not make sense for small statistics
0581   if (numResiduals() < 25)
0582     return;
0583 
0584   int nbefore = numResiduals();
0585 
0586   //just to be sure (can't see why it might ever be more then 10)
0587   assert(nvar <= 10);
0588 
0589   // estimate nvar-D ellipsoid center & axes
0590   for (int v = 0; v < nvar; v++) {
0591     int which = vars[v];
0592     histogramChi2GaussianFit(which, m_center[which], m_radii[which]);
0593     m_radii[which] = nsigma * m_radii[which];
0594   }
0595 
0596   // filter out residuals that don't fit into the ellipsoid
0597   std::vector<double *>::iterator r = m_residuals.begin();
0598   while (r != m_residuals.end()) {
0599     double ellipsoid_sum = 0;
0600     for (int v = 0; v < nvar; v++) {
0601       int which = vars[v];
0602       if (m_radii[which] == 0.)
0603         continue;
0604       ellipsoid_sum += pow(((*r)[which] - m_center[which]) / m_radii[which], 2);
0605     }
0606     if (ellipsoid_sum <= 1.)
0607       ++r;
0608     else {
0609       m_residuals_ok[r - m_residuals.begin()] = false;
0610       ++r;
0611       // delete [] (*r);
0612       // r = m_residuals.erase(r);
0613     }
0614   }
0615   std::cout << " N residuals " << nbefore << " -> " << numResiduals() << std::endl;
0616 }
0617 
0618 // pre-selection using robust covariance estimator
0619 // THIS WAS selectPeakResiduals but I changed it to use OTHER simple function as default
0620 void MuonResidualsFitter::selectPeakResiduals_simple(double nsigma, int nvar, int *vars) {
0621   //std::cout<<"doing selectpeakresiduals: nsig="<<nsigma<<" nvar="<<nvar<<" vars=";
0622   for (int i = 0; i < nvar; ++i)
0623     std::cout << vars[i] << " ";
0624   std::cout << std::endl;
0625 
0626   // does not make sense for small statistics set to 50
0627   // YP changed it to 10 for test
0628   if (numResiduals() < 10)
0629     return;
0630   // if (numResiduals()<50) return;
0631 
0632   size_t nbefore = numResiduals();
0633   std::cout << " N residuals " << nbefore << " ~ "
0634             << (size_t)std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true) << std::endl;
0635   //just to be sure (can't see why it might ever be more then 10)
0636   assert(nvar <= 10 && nvar > 0);
0637 
0638   std::vector<double *>::iterator r = m_residuals.begin();
0639 
0640   // it's awkward, but the 1D case has to be handled separately
0641   if (nvar == 1) {
0642     std::cout << "1D case" << std::endl;
0643     // get robust estimates for the peak and sigma
0644     double *data = new double[nbefore];
0645     for (size_t i = 0; i < nbefore; i++)
0646       data[i] = m_residuals[i][vars[0]];
0647     double peak, sigma;
0648     TRobustEstimator re;
0649     re.EvaluateUni(nbefore, data, peak, sigma);
0650 
0651     // filter out residuals that are more then nsigma away from the peak
0652     while (r != m_residuals.end()) {
0653       double distance = fabs(((*r)[vars[0]] - peak) / sigma);
0654       if (distance <= nsigma)
0655         ++r;
0656       else {
0657         m_residuals_ok[r - m_residuals.begin()] = false;
0658         ++r;
0659         //delete [] (*r);
0660         //r = m_residuals.erase(r);
0661       }
0662     }
0663     std::cout << " N residuals " << nbefore << " -> " << numResiduals() << std::endl;
0664     return;
0665   }  // end 1D case
0666 
0667   // initialize and run the robust estimator for D>1
0668   std::cout << "D>1 case" << std::endl;
0669   TRobustEstimator re(nbefore + 1, nvar);
0670   std::cout << "nbefore " << nbefore << " nvar " << nvar << std::endl;
0671   r = m_residuals.begin();
0672   std::cout << "+++++ JUST before loop while (r != m_residuals.end())" << std::endl;
0673   int counter1 = 0;
0674   while (r != m_residuals.end()) {
0675     double *row = new double[nvar];
0676     for (int v = 0; v < nvar; v++)
0677       row[v] = (*r)[vars[v]];
0678     re.AddRow(row);
0679     // delete[] row;
0680     ++r;
0681     counter1++;
0682   }
0683   std::cout << "counter1 " << counter1 << std::endl;
0684   std::cout << "+++++ JUST after loop while (r != m_residuals.end())" << std::endl;
0685   re.Evaluate();
0686   std::cout << "+++++ JUST after re.Evaluate()" << std::endl;
0687 
0688   // get nvar-dimensional ellipsoid center & covariance
0689   TVectorD M(nvar);
0690   re.GetMean(M);
0691   TMatrixDSym Cov(nvar);
0692   re.GetCovariance(Cov);
0693   Cov.Invert();
0694 
0695   // calculate the normalized radius for this nvar-dimensional ellipsoid from a 1D-Gaussian nsigma equivalent distance
0696   double conf_1d = TMath::Erf(nsigma / sqrt(2));
0697   double surf_radius = sqrt(TMath::ChisquareQuantile(conf_1d, nvar));
0698 
0699   // filter out residuals that are outside of the covariance ellipsoid with the above normalized radius
0700   r = m_residuals.begin();
0701   while (r != m_residuals.end()) {
0702     TVectorD res(nvar);
0703     for (int v = 0; v < nvar; v++)
0704       res[v] = (*r)[vars[v]];
0705     double distance = sqrt(Cov.Similarity(res - M));
0706     if (distance <= surf_radius)
0707       ++r;
0708     else {
0709       m_residuals_ok[r - m_residuals.begin()] = false;
0710       ++r;
0711       //delete [] (*r);
0712       //r = m_residuals.erase(r);
0713     }
0714   }
0715   std::cout << " N residuals " << nbefore << " -> "
0716             << (size_t)std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true) << std::endl;
0717 }
0718 
0719 void MuonResidualsFitter::fiducialCuts(double xMin, double xMax, double yMin, double yMax, bool fidcut1) {
0720   int iResidual = -1;
0721 
0722   int n_station = 9999;
0723   int n_wheel = 9999;
0724   int n_sector = 9999;
0725 
0726   double positionX = 9999.;
0727   double positionY = 9999.;
0728 
0729   double chambw = 9999.;
0730   double chambl = 9999.;
0731 
0732   for (std::vector<double *>::const_iterator r = residuals_begin(); r != residuals_end(); ++r) {
0733     iResidual++;
0734     if (!m_residuals_ok[iResidual])
0735       continue;
0736 
0737     if ((*r)[15] >
0738         0.0001) {  // this value is greater than zero (chamber width) for 6DOFs stations 1,2,3 better to change for type()!!!
0739       n_station = (*r)[12];
0740       n_wheel = (*r)[13];
0741       n_sector = (*r)[14];
0742       positionX = (*r)[4];
0743       positionY = (*r)[5];
0744       chambw = (*r)[15];
0745       chambl = (*r)[16];
0746     } else {  // in case of 5DOF residual the residual object index is different
0747       n_station = (*r)[10];
0748       n_wheel = (*r)[11];
0749       n_sector = (*r)[12];
0750       positionX = (*r)[2];
0751       positionY = (*r)[3];
0752       chambw = (*r)[13];
0753       chambl = (*r)[14];
0754     }
0755 
0756     if (fidcut1) {  // this is the standard fiducial cut used so far 80x80 cm in x,y
0757       if (positionX >= xMax || positionX <= xMin)
0758         m_residuals_ok[iResidual] = false;
0759       if (positionY >= yMax || positionY <= yMin)
0760         m_residuals_ok[iResidual] = false;
0761     }
0762 
0763     // Implementation of new fiducial cut
0764 
0765     double dtrkchamx = (chambw / 2.) - positionX;  // variables to cut tracks on the edge of the chambers
0766     double dtrkchamy = (chambl / 2.) - positionY;
0767 
0768     if (!fidcut1) {
0769       if (n_station == 4) {
0770         if ((n_wheel == -1 && n_sector == 3) ||
0771             (n_wheel == 1 &&
0772              n_sector == 4)) {  // FOR SHORT CHAMBER LENGTH IN:  WHEEL 1 SECTOR 4  AND  WHEEL -1 SECTOR 3
0773           if ((n_sector == 1 || n_sector == 2 || n_sector == 3 || n_sector == 5 || n_sector == 6 || n_sector == 7 ||
0774                n_sector == 8 || n_sector == 12) &&
0775               ((dtrkchamx < 40 || dtrkchamx > 380) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
0776             m_residuals_ok[iResidual] = false;
0777           if ((n_sector == 4 || n_sector == 13) &&
0778               ((dtrkchamx < 40 || dtrkchamx > 280) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
0779             m_residuals_ok[iResidual] = false;
0780           if ((n_sector == 9 || n_sector == 11) &&
0781               ((dtrkchamx < 40 || dtrkchamx > 180) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
0782             m_residuals_ok[iResidual] = false;
0783           if ((n_sector == 10 || n_sector == 14) &&
0784               ((dtrkchamx < 40 || dtrkchamx > 220) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
0785             m_residuals_ok[iResidual] = false;
0786         } else {
0787           if ((n_sector == 1 || n_sector == 2 || n_sector == 3 || n_sector == 5 || n_sector == 6 || n_sector == 7 ||
0788                n_sector == 8 || n_sector == 12) &&
0789               ((dtrkchamx < 40 || dtrkchamx > 380) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
0790             m_residuals_ok[iResidual] = false;
0791           if ((n_sector == 4 || n_sector == 13) &&
0792               ((dtrkchamx < 40 || dtrkchamx > 280) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
0793             m_residuals_ok[iResidual] = false;
0794           if ((n_sector == 9 || n_sector == 11) &&
0795               ((dtrkchamx < 40 || dtrkchamx > 180) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
0796             m_residuals_ok[iResidual] = false;
0797           if ((n_sector == 10 || n_sector == 14) &&
0798               ((dtrkchamx < 40 || dtrkchamx > 220) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
0799             m_residuals_ok[iResidual] = false;
0800         }
0801       } else {
0802         if ((n_wheel == -1 && n_sector == 3) || (n_wheel == 1 && n_sector == 4)) {
0803           if (n_station == 1 && ((dtrkchamx < 30.0 || dtrkchamx > 190.0) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
0804             m_residuals_ok[iResidual] = false;
0805           if (n_station == 2 && ((dtrkchamx < 30.0 || dtrkchamx > 240.0) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
0806             m_residuals_ok[iResidual] = false;
0807           if (n_station == 3 && ((dtrkchamx < 30.0 || dtrkchamx > 280.0) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
0808             m_residuals_ok[iResidual] = false;
0809         } else {
0810           if (n_station == 1 && ((dtrkchamx < 30.0 || dtrkchamx > 190.0) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
0811             m_residuals_ok[iResidual] = false;
0812           if (n_station == 2 && ((dtrkchamx < 30.0 || dtrkchamx > 240.0) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
0813             m_residuals_ok[iResidual] = false;
0814           if (n_station == 3 && ((dtrkchamx < 30.0 || dtrkchamx > 280.0) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
0815             m_residuals_ok[iResidual] = false;
0816         }
0817       }
0818     }
0819   }
0820 }
0821 
0822 void MuonResidualsFitter::correctBField(int idx_momentum, int idx_q) {
0823   const int Nbin = 17;
0824   // find max 1/pt and bin width
0825   double min_pt = 9999.;
0826   for (std::vector<double *>::const_iterator r = residuals_begin(); r != residuals_end(); ++r) {
0827     double pt = fabs((*r)[idx_momentum]);
0828     if (pt < min_pt)
0829       min_pt = pt;
0830   }
0831   min_pt -= 0.01;  // to prevent bin # overflow
0832   const double bin_width = 1. / min_pt / Nbin;
0833 
0834   // fill indices of positive and negative charge residuals in each bin
0835   std::vector<size_t> pos[Nbin], neg[Nbin], to_erase;
0836   for (size_t i = 0; i < m_residuals.size(); i++) {
0837     if (!m_residuals_ok[i])
0838       continue;
0839     int bin = (int)floor(1. / fabs(m_residuals[i][idx_momentum]) / bin_width);
0840     if (m_residuals[i][idx_q] > 0)
0841       pos[bin].push_back(i);
0842     else
0843       neg[bin].push_back(i);
0844   }
0845 
0846   // equalize pos and neg in each bin
0847   for (int j = 0; j < Nbin; j++) {
0848     size_t psize = pos[j].size();
0849     size_t nsize = neg[j].size();
0850     if (psize == nsize)
0851       continue;
0852 
0853     std::set<int> idx_set;  // use a set to collect certain number of unique random indices to erase
0854     if (psize > nsize) {
0855       while (idx_set.size() < psize - nsize)
0856         idx_set.insert(gRandom->Integer(psize));
0857       for (std::set<int>::iterator it = idx_set.begin(); it != idx_set.end(); it++)
0858         to_erase.push_back(pos[j][*it]);
0859     } else {
0860       while (idx_set.size() < nsize - psize)
0861         idx_set.insert(gRandom->Integer(nsize));
0862       for (std::set<int>::iterator it = idx_set.begin(); it != idx_set.end(); it++)
0863         to_erase.push_back(neg[j][*it]);
0864     }
0865   }
0866   // sort in descending order, so we safely go from higher to lower indices:
0867   std::sort(to_erase.begin(), to_erase.end(), std::greater<int>());
0868   for (std::vector<size_t>::const_iterator e = to_erase.begin(); e != to_erase.end(); ++e) {
0869     m_residuals_ok[*e] = false;
0870     //delete[] *(m_residuals.begin() + *e);
0871     //m_residuals.erase(m_residuals.begin() + *e);
0872   }
0873 
0874   std::vector<size_t> apos[Nbin], aneg[Nbin];
0875   for (size_t i = 0; i < m_residuals.size(); i++) {
0876     if (!m_residuals_ok[i])
0877       continue;
0878     int bin = (int)floor(1. / fabs(m_residuals[i][idx_momentum]) / bin_width);
0879     if (m_residuals[i][idx_q] > 0)
0880       apos[bin].push_back(i);
0881     else
0882       aneg[bin].push_back(i);
0883   }
0884   for (int j = 0; j < Nbin; j++)
0885     std::cout << "bin " << j << ": [pos,neg] sizes before & after: [" << pos[j].size() << "," << neg[j].size()
0886               << "] -> [" << apos[j].size() << "," << aneg[j].size() << "]" << std::endl;
0887   std::cout << " N residuals " << m_residuals.size() << " -> "
0888             << (size_t)std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true) << std::endl;
0889 }
0890 
0891 void MuonResidualsFitter::eraseNotSelectedResiduals() {
0892   // it should probably be faster then doing erase
0893   size_t n_ok = (size_t)std::count(m_residuals_ok.begin(), m_residuals_ok.end(), true);
0894   std::vector<double *> tmp(n_ok, nullptr);
0895   std::cout << "residuals sizes: all=" << m_residuals.size() << " good=" << n_ok << std::endl;
0896   int iok = 0;
0897   for (size_t i = 0; i < m_residuals.size(); i++) {
0898     if (!m_residuals_ok[i])
0899       continue;
0900     tmp[iok++] = m_residuals[i];
0901   }
0902   m_residuals.swap(tmp);
0903 
0904   std::vector<bool> tmp_ok(n_ok, true);
0905   m_residuals_ok.swap(tmp_ok);
0906 
0907   std::cout << "residuals size after eraseNotSelectedResiduals =" << m_residuals.size()
0908             << "  ok size=" << m_residuals_ok.size() << std::endl;
0909 }