CMS 3D CMS Logo

MuonResiduals5DOFFitter.cc
Go to the documentation of this file.
1 // $Id:$
2 
3 #ifdef STANDALONE_FITTER
5 #else
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_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 alpha, double resslope)
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  + resslope * alpha;
33  }
34 
35  double residual_dxdz(double delta_x, 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  {
39  return -track_dxdz * track_dydz * delta_phix
40  + (1. + track_dxdz * track_dxdz) * delta_phiy
41  - track_dydz * delta_phiz;
42  }
43 
44  Double_t residual_x_trackx_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9], p[10], p[11]); }
45  Double_t residual_x_tracky_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9], p[10], p[11]); }
46  Double_t residual_x_trackdxdz_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9], p[10], p[11]); }
47  Double_t residual_x_trackdydz_TF1(Double_t *xx, Double_t *p) { return 10.*residual_x(p[0], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0], p[10], p[11]); }
48 
49  Double_t residual_dxdz_trackx_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9]); }
50  Double_t residual_dxdz_tracky_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9]); }
51  Double_t residual_dxdz_trackdxdz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9]); }
52  Double_t residual_dxdz_trackdydz_TF1(Double_t *xx, Double_t *p) { return 1000.*residual_dxdz(p[0], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0]); }
53 }
54 
55 
56 void MuonResiduals5DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
57 {
58  MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo*)(minuit->GetObjectFit());
59  MuonResidualsFitter *fitter = fitinfo->fitter();
60 
61  fval = 0.;
62  for (std::vector<double*>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end(); ++resiter)
63  {
64  const double residual = (*resiter)[MuonResiduals5DOFFitter::kResid];
65  const double resslope = (*resiter)[MuonResiduals5DOFFitter::kResSlope];
66  const double positionX = (*resiter)[MuonResiduals5DOFFitter::kPositionX];
67  const double positionY = (*resiter)[MuonResiduals5DOFFitter::kPositionY];
68  const double angleX = (*resiter)[MuonResiduals5DOFFitter::kAngleX];
69  const double angleY = (*resiter)[MuonResiduals5DOFFitter::kAngleY];
70  const double redchi2 = (*resiter)[MuonResiduals5DOFFitter::kRedChi2];
71 
72  const double alignx = par[MuonResiduals5DOFFitter::kAlignX];
73  const double alignz = par[MuonResiduals5DOFFitter::kAlignZ];
74  const double alignphix = par[MuonResiduals5DOFFitter::kAlignPhiX];
75  const double alignphiy = par[MuonResiduals5DOFFitter::kAlignPhiY];
76  const double alignphiz = par[MuonResiduals5DOFFitter::kAlignPhiZ];
77  const double residsigma = par[MuonResiduals5DOFFitter::kResidSigma];
78  const double resslopesigma = par[MuonResiduals5DOFFitter::kResSlopeSigma];
79  const double alpha = par[MuonResiduals5DOFFitter::kAlpha];
80  const double residgamma = par[MuonResiduals5DOFFitter::kResidGamma];
81  const double resslopegamma = par[MuonResiduals5DOFFitter::kResSlopeGamma];
82 
83  double coeff = alpha;
86  double residpeak = residual_x(alignx, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coeff, resslope);
87  double resslopepeak = residual_dxdz(alignx, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
88 
89  double weight = (1./redchi2) * number_of_hits / sum_of_weights;
90  if (!weight_alignment) weight = 1.;
91 
92  if (!weight_alignment || TMath::Prob(redchi2*8, 8) < 0.99) // no spikes allowed
93  {
96  fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
97  fval += -weight * MuonResidualsFitter_logPureGaussian(resslope, resslopepeak, resslopesigma);
98  }
99  else if (fitter->useRes() == MuonResidualsFitter::k1100) {
100  fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
101  }
102  else if (fitter->useRes() == MuonResidualsFitter::k0010) {
103  fval += -weight * MuonResidualsFitter_logPureGaussian(resslope, resslopepeak, resslopesigma);
104  }
105  }
107  if (fitter->useRes() == MuonResidualsFitter::k1111 || fitter->useRes() == MuonResidualsFitter::k1110 || fitter->useRes() == MuonResidualsFitter::k1010) {
108  fval += -weight * MuonResidualsFitter_logPureGaussian2D(residual, resslope, residpeak, resslopepeak, residsigma, resslopesigma, alpha);
109  }
110  else if (fitter->useRes() == MuonResidualsFitter::k1100) {
111  fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
112  }
113  else if (fitter->useRes() == MuonResidualsFitter::k0010) {
114  fval += -weight * MuonResidualsFitter_logPureGaussian(resslope, resslopepeak, resslopesigma);
115  }
116  }
117  else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
118  fval += -weight * MuonResidualsFitter_logPowerLawTails(residual, residpeak, residsigma, residgamma);
119  fval += -weight * MuonResidualsFitter_logPowerLawTails(resslope, resslopepeak, resslopesigma, resslopegamma);
120  }
121  else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
122  fval += -weight * MuonResidualsFitter_logROOTVoigt(residual, residpeak, residsigma, residgamma);
123  fval += -weight * MuonResidualsFitter_logROOTVoigt(resslope, resslopepeak, resslopesigma, resslopegamma);
124  }
126  fval += -weight * MuonResidualsFitter_logGaussPowerTails(residual, residpeak, residsigma);
127  fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslope, resslopepeak, resslopesigma);
128  }
129  else { assert(false); }
130  }
131  }
132 }
133 
134 
136 {
138 }
139 
140 
141 void MuonResiduals5DOFFitter::inform(TMinuit *tMinuit)
142 {
143  minuit = tMinuit;
144 }
145 
146 
148 {
149  sum_of_weights = 0.;
150  number_of_hits = 0.;
151  weight_alignment = m_weightAlignment;
152  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
153  if (m_weightAlignment) {
154  double redchi2 = (*resiter)[MuonResiduals5DOFFitter::kRedChi2];
155  if (TMath::Prob(redchi2*8, 8) < 0.99) {
156  sum_of_weights += 1./redchi2;
157  number_of_hits += 1.;
158  }
159  }
160  else {
161  sum_of_weights += 1.;
162  number_of_hits += 1.;
163  }
164  }
165  return sum_of_weights;
166 }
167 
168 
170 {
171  initialize_table(); // if not already initialized
172  sumofweights();
173 
174  double res_std = 0.5;
175  double resslope_std = 0.005;
176 
178  std::string names[10] = {"AlignX","AlignZ","AlignPhiX","AlignPhiY","AlignPhiZ", "ResidSigma","ResSlopeSigma", "Alpha", "ResidGamma","ResSlopeGamma"};
179  double starts[10] = {0., 0., 0., 0., 0., res_std, resslope_std, 0., 0.1*res_std, 0.1*resslope_std};
180  double steps[10] = {0.1, 0.1, 0.001, 0.001, 0.001, 0.001*res_std, 0.001*resslope_std, 0.001, 0.01*res_std, 0.01*resslope_std};
181  double lows[10] = {0., 0., 0., 0., 0., 0., 0., -1., 0., 0.};
182  double highs[10] = {0., 0., 0., 0., 0., 10., 0.1, 1., 0., 0.};
183 
184  std::vector<int> num(nums, nums+5);
185  std::vector<std::string> name(names, names+5);
186  std::vector<double> start(starts, starts+5);
187  std::vector<double> step(steps, steps+5);
188  std::vector<double> low(lows, lows+5);
189  std::vector<double> high(highs, highs+5);
190 
191  bool add_alpha = ( residualsModel() == kPureGaussian2D);
192  bool add_gamma = ( residualsModel() == kROOTVoigt || residualsModel() == kPowerLawTails);
193 
194  int idx[4], ni = 0;
195  if (useRes() == k1111 || useRes() == k1110 || useRes() == k1010) {
196  for(ni=0; ni<2; ni++) idx[ni] = ni+5;
197  if (add_alpha) idx[ni++] = 7;
198  else if (add_gamma) for(; ni<4; ni++) idx[ni] = ni+6;
199  if (!add_alpha) fix(kAlpha);
200  }
201  else if (useRes() == k1100) {
202  idx[ni++] = 5;
203  if (add_gamma) idx[ni++] = 8;
205  fix(kAlpha);
206  }
207  else if (useRes() == k0010) {
208  idx[ni++] = 6;
209  if (add_gamma) idx[ni++] = 9;
210  fix(kResidSigma);
211  fix(kAlpha);
212  }
213  for (int i=0; i<ni; i++){
214  num.push_back(nums[idx[i]]);
215  name.push_back(names[idx[i]]);
216  start.push_back(starts[idx[i]]);
217  step.push_back(steps[idx[i]]);
218  low.push_back(lows[idx[i]]);
219  high.push_back(highs[idx[i]]);
220  }
221 
222  return dofit(&MuonResiduals5DOFFitter_FCN, num, name, start, step, low, high);
223 }
224 
225 
227 {
228  sumofweights();
229 
230  double mean_residual = 0., mean_resslope = 0.;
231  double mean_trackx = 0., mean_tracky = 0., mean_trackdxdz = 0., mean_trackdydz = 0.;
232  double sum_w = 0.;
233 
234  for (std::vector<double*>::const_iterator rit = residuals_begin(); rit != residuals_end(); ++rit)
235  {
236  const double redchi2 = (*rit)[kRedChi2];
237  double weight = 1./redchi2;
238  if (!m_weightAlignment) weight = 1.;
239 
240  if (!m_weightAlignment || TMath::Prob(redchi2*6, 6) < 0.99) // no spikes allowed
241  {
242  double factor_w = 1./(sum_w + weight);
243  mean_residual = factor_w * (sum_w * mean_residual + weight * (*rit)[kResid]);
244  mean_resslope = factor_w * (sum_w * mean_resslope + weight * (*rit)[kResSlope]);
245  mean_trackx = factor_w * (sum_w * mean_trackx + weight * (*rit)[kPositionX]);
246  mean_tracky = factor_w * (sum_w * mean_tracky + weight * (*rit)[kPositionY]);
247  mean_trackdxdz = factor_w * (sum_w * mean_trackdxdz + weight * (*rit)[kAngleX]);
248  mean_trackdydz = factor_w * (sum_w * mean_trackdydz + weight * (*rit)[kAngleY]);
249  sum_w += weight;
250  }
251  }
252 
253  std::string name_residual, name_resslope, name_residual_raw, name_resslope_raw, name_residual_cut, name_alpha;
254  std::string name_residual_trackx, name_resslope_trackx;
255  std::string name_residual_tracky, name_resslope_tracky;
256  std::string name_residual_trackdxdz, name_resslope_trackdxdz;
257  std::string name_residual_trackdydz, name_resslope_trackdydz;
258 
259  name_residual = name + "_residual";
260  name_resslope = name + "_resslope";
261  name_residual_raw = name + "_residual_raw";
262  name_resslope_raw = name + "_resslope_raw";
263  name_residual_cut = name + "_residual_cut";
264  name_alpha = name + "_alpha";
265  name_residual_trackx = name + "_residual_trackx";
266  name_resslope_trackx = name + "_resslope_trackx";
267  name_residual_tracky = name + "_residual_tracky";
268  name_resslope_tracky = name + "_resslope_tracky";
269  name_residual_trackdxdz = name + "_residual_trackdxdz";
270  name_resslope_trackdxdz = name + "_resslope_trackdxdz";
271  name_residual_trackdydz = name + "_residual_trackdydz";
272  name_resslope_trackdydz = name + "_resslope_trackdydz";
273 
274  double width = ali->surface().width();
275  double length = ali->surface().length();
276  int bins_residual = 150, bins_resslope = 100;
277  double min_residual = -75., max_residual = 75.;
278  double min_resslope = -50., max_resslope = 50.;
279  double min_trackx = -width/2., max_trackx = width/2.;
280  double min_tracky = -length/2., max_tracky = length/2.;
281  double min_trackdxdz = -1.5, max_trackdxdz = 1.5;
282  double min_trackdydz = -1.5, max_trackdydz = 1.5;
283 
284  TH1F *hist_residual = dir->make<TH1F>(name_residual.c_str(), "", bins_residual, min_residual, max_residual);
285  TH1F *hist_resslope = dir->make<TH1F>(name_resslope.c_str(), "", bins_resslope, min_resslope, max_resslope);
286  TH1F *hist_residual_raw = dir->make<TH1F>(name_residual_raw.c_str(), "", bins_residual, min_residual, max_residual);
287  TH1F *hist_resslope_raw = dir->make<TH1F>(name_resslope_raw.c_str(), "", bins_resslope, min_resslope, max_resslope);
288  TH1F *hist_residual_cut = dir->make<TH1F>(name_residual_cut.c_str(), "", bins_residual, min_residual, max_residual);
289  TH2F *hist_alpha = dir->make<TH2F>(name_alpha.c_str(), "", 50, min_resslope, max_resslope, 50, -50., 50.);
290 
291  TProfile *hist_residual_trackx = dir->make<TProfile>(name_residual_trackx.c_str(), "", 50, min_trackx, max_trackx, min_residual, max_residual);
292  TProfile *hist_resslope_trackx = dir->make<TProfile>(name_resslope_trackx.c_str(), "", 50, min_trackx, max_trackx, min_resslope, max_resslope);
293  TProfile *hist_residual_tracky = dir->make<TProfile>(name_residual_tracky.c_str(), "", 50, min_tracky, max_tracky, min_residual, max_residual);
294  TProfile *hist_resslope_tracky = dir->make<TProfile>(name_resslope_tracky.c_str(), "", 50, min_tracky, max_tracky, min_resslope, max_resslope);
295  TProfile *hist_residual_trackdxdz = dir->make<TProfile>(name_residual_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_residual, max_residual);
296  TProfile *hist_resslope_trackdxdz = dir->make<TProfile>(name_resslope_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_resslope, max_resslope);
297  TProfile *hist_residual_trackdydz = dir->make<TProfile>(name_residual_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_residual, max_residual);
298  TProfile *hist_resslope_trackdydz = dir->make<TProfile>(name_resslope_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_resslope, max_resslope);
299 
300  hist_residual_trackx->SetAxisRange(-10., 10., "Y");
301  hist_resslope_trackx->SetAxisRange(-10., 10., "Y");
302  hist_residual_tracky->SetAxisRange(-10., 10., "Y");
303  hist_resslope_tracky->SetAxisRange(-10., 10., "Y");
304  hist_residual_trackdxdz->SetAxisRange(-10., 10., "Y");
305  hist_resslope_trackdxdz->SetAxisRange(-10., 10., "Y");
306  hist_residual_trackdydz->SetAxisRange(-10., 10., "Y");
307  hist_resslope_trackdydz->SetAxisRange(-10., 10., "Y");
308 
309  name_residual += "_fit";
310  name_resslope += "_fit";
311  name_alpha += "_fit";
312  name_residual_trackx += "_fit";
313  name_resslope_trackx += "_fit";
314  name_residual_tracky += "_fit";
315  name_resslope_tracky += "_fit";
316  name_residual_trackdxdz += "_fit";
317  name_resslope_trackdxdz += "_fit";
318  name_residual_trackdydz += "_fit";
319  name_resslope_trackdydz += "_fit";
320 
321  TF1 *fit_residual = NULL;
322  TF1 *fit_resslope = NULL;
324  fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_residual, max_residual, 3);
325  fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual)/bins_residual, 10.*value(kAlignX), 10.*value(kResidSigma));
326  const double er_res[3] = {0., 10.*errorerror(kAlignX), 10.*errorerror(kResidSigma)};
327  fit_residual->SetParErrors(er_res);
328  fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_resslope, max_resslope, 3);
329  fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope)/bins_resslope, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma));
330  const double er_resslope[3] = {0., 1000.*errorerror(kAlignPhiY), 1000.*errorerror(kResSlopeSigma)};
331  fit_resslope->SetParErrors(er_resslope);
332  }
333  else if (residualsModel() == kPowerLawTails) {
334  fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_residual, max_residual, 4);
335  fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual)/bins_residual, 10.*value(kAlignX), 10.*value(kResidSigma), 10.*value(kResidGamma));
336  fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_resslope, max_resslope, 4);
337  fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope)/bins_resslope, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma), 1000.*value(kResSlopeGamma));
338  }
339  else if (residualsModel() == kROOTVoigt) {
340  fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_residual, max_residual, 4);
341  fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual)/bins_residual, 10.*value(kAlignX), 10.*value(kResidSigma), 10.*value(kResidGamma));
342  fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_resslope, max_resslope, 4);
343  fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope)/bins_resslope, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma), 1000.*value(kResSlopeGamma));
344  }
345  else if (residualsModel() == kGaussPowerTails) {
346  fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_residual, max_residual, 3);
347  fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual)/bins_residual, 10.*value(kAlignX), 10.*value(kResidSigma));
348  fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_resslope, max_resslope, 3);
349  fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope)/bins_resslope, 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma));
350  }
351  else { assert(false); }
352 
353  fit_residual->SetLineColor(2); fit_residual->SetLineWidth(2); fit_residual->Write();
354  fit_resslope->SetLineColor(2); fit_resslope->SetLineWidth(2); fit_resslope->Write();
355 
356  TF1 *fit_alpha = new TF1(name_alpha.c_str(), "[0] + x*[1]", min_resslope, max_resslope);
357  double a = 10.*value(kAlignX), b = 10.*value(kAlpha)/1000.;
359  {
360  double sx = 10.*value(kResidSigma), sy = 1000.*value(kResSlopeSigma), r = value(kAlpha);
361  a = mean_residual;
362  b = 0.;
363  if ( sx != 0. )
364  {
365  b = 1./(sy/sx*r);
366  a = - b * mean_resslope;
367  }
368  }
369  fit_alpha->SetParameters(a, b);
370  fit_alpha->SetLineColor(2); fit_alpha->SetLineWidth(2); fit_alpha->Write();
371 
372  TProfile *fit_residual_trackx = dir->make<TProfile>(name_residual_trackx.c_str(), "", 50, min_trackx, max_trackx);
373  TProfile *fit_resslope_trackx = dir->make<TProfile>(name_resslope_trackx.c_str(), "", 50, min_trackx, max_trackx);
374  TProfile *fit_residual_tracky = dir->make<TProfile>(name_residual_tracky.c_str(), "", 50, min_tracky, max_tracky);
375  TProfile *fit_resslope_tracky = dir->make<TProfile>(name_resslope_tracky.c_str(), "", 50, min_tracky, max_tracky);
376  TProfile *fit_residual_trackdxdz = dir->make<TProfile>(name_residual_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz);
377  TProfile *fit_resslope_trackdxdz = dir->make<TProfile>(name_resslope_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz);
378  TProfile *fit_residual_trackdydz = dir->make<TProfile>(name_residual_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz);
379  TProfile *fit_resslope_trackdydz = dir->make<TProfile>(name_resslope_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz);
380 
381  fit_residual_trackx->SetLineColor(2); fit_residual_trackx->SetLineWidth(2);
382  fit_resslope_trackx->SetLineColor(2); fit_resslope_trackx->SetLineWidth(2);
383  fit_residual_tracky->SetLineColor(2); fit_residual_tracky->SetLineWidth(2);
384  fit_resslope_tracky->SetLineColor(2); fit_resslope_tracky->SetLineWidth(2);
385  fit_residual_trackdxdz->SetLineColor(2); fit_residual_trackdxdz->SetLineWidth(2);
386  fit_resslope_trackdxdz->SetLineColor(2); fit_resslope_trackdxdz->SetLineWidth(2);
387  fit_residual_trackdydz->SetLineColor(2); fit_residual_trackdydz->SetLineWidth(2);
388  fit_resslope_trackdydz->SetLineColor(2); fit_resslope_trackdydz->SetLineWidth(2);
389 
390  name_residual_trackx += "line";
391  name_resslope_trackx += "line";
392  name_residual_tracky += "line";
393  name_resslope_tracky += "line";
394  name_residual_trackdxdz += "line";
395  name_resslope_trackdxdz += "line";
396  name_residual_trackdydz += "line";
397  name_resslope_trackdydz += "line";
398 
399  TF1 *fitline_residual_trackx = new TF1(name_residual_trackx.c_str(), residual_x_trackx_TF1, min_trackx, max_trackx, 12);
400  TF1 *fitline_resslope_trackx = new TF1(name_resslope_trackx.c_str(), residual_dxdz_trackx_TF1, min_trackx, max_trackx, 12);
401  TF1 *fitline_residual_tracky = new TF1(name_residual_tracky.c_str(), residual_x_tracky_TF1, min_tracky, max_tracky, 12);
402  TF1 *fitline_resslope_tracky = new TF1(name_resslope_tracky.c_str(), residual_dxdz_tracky_TF1, min_tracky, max_tracky, 12);
403  TF1 *fitline_residual_trackdxdz = new TF1(name_residual_trackdxdz.c_str(), residual_x_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
404  TF1 *fitline_resslope_trackdxdz = new TF1(name_resslope_trackdxdz.c_str(), residual_dxdz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
405  TF1 *fitline_residual_trackdydz = new TF1(name_residual_trackdydz.c_str(), residual_x_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);
406  TF1 *fitline_resslope_trackdydz = new TF1(name_resslope_trackdydz.c_str(), residual_dxdz_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);
407 
408  std::vector<TF1*> fitlines;
409  fitlines.push_back(fitline_residual_trackx);
410  fitlines.push_back(fitline_resslope_trackx);
411  fitlines.push_back(fitline_residual_tracky);
412  fitlines.push_back(fitline_resslope_tracky);
413  fitlines.push_back(fitline_residual_trackdxdz);
414  fitlines.push_back(fitline_resslope_trackdxdz);
415  fitlines.push_back(fitline_residual_trackdydz);
416  fitlines.push_back(fitline_resslope_trackdydz);
417 
418  double fitparameters[12] = {value(kAlignX), 0., value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ),
419  mean_trackx, mean_tracky, mean_trackdxdz, mean_trackdydz, value(kAlpha), mean_resslope};
420  if (residualsModel() == kPureGaussian2D) fitparameters[10] = 0.;
421 
422  for(std::vector<TF1*>::const_iterator itr = fitlines.begin(); itr != fitlines.end(); itr++)
423  {
424  (*itr)->SetParameters(fitparameters);
425  (*itr)->SetLineColor(2);
426  (*itr)->SetLineWidth(2);
427  (*itr)->Write();
428  }
429 
430  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
431  const double resid = (*resiter)[kResid];
432  const double resslope = (*resiter)[kResSlope];
433  const double positionX = (*resiter)[kPositionX];
434  const double positionY = (*resiter)[kPositionY];
435  const double angleX = (*resiter)[kAngleX];
436  const double angleY = (*resiter)[kAngleY];
437  const double redchi2 = (*resiter)[kRedChi2];
438  double weight = 1./redchi2;
439  if (!m_weightAlignment) weight = 1.;
440 
441  if (!m_weightAlignment || TMath::Prob(redchi2*8, 8) < 0.99) { // no spikes allowed
442  hist_alpha->Fill(1000.*resslope, 10.*resid);
443 
444  double coeff = value(kAlpha);
445  if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) coeff = 0.;
446  double geom_resid = residual_x(value(kAlignX), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY, coeff, resslope);
447  hist_residual->Fill(10.*(resid - geom_resid + value(kAlignX)), weight);
448  hist_residual_trackx->Fill(positionX, 10.*resid, weight);
449  hist_residual_tracky->Fill(positionY, 10.*resid, weight);
450  hist_residual_trackdxdz->Fill(angleX, 10.*resid, weight);
451  hist_residual_trackdydz->Fill(angleY, 10.*resid, weight);
452  fit_residual_trackx->Fill(positionX, 10.*geom_resid, weight);
453  fit_residual_tracky->Fill(positionY, 10.*geom_resid, weight);
454  fit_residual_trackdxdz->Fill(angleX, 10.*geom_resid, weight);
455  fit_residual_trackdydz->Fill(angleY, 10.*geom_resid, weight);
456 
457  double geom_resslope = residual_dxdz(value(kAlignX), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY);
458  hist_resslope->Fill(1000.*(resslope - geom_resslope + value(kAlignPhiY)), weight);
459  hist_resslope_trackx->Fill(positionX, 1000.*resslope, weight);
460  hist_resslope_tracky->Fill(positionY, 1000.*resslope, weight);
461  hist_resslope_trackdxdz->Fill(angleX, 1000.*resslope, weight);
462  hist_resslope_trackdydz->Fill(angleY, 1000.*resslope, weight);
463  fit_resslope_trackx->Fill(positionX, 1000.*geom_resslope, weight);
464  fit_resslope_tracky->Fill(positionY, 1000.*geom_resslope, weight);
465  fit_resslope_trackdxdz->Fill(angleX, 1000.*geom_resslope, weight);
466  fit_resslope_trackdydz->Fill(angleY, 1000.*geom_resslope, weight);
467  }
468 
469  hist_residual_raw->Fill(10.*resid);
470  hist_resslope_raw->Fill(1000.*resslope);
471  if (fabs(resslope) < 0.005) hist_residual_cut->Fill(10.*resid);
472  }
473 
474  double chi2 = 0.;
475  double ndof = 0.;
476  for (int i = 1; i <= hist_residual->GetNbinsX(); i++) {
477  double xi = hist_residual->GetBinCenter(i);
478  double yi = hist_residual->GetBinContent(i);
479  double yerri = hist_residual->GetBinError(i);
480  double yth = fit_residual->Eval(xi);
481  if (yerri > 0.) {
482  chi2 += pow((yth - yi)/yerri, 2);
483  ndof += 1.;
484  }
485  }
486  for (int i = 1; i <= hist_resslope->GetNbinsX(); i++) {
487  double xi = hist_resslope->GetBinCenter(i);
488  double yi = hist_resslope->GetBinContent(i);
489  double yerri = hist_resslope->GetBinError(i);
490  double yth = fit_resslope->Eval(xi);
491  if (yerri > 0.) {
492  chi2 += pow((yth - yi)/yerri, 2);
493  ndof += 1.;
494  }
495  }
496  ndof -= npar();
497 
498  return (ndof > 0. ? chi2 / ndof : -1.);
499 }
500 
501 
502 TTree * MuonResiduals5DOFFitter::readNtuple(std::string fname, unsigned int wheel, unsigned int station, unsigned int sector, unsigned int preselected)
503 {
504  TFile *f = new TFile(fname.c_str());
505  TTree *t = (TTree*)f->Get("mual_ttree");
506 
507  // Create new temporary file
508  TFile *tmpf = new TFile("small_tree_fit.root","recreate");
509  assert(tmpf!=0);
510 
511  // filter the tree (temporarily lives in the new file)
512  TTree *tt = t->CopyTree(Form("is_dt && ring_wheel==%d && station==%d && sector==%d && select==%d", wheel, station, sector, (bool)preselected));
513 
515  tt->SetBranchAddress("res_x", &r.res_x);
516  tt->SetBranchAddress("res_slope_x", &r.res_slope_x);
517  tt->SetBranchAddress("pos_x", &r.pos_x);
518  tt->SetBranchAddress("pos_y", &r.pos_y);
519  tt->SetBranchAddress("angle_x", &r.angle_x);
520  tt->SetBranchAddress("angle_y", &r.angle_y);
521  tt->SetBranchAddress("pz", &r.pz);
522  tt->SetBranchAddress("pt", &r.pt);
523  tt->SetBranchAddress("q", &r.q);
524 
525  Long64_t nentries = tt->GetEntries();
526  for (Long64_t i=0;i<nentries; i++)
527  {
528  tt->GetEntry(i);
529  double *rdata = new double[MuonResiduals5DOFFitter::kNData];
530  rdata[kResid] = r.res_x;
531  rdata[kResSlope] = r.res_slope_x;
532  rdata[kPositionX] = r.pos_x;
533  rdata[kPositionY] = r.pos_y;
534  rdata[kAngleX] = r.angle_x;
535  rdata[kAngleY] = r.angle_y;
536  rdata[kRedChi2] = 0.1;
537  rdata[kPz] = r.pz;
538  rdata[kPt] = r.pt;
539  rdata[kCharge] = r.q;
540  fill(rdata);
541  }
542  delete f;
543  //delete tmpf;
544  return tt;
545 }
align::Scalar width() const
double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma)
int i
Definition: DBlmapReader.cc:9
float alpha
Definition: AMPTWrapper.h:95
Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par)
static const HistoName names[]
double MuonResidualsFitter_logPureGaussian2D(double x, double y, double x0, double y0, double sx, double sy, double r)
void fill(double *residual)
MuonResidualsFitter * fitter()
#define NULL
Definition: scimark2.h:8
Definition: weight.py:1
double value(int parNum)
TTree * readNtuple(std::string fname, unsigned int wheel, unsigned int station, unsigned int sector, unsigned int preselected=1)
bool dofit(void(*fcn)(int &, double *, double &, double *, int), std::vector< int > &parNum, std::vector< std::string > &parName, std::vector< double > &start, std::vector< double > &step, std::vector< double > &low, std::vector< double > &high)
double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma)
void MuonResiduals5DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma)
double f[11][100]
Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par)
void fix(int parNum, bool dofix=true)
T * make(const Args &...args) const
make new ROOT object
virtual void correctBField()=0
const AlignableSurface & surface() const
Return the Surface (global position and orientation) of the object.
Definition: Alignable.h:135
Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par)
align::Scalar length() const
double plot(std::string name, TFileDirectory *dir, Alignable *ali)
double errorerror(int parNum)
double b
Definition: hdecay.h:120
std::vector< double * >::const_iterator residuals_end() const
int useRes(int pattern=-1)
string fname
main script
double a
Definition: hdecay.h:121
step
dbl *** dir
Definition: mlp_gen.cc:35
void inform(TMinuit *tMinuit)
Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par)
std::vector< double * >::const_iterator residuals_begin() const
double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40