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 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFFitter.h"
2 #include "TH2F.h"
3 // #include "TTree.h"
4 #include "TMath.h"
5 
6 static TMinuit *MuonResiduals6DOFFitter_TMinuit;
7 static double MuonResiduals6DOFFitter_sum_of_weights;
8 static double MuonResiduals6DOFFitter_number_of_hits;
9 static bool MuonResiduals6DOFFitter_weightAlignment;
10 
11 void MuonResiduals6DOFFitter::inform(TMinuit *tMinuit) {
12  MuonResiduals6DOFFitter_TMinuit = tMinuit;
13 }
14 
15 double MuonResiduals6DOFFitter_x(double delta_x, double delta_y, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz, double alphax, double residual_dxdz) {
16  return delta_x - track_dxdz * delta_z - track_y * track_dxdz * delta_phix + track_x * track_dxdz * delta_phiy - track_y * delta_phiz + residual_dxdz * alphax;
17 }
18 
19 double MuonResiduals6DOFFitter_y(double delta_x, double delta_y, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz, double alphay, double residual_dydz) {
20  return delta_y - track_dydz * delta_z - track_y * track_dydz * delta_phix + track_x * track_dydz * delta_phiy + track_x * delta_phiz + residual_dydz * alphay;
21 }
22 
23 double MuonResiduals6DOFFitter_dxdz(double delta_x, double delta_y, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz) {
24  return -track_dxdz * track_dydz * delta_phix + (1. + track_dxdz * track_dxdz) * delta_phiy - track_dydz * delta_phiz;
25 }
26 
27 double MuonResiduals6DOFFitter_dydz(double delta_x, double delta_y, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz) {
28  return -(1. + track_dydz * track_dydz) * delta_phix + track_dxdz * track_dydz * delta_phiy + track_dxdz * delta_phiz;
29 }
30 
31 Double_t MuonResiduals6DOFFitter_x_trackx_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFFitter_x(par[0], par[1], par[2], par[3], par[4], par[5], xvec[0], par[7], par[8], par[9], par[10], par[11]); }
32 Double_t MuonResiduals6DOFFitter_x_tracky_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFFitter_x(par[0], par[1], par[2], par[3], par[4], par[5], par[6], xvec[0], par[8], par[9], par[10], par[11]); }
33 Double_t MuonResiduals6DOFFitter_x_trackdxdz_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFFitter_x(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], xvec[0], par[9], par[10], par[11]); }
34 Double_t MuonResiduals6DOFFitter_x_trackdydz_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFFitter_x(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], xvec[0], par[10], par[11]); }
35 
36 Double_t MuonResiduals6DOFFitter_y_trackx_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFFitter_y(par[0], par[1], par[2], par[3], par[4], par[5], xvec[0], par[7], par[8], par[9], par[12], par[13]); }
37 Double_t MuonResiduals6DOFFitter_y_tracky_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFFitter_y(par[0], par[1], par[2], par[3], par[4], par[5], par[6], xvec[0], par[8], par[9], par[12], par[13]); }
38 Double_t MuonResiduals6DOFFitter_y_trackdxdz_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFFitter_y(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], xvec[0], par[9], par[12], par[13]); }
39 Double_t MuonResiduals6DOFFitter_y_trackdydz_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFFitter_y(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], xvec[0], par[12], par[13]); }
40 
41 Double_t MuonResiduals6DOFFitter_dxdz_trackx_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFFitter_dxdz(par[0], par[1], par[2], par[3], par[4], par[5], xvec[0], par[7], par[8], par[9]); }
42 Double_t MuonResiduals6DOFFitter_dxdz_tracky_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFFitter_dxdz(par[0], par[1], par[2], par[3], par[4], par[5], par[6], xvec[0], par[8], par[9]); }
43 Double_t MuonResiduals6DOFFitter_dxdz_trackdxdz_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFFitter_dxdz(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], xvec[0], par[9]); }
44 Double_t MuonResiduals6DOFFitter_dxdz_trackdydz_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFFitter_dxdz(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], xvec[0]); }
45 
46 Double_t MuonResiduals6DOFFitter_dydz_trackx_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFFitter_dydz(par[0], par[1], par[2], par[3], par[4], par[5], xvec[0], par[7], par[8], par[9]); }
47 Double_t MuonResiduals6DOFFitter_dydz_tracky_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFFitter_dydz(par[0], par[1], par[2], par[3], par[4], par[5], par[6], xvec[0], par[8], par[9]); }
48 Double_t MuonResiduals6DOFFitter_dydz_trackdxdz_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFFitter_dydz(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], xvec[0], par[9]); }
49 Double_t MuonResiduals6DOFFitter_dydz_trackdydz_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFFitter_dydz(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], xvec[0]); }
50 
51 void MuonResiduals6DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
52  MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo*)(MuonResiduals6DOFFitter_TMinuit->GetObjectFit());
53  MuonResidualsFitter *fitter = fitinfo->fitter();
54 
55  fval = 0.;
56  for (std::vector<double*>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end(); ++resiter) {
57  const double residX = (*resiter)[MuonResiduals6DOFFitter::kResidX];
58  const double residY = (*resiter)[MuonResiduals6DOFFitter::kResidY];
59  const double resslopeX = (*resiter)[MuonResiduals6DOFFitter::kResSlopeX];
60  const double resslopeY = (*resiter)[MuonResiduals6DOFFitter::kResSlopeY];
61  const double positionX = (*resiter)[MuonResiduals6DOFFitter::kPositionX];
62  const double positionY = (*resiter)[MuonResiduals6DOFFitter::kPositionY];
63  const double angleX = (*resiter)[MuonResiduals6DOFFitter::kAngleX];
64  const double angleY = (*resiter)[MuonResiduals6DOFFitter::kAngleY];
65  const double redchi2 = (*resiter)[MuonResiduals6DOFFitter::kRedChi2];
66 
67  const double alignx = par[MuonResiduals6DOFFitter::kAlignX];
68  const double aligny = par[MuonResiduals6DOFFitter::kAlignY];
69  const double alignz = par[MuonResiduals6DOFFitter::kAlignZ];
70  const double alignphix = par[MuonResiduals6DOFFitter::kAlignPhiX];
71  const double alignphiy = par[MuonResiduals6DOFFitter::kAlignPhiY];
72  const double alignphiz = par[MuonResiduals6DOFFitter::kAlignPhiZ];
73  const double resXsigma = par[MuonResiduals6DOFFitter::kResidXSigma];
74  const double resYsigma = par[MuonResiduals6DOFFitter::kResidYSigma];
75  const double slopeXsigma = par[MuonResiduals6DOFFitter::kResSlopeXSigma];
76  const double slopeYsigma = par[MuonResiduals6DOFFitter::kResSlopeYSigma];
77  const double alphax = par[MuonResiduals6DOFFitter::kAlphaX];
78  const double alphay = par[MuonResiduals6DOFFitter::kAlphaY];
79  const double resXgamma = par[MuonResiduals6DOFFitter::kResidXGamma];
80  const double resYgamma = par[MuonResiduals6DOFFitter::kResidYGamma];
81  const double slopeXgamma = par[MuonResiduals6DOFFitter::kResSlopeXGamma];
82  const double slopeYgamma = par[MuonResiduals6DOFFitter::kResSlopeYGamma];
83 
84  double residXpeak = MuonResiduals6DOFFitter_x(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, alphax, resslopeX);
85  double residYpeak = MuonResiduals6DOFFitter_y(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, alphay, resslopeY);
86  double slopeXpeak = MuonResiduals6DOFFitter_dxdz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
87  double slopeYpeak = MuonResiduals6DOFFitter_dydz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
88 
89  double weight = (1./redchi2) * MuonResiduals6DOFFitter_number_of_hits / MuonResiduals6DOFFitter_sum_of_weights;
90  if (!MuonResiduals6DOFFitter_weightAlignment) weight = 1.;
91 
92  if (!MuonResiduals6DOFFitter_weightAlignment || TMath::Prob(redchi2*12, 12) < 0.99) { // no spikes allowed
93 
94  if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
95  fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
96  fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
97  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
98  fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeY, slopeYpeak, slopeYsigma);
99  }
100  else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
101  fval += -weight * MuonResidualsFitter_logPowerLawTails(residX, residXpeak, resXsigma, resXgamma);
102  fval += -weight * MuonResidualsFitter_logPowerLawTails(residY, residYpeak, resYsigma, resYgamma);
103  fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
104  fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
105  }
106  else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
107  fval += -weight * MuonResidualsFitter_logROOTVoigt(residX, residXpeak, resXsigma, resXgamma);
108  fval += -weight * MuonResidualsFitter_logROOTVoigt(residY, residYpeak, resYsigma, resYgamma);
109  fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
110  fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
111  }
112  else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
113  fval += -weight * MuonResidualsFitter_logGaussPowerTails(residX, residXpeak, resXsigma);
114  fval += -weight * MuonResidualsFitter_logGaussPowerTails(residY, residYpeak, resYsigma);
115  fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeX, slopeXpeak, slopeXsigma);
116  fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeY, slopeYpeak, slopeYsigma);
117  }
118  else { assert(false); }
119 
120  }
121  }
122 }
123 
124 double MuonResiduals6DOFFitter::sumofweights() {
125  MuonResiduals6DOFFitter_sum_of_weights = 0.;
126  MuonResiduals6DOFFitter_number_of_hits = 0.;
127  MuonResiduals6DOFFitter_weightAlignment = m_weightAlignment;
128  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
129  if (m_weightAlignment) {
130  double redchi2 = (*resiter)[MuonResiduals6DOFFitter::kRedChi2];
131  if (TMath::Prob(redchi2*12, 12) < 0.99) { // no spikes allowed
132  MuonResiduals6DOFFitter_sum_of_weights += 1./redchi2;
133  MuonResiduals6DOFFitter_number_of_hits += 1.;
134  }
135  }
136  else {
137  MuonResiduals6DOFFitter_sum_of_weights += 1.;
138  MuonResiduals6DOFFitter_number_of_hits += 1.;
139  }
140  }
141  return MuonResiduals6DOFFitter_sum_of_weights;
142 }
143 
144 bool MuonResiduals6DOFFitter::fit(Alignable *ali) {
145  initialize_table(); // if not already initialized
146  sumofweights();
147 
148  double residx_mean = 0;
149  double residy_mean = 0;
150  double resslopex_mean = 0;
151  double resslopey_mean = 0;
152  double residx_stdev = 0.5;
153  double residy_stdev = 3.0;
154  double resslopex_stdev = 0.005;
155  double resslopey_stdev = 0.03;
156  double alphax_estimate = 0;
157  double alphay_estimate = 0;
158 
159  std::vector<int> num;
160  std::vector<std::string> name;
161  std::vector<double> start;
162  std::vector<double> step;
163  std::vector<double> low;
164  std::vector<double> high;
165 
166  if (fixed(kAlignX)) {
167  num.push_back(kAlignX); name.push_back(std::string("AlignX")); start.push_back(0.); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
168  } else {
169  num.push_back(kAlignX); name.push_back(std::string("AlignX")); start.push_back(residx_mean); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
170  }
171  if (fixed(kAlignY)) {
172  num.push_back(kAlignY); name.push_back(std::string("AlignY")); start.push_back(0.); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
173  } else {
174  num.push_back(kAlignY); name.push_back(std::string("AlignY")); start.push_back(residy_mean); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
175  }
176  num.push_back(kAlignZ); name.push_back(std::string("AlignZ")); start.push_back(0.); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
177  if (fixed(kAlignPhiX)) {
178  num.push_back(kAlignPhiX); name.push_back(std::string("AlignPhiX")); start.push_back(0.); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
179  } else {
180  num.push_back(kAlignPhiX); name.push_back(std::string("AlignPhiX")); start.push_back(-resslopey_mean); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
181  }
182  if (fixed(kAlignPhiY)) {
183  num.push_back(kAlignPhiY); name.push_back(std::string("AlignPhiY")); start.push_back(0.); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
184  } else {
185  num.push_back(kAlignPhiY); name.push_back(std::string("AlignPhiY")); start.push_back(resslopex_mean); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
186  }
187  num.push_back(kAlignPhiZ); name.push_back(std::string("AlignPhiZ")); start.push_back(0.); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
188  num.push_back(kResidXSigma); name.push_back(std::string("ResidXSigma")); start.push_back(residx_stdev); step.push_back(0.01*residx_stdev); low.push_back(0.); high.push_back(0.);
189  num.push_back(kResidYSigma); name.push_back(std::string("ResidYSigma")); start.push_back(residy_stdev); step.push_back(0.01*residy_stdev); low.push_back(0.); high.push_back(0.);
190  num.push_back(kResSlopeXSigma); name.push_back(std::string("ResSlopeXSigma")); start.push_back(resslopex_stdev); step.push_back(0.01*resslopex_stdev); low.push_back(0.); high.push_back(0.);
191  num.push_back(kResSlopeYSigma); name.push_back(std::string("ResSlopeYSigma")); start.push_back(resslopey_stdev); step.push_back(0.01*resslopey_stdev); low.push_back(0.); high.push_back(0.);
192  num.push_back(kAlphaX); name.push_back(std::string("AlphaX")); start.push_back(alphax_estimate); step.push_back(0.01*resslopex_stdev); low.push_back(0.); high.push_back(0.);
193  num.push_back(kAlphaY); name.push_back(std::string("AlphaY")); start.push_back(alphay_estimate); step.push_back(0.01*resslopey_stdev); low.push_back(0.); high.push_back(0.);
194  if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
195  num.push_back(kResidXGamma); name.push_back(std::string("ResidXGamma")); start.push_back(0.1*residx_stdev); step.push_back(0.01*residx_stdev); low.push_back(0.); high.push_back(0.);
196  num.push_back(kResidYGamma); name.push_back(std::string("ResidYGamma")); start.push_back(0.1*residy_stdev); step.push_back(0.01*residy_stdev); low.push_back(0.); high.push_back(0.);
197  num.push_back(kResSlopeXGamma); name.push_back(std::string("ResSlopeXGamma")); start.push_back(0.1*resslopex_stdev); step.push_back(0.01*resslopex_stdev); low.push_back(0.); high.push_back(0.);
198  num.push_back(kResSlopeYGamma); name.push_back(std::string("ResSlopeYGamma")); start.push_back(0.1*resslopey_stdev); step.push_back(0.01*resslopey_stdev); low.push_back(0.); high.push_back(0.);
199  }
200 
201  return dofit(&MuonResiduals6DOFFitter_FCN, num, name, start, step, low, high);
202 }
203 
204 double MuonResiduals6DOFFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
205  sumofweights();
206 
207 // std::stringstream name_ntuple;
208 // name_ntuple << name << "_ntuple";
209 // TTree *ntuple = dir->make<TTree>(name_ntuple.str().c_str(), "");
210 // Float_t ntuple_residX;
211 // Float_t ntuple_residY;
212 // Float_t ntuple_resslopeX;
213 // Float_t ntuple_resslopeY;
214 // Float_t ntuple_positionX;
215 // Float_t ntuple_positionY;
216 // Float_t ntuple_angleX;
217 // Float_t ntuple_angleY;
218 // Float_t ntuple_redchi2;
219 // Float_t ntuple_prob;
220 
221 // ntuple->Branch("residX", &ntuple_residX, "residX/F");
222 // ntuple->Branch("residY", &ntuple_residY, "residY/F");
223 // ntuple->Branch("resslopeX", &ntuple_resslopeX, "resslopeX/F");
224 // ntuple->Branch("resslopeY", &ntuple_resslopeY, "resslopeY/F");
225 // ntuple->Branch("positionX", &ntuple_positionX, "positionX/F");
226 // ntuple->Branch("positionY", &ntuple_positionY, "positionY/F");
227 // ntuple->Branch("angleX", &ntuple_angleX, "angleX/F");
228 // ntuple->Branch("angleY", &ntuple_angleY, "angleY/F");
229 // ntuple->Branch("redchi2", &ntuple_redchi2, "redchi2/F");
230 // ntuple->Branch("prob", &ntuple_prob, "prob/F");
231 
232 // for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
233 // ntuple_residX = (*resiter)[MuonResiduals6DOFFitter::kResidX];
234 // ntuple_residY = (*resiter)[MuonResiduals6DOFFitter::kResidY];
235 // ntuple_resslopeX = (*resiter)[MuonResiduals6DOFFitter::kResSlopeX];
236 // ntuple_resslopeY = (*resiter)[MuonResiduals6DOFFitter::kResSlopeY];
237 // ntuple_positionX = (*resiter)[MuonResiduals6DOFFitter::kPositionX];
238 // ntuple_positionY = (*resiter)[MuonResiduals6DOFFitter::kPositionY];
239 // ntuple_angleX = (*resiter)[MuonResiduals6DOFFitter::kAngleX];
240 // ntuple_angleY = (*resiter)[MuonResiduals6DOFFitter::kAngleY];
241 // ntuple_redchi2 = (*resiter)[MuonResiduals6DOFFitter::kRedChi2];
242 // ntuple_prob = TMath::Prob((*resiter)[MuonResiduals6DOFFitter::kRedChi2]);
243 
244 // ntuple->Fill();
245 // }
246 
247  std::stringstream 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;
248  std::stringstream name_x_trackx, name_y_trackx, name_dxdz_trackx, name_dydz_trackx;
249  std::stringstream name_x_tracky, name_y_tracky, name_dxdz_tracky, name_dydz_tracky;
250  std::stringstream name_x_trackdxdz, name_y_trackdxdz, name_dxdz_trackdxdz, name_dydz_trackdxdz;
251  std::stringstream name_x_trackdydz, name_y_trackdydz, name_dxdz_trackdydz, name_dydz_trackdydz;
252 
253  name_x << name << "_x";
254  name_y << name << "_y";
255  name_dxdz << name << "_dxdz";
256  name_dydz << name << "_dydz";
257  name_x_raw << name << "_x_raw";
258  name_y_raw << name << "_y_raw";
259  name_dxdz_raw << name << "_dxdz_raw";
260  name_dydz_raw << name << "_dydz_raw";
261  name_x_cut << name << "_x_cut";
262  name_y_cut << name << "_y_cut";
263  name_alphax << name << "_alphax";
264  name_alphay << name << "_alphay";
265  name_x_trackx << name << "_x_trackx";
266  name_y_trackx << name << "_y_trackx";
267  name_dxdz_trackx << name << "_dxdz_trackx";
268  name_dydz_trackx << name << "_dydz_trackx";
269  name_x_tracky << name << "_x_tracky";
270  name_y_tracky << name << "_y_tracky";
271  name_dxdz_tracky << name << "_dxdz_tracky";
272  name_dydz_tracky << name << "_dydz_tracky";
273  name_x_trackdxdz << name << "_x_trackdxdz";
274  name_y_trackdxdz << name << "_y_trackdxdz";
275  name_dxdz_trackdxdz << name << "_dxdz_trackdxdz";
276  name_dydz_trackdxdz << name << "_dydz_trackdxdz";
277  name_x_trackdydz << name << "_x_trackdydz";
278  name_y_trackdydz << name << "_y_trackdydz";
279  name_dxdz_trackdydz << name << "_dxdz_trackdydz";
280  name_dydz_trackdydz << name << "_dydz_trackdydz";
281 
282  double width = ali->surface().width();
283  double length = ali->surface().length();
284  double min_x = -100.; double max_x = 100.;
285  double min_y = -200.; double max_y = 200.;
286  double min_dxdz = -100.; double max_dxdz = 100.;
287  double min_dydz = -200.; double max_dydz = 200.;
288  double min_trackx = -width/2.; double max_trackx = width/2.;
289  double min_tracky = -length/2.; double max_tracky = length/2.;
290  double min_trackdxdz = -1.5; double max_trackdxdz = 1.5;
291  double min_trackdydz = -1.5; double max_trackdydz = 1.5;
292 
293  TH1F *hist_x = dir->make<TH1F>(name_x.str().c_str(), "", 100, min_x, max_x);
294  TH1F *hist_y = dir->make<TH1F>(name_y.str().c_str(), "", 100, min_y, max_y);
295  TH1F *hist_dxdz = dir->make<TH1F>(name_dxdz.str().c_str(), "", 100, min_dxdz, max_dxdz);
296  TH1F *hist_dydz = dir->make<TH1F>(name_dydz.str().c_str(), "", 100, min_dydz, max_dydz);
297  TH1F *hist_x_raw = dir->make<TH1F>(name_x_raw.str().c_str(), "", 100, min_x, max_x);
298  TH1F *hist_y_raw = dir->make<TH1F>(name_y_raw.str().c_str(), "", 100, min_y, max_y);
299  TH1F *hist_dxdz_raw = dir->make<TH1F>(name_dxdz_raw.str().c_str(), "", 100, min_dxdz, max_dxdz);
300  TH1F *hist_dydz_raw = dir->make<TH1F>(name_dydz_raw.str().c_str(), "", 100, min_dydz, max_dydz);
301  TH1F *hist_x_cut = dir->make<TH1F>(name_x_cut.str().c_str(), "", 100, min_x, max_x);
302  TH1F *hist_y_cut = dir->make<TH1F>(name_y_cut.str().c_str(), "", 100, min_y, max_y);
303  TH2F *hist_alphax = dir->make<TH2F>(name_alphax.str().c_str(), "", 40, min_dxdz, max_dxdz, 40, -20., 20.);
304  TH2F *hist_alphay = dir->make<TH2F>(name_alphay.str().c_str(), "", 40, min_dydz, max_dydz, 40, -100., 100.);
305  TProfile *hist_x_trackx = dir->make<TProfile>(name_x_trackx.str().c_str(), "", 100, min_trackx, max_trackx, min_x, max_x);
306  TProfile *hist_y_trackx = dir->make<TProfile>(name_y_trackx.str().c_str(), "", 100, min_trackx, max_trackx, min_y, max_y);
307  TProfile *hist_dxdz_trackx = dir->make<TProfile>(name_dxdz_trackx.str().c_str(), "", 100, min_trackx, max_trackx, min_dxdz, max_dxdz);
308  TProfile *hist_dydz_trackx = dir->make<TProfile>(name_dydz_trackx.str().c_str(), "", 100, min_trackx, max_trackx, min_dydz, max_dydz);
309  TProfile *hist_x_tracky = dir->make<TProfile>(name_x_tracky.str().c_str(), "", 100, min_tracky, max_tracky, min_x, max_x);
310  TProfile *hist_y_tracky = dir->make<TProfile>(name_y_tracky.str().c_str(), "", 100, min_tracky, max_tracky, min_y, max_y);
311  TProfile *hist_dxdz_tracky = dir->make<TProfile>(name_dxdz_tracky.str().c_str(), "", 100, min_tracky, max_tracky, min_dxdz, max_dxdz);
312  TProfile *hist_dydz_tracky = dir->make<TProfile>(name_dydz_tracky.str().c_str(), "", 100, min_tracky, max_tracky, min_dydz, max_dydz);
313  TProfile *hist_x_trackdxdz = dir->make<TProfile>(name_x_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_x, max_x);
314  TProfile *hist_y_trackdxdz = dir->make<TProfile>(name_y_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_y, max_y);
315  TProfile *hist_dxdz_trackdxdz = dir->make<TProfile>(name_dxdz_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_dxdz, max_dxdz);
316  TProfile *hist_dydz_trackdxdz = dir->make<TProfile>(name_dydz_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_dydz, max_dydz);
317  TProfile *hist_x_trackdydz = dir->make<TProfile>(name_x_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz, min_x, max_x);
318  TProfile *hist_y_trackdydz = dir->make<TProfile>(name_y_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz, min_y, max_y);
319  TProfile *hist_dxdz_trackdydz = dir->make<TProfile>(name_dxdz_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz, min_dxdz, max_dxdz);
320  TProfile *hist_dydz_trackdydz = dir->make<TProfile>(name_dydz_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz, min_dydz, max_dydz);
321 
322  hist_x_trackx->SetAxisRange(-10., 10., "Y");
323  hist_y_trackx->SetAxisRange(-20., 20., "Y");
324  hist_dxdz_trackx->SetAxisRange(-10., 10., "Y");
325  hist_dydz_trackx->SetAxisRange(-20., 20., "Y");
326  hist_x_tracky->SetAxisRange(-10., 10., "Y");
327  hist_y_tracky->SetAxisRange(-20., 20., "Y");
328  hist_dxdz_tracky->SetAxisRange(-10., 10., "Y");
329  hist_dydz_tracky->SetAxisRange(-20., 20., "Y");
330  hist_x_trackdxdz->SetAxisRange(-10., 10., "Y");
331  hist_y_trackdxdz->SetAxisRange(-20., 20., "Y");
332  hist_dxdz_trackdxdz->SetAxisRange(-10., 10., "Y");
333  hist_dydz_trackdxdz->SetAxisRange(-20., 20., "Y");
334  hist_x_trackdydz->SetAxisRange(-10., 10., "Y");
335  hist_y_trackdydz->SetAxisRange(-20., 20., "Y");
336  hist_dxdz_trackdydz->SetAxisRange(-10., 10., "Y");
337  hist_dydz_trackdydz->SetAxisRange(-20., 20., "Y");
338 
339  name_x << "_fit";
340  name_y << "_fit";
341  name_dxdz << "_fit";
342  name_dydz << "_fit";
343  name_alphax << "_fit";
344  name_alphay << "_fit";
345  name_x_trackx << "_fit";
346  name_y_trackx << "_fit";
347  name_dxdz_trackx << "_fit";
348  name_dydz_trackx << "_fit";
349  name_x_tracky << "_fit";
350  name_y_tracky << "_fit";
351  name_dxdz_tracky << "_fit";
352  name_dydz_tracky << "_fit";
353  name_x_trackdxdz << "_fit";
354  name_y_trackdxdz << "_fit";
355  name_dxdz_trackdxdz << "_fit";
356  name_dydz_trackdxdz << "_fit";
357  name_x_trackdydz << "_fit";
358  name_y_trackdydz << "_fit";
359  name_dxdz_trackdydz << "_fit";
360  name_dydz_trackdydz << "_fit";
361 
362  TF1 *fit_x = NULL;
363  TF1 *fit_y = NULL;
364  TF1 *fit_dxdz = NULL;
365  TF1 *fit_dydz = NULL;
366  if (residualsModel() == kPureGaussian) {
367  fit_x = new TF1(name_x.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, min_x, max_x, 3);
368  fit_x->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_x - min_x)/100., 10.*value(kAlignX), 10.*value(kResidXSigma));
369  fit_y = new TF1(name_y.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, min_y, max_y, 3);
370  fit_y->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_y - min_y)/100., 10.*value(kAlignY), 10.*value(kResidYSigma));
371  fit_dxdz = new TF1(name_dxdz.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dxdz, max_dxdz, 3);
372  fit_dxdz->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_dxdz - min_dxdz)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma));
373  fit_dydz = new TF1(name_dydz.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dydz, max_dydz, 3);
374  fit_dydz->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_dydz - min_dydz)/100., -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma));
375  }
376  else if (residualsModel() == kPowerLawTails) {
377  fit_x = new TF1(name_x.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, min_x, max_x, 4);
378  fit_x->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_x - min_x)/100., 10.*value(kAlignX), 10.*value(kResidXSigma), 10.*value(kResidXGamma));
379  fit_y = new TF1(name_y.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, min_y, max_y, 4);
380  fit_y->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_y - min_y)/100., 10.*value(kAlignY), 10.*value(kResidYSigma), 10.*value(kResidYGamma));
381  fit_dxdz = new TF1(name_dxdz.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dxdz, max_dxdz, 4);
382  fit_dxdz->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_dxdz - min_dxdz)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma), 1000.*value(kResSlopeXGamma));
383  fit_dydz = new TF1(name_dydz.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dydz, max_dydz, 4);
384  fit_dydz->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_dydz - min_dydz)/100., -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma), 1000.*value(kResSlopeYGamma));
385  }
386  else if (residualsModel() == kROOTVoigt) {
387  fit_x = new TF1(name_x.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_x, max_x, 4);
388  fit_x->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_x - min_x)/100., 10.*value(kAlignX), 10.*value(kResidXSigma), 10.*value(kResidXGamma));
389  fit_y = new TF1(name_y.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_y, max_y, 4);
390  fit_y->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_y - min_y)/100., 10.*value(kAlignY), 10.*value(kResidYSigma), 10.*value(kResidYGamma));
391  fit_dxdz = new TF1(name_dxdz.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dxdz, max_dxdz, 4);
392  fit_dxdz->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_dxdz - min_dxdz)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma), 1000.*value(kResSlopeXGamma));
393  fit_dydz = new TF1(name_dydz.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dydz, max_dydz, 4);
394  fit_dydz->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_dydz - min_dydz)/100., -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma), 1000.*value(kResSlopeYGamma));
395  }
396  else if (residualsModel() == kGaussPowerTails) {
397  fit_x = new TF1(name_x.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_x, max_x, 3);
398  fit_x->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_x - min_x)/100., 10.*value(kAlignX), 10.*value(kResidXSigma));
399  fit_y = new TF1(name_y.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_y, max_y, 3);
400  fit_y->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_y - min_y)/100., 10.*value(kAlignY), 10.*value(kResidYSigma));
401  fit_dxdz = new TF1(name_dxdz.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dxdz, max_dxdz, 3);
402  fit_dxdz->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_dxdz - min_dxdz)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeXSigma));
403  fit_dydz = new TF1(name_dydz.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dydz, max_dydz, 3);
404  fit_dydz->SetParameters(MuonResiduals6DOFFitter_sum_of_weights * (max_dydz - min_dydz)/100., -1000.*value(kAlignPhiX), 1000.*value(kResSlopeYSigma));
405  }
406  else { assert(false); }
407 
408  fit_x->SetLineColor(2); fit_x->SetLineWidth(2);
409  fit_y->SetLineColor(2); fit_y->SetLineWidth(2);
410  fit_dxdz->SetLineColor(2); fit_dxdz->SetLineWidth(2);
411  fit_dydz->SetLineColor(2); fit_dydz->SetLineWidth(2);
412  fit_x->Write();
413  fit_y->Write();
414  fit_dxdz->Write();
415  fit_dydz->Write();
416 
417  TF1 *fit_alphax = new TF1(name_alphax.str().c_str(), "[0] + x*[1]", min_dxdz, max_dxdz);
418  fit_alphax->SetParameters(10.*value(kAlignX), 10.*value(kAlphaX)/1000.);
419  TF1 *fit_alphay = new TF1(name_alphay.str().c_str(), "[0] + x*[1]", min_dydz, max_dydz);
420  fit_alphay->SetParameters(10.*value(kAlignY), 10.*value(kAlphaY)/1000.);
421 
422  fit_alphax->SetLineColor(2); fit_alphax->SetLineWidth(2);
423  fit_alphay->SetLineColor(2); fit_alphay->SetLineWidth(2);
424  fit_alphax->Write();
425  fit_alphay->Write();
426 
427  TProfile *fit_x_trackx = dir->make<TProfile>(name_x_trackx.str().c_str(), "", 100, min_trackx, max_trackx);
428  TProfile *fit_y_trackx = dir->make<TProfile>(name_y_trackx.str().c_str(), "", 100, min_trackx, max_trackx);
429  TProfile *fit_dxdz_trackx = dir->make<TProfile>(name_dxdz_trackx.str().c_str(), "", 100, min_trackx, max_trackx);
430  TProfile *fit_dydz_trackx = dir->make<TProfile>(name_dydz_trackx.str().c_str(), "", 100, min_trackx, max_trackx);
431  TProfile *fit_x_tracky = dir->make<TProfile>(name_x_tracky.str().c_str(), "", 100, min_tracky, max_tracky);
432  TProfile *fit_y_tracky = dir->make<TProfile>(name_y_tracky.str().c_str(), "", 100, min_tracky, max_tracky);
433  TProfile *fit_dxdz_tracky = dir->make<TProfile>(name_dxdz_tracky.str().c_str(), "", 100, min_tracky, max_tracky);
434  TProfile *fit_dydz_tracky = dir->make<TProfile>(name_dydz_tracky.str().c_str(), "", 100, min_tracky, max_tracky);
435  TProfile *fit_x_trackdxdz = dir->make<TProfile>(name_x_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz);
436  TProfile *fit_y_trackdxdz = dir->make<TProfile>(name_y_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz);
437  TProfile *fit_dxdz_trackdxdz = dir->make<TProfile>(name_dxdz_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz);
438  TProfile *fit_dydz_trackdxdz = dir->make<TProfile>(name_dydz_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz);
439  TProfile *fit_x_trackdydz = dir->make<TProfile>(name_x_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz);
440  TProfile *fit_y_trackdydz = dir->make<TProfile>(name_y_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz);
441  TProfile *fit_dxdz_trackdydz = dir->make<TProfile>(name_dxdz_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz);
442  TProfile *fit_dydz_trackdydz = dir->make<TProfile>(name_dydz_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz);
443 
444  fit_x_trackx->SetLineColor(2); fit_x_trackx->SetLineWidth(2);
445  fit_y_trackx->SetLineColor(2); fit_y_trackx->SetLineWidth(2);
446  fit_dxdz_trackx->SetLineColor(2); fit_dxdz_trackx->SetLineWidth(2);
447  fit_dydz_trackx->SetLineColor(2); fit_dydz_trackx->SetLineWidth(2);
448  fit_x_tracky->SetLineColor(2); fit_x_tracky->SetLineWidth(2);
449  fit_y_tracky->SetLineColor(2); fit_y_tracky->SetLineWidth(2);
450  fit_dxdz_tracky->SetLineColor(2); fit_dxdz_tracky->SetLineWidth(2);
451  fit_dydz_tracky->SetLineColor(2); fit_dydz_tracky->SetLineWidth(2);
452  fit_x_trackdxdz->SetLineColor(2); fit_x_trackdxdz->SetLineWidth(2);
453  fit_y_trackdxdz->SetLineColor(2); fit_y_trackdxdz->SetLineWidth(2);
454  fit_dxdz_trackdxdz->SetLineColor(2); fit_dxdz_trackdxdz->SetLineWidth(2);
455  fit_dydz_trackdxdz->SetLineColor(2); fit_dydz_trackdxdz->SetLineWidth(2);
456  fit_x_trackdydz->SetLineColor(2); fit_x_trackdydz->SetLineWidth(2);
457  fit_y_trackdydz->SetLineColor(2); fit_y_trackdydz->SetLineWidth(2);
458  fit_dxdz_trackdydz->SetLineColor(2); fit_dxdz_trackdydz->SetLineWidth(2);
459  fit_dydz_trackdydz->SetLineColor(2); fit_dydz_trackdydz->SetLineWidth(2);
460 
461  name_x_trackx << "line";
462  name_y_trackx << "line";
463  name_dxdz_trackx << "line";
464  name_dydz_trackx << "line";
465  name_x_tracky << "line";
466  name_y_tracky << "line";
467  name_dxdz_tracky << "line";
468  name_dydz_tracky << "line";
469  name_x_trackdxdz << "line";
470  name_y_trackdxdz << "line";
471  name_dxdz_trackdxdz << "line";
472  name_dydz_trackdxdz << "line";
473  name_x_trackdydz << "line";
474  name_y_trackdydz << "line";
475  name_dxdz_trackdydz << "line";
476  name_dydz_trackdydz << "line";
477 
478  TF1 *fitline_x_trackx = new TF1(name_x_trackx.str().c_str(), MuonResiduals6DOFFitter_x_trackx_TF1, min_trackx, max_trackx, 14);
479  TF1 *fitline_y_trackx = new TF1(name_y_trackx.str().c_str(), MuonResiduals6DOFFitter_y_trackx_TF1, min_trackx, max_trackx, 14);
480  TF1 *fitline_dxdz_trackx = new TF1(name_dxdz_trackx.str().c_str(), MuonResiduals6DOFFitter_dxdz_trackx_TF1, min_trackx, max_trackx, 14);
481  TF1 *fitline_dydz_trackx = new TF1(name_dydz_trackx.str().c_str(), MuonResiduals6DOFFitter_dydz_trackx_TF1, min_trackx, max_trackx, 14);
482  TF1 *fitline_x_tracky = new TF1(name_x_tracky.str().c_str(), MuonResiduals6DOFFitter_x_tracky_TF1, min_tracky, max_tracky, 14);
483  TF1 *fitline_y_tracky = new TF1(name_y_tracky.str().c_str(), MuonResiduals6DOFFitter_y_tracky_TF1, min_tracky, max_tracky, 14);
484  TF1 *fitline_dxdz_tracky = new TF1(name_dxdz_tracky.str().c_str(), MuonResiduals6DOFFitter_dxdz_tracky_TF1, min_tracky, max_tracky, 14);
485  TF1 *fitline_dydz_tracky = new TF1(name_dydz_tracky.str().c_str(), MuonResiduals6DOFFitter_dydz_tracky_TF1, min_tracky, max_tracky, 14);
486  TF1 *fitline_x_trackdxdz = new TF1(name_x_trackdxdz.str().c_str(), MuonResiduals6DOFFitter_x_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
487  TF1 *fitline_y_trackdxdz = new TF1(name_y_trackdxdz.str().c_str(), MuonResiduals6DOFFitter_y_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
488  TF1 *fitline_dxdz_trackdxdz = new TF1(name_dxdz_trackdxdz.str().c_str(), MuonResiduals6DOFFitter_dxdz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
489  TF1 *fitline_dydz_trackdxdz = new TF1(name_dydz_trackdxdz.str().c_str(), MuonResiduals6DOFFitter_dydz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
490  TF1 *fitline_x_trackdydz = new TF1(name_x_trackdydz.str().c_str(), MuonResiduals6DOFFitter_x_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
491  TF1 *fitline_y_trackdydz = new TF1(name_y_trackdydz.str().c_str(), MuonResiduals6DOFFitter_y_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
492  TF1 *fitline_dxdz_trackdydz = new TF1(name_dxdz_trackdydz.str().c_str(), MuonResiduals6DOFFitter_dxdz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
493  TF1 *fitline_dydz_trackdydz = new TF1(name_dydz_trackdydz.str().c_str(), MuonResiduals6DOFFitter_dydz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
494 
495  double sum_resslopex = 0.;
496  double sum_resslopey = 0.;
497  double sum_trackx = 0.;
498  double sum_tracky = 0.;
499  double sum_trackdxdz = 0.;
500  double sum_trackdydz = 0.;
501  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
502  const double resslopeX = (*resiter)[MuonResiduals6DOFFitter::kResSlopeX];
503  const double resslopeY = (*resiter)[MuonResiduals6DOFFitter::kResSlopeY];
504  const double positionX = (*resiter)[MuonResiduals6DOFFitter::kPositionX];
505  const double positionY = (*resiter)[MuonResiduals6DOFFitter::kPositionY];
506  const double angleX = (*resiter)[MuonResiduals6DOFFitter::kAngleX];
507  const double angleY = (*resiter)[MuonResiduals6DOFFitter::kAngleY];
508  const double redchi2 = (*resiter)[MuonResiduals6DOFFitter::kRedChi2];
509  double weight = 1./redchi2;
510  if (!m_weightAlignment) weight = 1.;
511 
512  if (!m_weightAlignment || TMath::Prob(redchi2*12, 12) < 0.99) { // no spikes allowed
513  sum_resslopex += weight * resslopeX;
514  sum_resslopey += weight * resslopeY;
515  sum_trackx += weight * positionX;
516  sum_tracky += weight * positionY;
517  sum_trackdxdz += weight * angleX;
518  sum_trackdydz += weight * angleY;
519  }
520  }
521  double mean_resslopex = sum_resslopex / MuonResiduals6DOFFitter_sum_of_weights;
522  double mean_resslopey = sum_resslopey / MuonResiduals6DOFFitter_sum_of_weights;
523  double mean_trackx = sum_trackx / MuonResiduals6DOFFitter_sum_of_weights;
524  double mean_tracky = sum_tracky / MuonResiduals6DOFFitter_sum_of_weights;
525  double mean_trackdxdz = sum_trackdxdz / MuonResiduals6DOFFitter_sum_of_weights;
526  double mean_trackdydz = sum_trackdydz / MuonResiduals6DOFFitter_sum_of_weights;
527 
528  double fitparameters[14];
529  fitparameters[0] = value(kAlignX);
530  fitparameters[1] = value(kAlignY);
531  fitparameters[2] = value(kAlignZ);
532  fitparameters[3] = value(kAlignPhiX);
533  fitparameters[4] = value(kAlignPhiY);
534  fitparameters[5] = value(kAlignPhiZ);
535  fitparameters[6] = mean_trackx;
536  fitparameters[7] = mean_tracky;
537  fitparameters[8] = mean_trackdxdz;
538  fitparameters[9] = mean_trackdydz;
539  fitparameters[10] = value(kAlphaX);
540  fitparameters[11] = mean_resslopex;
541  fitparameters[12] = value(kAlphaY);
542  fitparameters[13] = mean_resslopey;
543 
544  fitline_x_trackx->SetParameters(fitparameters);
545  fitline_y_trackx->SetParameters(fitparameters);
546  fitline_dxdz_trackx->SetParameters(fitparameters);
547  fitline_dydz_trackx->SetParameters(fitparameters);
548  fitline_x_tracky->SetParameters(fitparameters);
549  fitline_y_tracky->SetParameters(fitparameters);
550  fitline_dxdz_tracky->SetParameters(fitparameters);
551  fitline_dydz_tracky->SetParameters(fitparameters);
552  fitline_x_trackdxdz->SetParameters(fitparameters);
553  fitline_y_trackdxdz->SetParameters(fitparameters);
554  fitline_dxdz_trackdxdz->SetParameters(fitparameters);
555  fitline_dydz_trackdxdz->SetParameters(fitparameters);
556  fitline_x_trackdydz->SetParameters(fitparameters);
557  fitline_y_trackdydz->SetParameters(fitparameters);
558  fitline_dxdz_trackdydz->SetParameters(fitparameters);
559  fitline_dydz_trackdydz->SetParameters(fitparameters);
560 
561  fitline_x_trackx->SetLineColor(2); fitline_x_trackx->SetLineWidth(2);
562  fitline_y_trackx->SetLineColor(2); fitline_y_trackx->SetLineWidth(2);
563  fitline_dxdz_trackx->SetLineColor(2); fitline_dxdz_trackx->SetLineWidth(2);
564  fitline_dydz_trackx->SetLineColor(2); fitline_dydz_trackx->SetLineWidth(2);
565  fitline_x_tracky->SetLineColor(2); fitline_x_tracky->SetLineWidth(2);
566  fitline_y_tracky->SetLineColor(2); fitline_y_tracky->SetLineWidth(2);
567  fitline_dxdz_tracky->SetLineColor(2); fitline_dxdz_tracky->SetLineWidth(2);
568  fitline_dydz_tracky->SetLineColor(2); fitline_dydz_tracky->SetLineWidth(2);
569  fitline_x_trackdxdz->SetLineColor(2); fitline_x_trackdxdz->SetLineWidth(2);
570  fitline_y_trackdxdz->SetLineColor(2); fitline_y_trackdxdz->SetLineWidth(2);
571  fitline_dxdz_trackdxdz->SetLineColor(2); fitline_dxdz_trackdxdz->SetLineWidth(2);
572  fitline_dydz_trackdxdz->SetLineColor(2); fitline_dydz_trackdxdz->SetLineWidth(2);
573  fitline_x_trackdydz->SetLineColor(2); fitline_x_trackdydz->SetLineWidth(2);
574  fitline_y_trackdydz->SetLineColor(2); fitline_y_trackdydz->SetLineWidth(2);
575  fitline_dxdz_trackdydz->SetLineColor(2); fitline_dxdz_trackdydz->SetLineWidth(2);
576  fitline_dydz_trackdydz->SetLineColor(2); fitline_dydz_trackdydz->SetLineWidth(2);
577 
578  fitline_x_trackx->Write();
579  fitline_y_trackx->Write();
580  fitline_dxdz_trackx->Write();
581  fitline_dydz_trackx->Write();
582  fitline_x_tracky->Write();
583  fitline_y_tracky->Write();
584  fitline_dxdz_tracky->Write();
585  fitline_dydz_tracky->Write();
586  fitline_x_trackdxdz->Write();
587  fitline_y_trackdxdz->Write();
588  fitline_dxdz_trackdxdz->Write();
589  fitline_dydz_trackdxdz->Write();
590  fitline_x_trackdydz->Write();
591  fitline_y_trackdydz->Write();
592  fitline_dxdz_trackdydz->Write();
593  fitline_dydz_trackdydz->Write();
594 
595  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
596  const double residX = (*resiter)[MuonResiduals6DOFFitter::kResidX];
597  const double residY = (*resiter)[MuonResiduals6DOFFitter::kResidY];
598  const double resslopeX = (*resiter)[MuonResiduals6DOFFitter::kResSlopeX];
599  const double resslopeY = (*resiter)[MuonResiduals6DOFFitter::kResSlopeY];
600  const double positionX = (*resiter)[MuonResiduals6DOFFitter::kPositionX];
601  const double positionY = (*resiter)[MuonResiduals6DOFFitter::kPositionY];
602  const double angleX = (*resiter)[MuonResiduals6DOFFitter::kAngleX];
603  const double angleY = (*resiter)[MuonResiduals6DOFFitter::kAngleY];
604  const double redchi2 = (*resiter)[MuonResiduals6DOFFitter::kRedChi2];
605  double weight = 1./redchi2;
606  if (!m_weightAlignment) weight = 1.;
607 
608  if (!m_weightAlignment || TMath::Prob(redchi2*12, 12) < 0.99) { // no spikes allowed
609 
610  hist_alphax->Fill(1000.*resslopeX, 10.*residX);
611  hist_alphay->Fill(1000.*resslopeY, 10.*residY);
612 
613  double geom_residX = MuonResiduals6DOFFitter_x(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY, value(kAlphaX), resslopeX);
614  hist_x->Fill(10.*(residX - geom_residX + value(kAlignX)), weight);
615  hist_x_trackx->Fill(positionX, 10.*residX, weight);
616  hist_x_tracky->Fill(positionY, 10.*residX, weight);
617  hist_x_trackdxdz->Fill(angleX, 10.*residX, weight);
618  hist_x_trackdydz->Fill(angleY, 10.*residX, weight);
619  fit_x_trackx->Fill(positionX, 10.*geom_residX, weight);
620  fit_x_tracky->Fill(positionY, 10.*geom_residX, weight);
621  fit_x_trackdxdz->Fill(angleX, 10.*geom_residX, weight);
622  fit_x_trackdydz->Fill(angleY, 10.*geom_residX, weight);
623 
624  double geom_residY = MuonResiduals6DOFFitter_y(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY, value(kAlphaY), resslopeY);
625  hist_y->Fill(10.*(residY - geom_residY + value(kAlignY)), weight);
626  hist_y_trackx->Fill(positionX, 10.*residY, weight);
627  hist_y_tracky->Fill(positionY, 10.*residY, weight);
628  hist_y_trackdxdz->Fill(angleX, 10.*residY, weight);
629  hist_y_trackdydz->Fill(angleY, 10.*residY, weight);
630  fit_y_trackx->Fill(positionX, 10.*geom_residY, weight);
631  fit_y_tracky->Fill(positionY, 10.*geom_residY, weight);
632  fit_y_trackdxdz->Fill(angleX, 10.*geom_residY, weight);
633  fit_y_trackdydz->Fill(angleY, 10.*geom_residY, weight);
634 
635  double geom_resslopeX = MuonResiduals6DOFFitter_dxdz(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY);
636  hist_dxdz->Fill(1000.*(resslopeX - geom_resslopeX + value(kAlignPhiY)), weight);
637  hist_dxdz_trackx->Fill(positionX, 1000.*resslopeX, weight);
638  hist_dxdz_tracky->Fill(positionY, 1000.*resslopeX, weight);
639  hist_dxdz_trackdxdz->Fill(angleX, 1000.*resslopeX, weight);
640  hist_dxdz_trackdydz->Fill(angleY, 1000.*resslopeX, weight);
641  fit_dxdz_trackx->Fill(positionX, 1000.*geom_resslopeX, weight);
642  fit_dxdz_tracky->Fill(positionY, 1000.*geom_resslopeX, weight);
643  fit_dxdz_trackdxdz->Fill(angleX, 1000.*geom_resslopeX, weight);
644  fit_dxdz_trackdydz->Fill(angleY, 1000.*geom_resslopeX, weight);
645 
646  double geom_resslopeY = MuonResiduals6DOFFitter_dydz(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY);
647  hist_dydz->Fill(1000.*(resslopeY - geom_resslopeY - value(kAlignPhiX)), weight);
648  hist_dydz_trackx->Fill(positionX, 1000.*resslopeY, weight);
649  hist_dydz_tracky->Fill(positionY, 1000.*resslopeY, weight);
650  hist_dydz_trackdxdz->Fill(angleX, 1000.*resslopeY, weight);
651  hist_dydz_trackdydz->Fill(angleY, 1000.*resslopeY, weight);
652  fit_dydz_trackx->Fill(positionX, 1000.*geom_resslopeY, weight);
653  fit_dydz_tracky->Fill(positionY, 1000.*geom_resslopeY, weight);
654  fit_dydz_trackdxdz->Fill(angleX, 1000.*geom_resslopeY, weight);
655  fit_dydz_trackdydz->Fill(angleY, 1000.*geom_resslopeY, weight);
656  }
657 
658  hist_x_raw->Fill(10.*residX);
659  hist_y_raw->Fill(10.*residY);
660  hist_dxdz_raw->Fill(1000.*resslopeX);
661  hist_dydz_raw->Fill(1000.*resslopeY);
662  if (fabs(resslopeX) < 0.005) hist_x_cut->Fill(10.*residX);
663  if (fabs(resslopeY) < 0.030) hist_y_cut->Fill(10.*residY);
664  }
665 
666  double chi2 = 0.;
667  double ndof = 0.;
668  for (int i = 1; i <= hist_x->GetNbinsX(); i++) {
669  double xi = hist_x->GetBinCenter(i);
670  double yi = hist_x->GetBinContent(i);
671  double yerri = hist_x->GetBinError(i);
672  double yth = fit_x->Eval(xi);
673  if (yerri > 0.) {
674  chi2 += pow((yth - yi)/yerri, 2);
675  ndof += 1.;
676  }
677  }
678  for (int i = 1; i <= hist_y->GetNbinsX(); i++) {
679  double xi = hist_y->GetBinCenter(i);
680  double yi = hist_y->GetBinContent(i);
681  double yerri = hist_y->GetBinError(i);
682  double yth = fit_y->Eval(xi);
683  if (yerri > 0.) {
684  chi2 += pow((yth - yi)/yerri, 2);
685  ndof += 1.;
686  }
687  }
688  for (int i = 1; i <= hist_dxdz->GetNbinsX(); i++) {
689  double xi = hist_dxdz->GetBinCenter(i);
690  double yi = hist_dxdz->GetBinContent(i);
691  double yerri = hist_dxdz->GetBinError(i);
692  double yth = fit_dxdz->Eval(xi);
693  if (yerri > 0.) {
694  chi2 += pow((yth - yi)/yerri, 2);
695  ndof += 1.;
696  }
697  }
698  for (int i = 1; i <= hist_dydz->GetNbinsX(); i++) {
699  double xi = hist_dydz->GetBinCenter(i);
700  double yi = hist_dydz->GetBinContent(i);
701  double yerri = hist_dydz->GetBinError(i);
702  double yth = fit_dydz->Eval(xi);
703  if (yerri > 0.) {
704  chi2 += pow((yth - yi)/yerri, 2);
705  ndof += 1.;
706  }
707  }
708  ndof -= npar();
709 
710  return (ndof > 0. ? chi2 / ndof : -1.);
711 }
MuonResiduals6DOFFitter_y
double MuonResiduals6DOFFitter_y(double delta_x, double delta_y, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz, double alphay, double residual_dydz)
Definition: MuonResiduals6DOFFitter.cc:19
AlignableSurface::width
align::Scalar width() const
Definition: AlignableSurface.h:36
MuonResidualsFitter_logGaussPowerTails
double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma)
Definition: MuonResidualsFitter.cc:101
i
int i
Definition: DBlmapReader.cc:9
MuonResiduals6DOFFitter_y_trackx_TF1
Double_t MuonResiduals6DOFFitter_y_trackx_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:36
MuonResidualsFitter_powerLawTails_TF1
Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:89
MuonResiduals6DOFFitter_x_tracky_TF1
Double_t MuonResiduals6DOFFitter_x_tracky_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:32
launcher.step
list step
Definition: launcher.py:15
MuonResiduals6DOFFitter_y_tracky_TF1
Double_t MuonResiduals6DOFFitter_y_tracky_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:37
MuonResidualsFitterFitInfo::fitter
MuonResidualsFitter * fitter()
Definition: MuonResidualsFitter.h:141
NULL
#define NULL
Definition: scimark2.h:8
MuonResiduals6DOFFitter_x_trackx_TF1
Double_t MuonResiduals6DOFFitter_x_trackx_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:31
MuonResidualsFitter::value
double value(int parNum)
Definition: MuonResidualsFitter.h:92
MuonResiduals6DOFFitter_sum_of_weights
static double MuonResiduals6DOFFitter_sum_of_weights
Definition: MuonResiduals6DOFFitter.cc:7
MuonResiduals6DOFFitter::plot
double plot(std::string name, TFileDirectory *dir, Alignable *ali)
Definition: MuonResiduals6DOFFitter.cc:204
MuonResiduals6DOFFitter_y_trackdydz_TF1
Double_t MuonResiduals6DOFFitter_y_trackdydz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:39
MuonResiduals6DOFFitter_dxdz_trackx_TF1
Double_t MuonResiduals6DOFFitter_dxdz_trackx_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:41
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:192
MuonResidualsFitter_logROOTVoigt
double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma)
Definition: MuonResidualsFitter.cc:93
MuonResiduals6DOFFitter_dydz_trackdydz_TF1
Double_t MuonResiduals6DOFFitter_dydz_trackdydz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:49
MuonResiduals6DOFFitter_dydz_trackdxdz_TF1
Double_t MuonResiduals6DOFFitter_dydz_trackdxdz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:48
MuonResidualsFitter::fixed
bool fixed(int parNum)
Definition: MuonResidualsFitter.h:74
MuonResiduals6DOFFitter_dydz_trackx_TF1
Double_t MuonResiduals6DOFFitter_dydz_trackx_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:46
MuonResiduals6DOFFitter_dydz
double MuonResiduals6DOFFitter_dydz(double delta_x, double delta_y, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz)
Definition: MuonResiduals6DOFFitter.cc:27
MuonResidualsFitter_logPowerLawTails
double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma)
Definition: MuonResidualsFitter.cc:57
MuonResiduals6DOFFitter_dydz_tracky_TF1
Double_t MuonResiduals6DOFFitter_dydz_tracky_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:47
MuonResidualsFitter_GaussPowerTails_TF1
Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:112
Alignable::surface
const AlignableSurface & surface() const
Return the Surface (global position and orientation) of the object.
Definition: Alignable.h:126
MuonResiduals6DOFFitter_dxdz_trackdxdz_TF1
Double_t MuonResiduals6DOFFitter_dxdz_trackdxdz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:43
MuonResiduals6DOFFitter_dxdz
double MuonResiduals6DOFFitter_dxdz(double delta_x, double delta_y, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz)
Definition: MuonResiduals6DOFFitter.cc:23
MuonResidualsFitter_pureGaussian_TF1
Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:26
MuonResiduals6DOFFitter_x_trackdxdz_TF1
Double_t MuonResiduals6DOFFitter_x_trackdxdz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:33
MuonResiduals6DOFFitter::inform
void inform(TMinuit *tMinuit)
Definition: MuonResiduals6DOFFitter.cc:11
AlignableSurface::length
align::Scalar length() const
Definition: AlignableSurface.h:38
evf::utils::num
long long int num
Definition: procUtils.cc:71
MuonResidualsFitter::residuals_end
std::vector< double * >::const_iterator residuals_end() const
Definition: MuonResidualsFitter.h:113
MuonResiduals6DOFFitter_number_of_hits
static double MuonResiduals6DOFFitter_number_of_hits
Definition: MuonResiduals6DOFFitter.cc:8
MuonResiduals6DOFFitter_x_trackdydz_TF1
Double_t MuonResiduals6DOFFitter_x_trackdydz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:34
MuonResiduals6DOFFitter_x
double MuonResiduals6DOFFitter_x(double delta_x, double delta_y, double delta_z, double delta_phix, double delta_phiy, double delta_phiz, double track_x, double track_y, double track_dxdz, double track_dydz, double alphax, double residual_dxdz)
Definition: MuonResiduals6DOFFitter.cc:15
TFileDirectory::make
T * make() const
make new ROOT object
Definition: TFileDirectory.h:55
MuonResiduals6DOFFitter_weightAlignment
static bool MuonResiduals6DOFFitter_weightAlignment
Definition: MuonResiduals6DOFFitter.cc:9
MuonResiduals6DOFFitter_dxdz_trackdydz_TF1
Double_t MuonResiduals6DOFFitter_dxdz_trackdydz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:44
dir
dbl *** dir
Definition: mlp_gen.cc:35
MuonResiduals6DOFFitter_TMinuit
static TMinuit * MuonResiduals6DOFFitter_TMinuit
Definition: MuonResiduals6DOFFitter.cc:6
MuonResidualsFitter_ROOTVoigt_TF1
Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:97
MuonResidualsFitter::residuals_begin
std::vector< double * >::const_iterator residuals_begin() const
Definition: MuonResidualsFitter.h:112
MuonResiduals6DOFFitter_y_trackdxdz_TF1
Double_t MuonResiduals6DOFFitter_y_trackdxdz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:38
MuonResiduals6DOFFitter_FCN
void MuonResiduals6DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
Definition: MuonResiduals6DOFFitter.cc:51
MuonResidualsFitter_logPureGaussian
double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma)
Definition: MuonResidualsFitter.cc:20
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
MuonResiduals6DOFFitter_dxdz_tracky_TF1
Double_t MuonResiduals6DOFFitter_dxdz_tracky_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFFitter.cc:42
Gflash::par
const double par[8 *NPar][4]
Definition: GflashNameSpace.h:177
Generated for CMSSW Reference Manual by doxygen 1.8.5