CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonResiduals6DOFrphiFitter.cc
Go to the documentation of this file.
1 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFrphiFitter.h"
2 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
3 #include "TH2F.h"
4 #include "TMath.h"
5 
6 static TMinuit *MuonResiduals6DOFrphiFitter_TMinuit;
7 static double MuonResiduals6DOFrphiFitter_sum_of_weights;
8 static double MuonResiduals6DOFrphiFitter_number_of_hits;
9 static bool MuonResiduals6DOFrphiFitter_weightAlignment;
10 
11 const CSCGeometry *MuonResiduals6DOFrphiFitter_cscGeometry;
12 static CSCDetId MuonResiduals6DOFrphiFitter_cscDetId;
13 static double MuonResiduals6DOFrphiFitter_R;
14 
15 void MuonResiduals6DOFrphiFitter::inform(TMinuit *tMinuit) {
16  MuonResiduals6DOFrphiFitter_TMinuit = tMinuit;
17 }
18 
19 double MuonResiduals6DOFrphiFitter_residual(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 R, double alpha, double resslope) {
20  return delta_x - (track_x/R - 3.*pow(track_x/R, 3)) * delta_y - track_dxdz * delta_z - track_y * track_dxdz * delta_phix + track_x * track_dxdz * delta_phiy - track_y * delta_phiz + resslope * alpha;
21 }
22 
23 double MuonResiduals6DOFrphiFitter_resslope(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 R) {
24  return -track_dxdz/2./R * delta_y + (track_x/R - track_dxdz * track_dydz) * delta_phix + (1. + track_dxdz * track_dxdz) * delta_phiy - track_dydz * delta_phiz;
25 }
26 
27 double MuonResiduals6DOFrphiFitter_effectiveR(double x, double y) {
28 // CSCDetId id = MuonResiduals6DOFrphiFitter_cscDetId;
29 // CSCDetId layerId(id.endcap(), id.station(), id.ring(), id.chamber(), 3);
30 
31 // int strip = MuonResiduals6DOFrphiFitter_cscGeometry->layer(layerId)->geometry()->nearestStrip(LocalPoint(x, y, 0.));
32 // double angle = MuonResiduals6DOFrphiFitter_cscGeometry->layer(layerId)->geometry()->stripAngle(strip) - M_PI/2.;
33 
34 // if (fabs(angle) < 1e-10) return MuonResiduals6DOFrphiFitter_R;
35 // else return x/tan(angle) - y;
36 
37  return MuonResiduals6DOFrphiFitter_R;
38 }
39 
40 Double_t MuonResiduals6DOFrphiFitter_residual_trackx_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFrphiFitter_residual(par[0], par[1], par[2], par[3], par[4], par[5], xvec[0], par[7], par[8], par[9], MuonResiduals6DOFrphiFitter_effectiveR(xvec[0], par[7]), par[10], par[11]); }
41 Double_t MuonResiduals6DOFrphiFitter_residual_tracky_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFrphiFitter_residual(par[0], par[1], par[2], par[3], par[4], par[5], par[6], xvec[0], par[8], par[9], MuonResiduals6DOFrphiFitter_effectiveR(par[6], xvec[0]), par[10], par[11]); }
42 Double_t MuonResiduals6DOFrphiFitter_residual_trackdxdz_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFrphiFitter_residual(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], xvec[0], par[9], MuonResiduals6DOFrphiFitter_effectiveR(par[6], par[7]), par[10], par[11]); }
43 Double_t MuonResiduals6DOFrphiFitter_residual_trackdydz_TF1(Double_t *xvec, Double_t *par) { return 10.*MuonResiduals6DOFrphiFitter_residual(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], xvec[0], MuonResiduals6DOFrphiFitter_effectiveR(par[6], par[7]), par[10], par[11]); }
44 
45 Double_t MuonResiduals6DOFrphiFitter_resslope_trackx_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFrphiFitter_resslope(par[0], par[1], par[2], par[3], par[4], par[5], xvec[0], par[7], par[8], par[9], MuonResiduals6DOFrphiFitter_effectiveR(xvec[0], par[7])); }
46 Double_t MuonResiduals6DOFrphiFitter_resslope_tracky_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFrphiFitter_resslope(par[0], par[1], par[2], par[3], par[4], par[5], par[6], xvec[0], par[8], par[9], MuonResiduals6DOFrphiFitter_effectiveR(par[6], xvec[0])); }
47 Double_t MuonResiduals6DOFrphiFitter_resslope_trackdxdz_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFrphiFitter_resslope(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], xvec[0], par[9], MuonResiduals6DOFrphiFitter_effectiveR(par[6], par[7])); }
48 Double_t MuonResiduals6DOFrphiFitter_resslope_trackdydz_TF1(Double_t *xvec, Double_t *par) { return 1000.*MuonResiduals6DOFrphiFitter_resslope(par[0], par[1], par[2], par[3], par[4], par[5], par[6], par[7], par[8], xvec[0], MuonResiduals6DOFrphiFitter_effectiveR(par[6], par[7])); }
49 
50 void MuonResiduals6DOFrphiFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
51  MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo*)(MuonResiduals6DOFrphiFitter_TMinuit->GetObjectFit());
52  MuonResidualsFitter *fitter = fitinfo->fitter();
53 
54  fval = 0.;
55  for (std::vector<double*>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end(); ++resiter) {
56  const double residual = (*resiter)[MuonResiduals6DOFrphiFitter::kResid];
57  const double resslope = (*resiter)[MuonResiduals6DOFrphiFitter::kResSlope];
58  const double positionX = (*resiter)[MuonResiduals6DOFrphiFitter::kPositionX];
59  const double positionY = (*resiter)[MuonResiduals6DOFrphiFitter::kPositionY];
60  const double angleX = (*resiter)[MuonResiduals6DOFrphiFitter::kAngleX];
61  const double angleY = (*resiter)[MuonResiduals6DOFrphiFitter::kAngleY];
62  const double redchi2 = (*resiter)[MuonResiduals6DOFrphiFitter::kRedChi2];
63 
64  const double alignx = par[MuonResiduals6DOFrphiFitter::kAlignX];
65  const double aligny = par[MuonResiduals6DOFrphiFitter::kAlignY];
66  const double alignz = par[MuonResiduals6DOFrphiFitter::kAlignZ];
67  const double alignphix = par[MuonResiduals6DOFrphiFitter::kAlignPhiX];
68  const double alignphiy = par[MuonResiduals6DOFrphiFitter::kAlignPhiY];
69  const double alignphiz = par[MuonResiduals6DOFrphiFitter::kAlignPhiZ];
70  const double residsigma = par[MuonResiduals6DOFrphiFitter::kResidSigma];
71  const double resslopesigma = par[MuonResiduals6DOFrphiFitter::kResSlopeSigma];
72  const double alpha = par[MuonResiduals6DOFrphiFitter::kAlpha];
73  const double residgamma = par[MuonResiduals6DOFrphiFitter::kResidGamma];
74  const double resslopegamma = par[MuonResiduals6DOFrphiFitter::kResSlopeGamma];
75 
76  double residpeak = MuonResiduals6DOFrphiFitter_residual(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, MuonResiduals6DOFrphiFitter_effectiveR(positionX, positionY), alpha, resslope);
77  double resslopepeak = MuonResiduals6DOFrphiFitter_resslope(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, MuonResiduals6DOFrphiFitter_effectiveR(positionX, positionY));
78 
79  double weight = (1./redchi2) * MuonResiduals6DOFrphiFitter_number_of_hits / MuonResiduals6DOFrphiFitter_sum_of_weights;
80  if (!MuonResiduals6DOFrphiFitter_weightAlignment) weight = 1.;
81 
82  if (!MuonResiduals6DOFrphiFitter_weightAlignment || TMath::Prob(redchi2*6, 6) < 0.99) { // no spikes allowed
83 
84  if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
85  fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
86  fval += -weight * MuonResidualsFitter_logPureGaussian(resslope, resslopepeak, resslopesigma);
87  }
88  else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
89  fval += -weight * MuonResidualsFitter_logPowerLawTails(residual, residpeak, residsigma, residgamma);
90  fval += -weight * MuonResidualsFitter_logPowerLawTails(resslope, resslopepeak, resslopesigma, resslopegamma);
91  }
92  else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
93  fval += -weight * MuonResidualsFitter_logROOTVoigt(residual, residpeak, residsigma, residgamma);
94  fval += -weight * MuonResidualsFitter_logROOTVoigt(resslope, resslopepeak, resslopesigma, resslopegamma);
95  }
96  else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
97  fval += -weight * MuonResidualsFitter_logGaussPowerTails(residual, residpeak, residsigma);
98  fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslope, resslopepeak, resslopesigma);
99  }
100  else { assert(false); }
101 
102  }
103  }
104 }
105 
106 double MuonResiduals6DOFrphiFitter::sumofweights() {
107  MuonResiduals6DOFrphiFitter_sum_of_weights = 0.;
108  MuonResiduals6DOFrphiFitter_number_of_hits = 0.;
109  MuonResiduals6DOFrphiFitter_weightAlignment = m_weightAlignment;
110  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
111  if (m_weightAlignment) {
112  double redchi2 = (*resiter)[MuonResiduals6DOFrphiFitter::kRedChi2];
113  if (TMath::Prob(redchi2*6, 6) < 0.99) { // no spikes allowed
114  MuonResiduals6DOFrphiFitter_sum_of_weights += 1./redchi2;
115  MuonResiduals6DOFrphiFitter_number_of_hits += 1.;
116  }
117  }
118  else {
119  MuonResiduals6DOFrphiFitter_sum_of_weights += 1.;
120  MuonResiduals6DOFrphiFitter_number_of_hits += 1.;
121  }
122  }
123  return MuonResiduals6DOFrphiFitter_sum_of_weights;
124 }
125 
126 bool MuonResiduals6DOFrphiFitter::fit(Alignable *ali) {
127  initialize_table(); // if not already initialized
128  sumofweights();
129 
130  double resid_mean = 0;
131  double resslope_mean = 0;
132  double resid_stdev = 0.5;
133  double resslope_stdev = 0.005;
134  double alpha_estimate = 0;
135 
136  std::vector<int> num;
137  std::vector<std::string> name;
138  std::vector<double> start;
139  std::vector<double> step;
140  std::vector<double> low;
141  std::vector<double> high;
142 
143  if (fixed(kAlignX)) {
144  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.);
145  } else {
146  num.push_back(kAlignX); name.push_back(std::string("AlignX")); start.push_back(resid_mean); step.push_back(0.1); low.push_back(0.); high.push_back(0.);
147  }
148  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.);
149  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.);
150  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.);
151  if (fixed(kAlignPhiY)) {
152  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.);
153  } else {
154  num.push_back(kAlignPhiY); name.push_back(std::string("AlignPhiY")); start.push_back(resslope_mean); step.push_back(0.001); low.push_back(0.); high.push_back(0.);
155  }
156  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.);
157  num.push_back(kResidSigma); name.push_back(std::string("ResidSigma")); start.push_back(resid_stdev); step.push_back(0.01*resid_stdev); low.push_back(0.); high.push_back(0.);
158  num.push_back(kResSlopeSigma); name.push_back(std::string("ResSlopeSigma")); start.push_back(resslope_stdev); step.push_back(0.01*resslope_stdev); low.push_back(0.); high.push_back(0.);
159  num.push_back(kAlpha); name.push_back(std::string("Alpha")); start.push_back(alpha_estimate); step.push_back(0.01*resslope_stdev); low.push_back(0.); high.push_back(0.);
160  if (residualsModel() != kPureGaussian && residualsModel() != kGaussPowerTails) {
161  num.push_back(kResidGamma); name.push_back(std::string("ResidGamma")); start.push_back(0.1*resid_stdev); step.push_back(0.01*resid_stdev); low.push_back(0.); high.push_back(0.);
162  num.push_back(kResSlopeGamma); name.push_back(std::string("ResSlopeGamma")); start.push_back(0.1*resslope_stdev); step.push_back(0.01*resslope_stdev); low.push_back(0.); high.push_back(0.);
163  }
164 
165  return dofit(&MuonResiduals6DOFrphiFitter_FCN, num, name, start, step, low, high);
166 }
167 
168 double MuonResiduals6DOFrphiFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
169  sumofweights();
170 
171  MuonResiduals6DOFrphiFitter_R = sqrt(pow(ali->globalPosition().x(), 2) + pow(ali->globalPosition().y(), 2));
172  MuonResiduals6DOFrphiFitter_cscGeometry = m_cscGeometry;
173  MuonResiduals6DOFrphiFitter_cscDetId = CSCDetId(ali->geomDetId().rawId());
174 
175  std::stringstream name_residual, name_resslope, name_residual_raw, name_resslope_raw, name_residual_cut, name_alpha;
176  std::stringstream name_residual_trackx, name_resslope_trackx;
177  std::stringstream name_residual_tracky, name_resslope_tracky;
178  std::stringstream name_residual_trackdxdz, name_resslope_trackdxdz;
179  std::stringstream name_residual_trackdydz, name_resslope_trackdydz;
180 
181  name_residual << name << "_residual";
182  name_resslope << name << "_resslope";
183  name_residual_raw << name << "_residual_raw";
184  name_resslope_raw << name << "_resslope_raw";
185  name_residual_cut << name << "_residual_cut";
186  name_alpha << name << "_alpha";
187  name_residual_trackx << name << "_residual_trackx";
188  name_resslope_trackx << name << "_resslope_trackx";
189  name_residual_tracky << name << "_residual_tracky";
190  name_resslope_tracky << name << "_resslope_tracky";
191  name_residual_trackdxdz << name << "_residual_trackdxdz";
192  name_resslope_trackdxdz << name << "_resslope_trackdxdz";
193  name_residual_trackdydz << name << "_residual_trackdydz";
194  name_resslope_trackdydz << name << "_resslope_trackdydz";
195 
196  double width = ali->surface().width();
197  double length = ali->surface().length();
198  double min_residual = -100.; double max_residual = 100.;
199  double min_resslope = -100.; double max_resslope = 100.;
200  double min_trackx = -width/2.; double max_trackx = width/2.;
201  double min_tracky = -length/2.; double max_tracky = length/2.;
202  double min_trackdxdz = -1.5; double max_trackdxdz = 1.5;
203  double min_trackdydz = -1.5; double max_trackdydz = 1.5;
204 
205  TH1F *hist_residual = dir->make<TH1F>(name_residual.str().c_str(), "", 100, min_residual, max_residual);
206  TH1F *hist_resslope = dir->make<TH1F>(name_resslope.str().c_str(), "", 100, min_resslope, max_resslope);
207  TH1F *hist_residual_raw = dir->make<TH1F>(name_residual_raw.str().c_str(), "", 100, min_residual, max_residual);
208  TH1F *hist_resslope_raw = dir->make<TH1F>(name_resslope_raw.str().c_str(), "", 100, min_resslope, max_resslope);
209  TH1F *hist_residual_cut = dir->make<TH1F>(name_residual_cut.str().c_str(), "", 100, min_residual, max_residual);
210  TH2F *hist_alpha = dir->make<TH2F>(name_alpha.str().c_str(),"", 40, min_resslope, max_resslope, 40, -20., 20.);
211  TProfile *hist_residual_trackx = dir->make<TProfile>(name_residual_trackx.str().c_str(), "", 100, min_trackx, max_trackx, min_residual, max_residual);
212  TProfile *hist_resslope_trackx = dir->make<TProfile>(name_resslope_trackx.str().c_str(), "", 100, min_trackx, max_trackx, min_resslope, max_resslope);
213  TProfile *hist_residual_tracky = dir->make<TProfile>(name_residual_tracky.str().c_str(), "", 100, min_tracky, max_tracky, min_residual, max_residual);
214  TProfile *hist_resslope_tracky = dir->make<TProfile>(name_resslope_tracky.str().c_str(), "", 100, min_tracky, max_tracky, min_resslope, max_resslope);
215  TProfile *hist_residual_trackdxdz = dir->make<TProfile>(name_residual_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_residual, max_residual);
216  TProfile *hist_resslope_trackdxdz = dir->make<TProfile>(name_resslope_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_resslope, max_resslope);
217  TProfile *hist_residual_trackdydz = dir->make<TProfile>(name_residual_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz, min_residual, max_residual);
218  TProfile *hist_resslope_trackdydz = dir->make<TProfile>(name_resslope_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz, min_resslope, max_resslope);
219 
220  hist_residual_trackx->SetAxisRange(-10., 10., "Y");
221  hist_resslope_trackx->SetAxisRange(-10., 10., "Y");
222  hist_residual_tracky->SetAxisRange(-10., 10., "Y");
223  hist_resslope_tracky->SetAxisRange(-10., 10., "Y");
224  hist_residual_trackdxdz->SetAxisRange(-10., 10., "Y");
225  hist_resslope_trackdxdz->SetAxisRange(-10., 10., "Y");
226  hist_residual_trackdydz->SetAxisRange(-10., 10., "Y");
227  hist_resslope_trackdydz->SetAxisRange(-10., 10., "Y");
228 
229  name_residual << "_fit";
230  name_resslope << "_fit";
231  name_alpha << "_fit";
232  name_residual_trackx << "_fit";
233  name_resslope_trackx << "_fit";
234  name_residual_tracky << "_fit";
235  name_resslope_tracky << "_fit";
236  name_residual_trackdxdz << "_fit";
237  name_resslope_trackdxdz << "_fit";
238  name_residual_trackdydz << "_fit";
239  name_resslope_trackdydz << "_fit";
240 
241  TF1 *fit_residual = NULL;
242  TF1 *fit_resslope = NULL;
243  if (residualsModel() == kPureGaussian) {
244  fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, min_residual, max_residual, 3);
245  fit_residual->SetParameters(MuonResiduals6DOFrphiFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlignX), 10.*value(kResidSigma));
246  fit_resslope = new TF1(name_resslope.str().c_str(), MuonResidualsFitter_pureGaussian_TF1, min_resslope, max_resslope, 3);
247  fit_resslope->SetParameters(MuonResiduals6DOFrphiFitter_sum_of_weights * (max_resslope - min_resslope)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma));
248  }
249  else if (residualsModel() == kPowerLawTails) {
250  fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, min_residual, max_residual, 4);
251  fit_residual->SetParameters(MuonResiduals6DOFrphiFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlignX), 10.*value(kResidSigma), 10.*value(kResidGamma));
252  fit_resslope = new TF1(name_resslope.str().c_str(), MuonResidualsFitter_powerLawTails_TF1, min_resslope, max_resslope, 4);
253  fit_resslope->SetParameters(MuonResiduals6DOFrphiFitter_sum_of_weights * (max_resslope - min_resslope)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma), 1000.*value(kResSlopeGamma));
254  }
255  else if (residualsModel() == kROOTVoigt) {
256  fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_residual, max_residual, 4);
257  fit_residual->SetParameters(MuonResiduals6DOFrphiFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlignX), 10.*value(kResidSigma), 10.*value(kResidGamma));
258  fit_resslope = new TF1(name_resslope.str().c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_resslope, max_resslope, 4);
259  fit_resslope->SetParameters(MuonResiduals6DOFrphiFitter_sum_of_weights * (max_resslope - min_resslope)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma), 1000.*value(kResSlopeGamma));
260  }
261  else if (residualsModel() == kGaussPowerTails) {
262  fit_residual = new TF1(name_residual.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_residual, max_residual, 3);
263  fit_residual->SetParameters(MuonResiduals6DOFrphiFitter_sum_of_weights * (max_residual - min_residual)/100., 10.*value(kAlignX), 10.*value(kResidSigma));
264  fit_resslope = new TF1(name_resslope.str().c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_resslope, max_resslope, 3);
265  fit_resslope->SetParameters(MuonResiduals6DOFrphiFitter_sum_of_weights * (max_resslope - min_resslope)/100., 1000.*value(kAlignPhiY), 1000.*value(kResSlopeSigma));
266  }
267  else { assert(false); }
268 
269  fit_residual->SetLineColor(2); fit_residual->SetLineWidth(2);
270  fit_resslope->SetLineColor(2); fit_resslope->SetLineWidth(2);
271  fit_residual->Write();
272  fit_resslope->Write();
273 
274  TF1 *fit_alpha = new TF1(name_alpha.str().c_str(), "[0] + x*[1]", min_resslope, max_resslope);
275  fit_alpha->SetParameters(10.*value(kAlignX), 10.*value(kAlpha)/1000.);
276  fit_alpha->SetLineColor(2); fit_alpha->SetLineWidth(2);
277  fit_alpha->Write();
278 
279  TProfile *fit_residual_trackx = dir->make<TProfile>(name_residual_trackx.str().c_str(), "", 100, min_trackx, max_trackx);
280  TProfile *fit_resslope_trackx = dir->make<TProfile>(name_resslope_trackx.str().c_str(), "", 100, min_trackx, max_trackx);
281  TProfile *fit_residual_tracky = dir->make<TProfile>(name_residual_tracky.str().c_str(), "", 100, min_tracky, max_tracky);
282  TProfile *fit_resslope_tracky = dir->make<TProfile>(name_resslope_tracky.str().c_str(), "", 100, min_tracky, max_tracky);
283  TProfile *fit_residual_trackdxdz = dir->make<TProfile>(name_residual_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz);
284  TProfile *fit_resslope_trackdxdz = dir->make<TProfile>(name_resslope_trackdxdz.str().c_str(), "", 500, min_trackdxdz, max_trackdxdz);
285  TProfile *fit_residual_trackdydz = dir->make<TProfile>(name_residual_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz);
286  TProfile *fit_resslope_trackdydz = dir->make<TProfile>(name_resslope_trackdydz.str().c_str(), "", 500, min_trackdydz, max_trackdydz);
287 
288  fit_residual_trackx->SetLineColor(2); fit_residual_trackx->SetLineWidth(2);
289  fit_resslope_trackx->SetLineColor(2); fit_resslope_trackx->SetLineWidth(2);
290  fit_residual_tracky->SetLineColor(2); fit_residual_tracky->SetLineWidth(2);
291  fit_resslope_tracky->SetLineColor(2); fit_resslope_tracky->SetLineWidth(2);
292  fit_residual_trackdxdz->SetLineColor(2); fit_residual_trackdxdz->SetLineWidth(2);
293  fit_resslope_trackdxdz->SetLineColor(2); fit_resslope_trackdxdz->SetLineWidth(2);
294  fit_residual_trackdydz->SetLineColor(2); fit_residual_trackdydz->SetLineWidth(2);
295  fit_resslope_trackdydz->SetLineColor(2); fit_resslope_trackdydz->SetLineWidth(2);
296 
297  name_residual_trackx << "line";
298  name_resslope_trackx << "line";
299  name_residual_tracky << "line";
300  name_resslope_tracky << "line";
301  name_residual_trackdxdz << "line";
302  name_resslope_trackdxdz << "line";
303  name_residual_trackdydz << "line";
304  name_resslope_trackdydz << "line";
305 
306  TF1 *fitline_residual_trackx = new TF1(name_residual_trackx.str().c_str(), MuonResiduals6DOFrphiFitter_residual_trackx_TF1, min_trackx, max_trackx, 12);
307  TF1 *fitline_resslope_trackx = new TF1(name_resslope_trackx.str().c_str(), MuonResiduals6DOFrphiFitter_resslope_trackx_TF1, min_trackx, max_trackx, 12);
308  TF1 *fitline_residual_tracky = new TF1(name_residual_tracky.str().c_str(), MuonResiduals6DOFrphiFitter_residual_tracky_TF1, min_tracky, max_tracky, 12);
309  TF1 *fitline_resslope_tracky = new TF1(name_resslope_tracky.str().c_str(), MuonResiduals6DOFrphiFitter_resslope_tracky_TF1, min_tracky, max_tracky, 12);
310  TF1 *fitline_residual_trackdxdz = new TF1(name_residual_trackdxdz.str().c_str(), MuonResiduals6DOFrphiFitter_residual_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
311  TF1 *fitline_resslope_trackdxdz = new TF1(name_resslope_trackdxdz.str().c_str(), MuonResiduals6DOFrphiFitter_resslope_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
312  TF1 *fitline_residual_trackdydz = new TF1(name_residual_trackdydz.str().c_str(), MuonResiduals6DOFrphiFitter_residual_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);
313  TF1 *fitline_resslope_trackdydz = new TF1(name_resslope_trackdydz.str().c_str(), MuonResiduals6DOFrphiFitter_resslope_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);
314 
315  double sum_resslope = 0.;
316  double sum_trackx = 0.;
317  double sum_tracky = 0.;
318  double sum_trackdxdz = 0.;
319  double sum_trackdydz = 0.;
320  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
321  const double resslope = (*resiter)[MuonResiduals6DOFrphiFitter::kResSlope];
322  const double positionX = (*resiter)[MuonResiduals6DOFrphiFitter::kPositionX];
323  const double positionY = (*resiter)[MuonResiduals6DOFrphiFitter::kPositionY];
324  const double angleX = (*resiter)[MuonResiduals6DOFrphiFitter::kAngleX];
325  const double angleY = (*resiter)[MuonResiduals6DOFrphiFitter::kAngleY];
326  const double redchi2 = (*resiter)[MuonResiduals6DOFrphiFitter::kRedChi2];
327  double weight = 1./redchi2;
328  if (!m_weightAlignment) weight = 1.;
329 
330  if (!m_weightAlignment || TMath::Prob(redchi2*6, 6) < 0.99) { // no spikes allowed
331  sum_resslope += weight * resslope;
332  sum_trackx += weight * positionX;
333  sum_tracky += weight * positionY;
334  sum_trackdxdz += weight * angleX;
335  sum_trackdydz += weight * angleY;
336  }
337  }
338  double mean_resslope = sum_resslope / MuonResiduals6DOFrphiFitter_sum_of_weights;
339  double mean_trackx = sum_trackx / MuonResiduals6DOFrphiFitter_sum_of_weights;
340  double mean_tracky = sum_tracky / MuonResiduals6DOFrphiFitter_sum_of_weights;
341  double mean_trackdxdz = sum_trackdxdz / MuonResiduals6DOFrphiFitter_sum_of_weights;
342  double mean_trackdydz = sum_trackdydz / MuonResiduals6DOFrphiFitter_sum_of_weights;
343 
344  double fitparameters[12];
345  fitparameters[0] = value(kAlignX);
346  fitparameters[1] = value(kAlignY);
347  fitparameters[2] = value(kAlignZ);
348  fitparameters[3] = value(kAlignPhiX);
349  fitparameters[4] = value(kAlignPhiY);
350  fitparameters[5] = value(kAlignPhiZ);
351  fitparameters[6] = mean_trackx;
352  fitparameters[7] = mean_tracky;
353  fitparameters[8] = mean_trackdxdz;
354  fitparameters[9] = mean_trackdydz;
355  fitparameters[10] = value(kAlpha);
356  fitparameters[11] = mean_resslope;
357 
358  fitline_residual_trackx->SetParameters(fitparameters);
359  fitline_resslope_trackx->SetParameters(fitparameters);
360  fitline_residual_tracky->SetParameters(fitparameters);
361  fitline_resslope_tracky->SetParameters(fitparameters);
362  fitline_residual_trackdxdz->SetParameters(fitparameters);
363  fitline_resslope_trackdxdz->SetParameters(fitparameters);
364  fitline_residual_trackdydz->SetParameters(fitparameters);
365  fitline_resslope_trackdydz->SetParameters(fitparameters);
366 
367  fitline_residual_trackx->SetLineColor(2); fitline_residual_trackx->SetLineWidth(2);
368  fitline_resslope_trackx->SetLineColor(2); fitline_resslope_trackx->SetLineWidth(2);
369  fitline_residual_tracky->SetLineColor(2); fitline_residual_tracky->SetLineWidth(2);
370  fitline_resslope_tracky->SetLineColor(2); fitline_resslope_tracky->SetLineWidth(2);
371  fitline_residual_trackdxdz->SetLineColor(2); fitline_residual_trackdxdz->SetLineWidth(2);
372  fitline_resslope_trackdxdz->SetLineColor(2); fitline_resslope_trackdxdz->SetLineWidth(2);
373  fitline_residual_trackdydz->SetLineColor(2); fitline_residual_trackdydz->SetLineWidth(2);
374  fitline_resslope_trackdydz->SetLineColor(2); fitline_resslope_trackdydz->SetLineWidth(2);
375 
376  fitline_residual_trackx->Write();
377  fitline_resslope_trackx->Write();
378  fitline_residual_tracky->Write();
379  fitline_resslope_tracky->Write();
380  fitline_residual_trackdxdz->Write();
381  fitline_resslope_trackdxdz->Write();
382  fitline_residual_trackdydz->Write();
383  fitline_resslope_trackdydz->Write();
384 
385  for (std::vector<double*>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
386  const double resid = (*resiter)[MuonResiduals6DOFrphiFitter::kResid];
387  const double resslope = (*resiter)[MuonResiduals6DOFrphiFitter::kResSlope];
388  const double positionX = (*resiter)[MuonResiduals6DOFrphiFitter::kPositionX];
389  const double positionY = (*resiter)[MuonResiduals6DOFrphiFitter::kPositionY];
390  const double angleX = (*resiter)[MuonResiduals6DOFrphiFitter::kAngleX];
391  const double angleY = (*resiter)[MuonResiduals6DOFrphiFitter::kAngleY];
392  const double redchi2 = (*resiter)[MuonResiduals6DOFrphiFitter::kRedChi2];
393  double weight = 1./redchi2;
394  if (!m_weightAlignment) weight = 1.;
395 
396  if (!m_weightAlignment || TMath::Prob(redchi2*6, 6) < 0.99) { // no spikes allowed
397  hist_alpha->Fill(1000.*resslope, 10.*resid);
398 
399  double geom_resid = MuonResiduals6DOFrphiFitter_residual(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY, MuonResiduals6DOFrphiFitter_effectiveR(positionX, positionY), value(kAlpha), resslope);
400  hist_residual->Fill(10.*(resid - geom_resid + value(kAlignX)), weight);
401  hist_residual_trackx->Fill(positionX, 10.*resid, weight);
402  hist_residual_tracky->Fill(positionY, 10.*resid, weight);
403  hist_residual_trackdxdz->Fill(angleX, 10.*resid, weight);
404  hist_residual_trackdydz->Fill(angleY, 10.*resid, weight);
405  fit_residual_trackx->Fill(positionX, 10.*geom_resid, weight);
406  fit_residual_tracky->Fill(positionY, 10.*geom_resid, weight);
407  fit_residual_trackdxdz->Fill(angleX, 10.*geom_resid, weight);
408  fit_residual_trackdydz->Fill(angleY, 10.*geom_resid, weight);
409 
410  double geom_resslope = MuonResiduals6DOFrphiFitter_resslope(value(kAlignX), value(kAlignY), value(kAlignZ), value(kAlignPhiX), value(kAlignPhiY), value(kAlignPhiZ), positionX, positionY, angleX, angleY, MuonResiduals6DOFrphiFitter_effectiveR(positionX, positionY));
411  hist_resslope->Fill(1000.*(resslope - geom_resslope + value(kAlignPhiY)), weight);
412  hist_resslope_trackx->Fill(positionX, 1000.*resslope, weight);
413  hist_resslope_tracky->Fill(positionY, 1000.*resslope, weight);
414  hist_resslope_trackdxdz->Fill(angleX, 1000.*resslope, weight);
415  hist_resslope_trackdydz->Fill(angleY, 1000.*resslope, weight);
416  fit_resslope_trackx->Fill(positionX, 1000.*geom_resslope, weight);
417  fit_resslope_tracky->Fill(positionY, 1000.*geom_resslope, weight);
418  fit_resslope_trackdxdz->Fill(angleX, 1000.*geom_resslope, weight);
419  fit_resslope_trackdydz->Fill(angleY, 1000.*geom_resslope, weight);
420  }
421 
422  hist_residual_raw->Fill(10.*resid);
423  hist_resslope_raw->Fill(1000.*resslope);
424  if (fabs(resslope) < 0.005) hist_residual_cut->Fill(10.*resid);
425  }
426 
427  double chi2 = 0.;
428  double ndof = 0.;
429  for (int i = 1; i <= hist_residual->GetNbinsX(); i++) {
430  double xi = hist_residual->GetBinCenter(i);
431  double yi = hist_residual->GetBinContent(i);
432  double yerri = hist_residual->GetBinError(i);
433  double yth = fit_residual->Eval(xi);
434  if (yerri > 0.) {
435  chi2 += pow((yth - yi)/yerri, 2);
436  ndof += 1.;
437  }
438  }
439  for (int i = 1; i <= hist_resslope->GetNbinsX(); i++) {
440  double xi = hist_resslope->GetBinCenter(i);
441  double yi = hist_resslope->GetBinContent(i);
442  double yerri = hist_resslope->GetBinError(i);
443  double yth = fit_resslope->Eval(xi);
444  if (yerri > 0.) {
445  chi2 += pow((yth - yi)/yerri, 2);
446  ndof += 1.;
447  }
448  }
449  ndof -= npar();
450 
451  return (ndof > 0. ? chi2 / ndof : -1.);
452 }
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
alpha
float alpha
Definition: AMPTWrapper.h:95
MuonResiduals6DOFrphiFitter_resslope_trackdxdz_TF1
Double_t MuonResiduals6DOFrphiFitter_resslope_trackdxdz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFrphiFitter.cc:47
ExpressReco_HICollisions_FallBack.chi2
tuple chi2
Definition: ExpressReco_HICollisions_FallBack.py:26076
MuonResiduals6DOFrphiFitter_R
static double MuonResiduals6DOFrphiFitter_R
Definition: MuonResiduals6DOFrphiFitter.cc:13
MuonResidualsFitter_powerLawTails_TF1
Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:89
MuonResiduals6DOFrphiFitter_residual_trackdxdz_TF1
Double_t MuonResiduals6DOFrphiFitter_residual_trackdxdz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFrphiFitter.cc:42
PV3DBase::y
T y() const
Definition: PV3DBase.h:57
MuonResidualsFitterFitInfo::fitter
MuonResidualsFitter * fitter()
Definition: MuonResidualsFitter.h:141
MuonResiduals6DOFrphiFitter_residual_tracky_TF1
Double_t MuonResiduals6DOFrphiFitter_residual_tracky_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFrphiFitter.cc:41
NULL
#define NULL
Definition: scimark2.h:8
MuonResidualsFitter::value
double value(int parNum)
Definition: MuonResidualsFitter.h:92
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
MuonResiduals6DOFrphiFitter_resslope
double MuonResiduals6DOFrphiFitter_resslope(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 R)
Definition: MuonResiduals6DOFrphiFitter.cc:23
AlCaRecoCosmics_cfg.name
tuple name
Definition: AlCaRecoCosmics_cfg.py:24
MuonResiduals6DOFrphiFitter_effectiveR
double MuonResiduals6DOFrphiFitter_effectiveR(double x, double y)
Definition: MuonResiduals6DOFrphiFitter.cc:27
MuonResiduals6DOFrphiFitter_number_of_hits
static double MuonResiduals6DOFrphiFitter_number_of_hits
Definition: MuonResiduals6DOFrphiFitter.cc:8
DetId::rawId
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
MuonResiduals6DOFrphiFitter_resslope_trackdydz_TF1
Double_t MuonResiduals6DOFrphiFitter_resslope_trackdydz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFrphiFitter.cc:48
ExpressReco_HICollisions_FallBack.y
tuple y
Definition: ExpressReco_HICollisions_FallBack.py:4645
MuonResidualsFitter::fixed
bool fixed(int parNum)
Definition: MuonResidualsFitter.h:74
MuonResiduals6DOFrphiFitter_weightAlignment
static bool MuonResiduals6DOFrphiFitter_weightAlignment
Definition: MuonResiduals6DOFrphiFitter.cc:9
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:28
MuonResiduals6DOFrphiFitter_FCN
void MuonResiduals6DOFrphiFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag)
Definition: MuonResiduals6DOFrphiFitter.cc:50
MuonResiduals6DOFrphiFitter_residual
double MuonResiduals6DOFrphiFitter_residual(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 R, double alpha, double resslope)
Definition: MuonResiduals6DOFrphiFitter.cc:19
MuonResidualsFitter_logPowerLawTails
double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma)
Definition: MuonResidualsFitter.cc:57
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
ExpressReco_HICollisions_FallBack.step
tuple step
Definition: ExpressReco_HICollisions_FallBack.py:1252
MuonResiduals6DOFrphiFitter_TMinuit
static TMinuit * MuonResiduals6DOFrphiFitter_TMinuit
Definition: MuonResiduals6DOFrphiFitter.cc:6
MuonResidualsFitter_pureGaussian_TF1
Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResidualsFitter.cc:26
AlignableSurface::length
align::Scalar length() const
Definition: AlignableSurface.h:38
evf::utils::num
long long int num
Definition: procUtils.cc:71
MuonResiduals6DOFrphiFitter_sum_of_weights
static double MuonResiduals6DOFrphiFitter_sum_of_weights
Definition: MuonResiduals6DOFrphiFitter.cc:7
MuonResidualsFitter::residuals_end
std::vector< double * >::const_iterator residuals_end() const
Definition: MuonResidualsFitter.h:113
MuonResiduals6DOFrphiFitter_resslope_trackx_TF1
Double_t MuonResiduals6DOFrphiFitter_resslope_trackx_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFrphiFitter.cc:45
MuonResiduals6DOFrphiFitter_residual_trackdydz_TF1
Double_t MuonResiduals6DOFrphiFitter_residual_trackdydz_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFrphiFitter.cc:43
ExpressReco_HICollisions_FallBack.x
tuple x
Definition: ExpressReco_HICollisions_FallBack.py:4646
MuonResiduals6DOFrphiFitter_resslope_tracky_TF1
Double_t MuonResiduals6DOFrphiFitter_resslope_tracky_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFrphiFitter.cc:46
TFileDirectory::make
T * make() const
make new ROOT object
Definition: TFileDirectory.h:59
Alignable::globalPosition
const PositionType & globalPosition() const
Return the global position of the object.
Definition: Alignable.h:129
MuonResiduals6DOFrphiFitter::plot
double plot(std::string name, TFileDirectory *dir, Alignable *ali)
Definition: MuonResiduals6DOFrphiFitter.cc:168
MuonResiduals6DOFrphiFitter_cscGeometry
const CSCGeometry * MuonResiduals6DOFrphiFitter_cscGeometry
Definition: MuonResiduals6DOFrphiFitter.cc:11
dir
dbl *** dir
Definition: mlp_gen.cc:35
MuonResiduals6DOFrphiFitter_residual_trackx_TF1
Double_t MuonResiduals6DOFrphiFitter_residual_trackx_TF1(Double_t *xvec, Double_t *par)
Definition: MuonResiduals6DOFrphiFitter.cc:40
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
PV3DBase::x
T x() const
Definition: PV3DBase.h:56
MuonResidualsFitter_logPureGaussian
double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma)
Definition: MuonResidualsFitter.cc:20
Alignable::geomDetId
const DetId & geomDetId() const
Definition: Alignable.h:177
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
MuonResiduals6DOFrphiFitter_cscDetId
static CSCDetId MuonResiduals6DOFrphiFitter_cscDetId
Definition: MuonResiduals6DOFrphiFitter.cc:12
Gflash::par
const double par[8 *NPar][4]
Definition: GflashNameSpace.h:177
Generated for CMSSW Reference Manual by doxygen 1.8.5