CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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 {
16  TMinuit *minuit;
17 
18  double sum_of_weights;
19  double number_of_hits;
20  bool weight_alignment;
21 
22  double residual_x(double delta_x, double delta_y, double delta_z,
23  double delta_phix, double delta_phiy, double delta_phiz,
24  double track_x, double track_y, double track_dxdz, double track_dydz,
25  double alphax, double residual_dxdz)
26  {
27  return delta_x
28  - track_dxdz * delta_z
29  - track_y * track_dxdz * delta_phix
30  + track_x * track_dxdz * delta_phiy
31  - track_y * delta_phiz
32  + residual_dxdz * alphax;
33  }
34 
35  double residual_y(double delta_x, double delta_y, double delta_z,
36  double delta_phix, double delta_phiy, double delta_phiz,
37  double track_x, double track_y, double track_dxdz, double track_dydz,
38  double alphay, double residual_dydz)
39  {
40  return delta_y
41  - track_dydz * delta_z
42  - track_y * track_dydz * delta_phix
43  + track_x * track_dydz * delta_phiy
44  + track_x * delta_phiz
45  + residual_dydz * alphay;
46  }
47 
48  double residual_dxdz(double delta_x, double delta_y, double delta_z,
49  double delta_phix, double delta_phiy, double delta_phiz,
50  double track_x, double track_y, double track_dxdz, double track_dydz)
51  {
52  return -track_dxdz * track_dydz * delta_phix
53  + (1. + track_dxdz * track_dxdz) * delta_phiy
54  - track_dydz * delta_phiz;
55  }
56 
57  double residual_dydz(double delta_x, double delta_y, double delta_z,
58  double delta_phix, double delta_phiy, double delta_phiz,
59  double track_x, double track_y, double track_dxdz, double track_dydz)
60  {
61  return -(1. + track_dydz * track_dydz) * delta_phix
62  + track_dxdz * track_dydz * delta_phiy
63  + track_dxdz * delta_phiz;
64  }
65 
66  Double_t residual_x_trackx_TF1(Double_t *xx, Double_t *p) { 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]); }
67  Double_t residual_x_tracky_TF1(Double_t *xx, Double_t *p) { 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]); }
68  Double_t residual_x_trackdxdz_TF1(Double_t *xx, Double_t *p) { 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]); }
69  Double_t residual_x_trackdydz_TF1(Double_t *xx, Double_t *p) { 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]); }
70 
71  Double_t residual_y_trackx_TF1(Double_t *xx, Double_t *p) { 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]); }
72  Double_t residual_y_tracky_TF1(Double_t *xx, Double_t *p) { 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]); }
73  Double_t residual_y_trackdxdz_TF1(Double_t *xx, Double_t *p) { 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]); }
74  Double_t residual_y_trackdydz_TF1(Double_t *xx, Double_t *p) { 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]); }
75 
76  Double_t residual_dxdz_trackx_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9]); }
77  Double_t residual_dxdz_tracky_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9]); }
78  Double_t residual_dxdz_trackdxdz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9]); }
79  Double_t residual_dxdz_trackdydz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0]); }
80 
81  Double_t residual_dydz_trackx_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9]); }
82  Double_t residual_dydz_tracky_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9]); }
83  Double_t residual_dydz_trackdxdz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9]); }
84  Double_t residual_dydz_trackdydz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0]); }
85 }
86 
87 
88 void MuonResiduals6DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
89 {
90  MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo*)(minuit->GetObjectFit());
91  MuonResidualsFitter *fitter = fitinfo->fitter();
92 
93  fval = 0.;
94  for (std::vector<double*>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end(); ++resiter)
95  {
96  const double residX = (*resiter)[MuonResiduals6DOFFitter::kResidX];
97  const double residY = (*resiter)[MuonResiduals6DOFFitter::kResidY];
98  const double resslopeX = (*resiter)[MuonResiduals6DOFFitter::kResSlopeX];
99  const double resslopeY = (*resiter)[MuonResiduals6DOFFitter::kResSlopeY];
100  const double positionX = (*resiter)[MuonResiduals6DOFFitter::kPositionX];
101  const double positionY = (*resiter)[MuonResiduals6DOFFitter::kPositionY];
102  const double angleX = (*resiter)[MuonResiduals6DOFFitter::kAngleX];
103  const double angleY = (*resiter)[MuonResiduals6DOFFitter::kAngleY];
104  const double redchi2 = (*resiter)[MuonResiduals6DOFFitter::kRedChi2];
105 
106  const double alignx = par[MuonResiduals6DOFFitter::kAlignX];
107  const double aligny = par[MuonResiduals6DOFFitter::kAlignY];
108  const double alignz = par[MuonResiduals6DOFFitter::kAlignZ];
109  const double alignphix = par[MuonResiduals6DOFFitter::kAlignPhiX];
110  const double alignphiy = par[MuonResiduals6DOFFitter::kAlignPhiY];
111  const double alignphiz = par[MuonResiduals6DOFFitter::kAlignPhiZ];
112  const double resXsigma = par[MuonResiduals6DOFFitter::kResidXSigma];
113  const double resYsigma = par[MuonResiduals6DOFFitter::kResidYSigma];
114  const double slopeXsigma = par[MuonResiduals6DOFFitter::kResSlopeXSigma];
115  const double slopeYsigma = par[MuonResiduals6DOFFitter::kResSlopeYSigma];
116  const double alphax = par[MuonResiduals6DOFFitter::kAlphaX];
117  const double alphay = par[MuonResiduals6DOFFitter::kAlphaY];
118  const double resXgamma = par[MuonResiduals6DOFFitter::kResidXGamma];
119  const double resYgamma = par[MuonResiduals6DOFFitter::kResidYGamma];
120  const double slopeXgamma = par[MuonResiduals6DOFFitter::kResSlopeXGamma];
121  const double slopeYgamma = par[MuonResiduals6DOFFitter::kResSlopeYGamma];
122 
123  double coefX = alphax, coefY = alphay;
124  if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian ||
125  fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D) coefX = coefY = 0.;
126  double residXpeak = residual_x(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coefX, resslopeX);
127  double residYpeak = residual_y(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coefY, resslopeY);
128  double slopeXpeak = residual_dxdz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
129  double slopeYpeak = residual_dydz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
130 
131  double weight = (1./redchi2) * number_of_hits / sum_of_weights;
132  if (!weight_alignment) weight = 1.;
133 
134  if (!weight_alignment || TMath::Prob(redchi2*12, 12) < 0.99) // no spikes allowed
135  {
136  if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
137  if (fitter->useRes() == MuonResidualsFitter::k1111) {
138  fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
139  fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
140  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
141  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeY, slopeYpeak, slopeYsigma);
142  }
143  else if (fitter->useRes() == MuonResidualsFitter::k1110) {
144  fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
145  fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
146  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
147  }
148  else if (fitter->useRes() == MuonResidualsFitter::k1100) {
149  fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
150  fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
151  }
152  else if (fitter->useRes() == MuonResidualsFitter::k1010) {
153  fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
154  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
155  }
156  else if (fitter->useRes() == MuonResidualsFitter::k0010) {
157  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
158  }
159  //std::cout<<"FCNx("<<residX<<","<<residXpeak<<","<<resXsigma<<") = "<<MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma)<<std::endl;
160  //std::cout<<"FCNy("<<residY<<","<<residYpeak<<","<<resYsigma<<") = "<<MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma)<<std::endl;
161  //std::cout<<"FCNsx("<<resslopeX<<","<<slopeXpeak<<","<<slopeXsigma<<") = "<<MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma)<<std::endl;
162  //std::cout<<"FCNsy("<<resslopeY<<","<<slopeYpeak<<","<<slopeYsigma<<") = "<<MuonResidualsFitter_logPureGaussian(resslopeY, slopeYpeak, slopeYsigma)<<std::endl;
163  }
164  else if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D) {
165  if (fitter->useRes() == MuonResidualsFitter::k1111) {
166  fval += -weight * MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
167  fval += -weight * MuonResidualsFitter_logPureGaussian2D(residY, resslopeY, residYpeak, slopeYpeak, resYsigma, slopeYsigma, alphay);
168  }
169  else if (fitter->useRes() == MuonResidualsFitter::k1110) {
170  fval += -weight * MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
171  fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
172  }
173  else if (fitter->useRes() == MuonResidualsFitter::k1100) {
174  fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
175  fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
176  }
177  else if (fitter->useRes() == MuonResidualsFitter::k1010) {
178  fval += -weight * MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
179  }
180  else if (fitter->useRes() == MuonResidualsFitter::k0010) {
181  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
182  }
183  //std::cout<<"FCNx("<<residX<<","<<resslopeX<<","<<residXpeak<<","<<slopeXpeak<<","<<resXsigma<<","<<slopeXsigma<<","<<alphax<<") = "<<MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax)<<std::endl;
184  //std::cout<<"FCNy("<<residY<<","<<resslopeY<<","<<residYpeak<<","<<slopeYpeak<<","<<resYsigma<<","<<slopeYsigma<<","<<alphay<<") = "<<MuonResidualsFitter_logPureGaussian2D(residY, resslopeY, residYpeak, slopeYpeak, resYsigma, slopeYsigma, alphay)<<std::endl;
185  }
186  else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
187  fval += -weight * MuonResidualsFitter_logPowerLawTails(residX, residXpeak, resXsigma, resXgamma);
188  fval += -weight * MuonResidualsFitter_logPowerLawTails(residY, residYpeak, resYsigma, resYgamma);
189  fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
190  fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
191  }
192  else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
193  fval += -weight * MuonResidualsFitter_logROOTVoigt(residX, residXpeak, resXsigma, resXgamma);
194  fval += -weight * MuonResidualsFitter_logROOTVoigt(residY, residYpeak, resYsigma, resYgamma);
195  fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
196  fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
197  }
198  else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
199  fval += -weight * MuonResidualsFitter_logGaussPowerTails(residX, residXpeak, resXsigma);
200  fval += -weight * MuonResidualsFitter_logGaussPowerTails(residY, residYpeak, resYsigma);
201  fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeX, slopeXpeak, slopeXsigma);
202  fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeY, slopeYpeak, slopeYsigma);
203  }
204  else { assert(false); }
205  }
206  }
207 /*
208  const double resXsigma = par[MuonResiduals6DOFFitter::kResidXSigma];
209  const double resYsigma = par[MuonResiduals6DOFFitter::kResidYSigma];
210  const double slopeXsigma = par[MuonResiduals6DOFFitter::kResSlopeXSigma];
211  const double slopeYsigma = par[MuonResiduals6DOFFitter::kResSlopeYSigma];
212  const double alphax = par[MuonResiduals6DOFFitter::kAlphaX];
213  const double alphay = par[MuonResiduals6DOFFitter::kAlphaY];
214  std::cout<<"fval="<<fval<<","<<resXsigma<<","<<slopeXsigma<<","<<alphax<<","<<resYsigma<<","<<slopeYsigma<<","<<alphay<<std::endl;
215 */
216 }
217 
218 
219 void MuonResiduals6DOFFitter::inform(TMinuit *tMinuit)
220 {
221  minuit = tMinuit;
222 }
223 
224 
225 void MuonResiduals6DOFFitter::correctBField()
226 {
227  MuonResidualsFitter::correctBField(kPt, kCharge);
228 }
229 
230 
231 double MuonResiduals6DOFFitter::sumofweights()
232 {
233  sum_of_weights = 0.;
234  number_of_hits = 0.;
235  weight_alignment = m_weightAlignment;
236  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
237  if (m_weightAlignment) {
238  double redchi2 = (*resiter)[kRedChi2];
239  if (TMath::Prob(redchi2*12, 12) < 0.99) { // no spikes allowed
240  sum_of_weights += 1./redchi2;
241  number_of_hits += 1.;
242  }
243  }
244  else {
245  sum_of_weights += 1.;
246  number_of_hits += 1.;
247  }
248  }
249  return sum_of_weights;
250 }
251 
252 
253 bool MuonResiduals6DOFFitter::fit(Alignable *ali)
254 {
255  initialize_table(); // if not already initialized
256  sumofweights();
257 
258  double resx_std = 0.5;
259  double resy_std = 3.0;
260  double resslopex_std = 0.002;
261  double resslopey_std = 0.005;
262 
263  int nums[16] = {kAlignX, kAlignY, kAlignZ, kAlignPhiX, kAlignPhiY, kAlignPhiZ, kResidXSigma, kResidYSigma, kResSlopeXSigma, kResSlopeYSigma,
264  kAlphaX, kAlphaY, kResidXGamma, kResidYGamma, kResSlopeXGamma, kResSlopeYGamma};
265  std::string names[16] = {"AlignX","AlignY","AlignZ","AlignPhiX","AlignPhiY","AlignPhiZ", "ResidXSigma","ResidYSigma","ResSlopeXSigma","ResSlopeYSigma",
266  "AlphaX","AlphaY", "ResidXGamma","ResidYGamma","ResSlopeXGamma","ResSlopeYGamma"};
267  double starts[16] = {0., 0., 0., 0., 0., 0., resx_std, resy_std, resslopex_std, resslopey_std,
268  0., 0., 0.1*resx_std, 0.1*resy_std, 0.1*resslopex_std, 0.1*resslopey_std};
269  double steps[16] = {0.1, 0.1, 0.1, 0.001, 0.001, 0.001, 0.001*resx_std, 0.001*resy_std, 0.001*resslopex_std, 0.001*resslopey_std,
270  0.001, 0.001, 0.01*resx_std, 0.01*resy_std, 0.01*resslopex_std, 0.01*resslopey_std};
271  double lows[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
272  -1., -1., 0., 0., 0., 0.};
273  double highs[16] = {0., 0., 0., 0., 0., 0., 10., 10., 0.1, 0.1,
274  1.,1., 0., 0., 0., 0.};
275 
276  std::vector<int> num(nums, nums+6);
277  std::vector<std::string> name(names, names+6);
278  std::vector<double> start(starts, starts+6);
279  std::vector<double> step(steps, steps+6);
280  std::vector<double> low(lows, lows+6);
281  std::vector<double> high(highs, highs+6);
282 
283  bool add_alpha = ( residualsModel() == kPureGaussian2D );
284  bool add_gamma = ( residualsModel() == kROOTVoigt || residualsModel() == kPowerLawTails );
285 
286  int idx[8], ni = 0;
287  if (useRes() == k1111) {
288  for(ni=0; ni<4; ni++) idx[ni] = ni+6;
289  if (add_alpha) for(; ni<6; ni++) idx[ni] = ni+6;
290  else if (add_gamma) for(; ni<8; ni++) idx[ni] = ni+8;
291  if (!add_alpha) fix(kAlphaX);
292  if (!add_alpha) fix(kAlphaY);
293  }
294  else if (useRes() == k1110) {
295  for(ni=0; ni<3; ni++) idx[ni] = ni+6;
296  if (add_alpha) idx[ni++] = 10;
297  else if (add_gamma) for(; ni<6; ni++) idx[ni] = ni+9;
298  fix(kResSlopeYSigma);
299  fix(kAlphaY);
300  if (!add_alpha) fix(kAlphaX);
301  }
302  else if (useRes() == k1100) {
303  for(ni=0; ni<2; ni++) idx[ni] = ni+6;
304  if (add_gamma) for(; ni<4; ni++) idx[ni] = ni+10;
305  fix(kResSlopeXSigma);
306  fix(kResSlopeYSigma);
307  fix(kAlphaX);
308  fix(kAlphaY);
309  }
310  else if (useRes() == k1010) {
311  idx[ni++] = 6; idx[ni++] = 8;
312  if (add_alpha) idx[ni++] = 10;
313  if (add_gamma) { idx[ni++] = 12; idx[ni++] = 14; }
314  fix(kResidYSigma);
315  fix(kResSlopeYSigma);
316  fix(kAlphaY);
317  if (!add_alpha) fix(kAlphaX);
318  }
319  else if (useRes() == k0010) {
320  idx[ni++] = 8;
321  if (add_gamma) idx[ni++] = 14;
322  fix(kResidXSigma);
323  fix(kResidYSigma);
324  fix(kResSlopeYSigma);
325  fix(kAlphaX);
326  fix(kAlphaY);
327  }
328  for (int i=0; i<ni; i++){
329  num.push_back(nums[idx[i]]);
330  name.push_back(names[idx[i]]);
331  start.push_back(starts[idx[i]]);
332  step.push_back(steps[idx[i]]);
333  low.push_back(lows[idx[i]]);
334  high.push_back(highs[idx[i]]);
335  }
336 
337  return dofit(&MuonResiduals6DOFFitter_FCN, num, name, start, step, low, high);
338 }
339 
340 
341 double MuonResiduals6DOFFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali)
342 {
343  sumofweights();
344 
345  double mean_residualx = 0., mean_residualy = 0., mean_resslopex = 0., mean_resslopey = 0.;
346  double mean_trackx = 0., mean_tracky = 0., mean_trackdxdz = 0., mean_trackdydz = 0.;
347  double sum_w = 0.;
348 
349  for (std::vector<double*>::const_iterator rit = residuals_begin(); rit != residuals_end(); ++rit)
350  {
351  const double redchi2 = (*rit)[kRedChi2];
352  double weight = 1./redchi2;
353  if (!m_weightAlignment) weight = 1.;
354 
355  if (!m_weightAlignment || TMath::Prob(redchi2*12, 12) < 0.99) // no spikes allowed
356  {
357  double factor_w = 1./(sum_w + weight);
358  mean_residualx = factor_w * (sum_w * mean_residualx + weight * (*rit)[kResidX]);
359  mean_residualy = factor_w * (sum_w * mean_residualy + weight * (*rit)[kResidY]);
360  mean_resslopex = factor_w * (sum_w * mean_resslopex + weight * (*rit)[kResSlopeX]);
361  mean_resslopey = factor_w * (sum_w * mean_resslopey + weight * (*rit)[kResSlopeY]);
362  mean_trackx = factor_w * (sum_w * mean_trackx + weight * (*rit)[kPositionX]);
363  mean_tracky = factor_w * (sum_w * mean_tracky + weight * (*rit)[kPositionY]);
364  mean_trackdxdz = factor_w * (sum_w * mean_trackdxdz + weight * (*rit)[kAngleX]);
365  mean_trackdydz = factor_w * (sum_w * mean_trackdydz + weight * (*rit)[kAngleY]);
366  sum_w += weight;
367  }
368  }
369 
370  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, name_y_cut, name_alphax, name_alphay;
371  std::string name_x_trackx, name_y_trackx, name_dxdz_trackx, name_dydz_trackx;
372  std::string name_x_tracky, name_y_tracky, name_dxdz_tracky, name_dydz_tracky;
373  std::string name_x_trackdxdz, name_y_trackdxdz, name_dxdz_trackdxdz, name_dydz_trackdxdz;
374  std::string name_x_trackdydz, name_y_trackdydz, name_dxdz_trackdydz, name_dydz_trackdydz;
375 
376  name_x = name + "_x";
377  name_y = name + "_y";
378  name_dxdz = name + "_dxdz";
379  name_dydz = name + "_dydz";
380  name_x_raw = name + "_x_raw";
381  name_y_raw = name + "_y_raw";
382  name_dxdz_raw = name + "_dxdz_raw";
383  name_dydz_raw = name + "_dydz_raw";
384  name_x_cut = name + "_x_cut";
385  name_y_cut = name + "_y_cut";
386  name_alphax = name + "_alphax";
387  name_alphay = name + "_alphay";
388  name_x_trackx = name + "_x_trackx";
389  name_y_trackx = name + "_y_trackx";
390  name_dxdz_trackx = name + "_dxdz_trackx";
391  name_dydz_trackx = name + "_dydz_trackx";
392  name_x_tracky = name + "_x_tracky";
393  name_y_tracky = name + "_y_tracky";
394  name_dxdz_tracky = name + "_dxdz_tracky";
395  name_dydz_tracky = name + "_dydz_tracky";
396  name_x_trackdxdz = name + "_x_trackdxdz";
397  name_y_trackdxdz = name + "_y_trackdxdz";
398  name_dxdz_trackdxdz = name + "_dxdz_trackdxdz";
399  name_dydz_trackdxdz = name + "_dydz_trackdxdz";
400  name_x_trackdydz = name + "_x_trackdydz";
401  name_y_trackdydz = name + "_y_trackdydz";
402  name_dxdz_trackdydz = name + "_dxdz_trackdydz";
403  name_dydz_trackdydz = name + "_dydz_trackdydz";
404 
405  double width = ali->surface().width();
406  double length = ali->surface().length();
407  int bins = 200;
408  double min_x = -100., max_x = 100.;
409  double min_y = -100., max_y = 100.;
410  double min_dxdz = -75., max_dxdz = 75.;
411  double min_dydz = -150., max_dydz = 150.;
412  double min_trackx = -width/2., max_trackx = width/2.;
413  double min_tracky = -length/2., max_tracky = length/2.;
414  double min_trackdxdz = -1.5, max_trackdxdz = 1.5;
415  double min_trackdydz = -1.5, max_trackdydz = 1.5;
416 
417  TH1F *hist_x = dir->make<TH1F>(name_x.c_str(), "", bins, min_x, max_x);
418  TH1F *hist_y = dir->make<TH1F>(name_y.c_str(), "", bins, min_y, max_y);
419  TH1F *hist_dxdz = dir->make<TH1F>(name_dxdz.c_str(), "", bins, min_dxdz, max_dxdz);
420  TH1F *hist_dydz = dir->make<TH1F>(name_dydz.c_str(), "", bins, min_dydz, max_dydz);
421  TH1F *hist_x_raw = dir->make<TH1F>(name_x_raw.c_str(), "", bins, min_x, max_x);
422  TH1F *hist_y_raw = dir->make<TH1F>(name_y_raw.c_str(), "", bins, min_y, max_y);
423  TH1F *hist_dxdz_raw = dir->make<TH1F>(name_dxdz_raw.c_str(), "", bins, min_dxdz, max_dxdz);
424  TH1F *hist_dydz_raw = dir->make<TH1F>(name_dydz_raw.c_str(), "", bins, min_dydz, max_dydz);
425  TH1F *hist_x_cut = dir->make<TH1F>(name_x_cut.c_str(), "", bins, min_x, max_x);
426  TH1F *hist_y_cut = dir->make<TH1F>(name_y_cut.c_str(), "", bins, min_y, max_y);
427  TH2F *hist_alphax = dir->make<TH2F>(name_alphax.c_str(), "", 50, 50, 50, 50, -50., 50.);
428  TH2F *hist_alphay = dir->make<TH2F>(name_alphay.c_str(), "", 75, 100, 100, 75, -75., 75.);
429 
430  TProfile *hist_x_trackx = dir->make<TProfile>(name_x_trackx.c_str(), "", 50, min_trackx, max_trackx, min_x, max_x);
431  TProfile *hist_y_trackx = dir->make<TProfile>(name_y_trackx.c_str(), "", 50, min_trackx, max_trackx, min_y, max_y);
432  TProfile *hist_dxdz_trackx = dir->make<TProfile>(name_dxdz_trackx.c_str(), "", 50, min_trackx, max_trackx, min_dxdz, max_dxdz);
433  TProfile *hist_dydz_trackx = dir->make<TProfile>(name_dydz_trackx.c_str(), "", 50, min_trackx, max_trackx, min_dydz, max_dydz);
434  TProfile *hist_x_tracky = dir->make<TProfile>(name_x_tracky.c_str(), "", 50, min_tracky, max_tracky, min_x, max_x);
435  TProfile *hist_y_tracky = dir->make<TProfile>(name_y_tracky.c_str(), "", 50, min_tracky, max_tracky, min_y, max_y);
436  TProfile *hist_dxdz_tracky = dir->make<TProfile>(name_dxdz_tracky.c_str(), "", 50, min_tracky, max_tracky, min_dxdz, max_dxdz);
437  TProfile *hist_dydz_tracky = dir->make<TProfile>(name_dydz_tracky.c_str(), "", 50, min_tracky, max_tracky, min_dydz, max_dydz);
438  TProfile *hist_x_trackdxdz = dir->make<TProfile>(name_x_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_x, max_x);
439  TProfile *hist_y_trackdxdz = dir->make<TProfile>(name_y_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_y, max_y);
440  TProfile *hist_dxdz_trackdxdz = dir->make<TProfile>(name_dxdz_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_dxdz, max_dxdz);
441  TProfile *hist_dydz_trackdxdz = dir->make<TProfile>(name_dydz_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_dydz, max_dydz);
442  TProfile *hist_x_trackdydz = dir->make<TProfile>(name_x_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_x, max_x);
443  TProfile *hist_y_trackdydz = dir->make<TProfile>(name_y_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_y, max_y);
444  TProfile *hist_dxdz_trackdydz = dir->make<TProfile>(name_dxdz_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_dxdz, max_dxdz);
445  TProfile *hist_dydz_trackdydz = dir->make<TProfile>(name_dydz_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_dydz, max_dydz);
446 
447  hist_x_trackx->SetAxisRange(-10., 10., "Y");
448  hist_y_trackx->SetAxisRange(-20., 20., "Y");
449  hist_dxdz_trackx->SetAxisRange(-10., 10., "Y");
450  hist_dydz_trackx->SetAxisRange(-20., 20., "Y");
451  hist_x_tracky->SetAxisRange(-10., 10., "Y");
452  hist_y_tracky->SetAxisRange(-20., 20., "Y");
453  hist_dxdz_tracky->SetAxisRange(-10., 10., "Y");
454  hist_dydz_tracky->SetAxisRange(-20., 20., "Y");
455  hist_x_trackdxdz->SetAxisRange(-10., 10., "Y");
456  hist_y_trackdxdz->SetAxisRange(-20., 20., "Y");
457  hist_dxdz_trackdxdz->SetAxisRange(-10., 10., "Y");
458  hist_dydz_trackdxdz->SetAxisRange(-20., 20., "Y");
459  hist_x_trackdydz->SetAxisRange(-10., 10., "Y");
460  hist_y_trackdydz->SetAxisRange(-20., 20., "Y");
461  hist_dxdz_trackdydz->SetAxisRange(-10., 10., "Y");
462  hist_dydz_trackdydz->SetAxisRange(-20., 20., "Y");
463 
464  name_x += "_fit";
465  name_y += "_fit";
466  name_dxdz += "_fit";
467  name_dydz += "_fit";
468  name_alphax += "_fit";
469  name_alphay += "_fit";
470  name_x_trackx += "_fit";
471  name_y_trackx += "_fit";
472  name_dxdz_trackx += "_fit";
473  name_dydz_trackx += "_fit";
474  name_x_tracky += "_fit";
475  name_y_tracky += "_fit";
476  name_dxdz_tracky += "_fit";
477  name_dydz_tracky += "_fit";
478  name_x_trackdxdz += "_fit";
479  name_y_trackdxdz += "_fit";
480  name_dxdz_trackdxdz += "_fit";
481  name_dydz_trackdxdz += "_fit";
482  name_x_trackdydz += "_fit";
483  name_y_trackdydz += "_fit";
484  name_dxdz_trackdydz += "_fit";
485  name_dydz_trackdydz += "_fit";
486 
487  TF1 *fit_x = NULL;
488  TF1 *fit_y = NULL;
489  TF1 *fit_dxdz = NULL;
490  TF1 *fit_dydz = NULL;
491  if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) {
492  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_x, max_x, 3);
493  fit_x->SetParameters(sum_of_weights * (max_x - min_x)/bins, 10.*value(kAlignX), 10.*value(kResidXSigma));
494  const double er_x[3] = {0., 10.*errorerror(kAlignX), 10.*errorerror(kResidXSigma)};
495  fit_x->SetParErrors(er_x);
496  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_y, max_y, 3);
497  fit_y->SetParameters(sum_of_weights * (max_y - min_y)/bins, 10.*value(kAlignY), 10.*value(kResidYSigma));
498  const double er_y[3] = {0., 10.*errorerror(kAlignY), 10.*errorerror(kResidYSigma)};
499  fit_y->SetParErrors(er_y);
500  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dxdz, max_dxdz, 3);
501  fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz)/bins, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma));
502  const double er_dxdz[3] = {0., 1000.*errorerror(kAlignPhiY), 1000.*errorerror(kResSlopeXSigma)};
503  fit_dxdz->SetParErrors(er_dxdz);
504  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dydz, max_dydz, 3);
505  fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz)/bins, -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma));
506  const double er_dydz[3] = {0., 1000.*errorerror(kAlignPhiX), 1000.*errorerror(kResSlopeYSigma)};
507  fit_dydz->SetParErrors(er_dydz);
508  }
509  else if (residualsModel() == kPowerLawTails) {
510  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_x, max_x, 4);
511  fit_x->SetParameters(sum_of_weights * (max_x - min_x)/bins, 10.*value(kAlignX), 10.*value(kResidXSigma), 10.*value(kResidXGamma));
512  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_y, max_y, 4);
513  fit_y->SetParameters(sum_of_weights * (max_y - min_y)/bins, 10.*value(kAlignY), 10.*value(kResidYSigma), 10.*value(kResidYGamma));
514  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dxdz, max_dxdz, 4);
515  fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz)/bins, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma), 1000.*value(kResSlopeXGamma));
516  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dydz, max_dydz, 4);
517  fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz)/bins, -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma), 1000.*value(kResSlopeYGamma));
518  }
519  else if (residualsModel() == kROOTVoigt) {
520  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_x, max_x, 4);
521  fit_x->SetParameters(sum_of_weights * (max_x - min_x)/bins, 10.*value(kAlignX), 10.*value(kResidXSigma), 10.*value(kResidXGamma));
522  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_y, max_y, 4);
523  fit_y->SetParameters(sum_of_weights * (max_y - min_y)/bins, 10.*value(kAlignY), 10.*value(kResidYSigma), 10.*value(kResidYGamma));
524  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dxdz, max_dxdz, 4);
525  fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz)/bins, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma), 1000.*value(kResSlopeXGamma));
526  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dydz, max_dydz, 4);
527  fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz)/bins, -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma), 1000.*value(kResSlopeYGamma));
528  }
529  else if (residualsModel() == kGaussPowerTails) {
530  fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_x, max_x, 3);
531  fit_x->SetParameters(sum_of_weights * (max_x - min_x)/bins, 10.*value(kAlignX), 10.*value(kResidXSigma));
532  fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_y, max_y, 3);
533  fit_y->SetParameters(sum_of_weights * (max_y - min_y)/bins, 10.*value(kAlignY), 10.*value(kResidYSigma));
534  fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dxdz, max_dxdz, 3);
535  fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz)/bins, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma));
536  fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dydz, max_dydz, 3);
537  fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz)/bins, -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma));
538  }
539  else { assert(false); }
540 
541  fit_x->SetLineColor(2); fit_x->SetLineWidth(2); fit_x->Write();
542  fit_y->SetLineColor(2); fit_y->SetLineWidth(2); fit_y->Write();
543  fit_dxdz->SetLineColor(2); fit_dxdz->SetLineWidth(2); fit_dxdz->Write();
544  fit_dydz->SetLineColor(2); fit_dydz->SetLineWidth(2); fit_dydz->Write();
545 
546  TF1 *fit_alphax = new TF1(name_alphax.c_str(), "[0] + x*[1]", min_dxdz, max_dxdz);
547  TF1 *fit_alphay = new TF1(name_alphay.c_str(), "[0] + x*[1]", min_dydz, max_dydz);
548  double aX = 10.*value(kAlignX), bX = 10.*value(kAlphaX)/1000.;
549  double aY = 10.*value(kAlignY), bY = 10.*value(kAlphaY)/1000.;
550  if (residualsModel() == kPureGaussian2D)
551  {
552  double sx = 10.*value(kResidXSigma), sy = 1000.*value(kResSlopeXSigma), r = value(kAlphaX);
553  aX = mean_residualx;
554  bX = 0.;
555  if ( sx != 0. )
556  {
557  bX = 1./(sy/sx*r);
558  aX = - bX * mean_resslopex;
559  }
560  sx = 10.*value(kResidYSigma); sy = 1000.*value(kResSlopeYSigma); r = value(kAlphaY);
561  aY = mean_residualx;
562  bY = 0.;
563  if ( sx != 0. )
564  {
565  bY = 1./(sy/sx*r);
566  aY = - bY * mean_resslopey;
567  }
568  }
569  fit_alphax->SetParameters(aX, bX);
570  fit_alphay->SetParameters(aY, bY);
571 
572  fit_alphax->SetLineColor(2); fit_alphax->SetLineWidth(2); fit_alphax->Write();
573  fit_alphay->SetLineColor(2); fit_alphay->SetLineWidth(2); fit_alphay->Write();
574 
575  TProfile *fit_x_trackx = dir->make<TProfile>(name_x_trackx.c_str(), "", 100, min_trackx, max_trackx);
576  TProfile *fit_y_trackx = dir->make<TProfile>(name_y_trackx.c_str(), "", 100, min_trackx, max_trackx);
577  TProfile *fit_dxdz_trackx = dir->make<TProfile>(name_dxdz_trackx.c_str(), "", 100, min_trackx, max_trackx);
578  TProfile *fit_dydz_trackx = dir->make<TProfile>(name_dydz_trackx.c_str(), "", 100, min_trackx, max_trackx);
579  TProfile *fit_x_tracky = dir->make<TProfile>(name_x_tracky.c_str(), "", 100, min_tracky, max_tracky);
580  TProfile *fit_y_tracky = dir->make<TProfile>(name_y_tracky.c_str(), "", 100, min_tracky, max_tracky);
581  TProfile *fit_dxdz_tracky = dir->make<TProfile>(name_dxdz_tracky.c_str(), "", 100, min_tracky, max_tracky);
582  TProfile *fit_dydz_tracky = dir->make<TProfile>(name_dydz_tracky.c_str(), "", 100, min_tracky, max_tracky);
583  TProfile *fit_x_trackdxdz = dir->make<TProfile>(name_x_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
584  TProfile *fit_y_trackdxdz = dir->make<TProfile>(name_y_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
585  TProfile *fit_dxdz_trackdxdz = dir->make<TProfile>(name_dxdz_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
586  TProfile *fit_dydz_trackdxdz = dir->make<TProfile>(name_dydz_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
587  TProfile *fit_x_trackdydz = dir->make<TProfile>(name_x_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
588  TProfile *fit_y_trackdydz = dir->make<TProfile>(name_y_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
589  TProfile *fit_dxdz_trackdydz = dir->make<TProfile>(name_dxdz_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
590  TProfile *fit_dydz_trackdydz = dir->make<TProfile>(name_dydz_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
591 
592  fit_x_trackx->SetLineColor(2); fit_x_trackx->SetLineWidth(2);
593  fit_y_trackx->SetLineColor(2); fit_y_trackx->SetLineWidth(2);
594  fit_dxdz_trackx->SetLineColor(2); fit_dxdz_trackx->SetLineWidth(2);
595  fit_dydz_trackx->SetLineColor(2); fit_dydz_trackx->SetLineWidth(2);
596  fit_x_tracky->SetLineColor(2); fit_x_tracky->SetLineWidth(2);
597  fit_y_tracky->SetLineColor(2); fit_y_tracky->SetLineWidth(2);
598  fit_dxdz_tracky->SetLineColor(2); fit_dxdz_tracky->SetLineWidth(2);
599  fit_dydz_tracky->SetLineColor(2); fit_dydz_tracky->SetLineWidth(2);
600  fit_x_trackdxdz->SetLineColor(2); fit_x_trackdxdz->SetLineWidth(2);
601  fit_y_trackdxdz->SetLineColor(2); fit_y_trackdxdz->SetLineWidth(2);
602  fit_dxdz_trackdxdz->SetLineColor(2); fit_dxdz_trackdxdz->SetLineWidth(2);
603  fit_dydz_trackdxdz->SetLineColor(2); fit_dydz_trackdxdz->SetLineWidth(2);
604  fit_x_trackdydz->SetLineColor(2); fit_x_trackdydz->SetLineWidth(2);
605  fit_y_trackdydz->SetLineColor(2); fit_y_trackdydz->SetLineWidth(2);
606  fit_dxdz_trackdydz->SetLineColor(2); fit_dxdz_trackdydz->SetLineWidth(2);
607  fit_dydz_trackdydz->SetLineColor(2); fit_dydz_trackdydz->SetLineWidth(2);
608 
609  name_x_trackx += "line";
610  name_y_trackx += "line";
611  name_dxdz_trackx += "line";
612  name_dydz_trackx += "line";
613  name_x_tracky += "line";
614  name_y_tracky += "line";
615  name_dxdz_tracky += "line";
616  name_dydz_tracky += "line";
617  name_x_trackdxdz += "line";
618  name_y_trackdxdz += "line";
619  name_dxdz_trackdxdz += "line";
620  name_dydz_trackdxdz += "line";
621  name_x_trackdydz += "line";
622  name_y_trackdydz += "line";
623  name_dxdz_trackdydz += "line";
624  name_dydz_trackdydz += "line";
625 
626  TF1 *fitline_x_trackx = new TF1(name_x_trackx.c_str(), residual_x_trackx_TF1, min_trackx, max_trackx, 14);
627  TF1 *fitline_y_trackx = new TF1(name_y_trackx.c_str(), residual_y_trackx_TF1, min_trackx, max_trackx, 14);
628  TF1 *fitline_dxdz_trackx = new TF1(name_dxdz_trackx.c_str(), residual_dxdz_trackx_TF1, min_trackx, max_trackx, 14);
629  TF1 *fitline_dydz_trackx = new TF1(name_dydz_trackx.c_str(), residual_dydz_trackx_TF1, min_trackx, max_trackx, 14);
630  TF1 *fitline_x_tracky = new TF1(name_x_tracky.c_str(), residual_x_tracky_TF1, min_tracky, max_tracky, 14);
631  TF1 *fitline_y_tracky = new TF1(name_y_tracky.c_str(), residual_y_tracky_TF1, min_tracky, max_tracky, 14);
632  TF1 *fitline_dxdz_tracky = new TF1(name_dxdz_tracky.c_str(), residual_dxdz_tracky_TF1, min_tracky, max_tracky, 14);
633  TF1 *fitline_dydz_tracky = new TF1(name_dydz_tracky.c_str(), residual_dydz_tracky_TF1, min_tracky, max_tracky, 14);
634  TF1 *fitline_x_trackdxdz = new TF1(name_x_trackdxdz.c_str(), residual_x_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
635  TF1 *fitline_y_trackdxdz = new TF1(name_y_trackdxdz.c_str(), residual_y_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
636  TF1 *fitline_dxdz_trackdxdz = new TF1(name_dxdz_trackdxdz.c_str(), residual_dxdz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
637  TF1 *fitline_dydz_trackdxdz = new TF1(name_dydz_trackdxdz.c_str(), residual_dydz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
638  TF1 *fitline_x_trackdydz = new TF1(name_x_trackdydz.c_str(), residual_x_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
639  TF1 *fitline_y_trackdydz = new TF1(name_y_trackdydz.c_str(), residual_y_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
640  TF1 *fitline_dxdz_trackdydz = new TF1(name_dxdz_trackdydz.c_str(), residual_dxdz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
641  TF1 *fitline_dydz_trackdydz = new TF1(name_dydz_trackdydz.c_str(), residual_dydz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
642 
643  std::vector<TF1*> fitlines;
644  fitlines.push_back(fitline_x_trackx);
645  fitlines.push_back(fitline_y_trackx);
646  fitlines.push_back(fitline_dxdz_trackx);
647  fitlines.push_back(fitline_dydz_trackx);
648  fitlines.push_back(fitline_x_tracky);
649  fitlines.push_back(fitline_y_tracky);
650  fitlines.push_back(fitline_dxdz_tracky);
651  fitlines.push_back(fitline_dydz_tracky);
652  fitlines.push_back(fitline_x_trackdxdz);
653  fitlines.push_back(fitline_y_trackdxdz);
654  fitlines.push_back(fitline_dxdz_trackdxdz);
655  fitlines.push_back(fitline_dydz_trackdxdz);
656  fitlines.push_back(fitline_x_trackdydz);
657  fitlines.push_back(fitline_y_trackdydz);
658  fitlines.push_back(fitline_dxdz_trackdydz);
659  fitlines.push_back(fitline_dydz_trackdydz);
660 
661  double fitparameters[14] = {value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ),
662  mean_trackx, mean_tracky, mean_trackdxdz, mean_trackdydz,
663  value(kAlphaX), mean_resslopex, value(kAlphaY), mean_resslopey};
664  if (residualsModel() == kPureGaussian2D) fitparameters[10] = fitparameters[12] = 0.;
665 
666  for(std::vector<TF1*>::const_iterator itr = fitlines.begin(); itr != fitlines.end(); itr++)
667  {
668  (*itr)->SetParameters(fitparameters);
669  (*itr)->SetLineColor(2);
670  (*itr)->SetLineWidth(2);
671  (*itr)->Write();
672  }
673 
674  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter)
675  {
676  const double residX = (*resiter)[kResidX];
677  const double residY = (*resiter)[kResidY];
678  const double resslopeX = (*resiter)[kResSlopeX];
679  const double resslopeY = (*resiter)[kResSlopeY];
680  const double positionX = (*resiter)[kPositionX];
681  const double positionY = (*resiter)[kPositionY];
682  const double angleX = (*resiter)[kAngleX];
683  const double angleY = (*resiter)[kAngleY];
684  const double redchi2 = (*resiter)[kRedChi2];
685  double weight = 1./redchi2;
686  if (!m_weightAlignment) weight = 1.;
687 
688  if (!m_weightAlignment || TMath::Prob(redchi2*12, 12) < 0.99) { // no spikes allowed
689 
690  hist_alphax->Fill(1000.*resslopeX, 10.*residX);
691  hist_alphay->Fill(1000.*resslopeY, 10.*residY);
692 
693  double coefX = value(kAlphaX), coefY = value(kAlphaY);
694  if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) coefX = coefY = 0.;
695  double geom_residX = residual_x(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY, coefX, resslopeX);
696  hist_x->Fill(10.*(residX - geom_residX + value(kAlignX)), weight);
697  hist_x_trackx->Fill(positionX, 10.*residX, weight);
698  hist_x_tracky->Fill(positionY, 10.*residX, weight);
699  hist_x_trackdxdz->Fill(angleX, 10.*residX, weight);
700  hist_x_trackdydz->Fill(angleY, 10.*residX, weight);
701  fit_x_trackx->Fill(positionX, 10.*geom_residX, weight);
702  fit_x_tracky->Fill(positionY, 10.*geom_residX, weight);
703  fit_x_trackdxdz->Fill(angleX, 10.*geom_residX, weight);
704  fit_x_trackdydz->Fill(angleY, 10.*geom_residX, weight);
705 
706  double geom_residY = residual_y(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY, coefY, resslopeY);
707  hist_y->Fill(10.*(residY - geom_residY + value(kAlignY)), weight);
708  hist_y_trackx->Fill(positionX, 10.*residY, weight);
709  hist_y_tracky->Fill(positionY, 10.*residY, weight);
710  hist_y_trackdxdz->Fill(angleX, 10.*residY, weight);
711  hist_y_trackdydz->Fill(angleY, 10.*residY, weight);
712  fit_y_trackx->Fill(positionX, 10.*geom_residY, weight);
713  fit_y_tracky->Fill(positionY, 10.*geom_residY, weight);
714  fit_y_trackdxdz->Fill(angleX, 10.*geom_residY, weight);
715  fit_y_trackdydz->Fill(angleY, 10.*geom_residY, weight);
716 
717  double geom_resslopeX = residual_dxdz(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY);
718  hist_dxdz->Fill(1000.*(resslopeX - geom_resslopeX + value(kAlignPhiY)), weight);
719  hist_dxdz_trackx->Fill(positionX, 1000.*resslopeX, weight);
720  hist_dxdz_tracky->Fill(positionY, 1000.*resslopeX, weight);
721  hist_dxdz_trackdxdz->Fill(angleX, 1000.*resslopeX, weight);
722  hist_dxdz_trackdydz->Fill(angleY, 1000.*resslopeX, weight);
723  fit_dxdz_trackx->Fill(positionX, 1000.*geom_resslopeX, weight);
724  fit_dxdz_tracky->Fill(positionY, 1000.*geom_resslopeX, weight);
725  fit_dxdz_trackdxdz->Fill(angleX, 1000.*geom_resslopeX, weight);
726  fit_dxdz_trackdydz->Fill(angleY, 1000.*geom_resslopeX, weight);
727 
728  double geom_resslopeY = residual_dydz(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY);
729  hist_dydz->Fill(1000.*(resslopeY - geom_resslopeY - value(kAlignPhiX)), weight);
730  hist_dydz_trackx->Fill(positionX, 1000.*resslopeY, weight);
731  hist_dydz_tracky->Fill(positionY, 1000.*resslopeY, weight);
732  hist_dydz_trackdxdz->Fill(angleX, 1000.*resslopeY, weight);
733  hist_dydz_trackdydz->Fill(angleY, 1000.*resslopeY, weight);
734  fit_dydz_trackx->Fill(positionX, 1000.*geom_resslopeY, weight);
735  fit_dydz_tracky->Fill(positionY, 1000.*geom_resslopeY, weight);
736  fit_dydz_trackdxdz->Fill(angleX, 1000.*geom_resslopeY, weight);
737  fit_dydz_trackdydz->Fill(angleY, 1000.*geom_resslopeY, weight);
738  }
739 
740  hist_x_raw->Fill(10.*residX);
741  hist_y_raw->Fill(10.*residY);
742  hist_dxdz_raw->Fill(1000.*resslopeX);
743  hist_dydz_raw->Fill(1000.*resslopeY);
744  if (fabs(resslopeX) < 0.005) hist_x_cut->Fill(10.*residX);
745  if (fabs(resslopeY) < 0.030) hist_y_cut->Fill(10.*residY);
746  }
747 
748  double chi2 = 0.;
749  double ndof = 0.;
750  for (int i = 1; i <= hist_x->GetNbinsX(); i++) {
751  double xi = hist_x->GetBinCenter(i);
752  double yi = hist_x->GetBinContent(i);
753  double yerri = hist_x->GetBinError(i);
754  double yth = fit_x->Eval(xi);
755  if (yerri > 0.) {
756  chi2 += pow((yth - yi)/yerri, 2);
757  ndof += 1.;
758  }
759  }
760  for (int i = 1; i <= hist_y->GetNbinsX(); i++) {
761  double xi = hist_y->GetBinCenter(i);
762  double yi = hist_y->GetBinContent(i);
763  double yerri = hist_y->GetBinError(i);
764  double yth = fit_y->Eval(xi);
765  if (yerri > 0.) {
766  chi2 += pow((yth - yi)/yerri, 2);
767  ndof += 1.;
768  }
769  }
770  for (int i = 1; i <= hist_dxdz->GetNbinsX(); i++) {
771  double xi = hist_dxdz->GetBinCenter(i);
772  double yi = hist_dxdz->GetBinContent(i);
773  double yerri = hist_dxdz->GetBinError(i);
774  double yth = fit_dxdz->Eval(xi);
775  if (yerri > 0.) {
776  chi2 += pow((yth - yi)/yerri, 2);
777  ndof += 1.;
778  }
779  }
780  for (int i = 1; i <= hist_dydz->GetNbinsX(); i++) {
781  double xi = hist_dydz->GetBinCenter(i);
782  double yi = hist_dydz->GetBinContent(i);
783  double yerri = hist_dydz->GetBinError(i);
784  double yth = fit_dydz->Eval(xi);
785  if (yerri > 0.) {
786  chi2 += pow((yth - yi)/yerri, 2);
787  ndof += 1.;
788  }
789  }
790  ndof -= npar();
791 
792  return (ndof > 0. ? chi2 / ndof : -1.);
793 }
794 
795 
796 TTree * MuonResiduals6DOFFitter::readNtuple(std::string fname, unsigned int wheel, unsigned int station, unsigned int sector, unsigned int preselected)
797 {
798  TFile *f = new TFile(fname.c_str());
799  TTree *t = (TTree*)f->Get("mual_ttree");
800 
801  // Create new temporary file
802  TFile *tmpf = new TFile("small_tree_fit.root","recreate");
803  assert(tmpf!=0);
804 
805  // filter the tree (temporarily lives in the new file)
806  TTree *tt = t->CopyTree(Form("is_dt && ring_wheel==%d && station==%d && sector==%d && select==%d", wheel, station, sector, (bool)preselected));
807 
808  MuonAlignmentTreeRow r;
809  tt->SetBranchAddress("res_x", &r.res_x);
810  tt->SetBranchAddress("res_slope_x", &r.res_slope_x);
811  tt->SetBranchAddress("res_y", &r.res_y);
812  tt->SetBranchAddress("res_slope_y", &r.res_slope_y);
813  tt->SetBranchAddress("pos_x", &r.pos_x);
814  tt->SetBranchAddress("pos_y", &r.pos_y);
815  tt->SetBranchAddress("angle_x", &r.angle_x);
816  tt->SetBranchAddress("angle_y", &r.angle_y);
817  tt->SetBranchAddress("pz", &r.pz);
818  tt->SetBranchAddress("pt", &r.pt);
819  tt->SetBranchAddress("q", &r.q);
820 
821  Long64_t nentries = tt->GetEntries();
822  for (Long64_t i=0;i<nentries; i++)
823  {
824  tt->GetEntry(i);
825  double *rdata = new double[MuonResiduals6DOFFitter::kNData];
826  rdata[kResidX] = r.res_x;
827  rdata[kResidY] = r.res_y;
828  rdata[kResSlopeX] = r.res_slope_x;
829  rdata[kResSlopeY] = r.res_slope_y;
830  rdata[kPositionX] = r.pos_x;
831  rdata[kPositionY] = r.pos_y;
832  rdata[kAngleX] = r.angle_x;
833  rdata[kAngleY] = r.angle_y;
834  rdata[kRedChi2] = 0.1;
835  rdata[kPz] = r.pz;
836  rdata[kPt] = r.pt;
837  rdata[kCharge] = r.q;
838  fill(rdata);
839  }
840  delete f;
841  //delete tmpf;
842  return tt;
843 }
AlignableSurface::width
align::Scalar width() const
Definition: AlignableSurface.h:36
MuonResidualsFitter_logGaussPowerTails
double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma)
Definition: MuonResidualsFitter.cc:140
i
int i
Definition: DBlmapReader.cc:9
dqm_diff.start
tuple start
Check for commandline option errors.
Definition: dqm_diff.py:58
create_public_lumi_plots.width
int width
Definition: create_public_lumi_plots.py:1009
MuonResidualsFitter_powerLawTails_TF1
Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:122
cscdqm::h::names
static const HistoName names[]
Definition: CSCDQM_HistoDef.h:831
MuonResidualsFitter_logPureGaussian2D
double MuonResidualsFitter_logPureGaussian2D(double x, double y, double x0, double y0, double sx, double sy, double r)
Definition: MuonResidualsFitter.cc:46
relval_machine.steps
tuple steps
Definition: relval_machine.py:8
assert
assert(m_qm.get())
puppiForMET_cff.weight
tuple weight
Definition: puppiForMET_cff.py:35
MuonResidualsFitter::fill
void fill(double *residual)
Definition: MuonResidualsFitter.cc:205
MuonResidualsFitterFitInfo::fitter
MuonResidualsFitter * fitter()
Definition: MuonResidualsFitter.h:252
NULL
#define NULL
Definition: scimark2.h:8
MuonResidualsFitter::value
double value(int parNum)
Definition: MuonResidualsFitter.h:161
MuonResiduals6DOFFitter::plot
double plot(std::string name, TFileDirectory *dir, Alignable *ali)
Definition: MuonResiduals6DOFFitter.cc:341
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:301
MuonResidualsFitter_logROOTVoigt
double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma)
Definition: MuonResidualsFitter.cc:128
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
MuonResiduals6DOFFitter::readNtuple
TTree * readNtuple(std::string fname, unsigned int wheel, unsigned int station, unsigned int sector, unsigned int preselected=1)
Definition: MuonResiduals6DOFFitter.cc:796
HLT_FULL_cff.chi2
tuple chi2
Definition: HLT_FULL_cff.py:61454
MuonResidualsFitter_logPowerLawTails
double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma)
Definition: MuonResidualsFitter.cc:86
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:153
MuonResidualsFitter::fix
void fix(int parNum, bool dofix=true)
Definition: MuonResidualsFitter.cc:189
TFileDirectory::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileDirectory.h:59
MuonResidualsFitter::correctBField
virtual void correctBField()=0
Alignable::surface
const AlignableSurface & surface() const
Return the Surface (global position and orientation) of the object.
Definition: Alignable.h:131
MuonResidualsFitter_pureGaussian_TF1
Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:39
MuonResiduals6DOFFitter::inform
void inform(TMinuit *tMinuit)
Definition: MuonResiduals6DOFFitter.cc:219
AlignableSurface::length
align::Scalar length() const
Definition: AlignableSurface.h:38
MuonResidualsFitter::errorerror
double errorerror(int parNum)
Definition: MuonResidualsFitter.h:162
customizeTrackingMonitorSeedNumber.idx
tuple idx
DEBUGGING if hasattr(process,&quot;trackMonIterativeTracking2012&quot;): print &quot;trackMonIterativeTracking2012 D...
Definition: customizeTrackingMonitorSeedNumber.py:15
MuonResidualsFitter::residuals_end
std::vector< double * >::const_iterator residuals_end() const
Definition: MuonResidualsFitter.h:199
MuonResidualsFitter::useRes
int useRes(int pattern=-1)
Definition: MuonResidualsFitter.h:139
alignmentValidation.fname
string fname
main script
Definition: alignmentValidation.py:957
relval_parameters_module.step
string step
Definition: relval_parameters_module.py:36
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
dir
dbl *** dir
Definition: mlp_gen.cc:35
MuonResidualsFitter_ROOTVoigt_TF1
Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:134
MuonResidualsFitter::residuals_begin
std::vector< double * >::const_iterator residuals_begin() const
Definition: MuonResidualsFitter.h:198
MuonResiduals6DOFFitter_FCN
void MuonResiduals6DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
Definition: MuonResiduals6DOFFitter.cc:88
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:40
Generated for CMSSW Reference Manual by doxygen 1.8.5