CMS 3D CMS Logo
MuonResiduals6DOFFitter.cc
Go to the documentation of this file.
1 // $Id:$
2 
3 #ifdef STANDALONE_FITTER
4 #include "MuonResiduals6DOFFitter.h"
5 #else
6 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFFitter.h"
7 #endif
8 
9 #include "TH2F.h"
10 #include "TMath.h"
11 #include "TTree.h"
12 #include "TFile.h"
13 
14 namespace {
15  TMinuit *minuit;
16 
17  double sum_of_weights;
18  double number_of_hits;
19  bool weight_alignment;
20 
21  double residual_x(double delta_x,
22  double delta_y,
23  double delta_z,
24  double delta_phix,
25  double delta_phiy,
26  double delta_phiz,
27  double track_x,
28  double track_y,
29  double track_dxdz,
30  double track_dydz,
31  double alphax,
32  double residual_dxdz) {
33  return delta_x - track_dxdz * delta_z - track_y * track_dxdz * delta_phix + track_x * track_dxdz * delta_phiy -
34  track_y * delta_phiz + residual_dxdz * alphax;
35  }
36 
37  double residual_y(double delta_x,
38  double delta_y,
39  double delta_z,
40  double delta_phix,
41  double delta_phiy,
42  double delta_phiz,
43  double track_x,
44  double track_y,
45  double track_dxdz,
46  double track_dydz,
47  double alphay,
48  double residual_dydz) {
49  return delta_y - track_dydz * delta_z - track_y * track_dydz * delta_phix + track_x * track_dydz * delta_phiy +
50  track_x * delta_phiz + residual_dydz * alphay;
51  }
52 
53  double residual_dxdz(double delta_x,
54  double delta_y,
55  double delta_z,
56  double delta_phix,
57  double delta_phiy,
58  double delta_phiz,
59  double track_x,
60  double track_y,
61  double track_dxdz,
62  double track_dydz) {
63  return -track_dxdz * track_dydz * delta_phix + (1. + track_dxdz * track_dxdz) * delta_phiy -
64  track_dydz * delta_phiz;
65  }
66 
67  double residual_dydz(double delta_x,
68  double delta_y,
69  double delta_z,
70  double delta_phix,
71  double delta_phiy,
72  double delta_phiz,
73  double track_x,
74  double track_y,
75  double track_dxdz,
76  double track_dydz) {
77  return -(1. + track_dydz * track_dydz) * delta_phix + track_dxdz * track_dydz * delta_phiy +
78  track_dxdz * delta_phiz;
79  }
80 
81  Double_t residual_x_trackx_TF1(Double_t *xx, Double_t *p) {
82  return 10. * residual_x(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9], p[10], p[11]);
83  }
84  Double_t residual_x_tracky_TF1(Double_t *xx, Double_t *p) {
85  return 10. * residual_x(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9], p[10], p[11]);
86  }
87  Double_t residual_x_trackdxdz_TF1(Double_t *xx, Double_t *p) {
88  return 10. * residual_x(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9], p[10], p[11]);
89  }
90  Double_t residual_x_trackdydz_TF1(Double_t *xx, Double_t *p) {
91  return 10. * residual_x(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0], p[10], p[11]);
92  }
93 
94  Double_t residual_y_trackx_TF1(Double_t *xx, Double_t *p) {
95  return 10. * residual_y(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9], p[12], p[13]);
96  }
97  Double_t residual_y_tracky_TF1(Double_t *xx, Double_t *p) {
98  return 10. * residual_y(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9], p[12], p[13]);
99  }
100  Double_t residual_y_trackdxdz_TF1(Double_t *xx, Double_t *p) {
101  return 10. * residual_y(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9], p[12], p[13]);
102  }
103  Double_t residual_y_trackdydz_TF1(Double_t *xx, Double_t *p) {
104  return 10. * residual_y(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0], p[12], p[13]);
105  }
106 
107  Double_t residual_dxdz_trackx_TF1(Double_t *xx, Double_t *p) {
108  return 1000. * residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9]);
109  }
110  Double_t residual_dxdz_tracky_TF1(Double_t *xx, Double_t *p) {
111  return 1000. * residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9]);
112  }
113  Double_t residual_dxdz_trackdxdz_TF1(Double_t *xx, Double_t *p) {
114  return 1000. * residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9]);
115  }
116  Double_t residual_dxdz_trackdydz_TF1(Double_t *xx, Double_t *p) {
117  return 1000. * residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0]);
118  }
119 
120  Double_t residual_dydz_trackx_TF1(Double_t *xx, Double_t *p) {
121  return 1000. * residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9]);
122  }
123  Double_t residual_dydz_tracky_TF1(Double_t *xx, Double_t *p) {
124  return 1000. * residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9]);
125  }
126  Double_t residual_dydz_trackdxdz_TF1(Double_t *xx, Double_t *p) {
127  return 1000. * residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9]);
128  }
129  Double_t residual_dydz_trackdydz_TF1(Double_t *xx, Double_t *p) {
130  return 1000. * residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0]);
131  }
132 } // namespace
133 
134 void MuonResiduals6DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
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 }
261 
262 void MuonResiduals6DOFFitter::inform(TMinuit *tMinuit) { minuit = tMinuit; }
263 
264 void MuonResiduals6DOFFitter::correctBField() { MuonResidualsFitter::correctBField(kPt, kCharge); }
265 
266 double MuonResiduals6DOFFitter::sumofweights() {
267  sum_of_weights = 0.;
268  number_of_hits = 0.;
269  weight_alignment = m_weightAlignment;
270  for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
271  if (m_weightAlignment) {
272  double redchi2 = (*resiter)[kRedChi2];
273  if (TMath::Prob(redchi2 * 12, 12) < 0.99) { // no spikes allowed
274  sum_of_weights += 1. / redchi2;
275  number_of_hits += 1.;
276  }
277  } else {
278  sum_of_weights += 1.;
279  number_of_hits += 1.;
280  }
281  }
282  return sum_of_weights;
283 }
284 
285 bool MuonResiduals6DOFFitter::fit(Alignable *ali) {
286  initialize_table(); // if not already initialized
287  sumofweights();
288 
289  double resx_std = 0.5;
290  double resy_std = 3.0;
291  double resslopex_std = 0.002;
292  double resslopey_std = 0.005;
293 
294  int nums[16] = {kAlignX,
295  kAlignY,
296  kAlignZ,
297  kAlignPhiX,
298  kAlignPhiY,
299  kAlignPhiZ,
300  kResidXSigma,
301  kResidYSigma,
302  kResSlopeXSigma,
303  kResSlopeYSigma,
304  kAlphaX,
305  kAlphaY,
306  kResidXGamma,
307  kResidYGamma,
308  kResSlopeXGamma,
309  kResSlopeYGamma};
310  std::string names[16] = {"AlignX",
311  "AlignY",
312  "AlignZ",
313  "AlignPhiX",
314  "AlignPhiY",
315  "AlignPhiZ",
316  "ResidXSigma",
317  "ResidYSigma",
318  "ResSlopeXSigma",
319  "ResSlopeYSigma",
320  "AlphaX",
321  "AlphaY",
322  "ResidXGamma",
323  "ResidYGamma",
324  "ResSlopeXGamma",
325  "ResSlopeYGamma"};
326  double starts[16] = {0.,
327  0.,
328  0.,
329  0.,
330  0.,
331  0.,
332  resx_std,
333  resy_std,
334  resslopex_std,
335  resslopey_std,
336  0.,
337  0.,
338  0.1 * resx_std,
339  0.1 * resy_std,
340  0.1 * resslopex_std,
341  0.1 * resslopey_std};
342  double steps[16] = {0.1,
343  0.1,
344  0.1,
345  0.001,
346  0.001,
347  0.001,
348  0.001 * resx_std,
349  0.001 * resy_std,
350  0.001 * resslopex_std,
351  0.001 * resslopey_std,
352  0.001,
353  0.001,
354  0.01 * resx_std,
355  0.01 * resy_std,
356  0.01 * resslopex_std,
357  0.01 * resslopey_std};
358  double lows[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., -1., 0., 0., 0., 0.};
359  double highs[16] = {0., 0., 0., 0., 0., 0., 10., 10., 0.1, 0.1, 1., 1., 0., 0., 0., 0.};
360 
361  std::vector<int> num(nums, nums + 6);
362  std::vector<std::string> name(names, names + 6);
363  std::vector<double> start(starts, starts + 6);
364  std::vector<double> step(steps, steps + 6);
365  std::vector<double> low(lows, lows + 6);
366  std::vector<double> high(highs, highs + 6);
367 
368  bool add_alpha = (residualsModel() == kPureGaussian2D);
369  bool add_gamma = (residualsModel() == kROOTVoigt || residualsModel() == kPowerLawTails);
370 
371  int idx[8], ni = 0;
372  if (useRes() == k1111) {
373  for (ni = 0; ni < 4; ni++)
374  idx[ni] = ni + 6;
375  if (add_alpha)
376  for (; ni < 6; ni++)
377  idx[ni] = ni + 6;
378  else if (add_gamma)
379  for (; ni < 8; ni++)
380  idx[ni] = ni + 8;
381  if (!add_alpha)
382  fix(kAlphaX);
383  if (!add_alpha)
384  fix(kAlphaY);
385  } else if (useRes() == k1110) {
386  for (ni = 0; ni < 3; ni++)
387  idx[ni] = ni + 6;
388  if (add_alpha)
389  idx[ni++] = 10;
390  else if (add_gamma)
391  for (; ni < 6; ni++)
392  idx[ni] = ni + 9;
393  fix(kResSlopeYSigma);
394  fix(kAlphaY);
395  if (!add_alpha)
396  fix(kAlphaX);
397  } else if (useRes() == k1100) {
398  for (ni = 0; ni < 2; ni++)
399  idx[ni] = ni + 6;
400  if (add_gamma)
401  for (; ni < 4; ni++)
402  idx[ni] = ni + 10;
403  fix(kResSlopeXSigma);
404  fix(kResSlopeYSigma);
405  fix(kAlphaX);
406  fix(kAlphaY);
407  } else if (useRes() == k1010) {
408  idx[ni++] = 6;
409  idx[ni++] = 8;
410  if (add_alpha)
411  idx[ni++] = 10;
412  if (add_gamma) {
413  idx[ni++] = 12;
414  idx[ni++] = 14;
415  }
416  fix(kResidYSigma);
417  fix(kResSlopeYSigma);
418  fix(kAlphaY);
419  if (!add_alpha)
420  fix(kAlphaX);
421  } else if (useRes() == k0010) {
422  idx[ni++] = 8;
423  if (add_gamma)
424  idx[ni++] = 14;
425  fix(kResidXSigma);
426  fix(kResidYSigma);
427  fix(kResSlopeYSigma);
428  fix(kAlphaX);
429  fix(kAlphaY);
430  }
431  for (int i = 0; i < ni; i++) {
432  num.push_back(nums[idx[i]]);
433  name.push_back(names[idx[i]]);
434  start.push_back(starts[idx[i]]);
435  step.push_back(steps[idx[i]]);
436  low.push_back(lows[idx[i]]);
437  high.push_back(highs[idx[i]]);
438  }
439 
440  return dofit(&MuonResiduals6DOFFitter_FCN, num, name, start, step, low, high);
441 }
442 
443 double MuonResiduals6DOFFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
444  sumofweights();
445 
446  double mean_residualx = 0., mean_residualy = 0., mean_resslopex = 0., mean_resslopey = 0.;
447  double mean_trackx = 0., mean_tracky = 0., mean_trackdxdz = 0., mean_trackdydz = 0.;
448  double sum_w = 0.;
449 
450  for (std::vector<double *>::const_iterator rit = residuals_begin(); rit != residuals_end(); ++rit) {
451  const double redchi2 = (*rit)[kRedChi2];
452  double weight = 1. / redchi2;
453  if (!m_weightAlignment)
454  weight = 1.;
455 
456  if (!m_weightAlignment || TMath::Prob(redchi2 * 12, 12) < 0.99) // no spikes allowed
457  {
458  double factor_w = 1. / (sum_w + weight);
459  mean_residualx = factor_w * (sum_w * mean_residualx + weight * (*rit)[kResidX]);
460  mean_residualy = factor_w * (sum_w * mean_residualy + weight * (*rit)[kResidY]);
461  mean_resslopex = factor_w * (sum_w * mean_resslopex + weight * (*rit)[kResSlopeX]);
462  mean_resslopey = factor_w * (sum_w * mean_resslopey + weight * (*rit)[kResSlopeY]);
463  mean_trackx = factor_w * (sum_w * mean_trackx + weight * (*rit)[kPositionX]);
464  mean_tracky = factor_w * (sum_w * mean_tracky + weight * (*rit)[kPositionY]);
465  mean_trackdxdz = factor_w * (sum_w * mean_trackdxdz + weight * (*rit)[kAngleX]);
466  mean_trackdydz = factor_w * (sum_w * mean_trackdydz + weight * (*rit)[kAngleY]);
467  sum_w += weight;
468  }
469  }
470 
471  std::string name_x, name_y, name_dxdz, name_dydz, name_x_raw, name_y_raw, name_dxdz_raw, name_dydz_raw, name_x_cut,
472  name_y_cut, name_alphax, name_alphay;
473  std::string name_x_trackx, name_y_trackx, name_dxdz_trackx, name_dydz_trackx;
474  std::string name_x_tracky, name_y_tracky, name_dxdz_tracky, name_dydz_tracky;
475  std::string name_x_trackdxdz, name_y_trackdxdz, name_dxdz_trackdxdz, name_dydz_trackdxdz;
476  std::string name_x_trackdydz, name_y_trackdydz, name_dxdz_trackdydz, name_dydz_trackdydz;
477 
478  name_x = name + "_x";
479  name_y = name + "_y";
480  name_dxdz = name + "_dxdz";
481  name_dydz = name + "_dydz";
482  name_x_raw = name + "_x_raw";
483  name_y_raw = name + "_y_raw";
484  name_dxdz_raw = name + "_dxdz_raw";
485  name_dydz_raw = name + "_dydz_raw";
486  name_x_cut = name + "_x_cut";
487  name_y_cut = name + "_y_cut";
488  name_alphax = name + "_alphax";
489  name_alphay = name + "_alphay";
490  name_x_trackx = name + "_x_trackx";
491  name_y_trackx = name + "_y_trackx";
492  name_dxdz_trackx = name + "_dxdz_trackx";
493  name_dydz_trackx = name + "_dydz_trackx";
494  name_x_tracky = name + "_x_tracky";
495  name_y_tracky = name + "_y_tracky";
496  name_dxdz_tracky = name + "_dxdz_tracky";
497  name_dydz_tracky = name + "_dydz_tracky";
498  name_x_trackdxdz = name + "_x_trackdxdz";
499  name_y_trackdxdz = name + "_y_trackdxdz";
500  name_dxdz_trackdxdz = name + "_dxdz_trackdxdz";
501  name_dydz_trackdxdz = name + "_dydz_trackdxdz";
502  name_x_trackdydz = name + "_x_trackdydz";
503  name_y_trackdydz = name + "_y_trackdydz";
504  name_dxdz_trackdydz = name + "_dxdz_trackdydz";
505  name_dydz_trackdydz = name + "_dydz_trackdydz";
506 
507  double width = ali->surface().width();
508  double length = ali->surface().length();
509  int bins = 200;
510  double min_x = -100., max_x = 100.;
511  double min_y = -100., max_y = 100.;
512  double min_dxdz = -75., max_dxdz = 75.;
513  double min_dydz = -150., max_dydz = 150.;
514  double min_trackx = -width / 2., max_trackx = width / 2.;
515  double min_tracky = -length / 2., max_tracky = length / 2.;
516  double min_trackdxdz = -1.5, max_trackdxdz = 1.5;
517  double min_trackdydz = -1.5, max_trackdydz = 1.5;
518 
519  TH1F *hist_x = dir->make<TH1F>(name_x.c_str(), "", bins, min_x, max_x);
520  TH1F *hist_y = dir->make<TH1F>(name_y.c_str(), "", bins, min_y, max_y);
521  TH1F *hist_dxdz = dir->make<TH1F>(name_dxdz.c_str(), "", bins, min_dxdz, max_dxdz);
522  TH1F *hist_dydz = dir->make<TH1F>(name_dydz.c_str(), "", bins, min_dydz, max_dydz);
523  TH1F *hist_x_raw = dir->make<TH1F>(name_x_raw.c_str(), "", bins, min_x, max_x);
524  TH1F *hist_y_raw = dir->make<TH1F>(name_y_raw.c_str(), "", bins, min_y, max_y);
525  TH1F *hist_dxdz_raw = dir->make<TH1F>(name_dxdz_raw.c_str(), "", bins, min_dxdz, max_dxdz);
526  TH1F *hist_dydz_raw = dir->make<TH1F>(name_dydz_raw.c_str(), "", bins, min_dydz, max_dydz);
527  TH1F *hist_x_cut = dir->make<TH1F>(name_x_cut.c_str(), "", bins, min_x, max_x);
528  TH1F *hist_y_cut = dir->make<TH1F>(name_y_cut.c_str(), "", bins, min_y, max_y);
529  TH2F *hist_alphax = dir->make<TH2F>(name_alphax.c_str(), "", 50, 50, 50, 50, -50., 50.);
530  TH2F *hist_alphay = dir->make<TH2F>(name_alphay.c_str(), "", 75, 100, 100, 75, -75., 75.);
531 
532  TProfile *hist_x_trackx = dir->make<TProfile>(name_x_trackx.c_str(), "", 50, min_trackx, max_trackx, min_x, max_x);
533  TProfile *hist_y_trackx = dir->make<TProfile>(name_y_trackx.c_str(), "", 50, min_trackx, max_trackx, min_y, max_y);
534  TProfile *hist_dxdz_trackx =
535  dir->make<TProfile>(name_dxdz_trackx.c_str(), "", 50, min_trackx, max_trackx, min_dxdz, max_dxdz);
536  TProfile *hist_dydz_trackx =
537  dir->make<TProfile>(name_dydz_trackx.c_str(), "", 50, min_trackx, max_trackx, min_dydz, max_dydz);
538  TProfile *hist_x_tracky = dir->make<TProfile>(name_x_tracky.c_str(), "", 50, min_tracky, max_tracky, min_x, max_x);
539  TProfile *hist_y_tracky = dir->make<TProfile>(name_y_tracky.c_str(), "", 50, min_tracky, max_tracky, min_y, max_y);
540  TProfile *hist_dxdz_tracky =
541  dir->make<TProfile>(name_dxdz_tracky.c_str(), "", 50, min_tracky, max_tracky, min_dxdz, max_dxdz);
542  TProfile *hist_dydz_tracky =
543  dir->make<TProfile>(name_dydz_tracky.c_str(), "", 50, min_tracky, max_tracky, min_dydz, max_dydz);
544  TProfile *hist_x_trackdxdz =
545  dir->make<TProfile>(name_x_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_x, max_x);
546  TProfile *hist_y_trackdxdz =
547  dir->make<TProfile>(name_y_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_y, max_y);
548  TProfile *hist_dxdz_trackdxdz =
549  dir->make<TProfile>(name_dxdz_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_dxdz, max_dxdz);
550  TProfile *hist_dydz_trackdxdz =
551  dir->make<TProfile>(name_dydz_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_dydz, max_dydz);
552  TProfile *hist_x_trackdydz =
553  dir->make<TProfile>(name_x_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_x, max_x);
554  TProfile *hist_y_trackdydz =
555  dir->make<TProfile>(name_y_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_y, max_y);
556  TProfile *hist_dxdz_trackdydz =
557  dir->make<TProfile>(name_dxdz_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_dxdz, max_dxdz);
558  TProfile *hist_dydz_trackdydz =
559  dir->make<TProfile>(name_dydz_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_dydz, max_dydz);
560 
561  hist_x_trackx->SetAxisRange(-10., 10., "Y");
562  hist_y_trackx->SetAxisRange(-20., 20., "Y");
563  hist_dxdz_trackx->SetAxisRange(-10., 10., "Y");
564  hist_dydz_trackx->SetAxisRange(-20., 20., "Y");
565  hist_x_tracky->SetAxisRange(-10., 10., "Y");
566  hist_y_tracky->SetAxisRange(-20., 20., "Y");
567  hist_dxdz_tracky->SetAxisRange(-10., 10., "Y");
568  hist_dydz_tracky->SetAxisRange(-20., 20., "Y");
569  hist_x_trackdxdz->SetAxisRange(-10., 10., "Y");
570  hist_y_trackdxdz->SetAxisRange(-20., 20., "Y");
571  hist_dxdz_trackdxdz->SetAxisRange(-10., 10., "Y");
572  hist_dydz_trackdxdz->SetAxisRange(-20., 20., "Y");
573  hist_x_trackdydz->SetAxisRange(-10., 10., "Y");
574  hist_y_trackdydz->SetAxisRange(-20., 20., "Y");
575  hist_dxdz_trackdydz->SetAxisRange(-10., 10., "Y");
576  hist_dydz_trackdydz->SetAxisRange(-20., 20., "Y");
577 
578  name_x += "_fit";
579  name_y += "_fit";
580  name_dxdz += "_fit";
581  name_dydz += "_fit";
582  name_alphax += "_fit";
583  name_alphay += "_fit";
584  name_x_trackx += "_fit";
585  name_y_trackx += "_fit";
586  name_dxdz_trackx += "_fit";
587  name_dydz_trackx += "_fit";
588  name_x_tracky += "_fit";
589  name_y_tracky += "_fit";
590  name_dxdz_tracky += "_fit";
591  name_dydz_tracky += "_fit";
592  name_x_trackdxdz += "_fit";
593  name_y_trackdxdz += "_fit";
594  name_dxdz_trackdxdz += "_fit";
595  name_dydz_trackdxdz += "_fit";
596  name_x_trackdydz += "_fit";
597  name_y_trackdydz += "_fit";
598  name_dxdz_trackdydz += "_fit";
599  name_dydz_trackdydz += "_fit";
600 
601  TF1 *fit_x = nullptr;
602  TF1 *fit_y = nullptr;
603  TF1 *fit_dxdz = nullptr;
604  TF1 *fit_dydz = nullptr;
605  if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) {
606  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_x, max_x, 3);
607  fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins, 10. * value(kAlignX), 10. * value(kResidXSigma));
608  const double er_x[3] = {0., 10. * errorerror(kAlignX), 10. * errorerror(kResidXSigma)};
609  fit_x->SetParErrors(er_x);
610  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_y, max_y, 3);
611  fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins, 10. * value(kAlignY), 10. * value(kResidYSigma));
612  const double er_y[3] = {0., 10. * errorerror(kAlignY), 10. * errorerror(kResidYSigma)};
613  fit_y->SetParErrors(er_y);
614  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dxdz, max_dxdz, 3);
615  fit_dxdz->SetParameters(
616  sum_of_weights * (max_dxdz - min_dxdz) / bins, 1000. * value(kAlignPhiY), 1000. * value(kResSlopeXSigma));
617  const double er_dxdz[3] = {0., 1000. * errorerror(kAlignPhiY), 1000. * errorerror(kResSlopeXSigma)};
618  fit_dxdz->SetParErrors(er_dxdz);
619  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dydz, max_dydz, 3);
620  fit_dydz->SetParameters(
621  sum_of_weights * (max_dydz - min_dydz) / bins, -1000. * value(kAlignPhiX), 1000. * value(kResSlopeYSigma));
622  const double er_dydz[3] = {0., 1000. * errorerror(kAlignPhiX), 1000. * errorerror(kResSlopeYSigma)};
623  fit_dydz->SetParErrors(er_dydz);
624  } else if (residualsModel() == kPowerLawTails) {
625  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_x, max_x, 4);
626  fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins,
627  10. * value(kAlignX),
628  10. * value(kResidXSigma),
629  10. * value(kResidXGamma));
630  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_y, max_y, 4);
631  fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins,
632  10. * value(kAlignY),
633  10. * value(kResidYSigma),
634  10. * value(kResidYGamma));
635  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dxdz, max_dxdz, 4);
636  fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz) / bins,
637  1000. * value(kAlignPhiY),
638  1000. * value(kResSlopeXSigma),
639  1000. * value(kResSlopeXGamma));
640  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dydz, max_dydz, 4);
641  fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz) / bins,
642  -1000. * value(kAlignPhiX),
643  1000. * value(kResSlopeYSigma),
644  1000. * value(kResSlopeYGamma));
645  } else if (residualsModel() == kROOTVoigt) {
646  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_x, max_x, 4);
647  fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins,
648  10. * value(kAlignX),
649  10. * value(kResidXSigma),
650  10. * value(kResidXGamma));
651  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_y, max_y, 4);
652  fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins,
653  10. * value(kAlignY),
654  10. * value(kResidYSigma),
655  10. * value(kResidYGamma));
656  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dxdz, max_dxdz, 4);
657  fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz) / bins,
658  1000. * value(kAlignPhiY),
659  1000. * value(kResSlopeXSigma),
660  1000. * value(kResSlopeXGamma));
661  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dydz, max_dydz, 4);
662  fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz) / bins,
663  -1000. * value(kAlignPhiX),
664  1000. * value(kResSlopeYSigma),
665  1000. * value(kResSlopeYGamma));
666  } else if (residualsModel() == kGaussPowerTails) {
667  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_x, max_x, 3);
668  fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins, 10. * value(kAlignX), 10. * value(kResidXSigma));
669  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_y, max_y, 3);
670  fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins, 10. * value(kAlignY), 10. * value(kResidYSigma));
671  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dxdz, max_dxdz, 3);
672  fit_dxdz->SetParameters(
673  sum_of_weights * (max_dxdz - min_dxdz) / bins, 1000. * value(kAlignPhiY), 1000. * value(kResSlopeXSigma));
674  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dydz, max_dydz, 3);
675  fit_dydz->SetParameters(
676  sum_of_weights * (max_dydz - min_dydz) / bins, -1000. * value(kAlignPhiX), 1000. * value(kResSlopeYSigma));
677  } else {
678  assert(false);
679  }
680 
681  fit_x->SetLineColor(2);
682  fit_x->SetLineWidth(2);
683  fit_x->Write();
684  fit_y->SetLineColor(2);
685  fit_y->SetLineWidth(2);
686  fit_y->Write();
687  fit_dxdz->SetLineColor(2);
688  fit_dxdz->SetLineWidth(2);
689  fit_dxdz->Write();
690  fit_dydz->SetLineColor(2);
691  fit_dydz->SetLineWidth(2);
692  fit_dydz->Write();
693 
694  TF1 *fit_alphax = new TF1(name_alphax.c_str(), "[0] + x*[1]", min_dxdz, max_dxdz);
695  TF1 *fit_alphay = new TF1(name_alphay.c_str(), "[0] + x*[1]", min_dydz, max_dydz);
696  double aX = 10. * value(kAlignX), bX = 10. * value(kAlphaX) / 1000.;
697  double aY = 10. * value(kAlignY), bY = 10. * value(kAlphaY) / 1000.;
698  if (residualsModel() == kPureGaussian2D) {
699  double sx = 10. * value(kResidXSigma), sy = 1000. * value(kResSlopeXSigma), r = value(kAlphaX);
700  aX = mean_residualx;
701  bX = 0.;
702  if (sx != 0.) {
703  bX = 1. / (sy / sx * r);
704  aX = -bX * mean_resslopex;
705  }
706  sx = 10. * value(kResidYSigma);
707  sy = 1000. * value(kResSlopeYSigma);
708  r = value(kAlphaY);
709  aY = mean_residualx;
710  bY = 0.;
711  if (sx != 0.) {
712  bY = 1. / (sy / sx * r);
713  aY = -bY * mean_resslopey;
714  }
715  }
716  fit_alphax->SetParameters(aX, bX);
717  fit_alphay->SetParameters(aY, bY);
718 
719  fit_alphax->SetLineColor(2);
720  fit_alphax->SetLineWidth(2);
721  fit_alphax->Write();
722  fit_alphay->SetLineColor(2);
723  fit_alphay->SetLineWidth(2);
724  fit_alphay->Write();
725 
726  TProfile *fit_x_trackx = dir->make<TProfile>(name_x_trackx.c_str(), "", 100, min_trackx, max_trackx);
727  TProfile *fit_y_trackx = dir->make<TProfile>(name_y_trackx.c_str(), "", 100, min_trackx, max_trackx);
728  TProfile *fit_dxdz_trackx = dir->make<TProfile>(name_dxdz_trackx.c_str(), "", 100, min_trackx, max_trackx);
729  TProfile *fit_dydz_trackx = dir->make<TProfile>(name_dydz_trackx.c_str(), "", 100, min_trackx, max_trackx);
730  TProfile *fit_x_tracky = dir->make<TProfile>(name_x_tracky.c_str(), "", 100, min_tracky, max_tracky);
731  TProfile *fit_y_tracky = dir->make<TProfile>(name_y_tracky.c_str(), "", 100, min_tracky, max_tracky);
732  TProfile *fit_dxdz_tracky = dir->make<TProfile>(name_dxdz_tracky.c_str(), "", 100, min_tracky, max_tracky);
733  TProfile *fit_dydz_tracky = dir->make<TProfile>(name_dydz_tracky.c_str(), "", 100, min_tracky, max_tracky);
734  TProfile *fit_x_trackdxdz = dir->make<TProfile>(name_x_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
735  TProfile *fit_y_trackdxdz = dir->make<TProfile>(name_y_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
736  TProfile *fit_dxdz_trackdxdz =
737  dir->make<TProfile>(name_dxdz_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
738  TProfile *fit_dydz_trackdxdz =
739  dir->make<TProfile>(name_dydz_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
740  TProfile *fit_x_trackdydz = dir->make<TProfile>(name_x_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
741  TProfile *fit_y_trackdydz = dir->make<TProfile>(name_y_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
742  TProfile *fit_dxdz_trackdydz =
743  dir->make<TProfile>(name_dxdz_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
744  TProfile *fit_dydz_trackdydz =
745  dir->make<TProfile>(name_dydz_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
746 
747  fit_x_trackx->SetLineColor(2);
748  fit_x_trackx->SetLineWidth(2);
749  fit_y_trackx->SetLineColor(2);
750  fit_y_trackx->SetLineWidth(2);
751  fit_dxdz_trackx->SetLineColor(2);
752  fit_dxdz_trackx->SetLineWidth(2);
753  fit_dydz_trackx->SetLineColor(2);
754  fit_dydz_trackx->SetLineWidth(2);
755  fit_x_tracky->SetLineColor(2);
756  fit_x_tracky->SetLineWidth(2);
757  fit_y_tracky->SetLineColor(2);
758  fit_y_tracky->SetLineWidth(2);
759  fit_dxdz_tracky->SetLineColor(2);
760  fit_dxdz_tracky->SetLineWidth(2);
761  fit_dydz_tracky->SetLineColor(2);
762  fit_dydz_tracky->SetLineWidth(2);
763  fit_x_trackdxdz->SetLineColor(2);
764  fit_x_trackdxdz->SetLineWidth(2);
765  fit_y_trackdxdz->SetLineColor(2);
766  fit_y_trackdxdz->SetLineWidth(2);
767  fit_dxdz_trackdxdz->SetLineColor(2);
768  fit_dxdz_trackdxdz->SetLineWidth(2);
769  fit_dydz_trackdxdz->SetLineColor(2);
770  fit_dydz_trackdxdz->SetLineWidth(2);
771  fit_x_trackdydz->SetLineColor(2);
772  fit_x_trackdydz->SetLineWidth(2);
773  fit_y_trackdydz->SetLineColor(2);
774  fit_y_trackdydz->SetLineWidth(2);
775  fit_dxdz_trackdydz->SetLineColor(2);
776  fit_dxdz_trackdydz->SetLineWidth(2);
777  fit_dydz_trackdydz->SetLineColor(2);
778  fit_dydz_trackdydz->SetLineWidth(2);
779 
780  name_x_trackx += "line";
781  name_y_trackx += "line";
782  name_dxdz_trackx += "line";
783  name_dydz_trackx += "line";
784  name_x_tracky += "line";
785  name_y_tracky += "line";
786  name_dxdz_tracky += "line";
787  name_dydz_tracky += "line";
788  name_x_trackdxdz += "line";
789  name_y_trackdxdz += "line";
790  name_dxdz_trackdxdz += "line";
791  name_dydz_trackdxdz += "line";
792  name_x_trackdydz += "line";
793  name_y_trackdydz += "line";
794  name_dxdz_trackdydz += "line";
795  name_dydz_trackdydz += "line";
796 
797  TF1 *fitline_x_trackx = new TF1(name_x_trackx.c_str(), residual_x_trackx_TF1, min_trackx, max_trackx, 14);
798  TF1 *fitline_y_trackx = new TF1(name_y_trackx.c_str(), residual_y_trackx_TF1, min_trackx, max_trackx, 14);
799  TF1 *fitline_dxdz_trackx = new TF1(name_dxdz_trackx.c_str(), residual_dxdz_trackx_TF1, min_trackx, max_trackx, 14);
800  TF1 *fitline_dydz_trackx = new TF1(name_dydz_trackx.c_str(), residual_dydz_trackx_TF1, min_trackx, max_trackx, 14);
801  TF1 *fitline_x_tracky = new TF1(name_x_tracky.c_str(), residual_x_tracky_TF1, min_tracky, max_tracky, 14);
802  TF1 *fitline_y_tracky = new TF1(name_y_tracky.c_str(), residual_y_tracky_TF1, min_tracky, max_tracky, 14);
803  TF1 *fitline_dxdz_tracky = new TF1(name_dxdz_tracky.c_str(), residual_dxdz_tracky_TF1, min_tracky, max_tracky, 14);
804  TF1 *fitline_dydz_tracky = new TF1(name_dydz_tracky.c_str(), residual_dydz_tracky_TF1, min_tracky, max_tracky, 14);
805  TF1 *fitline_x_trackdxdz =
806  new TF1(name_x_trackdxdz.c_str(), residual_x_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
807  TF1 *fitline_y_trackdxdz =
808  new TF1(name_y_trackdxdz.c_str(), residual_y_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
809  TF1 *fitline_dxdz_trackdxdz =
810  new TF1(name_dxdz_trackdxdz.c_str(), residual_dxdz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
811  TF1 *fitline_dydz_trackdxdz =
812  new TF1(name_dydz_trackdxdz.c_str(), residual_dydz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
813  TF1 *fitline_x_trackdydz =
814  new TF1(name_x_trackdydz.c_str(), residual_x_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
815  TF1 *fitline_y_trackdydz =
816  new TF1(name_y_trackdydz.c_str(), residual_y_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
817  TF1 *fitline_dxdz_trackdydz =
818  new TF1(name_dxdz_trackdydz.c_str(), residual_dxdz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
819  TF1 *fitline_dydz_trackdydz =
820  new TF1(name_dydz_trackdydz.c_str(), residual_dydz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
821 
822  std::vector<TF1 *> fitlines;
823  fitlines.push_back(fitline_x_trackx);
824  fitlines.push_back(fitline_y_trackx);
825  fitlines.push_back(fitline_dxdz_trackx);
826  fitlines.push_back(fitline_dydz_trackx);
827  fitlines.push_back(fitline_x_tracky);
828  fitlines.push_back(fitline_y_tracky);
829  fitlines.push_back(fitline_dxdz_tracky);
830  fitlines.push_back(fitline_dydz_tracky);
831  fitlines.push_back(fitline_x_trackdxdz);
832  fitlines.push_back(fitline_y_trackdxdz);
833  fitlines.push_back(fitline_dxdz_trackdxdz);
834  fitlines.push_back(fitline_dydz_trackdxdz);
835  fitlines.push_back(fitline_x_trackdydz);
836  fitlines.push_back(fitline_y_trackdydz);
837  fitlines.push_back(fitline_dxdz_trackdydz);
838  fitlines.push_back(fitline_dydz_trackdydz);
839 
840  double fitparameters[14] = {value(kAlignX),
841  value(kAlignY),
842  value(kAlignZ),
843  value(kAlignPhiX),
844  value(kAlignPhiY),
845  value(kAlignPhiZ),
846  mean_trackx,
847  mean_tracky,
848  mean_trackdxdz,
849  mean_trackdydz,
850  value(kAlphaX),
851  mean_resslopex,
852  value(kAlphaY),
853  mean_resslopey};
854  if (residualsModel() == kPureGaussian2D)
855  fitparameters[10] = fitparameters[12] = 0.;
856 
857  for (std::vector<TF1 *>::const_iterator itr = fitlines.begin(); itr != fitlines.end(); itr++) {
858  (*itr)->SetParameters(fitparameters);
859  (*itr)->SetLineColor(2);
860  (*itr)->SetLineWidth(2);
861  (*itr)->Write();
862  }
863 
864  for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
865  const double residX = (*resiter)[kResidX];
866  const double residY = (*resiter)[kResidY];
867  const double resslopeX = (*resiter)[kResSlopeX];
868  const double resslopeY = (*resiter)[kResSlopeY];
869  const double positionX = (*resiter)[kPositionX];
870  const double positionY = (*resiter)[kPositionY];
871  const double angleX = (*resiter)[kAngleX];
872  const double angleY = (*resiter)[kAngleY];
873  const double redchi2 = (*resiter)[kRedChi2];
874  double weight = 1. / redchi2;
875  if (!m_weightAlignment)
876  weight = 1.;
877 
878  if (!m_weightAlignment || TMath::Prob(redchi2 * 12, 12) < 0.99) { // no spikes allowed
879 
880  hist_alphax->Fill(1000. * resslopeX, 10. * residX);
881  hist_alphay->Fill(1000. * resslopeY, 10. * residY);
882 
883  double coefX = value(kAlphaX), coefY = value(kAlphaY);
884  if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D)
885  coefX = coefY = 0.;
886  double geom_residX = residual_x(value(kAlignX),
887  value(kAlignY),
888  value(kAlignZ),
889  value(kAlignPhiX),
890  value(kAlignPhiY),
891  value(kAlignPhiZ),
892  positionX,
893  positionY,
894  angleX,
895  angleY,
896  coefX,
897  resslopeX);
898  hist_x->Fill(10. * (residX - geom_residX + value(kAlignX)), weight);
899  hist_x_trackx->Fill(positionX, 10. * residX, weight);
900  hist_x_tracky->Fill(positionY, 10. * residX, weight);
901  hist_x_trackdxdz->Fill(angleX, 10. * residX, weight);
902  hist_x_trackdydz->Fill(angleY, 10. * residX, weight);
903  fit_x_trackx->Fill(positionX, 10. * geom_residX, weight);
904  fit_x_tracky->Fill(positionY, 10. * geom_residX, weight);
905  fit_x_trackdxdz->Fill(angleX, 10. * geom_residX, weight);
906  fit_x_trackdydz->Fill(angleY, 10. * geom_residX, weight);
907 
908  double geom_residY = residual_y(value(kAlignX),
909  value(kAlignY),
910  value(kAlignZ),
911  value(kAlignPhiX),
912  value(kAlignPhiY),
913  value(kAlignPhiZ),
914  positionX,
915  positionY,
916  angleX,
917  angleY,
918  coefY,
919  resslopeY);
920  hist_y->Fill(10. * (residY - geom_residY + value(kAlignY)), weight);
921  hist_y_trackx->Fill(positionX, 10. * residY, weight);
922  hist_y_tracky->Fill(positionY, 10. * residY, weight);
923  hist_y_trackdxdz->Fill(angleX, 10. * residY, weight);
924  hist_y_trackdydz->Fill(angleY, 10. * residY, weight);
925  fit_y_trackx->Fill(positionX, 10. * geom_residY, weight);
926  fit_y_tracky->Fill(positionY, 10. * geom_residY, weight);
927  fit_y_trackdxdz->Fill(angleX, 10. * geom_residY, weight);
928  fit_y_trackdydz->Fill(angleY, 10. * geom_residY, weight);
929 
930  double geom_resslopeX = residual_dxdz(value(kAlignX),
931  value(kAlignY),
932  value(kAlignZ),
933  value(kAlignPhiX),
934  value(kAlignPhiY),
935  value(kAlignPhiZ),
936  positionX,
937  positionY,
938  angleX,
939  angleY);
940  hist_dxdz->Fill(1000. * (resslopeX - geom_resslopeX + value(kAlignPhiY)), weight);
941  hist_dxdz_trackx->Fill(positionX, 1000. * resslopeX, weight);
942  hist_dxdz_tracky->Fill(positionY, 1000. * resslopeX, weight);
943  hist_dxdz_trackdxdz->Fill(angleX, 1000. * resslopeX, weight);
944  hist_dxdz_trackdydz->Fill(angleY, 1000. * resslopeX, weight);
945  fit_dxdz_trackx->Fill(positionX, 1000. * geom_resslopeX, weight);
946  fit_dxdz_tracky->Fill(positionY, 1000. * geom_resslopeX, weight);
947  fit_dxdz_trackdxdz->Fill(angleX, 1000. * geom_resslopeX, weight);
948  fit_dxdz_trackdydz->Fill(angleY, 1000. * geom_resslopeX, weight);
949 
950  double geom_resslopeY = residual_dydz(value(kAlignX),
951  value(kAlignY),
952  value(kAlignZ),
953  value(kAlignPhiX),
954  value(kAlignPhiY),
955  value(kAlignPhiZ),
956  positionX,
957  positionY,
958  angleX,
959  angleY);
960  hist_dydz->Fill(1000. * (resslopeY - geom_resslopeY - value(kAlignPhiX)), weight);
961  hist_dydz_trackx->Fill(positionX, 1000. * resslopeY, weight);
962  hist_dydz_tracky->Fill(positionY, 1000. * resslopeY, weight);
963  hist_dydz_trackdxdz->Fill(angleX, 1000. * resslopeY, weight);
964  hist_dydz_trackdydz->Fill(angleY, 1000. * resslopeY, weight);
965  fit_dydz_trackx->Fill(positionX, 1000. * geom_resslopeY, weight);
966  fit_dydz_tracky->Fill(positionY, 1000. * geom_resslopeY, weight);
967  fit_dydz_trackdxdz->Fill(angleX, 1000. * geom_resslopeY, weight);
968  fit_dydz_trackdydz->Fill(angleY, 1000. * geom_resslopeY, weight);
969  }
970 
971  hist_x_raw->Fill(10. * residX);
972  hist_y_raw->Fill(10. * residY);
973  hist_dxdz_raw->Fill(1000. * resslopeX);
974  hist_dydz_raw->Fill(1000. * resslopeY);
975  if (fabs(resslopeX) < 0.005)
976  hist_x_cut->Fill(10. * residX);
977  if (fabs(resslopeY) < 0.030)
978  hist_y_cut->Fill(10. * residY);
979  }
980 
981  double chi2 = 0.;
982  double ndof = 0.;
983  for (int i = 1; i <= hist_x->GetNbinsX(); i++) {
984  double xi = hist_x->GetBinCenter(i);
985  double yi = hist_x->GetBinContent(i);
986  double yerri = hist_x->GetBinError(i);
987  double yth = fit_x->Eval(xi);
988  if (yerri > 0.) {
989  chi2 += pow((yth - yi) / yerri, 2);
990  ndof += 1.;
991  }
992  }
993  for (int i = 1; i <= hist_y->GetNbinsX(); i++) {
994  double xi = hist_y->GetBinCenter(i);
995  double yi = hist_y->GetBinContent(i);
996  double yerri = hist_y->GetBinError(i);
997  double yth = fit_y->Eval(xi);
998  if (yerri > 0.) {
999  chi2 += pow((yth - yi) / yerri, 2);
1000  ndof += 1.;
1001  }
1002  }
1003  for (int i = 1; i <= hist_dxdz->GetNbinsX(); i++) {
1004  double xi = hist_dxdz->GetBinCenter(i);
1005  double yi = hist_dxdz->GetBinContent(i);
1006  double yerri = hist_dxdz->GetBinError(i);
1007  double yth = fit_dxdz->Eval(xi);
1008  if (yerri > 0.) {
1009  chi2 += pow((yth - yi) / yerri, 2);
1010  ndof += 1.;
1011  }
1012  }
1013  for (int i = 1; i <= hist_dydz->GetNbinsX(); i++) {
1014  double xi = hist_dydz->GetBinCenter(i);
1015  double yi = hist_dydz->GetBinContent(i);
1016  double yerri = hist_dydz->GetBinError(i);
1017  double yth = fit_dydz->Eval(xi);
1018  if (yerri > 0.) {
1019  chi2 += pow((yth - yi) / yerri, 2);
1020  ndof += 1.;
1021  }
1022  }
1023  ndof -= npar();
1024 
1025  return (ndof > 0. ? chi2 / ndof : -1.);
1026 }
1027 
1028 TTree *MuonResiduals6DOFFitter::readNtuple(
1029  std::string fname, unsigned int wheel, unsigned int station, unsigned int sector, unsigned int preselected) {
1030  TFile *f = new TFile(fname.c_str());
1031  TTree *t = (TTree *)f->Get("mual_ttree");
1032 
1033  // Create new temporary file
1034  TFile *tmpf = new TFile("small_tree_fit.root", "recreate");
1035  assert(tmpf != nullptr);
1036 
1037  // filter the tree (temporarily lives in the new file)
1038  TTree *tt = t->CopyTree(Form(
1039  "is_dt && ring_wheel==%d && station==%d && sector==%d && select==%d", wheel, station, sector, (bool)preselected));
1040 
1041  MuonAlignmentTreeRow r;
1042  tt->SetBranchAddress("res_x", &r.res_x);
1043  tt->SetBranchAddress("res_slope_x", &r.res_slope_x);
1044  tt->SetBranchAddress("res_y", &r.res_y);
1045  tt->SetBranchAddress("res_slope_y", &r.res_slope_y);
1046  tt->SetBranchAddress("pos_x", &r.pos_x);
1047  tt->SetBranchAddress("pos_y", &r.pos_y);
1048  tt->SetBranchAddress("angle_x", &r.angle_x);
1049  tt->SetBranchAddress("angle_y", &r.angle_y);
1050  tt->SetBranchAddress("pz", &r.pz);
1051  tt->SetBranchAddress("pt", &r.pt);
1052  tt->SetBranchAddress("q", &r.q);
1053 
1054  Long64_t nentries = tt->GetEntries();
1055  for (Long64_t i = 0; i < nentries; i++) {
1056  tt->GetEntry(i);
1057  double *rdata = new double[MuonResiduals6DOFFitter::kNData];
1058  rdata[kResidX] = r.res_x;
1059  rdata[kResidY] = r.res_y;
1060  rdata[kResSlopeX] = r.res_slope_x;
1061  rdata[kResSlopeY] = r.res_slope_y;
1062  rdata[kPositionX] = r.pos_x;
1063  rdata[kPositionY] = r.pos_y;
1064  rdata[kAngleX] = r.angle_x;
1065  rdata[kAngleY] = r.angle_y;
1066  rdata[kRedChi2] = 0.1;
1067  rdata[kPz] = r.pz;
1068  rdata[kPt] = r.pt;
1069  rdata[kCharge] = r.q;
1070  fill(rdata);
1071  }
1072  delete f;
1073  //delete tmpf;
1074  return tt;
1075 }
start
Definition: start.py:1
MuonResidualsFitter_logGaussPowerTails
double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma)
Definition: MuonResidualsFitter.cc:128
LaserClient_cfi.high
high
Definition: LaserClient_cfi.py:52
LaserClient_cfi.low
low
Definition: LaserClient_cfi.py:54
mps_fire.i
i
Definition: mps_fire.py:429
MuonResidualsFitter_powerLawTails_TF1
Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:116
Alignable::surface
const AlignableSurface & surface() const
Return the Surface (global position and orientation) of the object.
Definition: Alignable.h:132
MuonResidualsFitter_logPureGaussian2D
double MuonResidualsFitter_logPureGaussian2D(double x, double y, double x0, double y0, double sx, double sy, double r)
Definition: MuonResidualsFitter.cc:43
heavyIonCSV_trainingSettings.idx
idx
Definition: heavyIonCSV_trainingSettings.py:5
MuonResiduals6DOFFitter::fit
bool fit(Alignable *ali) override
Definition: MuonResiduals6DOFFitter.cc:285
nano_mu_local_reco_cff.chi2
chi2
Definition: nano_mu_local_reco_cff.py:193
MuonResidualsFitter::residuals_begin
std::vector< double * >::const_iterator residuals_begin() const
Definition: MuonResidualsFitter.h:190
MuonResiduals6DOFFitter::inform
void inform(TMinuit *tMinuit) override
Definition: MuonResiduals6DOFFitter.cc:262
MuonResidualsFitter::fill
void fill(double *residual)
Definition: MuonResidualsFitter.cc:186
AlignableSurface::width
align::Scalar width() const
Definition: AlignableSurface.h:30
MuonResidualsFitterFitInfo::fitter
MuonResidualsFitter * fitter()
Definition: MuonResidualsFitter.h:253
weight
Definition: weight.py:1
fftjetcommon_cfi.sy
sy
Definition: fftjetcommon_cfi.py:203
MuonResidualsFitter::value
double value(int parNum)
Definition: MuonResidualsFitter.h:145
MuonResidualsFitter::residuals_end
std::vector< double * >::const_iterator residuals_end() const
Definition: MuonResidualsFitter.h:191
protons_cff.xi
xi
Definition: protons_cff.py:21
cms::cuda::assert
assert(be >=bs)
names
const std::string names[nVars_]
Definition: PhotonIDValueMapProducer.cc:124
MuonResidualsFitter::dofit
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)
Definition: MuonResidualsFitter.cc:272
MuonResidualsFitter_logROOTVoigt
double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma)
Definition: MuonResidualsFitter.cc:120
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
DeadROC_duringRun.dir
dir
Definition: DeadROC_duringRun.py:23
trigObjTnPSource_cfi.bins
bins
Definition: trigObjTnPSource_cfi.py:20
MuonResiduals6DOFFitter::plot
double plot(std::string name, TFileDirectory *dir, Alignable *ali) override
Definition: MuonResiduals6DOFFitter.cc:443
command_line.start
start
Definition: command_line.py:165
tt
Definition: TTTypes.h:54
MuonResiduals6DOFFitter::readNtuple
TTree * readNtuple(std::string fname, unsigned int wheel, unsigned int station, unsigned int sector, unsigned int preselected=1)
Definition: MuonResiduals6DOFFitter.cc:1028
MuonResidualsFitter_logPowerLawTails
double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma)
Definition: MuonResidualsFitter.cc:85
f
double f[11][100]
Definition: MuScleFitUtils.cc:78
MuonResidualsFitter_GaussPowerTails_TF1
Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:142
MuonResidualsFitter::fix
void fix(int parNum, bool dofix=true)
Definition: MuonResidualsFitter.cc:171
MuonResidualsFitter::correctBField
virtual void correctBField()=0
MuonResidualsFitter_pureGaussian_TF1
Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:38
EgammaValidation_cff.num
num
Definition: EgammaValidation_cff.py:33
MuonResidualsFitter::errorerror
double errorerror(int parNum)
Definition: MuonResidualsFitter.h:149
ApeEstimator_cff.width
width
Definition: ApeEstimator_cff.py:24
fftjetcommon_cfi.sx
sx
Definition: fftjetcommon_cfi.py:202
AlignableSurface::length
align::Scalar length() const
Definition: AlignableSurface.h:32
MuonResidualsFitter::useRes
int useRes(int pattern=-1)
Definition: MuonResidualsFitter.h:119
alignmentValidation.fname
string fname
main script
Definition: alignmentValidation.py:959
nano_mu_digi_cff.sector
sector
Definition: nano_mu_digi_cff.py:21
step
step
Definition: StallMonitor.cc:83
MuonResidualsFitter_ROOTVoigt_TF1
Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:124
mps_merge.weight
weight
Definition: mps_merge.py:88
MuonResiduals6DOFFitter_FCN
void MuonResiduals6DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
Definition: MuonResiduals6DOFFitter.cc:134
MuonResidualsFitter_logPureGaussian
double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma)
Definition: MuonResidualsFitter.cc:31
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
AlCaHLTBitMon_ParallelJobs.p
def p
Definition: AlCaHLTBitMon_ParallelJobs.py:153
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
Generated for CMSSW Reference Manual by doxygen 1.8.14