14 double MuonResiduals5DOFFitter_residual(
double delta_x,
double delta_z,
double delta_phix,
double delta_phiy,
double delta_phiz,
double track_x,
double track_y,
double track_dxdz,
double track_dydz,
double alpha,
double resslope) {
15 return delta_x - track_dxdz * delta_z - track_y * track_dxdz * delta_phix + track_x * track_dxdz * delta_phiy - track_y * delta_phiz + resslope *
alpha;
18 double MuonResiduals5DOFFitter_resslope(
double delta_x,
double delta_z,
double delta_phix,
double delta_phiy,
double delta_phiz,
double track_x,
double track_y,
double track_dxdz,
double track_dydz) {
19 return -track_dxdz * track_dydz * delta_phix + (1. + track_dxdz * track_dxdz) * delta_phiy - track_dydz * delta_phiz;
22 Double_t
MuonResiduals5DOFFitter_residual_trackx_TF1(Double_t *xvec, Double_t *
par) {
return 10.*
MuonResiduals5DOFFitter_residual(par[0], par[2], par[3], par[4], par[5], xvec[0], par[7], par[8], par[9], par[10], par[11]); }
23 Double_t
MuonResiduals5DOFFitter_residual_tracky_TF1(Double_t *xvec, Double_t *
par) {
return 10.*
MuonResiduals5DOFFitter_residual(par[0], par[2], par[3], par[4], par[5], par[6], xvec[0], par[8], par[9], par[10], par[11]); }
24 Double_t
MuonResiduals5DOFFitter_residual_trackdxdz_TF1(Double_t *xvec, Double_t *
par) {
return 10.*
MuonResiduals5DOFFitter_residual(par[0], par[2], par[3], par[4], par[5], par[6], par[7], xvec[0], par[9], par[10], par[11]); }
25 Double_t
MuonResiduals5DOFFitter_residual_trackdydz_TF1(Double_t *xvec, Double_t *
par) {
return 10.*
MuonResiduals5DOFFitter_residual(par[0], par[2], par[3], par[4], par[5], par[6], par[7], par[8], xvec[0], par[10], par[11]); }
57 double residpeak =
MuonResiduals5DOFFitter_residual(alignx, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, alpha, resslope);
81 else { assert(
false); }
94 if (TMath::Prob(redchi2*8, 8) < 0.99) {
111 double resid_mean = 0;
112 double resslope_mean = 0;
113 double resid_stdev = 0.5;
114 double resslope_stdev = 0.005;
115 double alpha_estimate = 0;
117 std::vector<int>
num;
118 std::vector<std::string>
name;
119 std::vector<double> start;
120 std::vector<double>
step;
121 std::vector<double> low;
122 std::vector<double> high;
125 num.push_back(
kAlignX); name.push_back(std::string(
"AlignX")); start.push_back(0.); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
127 num.push_back(
kAlignX); name.push_back(std::string(
"AlignX")); start.push_back(resid_mean); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
129 num.push_back(
kAlignZ); name.push_back(std::string(
"AlignZ")); start.push_back(0.); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
130 num.push_back(
kAlignPhiX); name.push_back(std::string(
"AlignPhiX")); start.push_back(0.); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
132 num.push_back(
kAlignPhiY); name.push_back(std::string(
"AlignPhiY")); start.push_back(0.); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
134 num.push_back(
kAlignPhiY); name.push_back(std::string(
"AlignPhiY")); start.push_back(resslope_mean); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
136 num.push_back(
kAlignPhiZ); name.push_back(std::string(
"AlignPhiZ")); start.push_back(0.); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
137 num.push_back(
kResidSigma); name.push_back(std::string(
"ResidSigma")); start.push_back(resid_stdev); step.push_back(0.01*resid_stdev); low.push_back(0.); high.push_back(0.);
138 num.push_back(
kResSlopeSigma); name.push_back(std::string(
"ResSlopeSigma")); start.push_back(resslope_stdev); step.push_back(0.01*resslope_stdev); low.push_back(0.); high.push_back(0.);
139 num.push_back(
kAlpha); name.push_back(std::string(
"Alpha")); start.push_back(alpha_estimate); step.push_back(0.01*resslope_stdev); low.push_back(0.); high.push_back(0.);
141 num.push_back(
kResidGamma); name.push_back(std::string(
"ResidGamma")); start.push_back(0.1*resid_stdev); step.push_back(0.01*resid_stdev); low.push_back(0.); high.push_back(0.);
142 num.push_back(
kResSlopeGamma); name.push_back(std::string(
"ResSlopeGamma")); start.push_back(0.1*resslope_stdev); step.push_back(0.01*resslope_stdev); low.push_back(0.); high.push_back(0.);
151 std::stringstream name_residual, name_resslope, name_residual_raw, name_resslope_raw, name_residual_cut, name_alpha;
152 std::stringstream name_residual_trackx, name_resslope_trackx;
153 std::stringstream name_residual_tracky, name_resslope_tracky;
154 std::stringstream name_residual_trackdxdz, name_resslope_trackdxdz;
155 std::stringstream name_residual_trackdydz, name_resslope_trackdydz;
157 name_residual << name <<
"_residual";
158 name_resslope << name <<
"_resslope";
159 name_residual_raw << name <<
"_residual_raw";
160 name_resslope_raw << name <<
"_resslope_raw";
161 name_residual_cut << name <<
"_residual_cut";
162 name_alpha << name <<
"_alpha";
163 name_residual_trackx << name <<
"_residual_trackx";
164 name_resslope_trackx << name <<
"_resslope_trackx";
165 name_residual_tracky << name <<
"_residual_tracky";
166 name_resslope_tracky << name <<
"_resslope_tracky";
167 name_residual_trackdxdz << name <<
"_residual_trackdxdz";
168 name_resslope_trackdxdz << name <<
"_resslope_trackdxdz";
169 name_residual_trackdydz << name <<
"_residual_trackdydz";
170 name_resslope_trackdydz << name <<
"_resslope_trackdydz";
174 double min_residual = -100.;
double max_residual = 100.;
175 double min_resslope = -100.;
double max_resslope = 100.;
176 double min_trackx = -width/2.;
double max_trackx = width/2.;
177 double min_tracky = -length/2.;
double max_tracky = length/2.;
178 double min_trackdxdz = -1.5;
double max_trackdxdz = 1.5;
179 double min_trackdydz = -1.5;
double max_trackdydz = 1.5;
181 TH1F *hist_residual = dir->
make<TH1F>(name_residual.str().c_str(),
"", 100, min_residual, max_residual);
182 TH1F *hist_resslope = dir->
make<TH1F>(name_resslope.str().c_str(),
"", 100, min_resslope, max_resslope);
183 TH1F *hist_residual_raw = dir->
make<TH1F>(name_residual_raw.str().c_str(),
"", 100, min_residual, max_residual);
184 TH1F *hist_resslope_raw = dir->
make<TH1F>(name_resslope_raw.str().c_str(),
"", 100, min_resslope, max_resslope);
185 TH1F *hist_residual_cut = dir->
make<TH1F>(name_residual_cut.str().c_str(),
"", 100, min_residual, max_residual);
186 TH2F *hist_alpha = dir->
make<TH2F>(name_alpha.str().c_str(),
"", 40, min_resslope, max_resslope, 40, -20., 20.);
187 TProfile *hist_residual_trackx = dir->
make<TProfile>(name_residual_trackx.str().c_str(),
"", 100, min_trackx, max_trackx, min_residual, max_residual);
188 TProfile *hist_resslope_trackx = dir->
make<TProfile>(name_resslope_trackx.str().c_str(),
"", 100, min_trackx, max_trackx, min_resslope, max_resslope);
189 TProfile *hist_residual_tracky = dir->
make<TProfile>(name_residual_tracky.str().c_str(),
"", 100, min_tracky, max_tracky, min_residual, max_residual);
190 TProfile *hist_resslope_tracky = dir->
make<TProfile>(name_resslope_tracky.str().c_str(),
"", 100, min_tracky, max_tracky, min_resslope, max_resslope);
191 TProfile *hist_residual_trackdxdz = dir->
make<TProfile>(name_residual_trackdxdz.str().c_str(),
"", 500, min_trackdxdz, max_trackdxdz, min_residual, max_residual);
192 TProfile *hist_resslope_trackdxdz = dir->
make<TProfile>(name_resslope_trackdxdz.str().c_str(),
"", 500, min_trackdxdz, max_trackdxdz, min_resslope, max_resslope);
193 TProfile *hist_residual_trackdydz = dir->
make<TProfile>(name_residual_trackdydz.str().c_str(),
"", 500, min_trackdydz, max_trackdydz, min_residual, max_residual);
194 TProfile *hist_resslope_trackdydz = dir->
make<TProfile>(name_resslope_trackdydz.str().c_str(),
"", 500, min_trackdydz, max_trackdydz, min_resslope, max_resslope);
196 hist_residual_trackx->SetAxisRange(-10., 10.,
"Y");
197 hist_resslope_trackx->SetAxisRange(-10., 10.,
"Y");
198 hist_residual_tracky->SetAxisRange(-10., 10.,
"Y");
199 hist_resslope_tracky->SetAxisRange(-10., 10.,
"Y");
200 hist_residual_trackdxdz->SetAxisRange(-10., 10.,
"Y");
201 hist_resslope_trackdxdz->SetAxisRange(-10., 10.,
"Y");
202 hist_residual_trackdydz->SetAxisRange(-10., 10.,
"Y");
203 hist_resslope_trackdydz->SetAxisRange(-10., 10.,
"Y");
205 name_residual <<
"_fit";
206 name_resslope <<
"_fit";
207 name_alpha <<
"_fit";
208 name_residual_trackx <<
"_fit";
209 name_resslope_trackx <<
"_fit";
210 name_residual_tracky <<
"_fit";
211 name_resslope_tracky <<
"_fit";
212 name_residual_trackdxdz <<
"_fit";
213 name_resslope_trackdxdz <<
"_fit";
214 name_residual_trackdydz <<
"_fit";
215 name_resslope_trackdydz <<
"_fit";
217 TF1 *fit_residual =
NULL;
218 TF1 *fit_resslope =
NULL;
243 else { assert(
false); }
245 fit_residual->SetLineColor(2); fit_residual->SetLineWidth(2);
246 fit_resslope->SetLineColor(2); fit_resslope->SetLineWidth(2);
247 fit_residual->Write();
248 fit_resslope->Write();
250 TF1 *fit_alpha =
new TF1(name_alpha.str().c_str(),
"[0] + x*[1]", min_resslope, max_resslope);
252 fit_alpha->SetLineColor(2); fit_alpha->SetLineWidth(2);
255 TProfile *fit_residual_trackx = dir->
make<TProfile>(name_residual_trackx.str().c_str(),
"", 100, min_trackx, max_trackx);
256 TProfile *fit_resslope_trackx = dir->
make<TProfile>(name_resslope_trackx.str().c_str(),
"", 100, min_trackx, max_trackx);
257 TProfile *fit_residual_tracky = dir->
make<TProfile>(name_residual_tracky.str().c_str(),
"", 100, min_tracky, max_tracky);
258 TProfile *fit_resslope_tracky = dir->
make<TProfile>(name_resslope_tracky.str().c_str(),
"", 100, min_tracky, max_tracky);
259 TProfile *fit_residual_trackdxdz = dir->
make<TProfile>(name_residual_trackdxdz.str().c_str(),
"", 500, min_trackdxdz, max_trackdxdz);
260 TProfile *fit_resslope_trackdxdz = dir->
make<TProfile>(name_resslope_trackdxdz.str().c_str(),
"", 500, min_trackdxdz, max_trackdxdz);
261 TProfile *fit_residual_trackdydz = dir->
make<TProfile>(name_residual_trackdydz.str().c_str(),
"", 500, min_trackdydz, max_trackdydz);
262 TProfile *fit_resslope_trackdydz = dir->
make<TProfile>(name_resslope_trackdydz.str().c_str(),
"", 500, min_trackdydz, max_trackdydz);
264 fit_residual_trackx->SetLineColor(2); fit_residual_trackx->SetLineWidth(2);
265 fit_resslope_trackx->SetLineColor(2); fit_resslope_trackx->SetLineWidth(2);
266 fit_residual_tracky->SetLineColor(2); fit_residual_tracky->SetLineWidth(2);
267 fit_resslope_tracky->SetLineColor(2); fit_resslope_tracky->SetLineWidth(2);
268 fit_residual_trackdxdz->SetLineColor(2); fit_residual_trackdxdz->SetLineWidth(2);
269 fit_resslope_trackdxdz->SetLineColor(2); fit_resslope_trackdxdz->SetLineWidth(2);
270 fit_residual_trackdydz->SetLineColor(2); fit_residual_trackdydz->SetLineWidth(2);
271 fit_resslope_trackdydz->SetLineColor(2); fit_resslope_trackdydz->SetLineWidth(2);
273 name_residual_trackx <<
"line";
274 name_resslope_trackx <<
"line";
275 name_residual_tracky <<
"line";
276 name_resslope_tracky <<
"line";
277 name_residual_trackdxdz <<
"line";
278 name_resslope_trackdxdz <<
"line";
279 name_residual_trackdydz <<
"line";
280 name_resslope_trackdydz <<
"line";
291 double sum_resslope = 0.;
292 double sum_trackx = 0.;
293 double sum_tracky = 0.;
294 double sum_trackdxdz = 0.;
295 double sum_trackdydz = 0.;
303 double weight = 1./redchi2;
307 sum_resslope += weight * resslope;
308 sum_trackx += weight * positionX;
309 sum_tracky += weight * positionY;
310 sum_trackdxdz += weight * angleX;
311 sum_trackdydz += weight * angleY;
320 double fitparameters[12];
322 fitparameters[1] = 0.;
327 fitparameters[6] = mean_trackx;
328 fitparameters[7] = mean_tracky;
329 fitparameters[8] = mean_trackdxdz;
330 fitparameters[9] = mean_trackdydz;
332 fitparameters[11] = mean_resslope;
334 fitline_residual_trackx->SetParameters(fitparameters);
335 fitline_resslope_trackx->SetParameters(fitparameters);
336 fitline_residual_tracky->SetParameters(fitparameters);
337 fitline_resslope_tracky->SetParameters(fitparameters);
338 fitline_residual_trackdxdz->SetParameters(fitparameters);
339 fitline_resslope_trackdxdz->SetParameters(fitparameters);
340 fitline_residual_trackdydz->SetParameters(fitparameters);
341 fitline_resslope_trackdydz->SetParameters(fitparameters);
343 fitline_residual_trackx->SetLineColor(2); fitline_residual_trackx->SetLineWidth(2);
344 fitline_resslope_trackx->SetLineColor(2); fitline_resslope_trackx->SetLineWidth(2);
345 fitline_residual_tracky->SetLineColor(2); fitline_residual_tracky->SetLineWidth(2);
346 fitline_resslope_tracky->SetLineColor(2); fitline_resslope_tracky->SetLineWidth(2);
347 fitline_residual_trackdxdz->SetLineColor(2); fitline_residual_trackdxdz->SetLineWidth(2);
348 fitline_resslope_trackdxdz->SetLineColor(2); fitline_resslope_trackdxdz->SetLineWidth(2);
349 fitline_residual_trackdydz->SetLineColor(2); fitline_residual_trackdydz->SetLineWidth(2);
350 fitline_resslope_trackdydz->SetLineColor(2); fitline_resslope_trackdydz->SetLineWidth(2);
352 fitline_residual_trackx->Write();
353 fitline_resslope_trackx->Write();
354 fitline_residual_tracky->Write();
355 fitline_resslope_tracky->Write();
356 fitline_residual_trackdxdz->Write();
357 fitline_resslope_trackdxdz->Write();
358 fitline_residual_trackdydz->Write();
359 fitline_resslope_trackdydz->Write();
369 double weight = 1./redchi2;
373 hist_alpha->Fill(1000.*resslope, 10.*resid);
375 double geom_resid =
MuonResiduals5DOFFitter_residual(
value(
kAlignX),
value(
kAlignZ),
value(
kAlignPhiX),
value(
kAlignPhiY),
value(
kAlignPhiZ), positionX, positionY, angleX, angleY,
value(
kAlpha), resslope);
376 hist_residual->Fill(10.*(resid - geom_resid +
value(
kAlignX)), weight);
377 hist_residual_trackx->Fill(positionX, 10.*resid, weight);
378 hist_residual_tracky->Fill(positionY, 10.*resid, weight);
379 hist_residual_trackdxdz->Fill(angleX, 10.*resid, weight);
380 hist_residual_trackdydz->Fill(angleY, 10.*resid, weight);
381 fit_residual_trackx->Fill(positionX, 10.*geom_resid, weight);
382 fit_residual_tracky->Fill(positionY, 10.*geom_resid, weight);
383 fit_residual_trackdxdz->Fill(angleX, 10.*geom_resid, weight);
384 fit_residual_trackdydz->Fill(angleY, 10.*geom_resid, weight);
387 hist_resslope->Fill(1000.*(resslope - geom_resslope +
value(
kAlignPhiY)), weight);
388 hist_resslope_trackx->Fill(positionX, 1000.*resslope, weight);
389 hist_resslope_tracky->Fill(positionY, 1000.*resslope, weight);
390 hist_resslope_trackdxdz->Fill(angleX, 1000.*resslope, weight);
391 hist_resslope_trackdydz->Fill(angleY, 1000.*resslope, weight);
392 fit_resslope_trackx->Fill(positionX, 1000.*geom_resslope, weight);
393 fit_resslope_tracky->Fill(positionY, 1000.*geom_resslope, weight);
394 fit_resslope_trackdxdz->Fill(angleX, 1000.*geom_resslope, weight);
395 fit_resslope_trackdydz->Fill(angleY, 1000.*geom_resslope, weight);
398 hist_residual_raw->Fill(10.*resid);
399 hist_resslope_raw->Fill(1000.*resslope);
400 if (fabs(resslope) < 0.005) hist_residual_cut->Fill(10.*resid);
405 for (
int i = 1;
i <= hist_residual->GetNbinsX();
i++) {
406 double xi = hist_residual->GetBinCenter(
i);
407 double yi = hist_residual->GetBinContent(
i);
408 double yerri = hist_residual->GetBinError(
i);
409 double yth = fit_residual->Eval(xi);
411 chi2 +=
pow((yth - yi)/yerri, 2);
415 for (
int i = 1;
i <= hist_resslope->GetNbinsX();
i++) {
416 double xi = hist_resslope->GetBinCenter(
i);
417 double yi = hist_resslope->GetBinContent(
i);
418 double yerri = hist_resslope->GetBinError(
i);
419 double yth = fit_resslope->Eval(xi);
421 chi2 +=
pow((yth - yi)/yerri, 2);
427 return (ndof > 0. ? chi2 / ndof : -1.);
align::Scalar width() const
double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma)
Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par)
Double_t MuonResiduals5DOFFitter_resslope_trackdydz_TF1(Double_t *xvec, Double_t *par)
static bool MuonResiduals5DOFFitter_weightAlignment
MuonResidualsFitter * fitter()
int residualsModel() const
Double_t MuonResiduals5DOFFitter_residual_trackdxdz_TF1(Double_t *xvec, Double_t *par)
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_t MuonResiduals5DOFFitter_resslope_trackx_TF1(Double_t *xvec, Double_t *par)
double MuonResiduals5DOFFitter_residual(double delta_x, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz, double alpha, double resslope)
Double_t MuonResiduals5DOFFitter_resslope_tracky_TF1(Double_t *xvec, Double_t *par)
Double_t MuonResiduals5DOFFitter_resslope_trackdxdz_TF1(Double_t *xvec, Double_t *par)
void MuonResiduals5DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma)
Double_t MuonResiduals5DOFFitter_residual_tracky_TF1(Double_t *xvec, Double_t *par)
Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par)
static double MuonResiduals5DOFFitter_sum_of_weights
double MuonResiduals5DOFFitter_resslope(double delta_x, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz)
const AlignableSurface & surface() const
Return the Surface (global position and orientation) of the object.
Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par)
Double_t MuonResiduals5DOFFitter_residual_trackx_TF1(Double_t *xvec, Double_t *par)
align::Scalar length() const
Double_t MuonResiduals5DOFFitter_residual_trackdydz_TF1(Double_t *xvec, Double_t *par)
double plot(std::string name, TFileDirectory *dir, Alignable *ali)
std::vector< double * >::const_iterator residuals_end() const
static TMinuit * MuonResiduals5DOFFitter_TMinuit
static double MuonResiduals5DOFFitter_number_of_hits
T * make() const
make new ROOT object
void inform(TMinuit *tMinuit)
Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par)
std::vector< double * >::const_iterator residuals_begin() const
double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma)
Power< A, B >::type pow(const A &a, const B &b)
const double par[8 *NPar][4]