CMS 3D CMS Logo
Functions
MuonResiduals6DOFFitter.cc File Reference
#include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFFitter.h"
#include "TH2F.h"
#include "TMath.h"
#include "TTree.h"
#include "TFile.h"

Go to the source code of this file.

Functions

void MuonResiduals6DOFFitter_FCN (int &npar, double *gin, double &fval, double *par, int iflag)
 

Function Documentation

◆ MuonResiduals6DOFFitter_FCN()

void MuonResiduals6DOFFitter_FCN ( int &  npar,
double *  gin,
double &  fval,
double *  par,
int  iflag 
)

Definition at line 134 of file MuonResiduals6DOFFitter.cc.

References cms::cuda::assert(), MuonResidualsFitterFitInfo::fitter(), MuonResidualsFitter::k0010, MuonResidualsFitter::k1010, MuonResidualsFitter::k1100, MuonResidualsFitter::k1110, MuonResidualsFitter::k1111, MuonResiduals6DOFFitter::kAlignPhiX, MuonResiduals6DOFFitter::kAlignPhiY, MuonResiduals6DOFFitter::kAlignPhiZ, MuonResiduals6DOFFitter::kAlignX, MuonResiduals6DOFFitter::kAlignY, MuonResiduals6DOFFitter::kAlignZ, MuonResiduals6DOFFitter::kAlphaX, MuonResiduals6DOFFitter::kAlphaY, MuonResiduals6DOFFitter::kAngleX, MuonResiduals6DOFFitter::kAngleY, MuonResidualsFitter::kGaussPowerTails, MuonResiduals6DOFFitter::kPositionX, MuonResiduals6DOFFitter::kPositionY, MuonResidualsFitter::kPowerLawTails, MuonResidualsFitter::kPureGaussian, MuonResidualsFitter::kPureGaussian2D, MuonResiduals6DOFFitter::kRedChi2, MuonResiduals6DOFFitter::kResidX, MuonResiduals6DOFFitter::kResidXGamma, MuonResiduals6DOFFitter::kResidXSigma, MuonResiduals6DOFFitter::kResidY, MuonResiduals6DOFFitter::kResidYGamma, MuonResiduals6DOFFitter::kResidYSigma, MuonResiduals6DOFFitter::kResSlopeX, MuonResiduals6DOFFitter::kResSlopeXGamma, MuonResiduals6DOFFitter::kResSlopeXSigma, MuonResiduals6DOFFitter::kResSlopeY, MuonResiduals6DOFFitter::kResSlopeYGamma, MuonResiduals6DOFFitter::kResSlopeYSigma, MuonResidualsFitter::kROOTVoigt, MuonResidualsFitter_logGaussPowerTails(), MuonResidualsFitter_logPowerLawTails(), MuonResidualsFitter_logPureGaussian(), MuonResidualsFitter_logPureGaussian2D(), MuonResidualsFitter_logROOTVoigt(), MuonResidualsFitter::residuals_begin(), MuonResidualsFitter::residuals_end(), MuonResidualsFitter::residualsModel(), and MuonResidualsFitter::useRes().

Referenced by MuonResiduals6DOFFitter::fit().

