3 #ifdef STANDALONE_FITTER
14 #include "TRobustEstimator.h"
33 static const double cgaus = 0.5 *
log(2. *
M_PI);
34 return (-
pow(residual - center, 2) * 0.5 / sigma / sigma) - cgaus -
log(sigma);
46 static const double c2gaus =
log(2. *
M_PI);
47 double one_r2 = 1. - r *
r;
50 return -0.5 / one_r2 * (
pow(dx / sx, 2) +
pow(dy / sy, 2) - 2. * r * dx / sx * dy / sy) - c2gaus -
55 double toversigma,
double gammaoversigma,
double max,
double step,
double power) {
56 if (gammaoversigma == 0.)
57 return (-toversigma * toversigma / 2.) -
log(
sqrt(2 *
M_PI));
61 double integrandplus_last = 0.;
62 double integrandminus_last = 0.;
63 for (
double inc = 0.; uplus <
max; inc +=
step) {
64 double uplus_last = uplus;
65 uplus =
pow(inc, power);
67 double integrandplus =
exp(-
pow(uplus - toversigma, 2) / 2.) / (uplus * uplus / gammaoversigma + gammaoversigma) /
69 double integrandminus =
exp(-
pow(-uplus - toversigma, 2) / 2.) / (uplus * uplus / gammaoversigma + gammaoversigma) /
72 sum += integrandplus * (uplus - uplus_last);
73 sum += 0.5 * fabs(integrandplus_last - integrandplus) * (uplus - uplus_last);
75 sum += integrandminus * (uplus - uplus_last);
76 sum += 0.5 * fabs(integrandminus_last - integrandminus) * (uplus - uplus_last);
78 integrandplus_last = integrandplus;
79 integrandminus_last = integrandminus;
87 double gammaoversigma = fabs(gamma / sigma);
88 double toversigma = fabs((residual - center) / sigma);
98 if (gsisbad || tsisbad) {
99 return log(gammaoversigma /
M_PI) -
log(toversigma * toversigma + gammaoversigma * gammaoversigma) -
log(sigma);
111 return val -
log(sigma);
121 return log(TMath::Voigt(residual - center, fabs(sigma), fabs(gamma) * 2.));
125 return par[0] * TMath::Voigt(xvec[0] - par[1], fabs(par[2]), fabs(par[3]) * 2.);
129 double x = residual - center;
130 double s = fabs(sigma);
132 static const double e2 =
exp(-2.), sqerf =
sqrt(2. *
M_PI) * erf(
sqrt(2.));
133 double a =
pow(m, 4) * e2;
134 double n = sqerf * s + 2. * a *
pow(m, -3) / 3.;
137 return -x * x / (2. * s *
s) -
log(n);
139 return log(a) - 4. *
log(fabs(x)) -
log(n);
147 static const double isqr2 = 1. /
sqrt(2.);
148 return (erf((high + center) * isqr2 / sigma) - erf((low + center) * isqr2 / sigma)) *
exp(0.5 / sigma / sigma) * 0.5;
152 : m_residualsModel(residualsModel),
154 m_useResiduals(useResiduals),
155 m_weightAlignment(weightAlignment),
159 m_loglikelihood(0.) {
162 throw cms::Exception(
"MuonResidualsFitter") <<
"unrecognized residualsModel";
210 std::ifstream convolution_table(
"convolution_table.txt");
211 if (convolution_table.is_open()) {
214 double tsbinsize = 0.;
215 double gsbinsize = 0.;
217 convolution_table >> numgsbins >> numtsbins >> tsbinsize >> gsbinsize;
220 throw cms::Exception(
"MuonResidualsFitter") <<
"convolution_table.txt has the wrong bin width/bin size. Throw "
221 "it away and let the fitter re-create the file.\n";
230 convolution_table >> read_gsbin >> read_tsbin >>
val;
231 if (read_gsbin != gsbin || read_tsbin != tsbin) {
233 <<
"convolution_table.txt is out of order. Throw it away and let the fitter re-create the file.\n";
238 convolution_table.close();
240 std::ofstream convolution_table2(
"convolution_table.txt");
242 if (!convolution_table2.is_open())
243 throw cms::Exception(
"MuonResidualsFitter") <<
"Couldn't write to file convolution_table.txt\n";
248 std::cout <<
"Initializing convolution look-up table (takes a few minutes)..." << std::endl;
252 std::cout <<
" gsbin " << gsbin <<
"/" << MuonResidualsFitter_numgsbins << std::endl;
267 convolution_table2.close();
268 std::cout <<
"Initialization done!" << std::endl;
273 std::vector<int> &parNum,
274 std::vector<std::string> &parName,
275 std::vector<double> &
start,
276 std::vector<double> &
step,
277 std::vector<double> &low,
278 std::vector<double> &high) {
288 std::vector<int>::const_iterator iNum = parNum.begin();
289 std::vector<std::string>::const_iterator iName = parName.begin();
290 std::vector<double>::const_iterator istart = start.begin();
291 std::vector<double>::const_iterator
istep = step.begin();
292 std::vector<double>::const_iterator ilow = low.begin();
293 std::vector<double>::const_iterator ihigh = high.begin();
297 for (; iNum != parNum.end(); ++iNum, ++iName, ++istart, ++
istep, ++ilow, ++ihigh) {
308 for (
int i = 0;
i < 10;
i++)
321 for (
int i = 0;
i < 10;
i++)
332 bool try_again =
false;
335 for (
int i = 0;
i < 10;
i++)
345 for (
int i = 0;
i < 10;
i++)
351 Double_t fmin, fedm, errdef;
352 Int_t npari, nparx, istat;
356 for (
int i = 0;
i < 10;
i++)
367 for (
int i = 0;
i <
npar();
i++) {
386 int whichcopy =
which;
388 fwrite(&rows,
sizeof(
long), 1, file);
389 fwrite(&cols,
sizeof(
int), 1, file);
390 fwrite(&whichcopy,
sizeof(
int), 1, file);
392 double *likeAChecksum =
new double[cols];
393 double *likeAChecksum2 =
new double[cols];
394 for (
int i = 0;
i < cols;
i++) {
395 likeAChecksum[
i] = 0.;
396 likeAChecksum2[
i] = 0.;
400 fwrite((*residual),
sizeof(
double), cols, file);
401 for (
int i = 0;
i < cols;
i++) {
402 if (fabs((*residual)[
i]) > likeAChecksum[i])
403 likeAChecksum[
i] = fabs((*residual)[i]);
404 if (fabs((*residual)[i]) < likeAChecksum2[i])
405 likeAChecksum2[
i] = fabs((*residual)[i]);
411 fwrite(likeAChecksum,
sizeof(
double), cols, file);
412 fwrite(likeAChecksum2,
sizeof(
double), cols, file);
414 delete[] likeAChecksum;
415 delete[] likeAChecksum2;
421 int readwhich = -100;
423 fread(&rows,
sizeof(
long), 1, file);
424 fread(&cols,
sizeof(
int), 1, file);
425 fread(&readwhich,
sizeof(
int), 1, file);
427 if (cols !=
ndata() || rows < 0 || readwhich != which) {
429 <<
"temporary file is corrupted (which = " << which <<
" readwhich = " << readwhich <<
" rows = " << rows
430 <<
" cols = " << cols <<
")\n";
433 double *likeAChecksum =
new double[cols];
434 double *likeAChecksum2 =
new double[cols];
435 for (
int i = 0;
i < cols;
i++) {
436 likeAChecksum[
i] = 0.;
437 likeAChecksum2[
i] = 0.;
441 for (
long row = 0; row <
rows; row++) {
442 double *residual =
new double[cols];
443 fread(residual,
sizeof(
double), cols, file);
445 for (
int i = 0;
i < cols;
i++) {
446 if (fabs(residual[
i]) > likeAChecksum[i])
447 likeAChecksum[
i] = fabs(residual[i]);
448 if (fabs(residual[i]) < likeAChecksum2[i])
449 likeAChecksum2[
i] = fabs(residual[i]);
453 double *readChecksum =
new double[cols];
454 double *readChecksum2 =
new double[cols];
455 fread(readChecksum,
sizeof(
double), cols, file);
456 fread(readChecksum2,
sizeof(
double), cols, file);
458 for (
int i = 0;
i < cols;
i++) {
459 if (fabs(likeAChecksum[
i] - readChecksum[
i]) > 1
e-10 ||
460 fabs(1. / likeAChecksum2[i] - 1. / readChecksum2[i]) > 1e10) {
462 <<
"temporary file is corrupted (which = " << which <<
" rows = " << rows <<
" likeAChecksum "
463 << likeAChecksum[
i] <<
" != readChecksum " << readChecksum[
i] <<
" "
464 <<
" likeAChecksum2 " << likeAChecksum2[
i] <<
" != readChecksum2 " << readChecksum2[
i] <<
")\n";
470 delete[] likeAChecksum;
471 delete[] likeAChecksum2;
487 hist->Fill(multiplier * (*
r)[which]);
504 double weight = 1. / (*r)[whichredchi2];
505 if (TMath::Prob(1. / weight * 12, 12) < 0.99)
506 hist->Fill(multiplier * (*
r)[which], weight);
515 if (fabs((*
r)[which]) < 50.) {
521 const int n_quantiles = 7;
522 double probabilities[n_quantiles] = {0.00135, 0.02275, 0.25, 0.5, 0.75, 0.97725, 0.99865};
524 double quantiles[n_quantiles];
525 std::sort(data, data + n);
526 TMath::Quantiles(n, n_quantiles, data, quantiles, probabilities,
true,
nullptr, 7);
528 double iqr = quantiles[4] - quantiles[2];
531 double hbin = 2 * iqr /
pow(n, 1. / 3);
535 nbins = (int)((b - a) / hbin + 3.);
538 for (
int i = 0;
i < n_quantiles;
i++)
543 std::cout <<
" optimal h=" << hbin <<
" nbins=" << nbins << std::endl;
550 if (a == b || a > b) {
556 TH1D *
hist =
new TH1D(
"htmp",
"", nbins, a, b);
558 hist->Fill((*
r)[which]);
561 TF1 *
f1 =
new TF1(
"f1",
"gaus", a, b);
562 f1->SetParameter(0, hist->GetEntries());
563 f1->SetParameter(1, 0);
564 f1->SetParameter(2, hist->GetRMS());
565 hist->Fit(
"f1",
"RQ");
568 fit_mean = f1->GetParameter(1);
569 fit_sigma = f1->GetParameter(2);
570 std::cout <<
" h(" << nbins <<
"," << a <<
"," << b <<
") mu=" << fit_mean <<
" sig=" << fit_sigma << std::endl;
590 for (
int v = 0;
v < nvar;
v++) {
597 std::vector<double *>::iterator
r =
m_residuals.begin();
599 double ellipsoid_sum = 0;
600 for (
int v = 0;
v < nvar;
v++) {
606 if (ellipsoid_sum <= 1.)
622 for (
int i = 0;
i < nvar; ++
i)
633 std::cout <<
" N residuals " << nbefore <<
" ~ "
636 assert(nvar <= 10 && nvar > 0);
638 std::vector<double *>::iterator
r =
m_residuals.begin();
644 double *
data =
new double[nbefore];
645 for (
size_t i = 0;
i < nbefore;
i++)
649 re.EvaluateUni(nbefore, data, peak, sigma);
653 double distance = fabs(((*r)[vars[0]] - peak) / sigma);
654 if (distance <= nsigma)
669 TRobustEstimator re(nbefore + 1, nvar);
670 std::cout <<
"nbefore " << nbefore <<
" nvar " << nvar << std::endl;
672 std::cout <<
"+++++ JUST before loop while (r != m_residuals.end())" << std::endl;
675 double *row =
new double[nvar];
676 for (
int v = 0;
v < nvar;
v++)
677 row[
v] = (*r)[vars[
v]];
683 std::cout <<
"counter1 " << counter1 << std::endl;
684 std::cout <<
"+++++ JUST after loop while (r != m_residuals.end())" << std::endl;
686 std::cout <<
"+++++ JUST after re.Evaluate()" << std::endl;
691 TMatrixDSym Cov(nvar);
692 re.GetCovariance(Cov);
696 double conf_1d = TMath::Erf(nsigma /
sqrt(2));
697 double surf_radius =
sqrt(TMath::ChisquareQuantile(conf_1d, nvar));
703 for (
int v = 0;
v < nvar;
v++)
704 res[
v] = (*r)[vars[
v]];
706 if (distance <= surf_radius)
715 std::cout <<
" N residuals " << nbefore <<
" -> "
722 int n_station = 9999;
726 double positionX = 9999.;
727 double positionY = 9999.;
729 double chambw = 9999.;
730 double chambl = 9999.;
739 n_station = (*r)[12];
747 n_station = (*r)[10];
757 if (positionX >= xMax || positionX <= xMin)
759 if (positionY >= yMax || positionY <= yMin)
765 double dtrkchamx = (chambw / 2.) - positionX;
766 double dtrkchamy = (chambl / 2.) - positionY;
769 if (n_station == 4) {
770 if ((n_wheel == -1 && n_sector == 3) ||
773 if ((n_sector == 1 || n_sector == 2 || n_sector == 3 || n_sector == 5 || n_sector == 6 || n_sector == 7 ||
774 n_sector == 8 || n_sector == 12) &&
775 ((dtrkchamx < 40 || dtrkchamx > 380) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
777 if ((n_sector == 4 || n_sector == 13) &&
778 ((dtrkchamx < 40 || dtrkchamx > 280) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
780 if ((n_sector == 9 || n_sector == 11) &&
781 ((dtrkchamx < 40 || dtrkchamx > 180) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
783 if ((n_sector == 10 || n_sector == 14) &&
784 ((dtrkchamx < 40 || dtrkchamx > 220) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
787 if ((n_sector == 1 || n_sector == 2 || n_sector == 3 || n_sector == 5 || n_sector == 6 || n_sector == 7 ||
788 n_sector == 8 || n_sector == 12) &&
789 ((dtrkchamx < 40 || dtrkchamx > 380) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
791 if ((n_sector == 4 || n_sector == 13) &&
792 ((dtrkchamx < 40 || dtrkchamx > 280) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
794 if ((n_sector == 9 || n_sector == 11) &&
795 ((dtrkchamx < 40 || dtrkchamx > 180) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
797 if ((n_sector == 10 || n_sector == 14) &&
798 ((dtrkchamx < 40 || dtrkchamx > 220) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
802 if ((n_wheel == -1 && n_sector == 3) || (n_wheel == 1 && n_sector == 4)) {
803 if (n_station == 1 && ((dtrkchamx < 30.0 || dtrkchamx > 190.0) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
805 if (n_station == 2 && ((dtrkchamx < 30.0 || dtrkchamx > 240.0) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
807 if (n_station == 3 && ((dtrkchamx < 30.0 || dtrkchamx > 280.0) || (dtrkchamy < 40.0 || dtrkchamy > 170.0)))
810 if (n_station == 1 && ((dtrkchamx < 30.0 || dtrkchamx > 190.0) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
812 if (n_station == 2 && ((dtrkchamx < 30.0 || dtrkchamx > 240.0) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
814 if (n_station == 3 && ((dtrkchamx < 30.0 || dtrkchamx > 280.0) || (dtrkchamy < 40.0 || dtrkchamy > 210.0)))
825 double min_pt = 9999.;
827 double pt = fabs((*
r)[idx_momentum]);
832 const double bin_width = 1. / min_pt / Nbin;
835 std::vector<size_t> pos[Nbin], neg[Nbin], to_erase;
839 int bin = (int)floor(1. / fabs(
m_residuals[i][idx_momentum]) / bin_width);
841 pos[
bin].push_back(i);
843 neg[
bin].push_back(i);
847 for (
int j = 0;
j < Nbin;
j++) {
848 size_t psize = pos[
j].size();
849 size_t nsize = neg[
j].size();
853 std::set<int> idx_set;
855 while (idx_set.size() < psize - nsize)
856 idx_set.insert(gRandom->Integer(psize));
857 for (std::set<int>::iterator it = idx_set.begin(); it != idx_set.end(); it++)
858 to_erase.push_back(pos[
j][*it]);
860 while (idx_set.size() < nsize - psize)
861 idx_set.insert(gRandom->Integer(nsize));
862 for (std::set<int>::iterator it = idx_set.begin(); it != idx_set.end(); it++)
863 to_erase.push_back(neg[
j][*it]);
867 std::sort(to_erase.begin(), to_erase.end(), std::greater<int>());
868 for (std::vector<size_t>::const_iterator
e = to_erase.begin();
e != to_erase.end(); ++
e) {
874 std::vector<size_t> apos[Nbin], aneg[Nbin];
878 int bin = (int)floor(1. / fabs(
m_residuals[i][idx_momentum]) / bin_width);
880 apos[
bin].push_back(i);
882 aneg[
bin].push_back(i);
884 for (
int j = 0;
j < Nbin;
j++)
885 std::cout <<
"bin " <<
j <<
": [pos,neg] sizes before & after: [" << pos[
j].
size() <<
"," << neg[
j].size()
886 <<
"] -> [" << apos[
j].size() <<
"," << aneg[
j].size() <<
"]" << std::endl;
894 std::vector<double *>
tmp(n_ok,
nullptr);
904 std::vector<bool> tmp_ok(n_ok,
true);
constexpr int32_t ceil(float num)
void histogramChi2GaussianFit(int which, double &fit_mean, double &fit_sigma)
double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma)
double MuonResidualsFitter_integrate_pureGaussian(double low, double high, double center, double sigma)
static std::vector< std::string > checklist log
std::vector< double * >::const_iterator residuals_begin() const
Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par)
virtual void inform(TMinuit *tMinuit)=0
std::vector< double * > m_residuals
void computeHistogramRangeAndBinning(int which, int &nbins, double &a, double &b)
void selectPeakResiduals(double nsigma, int nvar, int *vars)
const int MuonResidualsFitter_numgsbins
double MuonResidualsFitter_logPureGaussian2D(double x, double y, double x0, double y0, double sx, double sy, double r)
void write(FILE *file, int which=0)
long numResiduals() const
double covarianceElement(int parNum1, int parNum2)
void eraseNotSelectedResiduals()
int parNum2parIdx(int parNum)
std::vector< double > m_error
void read(FILE *file, int which=0)
void fill(double *residual)
int residualsModel() const
Exp< T >::type exp(const T &t)
bool dofit(void(*fcn)(int &, double *, double &, double *, int), std::vector< int > &parNum, std::vector< std::string > &parName, std::vector< double > &start, std::vector< double > &step, std::vector< double > &low, std::vector< double > &high)
double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma)
double MuonResidualsFitter_compute_log_convolution(double toversigma, double gammaoversigma, double max=1000., double step=0.001, double power=4.)
std::vector< double > m_value
const int MuonResidualsFitter_numtsbins
double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma)
MuonResidualsFitter(int residualsModel, int minHits, int useResiduals, bool weightAlignment=true)
const double MuonResidualsFitter_gsbinsize
uint16_t const *__restrict__ x
std::vector< bool > m_residuals_ok
Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par)
void fix(int parNum, bool dofix=true)
T * make(const Args &...args) const
make new ROOT object
virtual void correctBField()=0
const double MuonResidualsFitter_tsbinsize
static TMinuit * MuonResidualsFitter_TMinuit
Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par)
void plotsimple(std::string name, TFileDirectory *dir, int which, double multiplier)
void fcn(int &, double *, double &, double *, int)
void fiducialCuts(double xMin=-80.0, double xMax=80.0, double yMin=-80.0, double yMax=80.0, bool fidcut1=false)
char data[epos_bytes_allocation]
void plotweighted(std::string name, TFileDirectory *dir, int which, int whichredchi2, double multiplier)
bool MuonResidualsFitter_table_initialized
double MuonResidualsFitter_lookup_table[MuonResidualsFitter_numgsbins][MuonResidualsFitter_numtsbins]
std::vector< double * >::const_iterator residuals_end() const
virtual ~MuonResidualsFitter()
void selectPeakResiduals_simple(double nsigma, int nvar, int *vars)
std::vector< bool > m_fixed
Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par)
double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma)
tuple size
Write out results.
Power< A, B >::type pow(const A &a, const B &b)