134  {
135  MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo *)(minuit->GetObjectFit());
136  MuonResidualsFitter *fitter = fitinfo->fitter();
137 
138  fval = 0.;
139  for (std::vector<double *>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end();
140  ++resiter) {
141  const double residX = (*resiter)[MuonResiduals6DOFFitter::kResidX];
142  const double residY = (*resiter)[MuonResiduals6DOFFitter::kResidY];
143  const double resslopeX = (*resiter)[MuonResiduals6DOFFitter::kResSlopeX];
144  const double resslopeY = (*resiter)[MuonResiduals6DOFFitter::kResSlopeY];
145  const double positionX = (*resiter)[MuonResiduals6DOFFitter::kPositionX];
146  const double positionY = (*resiter)[MuonResiduals6DOFFitter::kPositionY];
147  const double angleX = (*resiter)[MuonResiduals6DOFFitter::kAngleX];
148  const double angleY = (*resiter)[MuonResiduals6DOFFitter::kAngleY];
149  const double redchi2 = (*resiter)[MuonResiduals6DOFFitter::kRedChi2];
150 
151  const double alignx = par[MuonResiduals6DOFFitter::kAlignX];
152  const double aligny = par[MuonResiduals6DOFFitter::kAlignY];
153  const double alignz = par[MuonResiduals6DOFFitter::kAlignZ];
154  const double alignphix = par[MuonResiduals6DOFFitter::kAlignPhiX];
155  const double alignphiy = par[MuonResiduals6DOFFitter::kAlignPhiY];
156  const double alignphiz = par[MuonResiduals6DOFFitter::kAlignPhiZ];
157  const double resXsigma = par[MuonResiduals6DOFFitter::kResidXSigma];
158  const double resYsigma = par[MuonResiduals6DOFFitter::kResidYSigma];
159  const double slopeXsigma = par[MuonResiduals6DOFFitter::kResSlopeXSigma];
160  const double slopeYsigma = par[MuonResiduals6DOFFitter::kResSlopeYSigma];
161  const double alphax = par[MuonResiduals6DOFFitter::kAlphaX];
162  const double alphay = par[MuonResiduals6DOFFitter::kAlphaY];
163  const double resXgamma = par[MuonResiduals6DOFFitter::kResidXGamma];
164  const double resYgamma = par[MuonResiduals6DOFFitter::kResidYGamma];
165  const double slopeXgamma = par[MuonResiduals6DOFFitter::kResSlopeXGamma];
166  const double slopeYgamma = par[MuonResiduals6DOFFitter::kResSlopeYGamma];
167 
168  double coefX = alphax, coefY = alphay;
169  if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian ||
170  fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D)
171  coefX = coefY = 0.;
172  double residXpeak = residual_x(
173  alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coefX, resslopeX);
174  double residYpeak = residual_y(
175  alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coefY, resslopeY);
176  double slopeXpeak =
177  residual_dxdz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
178  double slopeYpeak =
179  residual_dydz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
180 
181  double weight = (1. / redchi2) * number_of_hits / sum_of_weights;
182  if (!weight_alignment)
183  weight = 1.;
184 
185  if (!weight_alignment || TMath::Prob(redchi2 * 12, 12) < 0.99) // no spikes allowed
186  {
187  if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
188  if (fitter->useRes() == MuonResidualsFitter::k1111) {
189  fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
190  fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
191  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
192  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeY, slopeYpeak, slopeYsigma);
193  } else if (fitter->useRes() == MuonResidualsFitter::k1110) {
194  fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
195  fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
196  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
197  } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
198  fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
199  fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
200  } else if (fitter->useRes() == MuonResidualsFitter::k1010) {
201  fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
202  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
203  } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
204  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
205  }
206  //std::cout<<"FCNx("<<residX<<","<<residXpeak<<","<<resXsigma<<") = "<<MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma)<<std::endl;
207  //std::cout<<"FCNy("<<residY<<","<<residYpeak<<","<<resYsigma<<") = "<<MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma)<<std::endl;
208  //std::cout<<"FCNsx("<<resslopeX<<","<<slopeXpeak<<","<<slopeXsigma<<") = "<<MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma)<<std::endl;
209  //std::cout<<"FCNsy("<<resslopeY<<","<<slopeYpeak<<","<<slopeYsigma<<") = "<<MuonResidualsFitter_logPureGaussian(resslopeY, slopeYpeak, slopeYsigma)<<std::endl;
210  } else if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D) {
211  if (fitter->useRes() == MuonResidualsFitter::k1111) {
212  fval += -weight * MuonResidualsFitter_logPureGaussian2D(
213  residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
214  fval += -weight * MuonResidualsFitter_logPureGaussian2D(
215  residY, resslopeY, residYpeak, slopeYpeak, resYsigma, slopeYsigma, alphay);
216  } else if (fitter->useRes() == MuonResidualsFitter::k1110) {
217  fval += -weight * MuonResidualsFitter_logPureGaussian2D(
218  residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
219  fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
220  } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
221  fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
222  fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
223  } else if (fitter->useRes() == MuonResidualsFitter::k1010) {
224  fval += -weight * MuonResidualsFitter_logPureGaussian2D(
225  residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
226  } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
227  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
228  }
229  //std::cout<<"FCNx("<<residX<<","<<resslopeX<<","<<residXpeak<<","<<slopeXpeak<<","<<resXsigma<<","<<slopeXsigma<<","<<alphax<<") = "<<MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax)<<std::endl;
230  //std::cout<<"FCNy("<<residY<<","<<resslopeY<<","<<residYpeak<<","<<slopeYpeak<<","<<resYsigma<<","<<slopeYsigma<<","<<alphay<<") = "<<MuonResidualsFitter_logPureGaussian2D(residY, resslopeY, residYpeak, slopeYpeak, resYsigma, slopeYsigma, alphay)<<std::endl;
231  } else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
232  fval += -weight * MuonResidualsFitter_logPowerLawTails(residX, residXpeak, resXsigma, resXgamma);
233  fval += -weight * MuonResidualsFitter_logPowerLawTails(residY, residYpeak, resYsigma, resYgamma);
234  fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
235  fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
236  } else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
237  fval += -weight * MuonResidualsFitter_logROOTVoigt(residX, residXpeak, resXsigma, resXgamma);
238  fval += -weight * MuonResidualsFitter_logROOTVoigt(residY, residYpeak, resYsigma, resYgamma);
239  fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
240  fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
241  } else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
242  fval += -weight * MuonResidualsFitter_logGaussPowerTails(residX, residXpeak, resXsigma);
243  fval += -weight * MuonResidualsFitter_logGaussPowerTails(residY, residYpeak, resYsigma);
244  fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeX, slopeXpeak, slopeXsigma);
245  fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeY, slopeYpeak, slopeYsigma);
246  } else {
247  assert(false);
248  }
249  }
250  }
251  /*
252  const double resXsigma = par[MuonResiduals6DOFFitter::kResidXSigma];
253  const double resYsigma = par[MuonResiduals6DOFFitter::kResidYSigma];
254  const double slopeXsigma = par[MuonResiduals6DOFFitter::kResSlopeXSigma];
255  const double slopeYsigma = par[MuonResiduals6DOFFitter::kResSlopeYSigma];
256  const double alphax = par[MuonResiduals6DOFFitter::kAlphaX];
257  const double alphay = par[MuonResiduals6DOFFitter::kAlphaY];
258  std::cout<<"fval="<<fval<<","<<resXsigma<<","<<slopeXsigma<<","<<alphax<<","<<resYsigma<<","<<slopeYsigma<<","<<alphay<<std::endl;
259 */
260 }
MuonResidualsFitter_logGaussPowerTails
double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma)
Definition: MuonResidualsFitter.cc:128
MuonResidualsFitter_logPureGaussian2D
double MuonResidualsFitter_logPureGaussian2D(double x, double y, double x0, double y0, double sx, double sy, double r)
Definition: MuonResidualsFitter.cc:43
MuonResidualsFitter::residuals_begin
std::vector< double * >::const_iterator residuals_begin() const
Definition: MuonResidualsFitter.h:190
MuonResidualsFitterFitInfo::fitter
MuonResidualsFitter * fitter()
Definition: MuonResidualsFitter.h:253
weight
Definition: weight.py:1
MuonResidualsFitter::residuals_end
std::vector< double * >::const_iterator residuals_end() const
Definition: MuonResidualsFitter.h:191
cms::cuda::assert
assert(be >=bs)
MuonResidualsFitter_logROOTVoigt
double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma)
Definition: MuonResidualsFitter.cc:120
MuonResidualsFitter_logPowerLawTails
double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma)
Definition: MuonResidualsFitter.cc:85
MuonResidualsFitter::useRes
int useRes(int pattern=-1)
Definition: MuonResidualsFitter.h:119
MuonResidualsFitter_logPureGaussian
double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma)
Definition: MuonResidualsFitter.cc:31
Generated for CMSSW Reference Manual by doxygen 1.8.14