File indexing completed on 2023-03-17 10:39:51
0001
0002
0003 #ifdef STANDALONE_FITTER
0004 #include "MuonResiduals6DOFrphiFitter.h"
0005 #else
0006 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFrphiFitter.h"
0007
0008 #endif
0009
0010 #include "TH2F.h"
0011 #include "TMath.h"
0012 #include "TTree.h"
0013 #include "TFile.h"
0014
0015 namespace {
0016 TMinuit *minuit;
0017
0018 double sum_of_weights;
0019 double number_of_hits;
0020 bool weight_alignment;
0021
0022
0023
0024 double csc_R;
0025
0026 double getResidual(double delta_x,
0027 double delta_y,
0028 double delta_z,
0029 double delta_phix,
0030 double delta_phiy,
0031 double delta_phiz,
0032 double track_x,
0033 double track_y,
0034 double track_dxdz,
0035 double track_dydz,
0036 double R,
0037 double alpha,
0038 double resslope) {
0039 return delta_x - (track_x / R - 3. * pow(track_x / R, 3)) * delta_y - track_dxdz * delta_z -
0040 track_y * track_dxdz * delta_phix + track_x * track_dxdz * delta_phiy - track_y * delta_phiz +
0041 resslope * alpha;
0042 }
0043
0044 double getResSlope(double delta_x,
0045 double delta_y,
0046 double delta_z,
0047 double delta_phix,
0048 double delta_phiy,
0049 double delta_phiz,
0050 double track_x,
0051 double track_y,
0052 double track_dxdz,
0053 double track_dydz,
0054 double R) {
0055 return -track_dxdz / 2. / R * delta_y + (track_x / R - track_dxdz * track_dydz) * delta_phix +
0056 (1. + track_dxdz * track_dxdz) * delta_phiy - track_dydz * delta_phiz;
0057 }
0058
0059 double effectiveR(double x, double y) {
0060
0061
0062
0063
0064
0065
0066 return csc_R;
0067 }
0068
0069 Double_t residual_trackx_TF1(Double_t *xx, Double_t *p) {
0070 return 10. *
0071 getResidual(
0072 p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9], effectiveR(xx[0], p[7]), p[10], p[11]);
0073 }
0074 Double_t residual_tracky_TF1(Double_t *xx, Double_t *p) {
0075 return 10. *
0076 getResidual(
0077 p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9], effectiveR(p[6], xx[0]), p[10], p[11]);
0078 }
0079 Double_t residual_trackdxdz_TF1(Double_t *xx, Double_t *p) {
0080 return 10. * getResidual(
0081 p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9], effectiveR(p[6], p[7]), p[10], p[11]);
0082 }
0083 Double_t residual_trackdydz_TF1(Double_t *xx, Double_t *p) {
0084 return 10. * getResidual(
0085 p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0], effectiveR(p[6], p[7]), p[10], p[11]);
0086 }
0087
0088 Double_t resslope_trackx_TF1(Double_t *xx, Double_t *p) {
0089 return 1000. * getResSlope(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9], effectiveR(xx[0], p[7]));
0090 }
0091 Double_t resslope_tracky_TF1(Double_t *xx, Double_t *p) {
0092 return 1000. * getResSlope(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9], effectiveR(p[6], xx[0]));
0093 }
0094 Double_t resslope_trackdxdz_TF1(Double_t *xx, Double_t *p) {
0095 return 1000. * getResSlope(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9], effectiveR(p[6], p[7]));
0096 }
0097 Double_t resslope_trackdydz_TF1(Double_t *xx, Double_t *p) {
0098 return 1000. * getResSlope(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0], effectiveR(p[6], p[7]));
0099 }
0100 }
0101
0102 void MuonResiduals6DOFrphiFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
0103 MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo *)(minuit->GetObjectFit());
0104 MuonResidualsFitter *fitter = fitinfo->fitter();
0105
0106 fval = 0.;
0107 for (std::vector<double *>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end();
0108 ++resiter) {
0109 const double residual = (*resiter)[MuonResiduals6DOFrphiFitter::kResid];
0110 const double resslope = (*resiter)[MuonResiduals6DOFrphiFitter::kResSlope];
0111 const double positionX = (*resiter)[MuonResiduals6DOFrphiFitter::kPositionX];
0112 const double positionY = (*resiter)[MuonResiduals6DOFrphiFitter::kPositionY];
0113 const double angleX = (*resiter)[MuonResiduals6DOFrphiFitter::kAngleX];
0114 const double angleY = (*resiter)[MuonResiduals6DOFrphiFitter::kAngleY];
0115 const double redchi2 = (*resiter)[MuonResiduals6DOFrphiFitter::kRedChi2];
0116
0117 const double alignx = par[MuonResiduals6DOFrphiFitter::kAlignX];
0118 const double aligny = par[MuonResiduals6DOFrphiFitter::kAlignY];
0119 const double alignz = par[MuonResiduals6DOFrphiFitter::kAlignZ];
0120 const double alignphix = par[MuonResiduals6DOFrphiFitter::kAlignPhiX];
0121 const double alignphiy = par[MuonResiduals6DOFrphiFitter::kAlignPhiY];
0122 const double alignphiz = par[MuonResiduals6DOFrphiFitter::kAlignPhiZ];
0123 const double residsigma = par[MuonResiduals6DOFrphiFitter::kResidSigma];
0124 const double resslopesigma = par[MuonResiduals6DOFrphiFitter::kResSlopeSigma];
0125 const double alpha = par[MuonResiduals6DOFrphiFitter::kAlpha];
0126 const double residgamma = par[MuonResiduals6DOFrphiFitter::kResidGamma];
0127 const double resslopegamma = par[MuonResiduals6DOFrphiFitter::kResSlopeGamma];
0128
0129 double coeff = alpha;
0130 if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian ||
0131 fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D)
0132 coeff = 0.;
0133 double effr = effectiveR(positionX, positionY);
0134 double residpeak = getResidual(alignx,
0135 aligny,
0136 alignz,
0137 alignphix,
0138 alignphiy,
0139 alignphiz,
0140 positionX,
0141 positionY,
0142 angleX,
0143 angleY,
0144 effr,
0145 coeff,
0146 resslope);
0147 double resslopepeak = getResSlope(
0148 alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, effr);
0149
0150 double weight = (1. / redchi2) * number_of_hits / sum_of_weights;
0151 if (!weight_alignment)
0152 weight = 1.;
0153
0154 if (!weight_alignment || TMath::Prob(redchi2 * 6, 6) < 0.99)
0155 {
0156 if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
0157 if (fitter->useRes() == MuonResidualsFitter::k1111 || fitter->useRes() == MuonResidualsFitter::k1110 ||
0158 fitter->useRes() == MuonResidualsFitter::k1010) {
0159 fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
0160 fval += -weight * MuonResidualsFitter_logPureGaussian(resslope, resslopepeak, resslopesigma);
0161 } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
0162 fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
0163 } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
0164 fval += -weight * MuonResidualsFitter_logPureGaussian(resslope, resslopepeak, resslopesigma);
0165 }
0166 } else if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D) {
0167 if (fitter->useRes() == MuonResidualsFitter::k1111 || fitter->useRes() == MuonResidualsFitter::k1110 ||
0168 fitter->useRes() == MuonResidualsFitter::k1010) {
0169 fval += -weight * MuonResidualsFitter_logPureGaussian2D(
0170 residual, resslope, residpeak, resslopepeak, residsigma, resslopesigma, alpha);
0171 } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
0172 fval += -weight * MuonResidualsFitter_logPureGaussian(residual, residpeak, residsigma);
0173 } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
0174 fval += -weight * MuonResidualsFitter_logPureGaussian(resslope, resslopepeak, resslopesigma);
0175 }
0176 } else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
0177 fval += -weight * MuonResidualsFitter_logPowerLawTails(residual, residpeak, residsigma, residgamma);
0178 fval += -weight * MuonResidualsFitter_logPowerLawTails(resslope, resslopepeak, resslopesigma, resslopegamma);
0179 } else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
0180 fval += -weight * MuonResidualsFitter_logROOTVoigt(residual, residpeak, residsigma, residgamma);
0181 fval += -weight * MuonResidualsFitter_logROOTVoigt(resslope, resslopepeak, resslopesigma, resslopegamma);
0182 } else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
0183 fval += -weight * MuonResidualsFitter_logGaussPowerTails(residual, residpeak, residsigma);
0184 fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslope, resslopepeak, resslopesigma);
0185 } else {
0186 assert(false);
0187 }
0188 }
0189 }
0190 }
0191
0192 void MuonResiduals6DOFrphiFitter::correctBField() { MuonResidualsFitter::correctBField(kPz, kCharge); }
0193
0194 void MuonResiduals6DOFrphiFitter::inform(TMinuit *tMinuit) { minuit = tMinuit; }
0195
0196 double MuonResiduals6DOFrphiFitter::sumofweights() {
0197 sum_of_weights = 0.;
0198 number_of_hits = 0.;
0199 weight_alignment = m_weightAlignment;
0200 for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0201 if (m_weightAlignment) {
0202 double redchi2 = (*resiter)[kRedChi2];
0203 if (TMath::Prob(redchi2 * 6, 6) < 0.99) {
0204 sum_of_weights += 1. / redchi2;
0205 number_of_hits += 1.;
0206 }
0207 } else {
0208 sum_of_weights += 1.;
0209 number_of_hits += 1.;
0210 }
0211 }
0212 return sum_of_weights;
0213 }
0214
0215 bool MuonResiduals6DOFrphiFitter::fit(Alignable *ali) {
0216
0217
0218
0219 #ifndef STANDALONE_FITTER
0220 csc_R = sqrt(pow(ali->globalPosition().x(), 2) + pow(ali->globalPosition().y(), 2));
0221 #else
0222 csc_R = 200;
0223 #endif
0224
0225 initialize_table();
0226 sumofweights();
0227
0228 double res_std = 0.5;
0229 double resslope_std = 0.002;
0230
0231 int nums[10] = {kAlignX,
0232 kAlignZ,
0233 kAlignPhiX,
0234 kAlignPhiY,
0235 kAlignPhiZ,
0236 kResidSigma,
0237 kResSlopeSigma,
0238 kAlpha,
0239 kResidGamma,
0240 kResSlopeGamma};
0241 std::string names[10] = {"AlignX",
0242 "AlignZ",
0243 "AlignPhiX",
0244 "AlignPhiY",
0245 "AlignPhiZ",
0246 "ResidSigma",
0247 "ResSlopeSigma",
0248 "Alpha",
0249 "ResidGamma",
0250 "ResSlopeGamma"};
0251 double starts[10] = {0., 0., 0., 0., 0., res_std, resslope_std, 0., 0.1 * res_std, 0.1 * resslope_std};
0252 double steps[10] = {
0253 0.1, 0.1, 0.001, 0.001, 0.001, 0.001 * res_std, 0.001 * resslope_std, 0.001, 0.01 * res_std, 0.01 * resslope_std};
0254 double lows[10] = {0., 0., 0., 0., 0., 0., 0., -1., 0., 0.};
0255 double highs[10] = {0., 0., 0., 0., 0., 10., 0.1, 1., 0., 0.};
0256
0257 std::vector<int> num(nums, nums + 5);
0258 std::vector<std::string> name(names, names + 5);
0259 std::vector<double> start(starts, starts + 5);
0260 std::vector<double> step(steps, steps + 5);
0261 std::vector<double> low(lows, lows + 5);
0262 std::vector<double> high(highs, highs + 5);
0263
0264 bool add_alpha = (residualsModel() == kPureGaussian2D);
0265 bool add_gamma = (residualsModel() == kROOTVoigt || residualsModel() == kPowerLawTails);
0266
0267 int idx[4], ni = 0;
0268 if (useRes() == k1111 || useRes() == k1110 || useRes() == k1010) {
0269 for (ni = 0; ni < 2; ni++)
0270 idx[ni] = ni + 5;
0271 if (add_alpha)
0272 idx[ni++] = 7;
0273 else if (add_gamma)
0274 for (; ni < 4; ni++)
0275 idx[ni] = ni + 6;
0276 if (!add_alpha)
0277 fix(kAlpha);
0278 } else if (useRes() == k1100) {
0279 idx[ni++] = 5;
0280 if (add_gamma)
0281 idx[ni++] = 8;
0282 fix(kResSlopeSigma);
0283 fix(kAlpha);
0284 } else if (useRes() == k0010) {
0285 idx[ni++] = 6;
0286 if (add_gamma)
0287 idx[ni++] = 9;
0288 fix(kResidSigma);
0289 fix(kAlpha);
0290 }
0291 for (int i = 0; i < ni; i++) {
0292 num.push_back(nums[idx[i]]);
0293 name.push_back(names[idx[i]]);
0294 start.push_back(starts[idx[i]]);
0295 step.push_back(steps[idx[i]]);
0296 low.push_back(lows[idx[i]]);
0297 high.push_back(highs[idx[i]]);
0298 }
0299
0300 return dofit(&MuonResiduals6DOFrphiFitter_FCN, num, name, start, step, low, high);
0301 }
0302
0303 double MuonResiduals6DOFrphiFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
0304 sumofweights();
0305
0306 double mean_residual = 0., mean_resslope = 0.;
0307 double mean_trackx = 0., mean_tracky = 0., mean_trackdxdz = 0., mean_trackdydz = 0.;
0308 double sum_w = 0.;
0309
0310 for (std::vector<double *>::const_iterator rit = residuals_begin(); rit != residuals_end(); ++rit) {
0311 const double redchi2 = (*rit)[kRedChi2];
0312 double weight = 1. / redchi2;
0313 if (!m_weightAlignment)
0314 weight = 1.;
0315
0316 if (!m_weightAlignment || TMath::Prob(redchi2 * 6, 6) < 0.99)
0317 {
0318 double factor_w = 1. / (sum_w + weight);
0319 mean_residual = factor_w * (sum_w * mean_residual + weight * (*rit)[kResid]);
0320 mean_resslope = factor_w * (sum_w * mean_resslope + weight * (*rit)[kResSlope]);
0321 mean_trackx = factor_w * (sum_w * mean_trackx + weight * (*rit)[kPositionX]);
0322 mean_tracky = factor_w * (sum_w * mean_tracky + weight * (*rit)[kPositionY]);
0323 mean_trackdxdz = factor_w * (sum_w * mean_trackdxdz + weight * (*rit)[kAngleX]);
0324 mean_trackdydz = factor_w * (sum_w * mean_trackdydz + weight * (*rit)[kAngleY]);
0325 sum_w += weight;
0326 }
0327 }
0328
0329 std::string name_residual, name_resslope, name_residual_raw, name_resslope_raw, name_residual_cut, name_alpha;
0330 std::string name_residual_trackx, name_resslope_trackx;
0331 std::string name_residual_tracky, name_resslope_tracky;
0332 std::string name_residual_trackdxdz, name_resslope_trackdxdz;
0333 std::string name_residual_trackdydz, name_resslope_trackdydz;
0334
0335 name_residual = name + "_residual";
0336 name_resslope = name + "_resslope";
0337 name_residual_raw = name + "_residual_raw";
0338 name_resslope_raw = name + "_resslope_raw";
0339 name_residual_cut = name + "_residual_cut";
0340 name_alpha = name + "_alpha";
0341 name_residual_trackx = name + "_residual_trackx";
0342 name_resslope_trackx = name + "_resslope_trackx";
0343 name_residual_tracky = name + "_residual_tracky";
0344 name_resslope_tracky = name + "_resslope_tracky";
0345 name_residual_trackdxdz = name + "_residual_trackdxdz";
0346 name_resslope_trackdxdz = name + "_resslope_trackdxdz";
0347 name_residual_trackdydz = name + "_residual_trackdydz";
0348 name_resslope_trackdydz = name + "_resslope_trackdydz";
0349
0350 double width = ali->surface().width();
0351 double length = ali->surface().length();
0352 int bins_residual = 100, bins_resslope = 100;
0353 double min_residual = -50., max_residual = 50.;
0354 double min_resslope = -50., max_resslope = 50.;
0355 double min_trackx = -width / 2., max_trackx = width / 2.;
0356 double min_tracky = -length / 2., max_tracky = length / 2.;
0357 double min_trackdxdz = -1.5, max_trackdxdz = 1.5;
0358 double min_trackdydz = -1.5, max_trackdydz = 1.5;
0359
0360 TH1F *hist_residual = dir->make<TH1F>(name_residual.c_str(), "", bins_residual, min_residual, max_residual);
0361 TH1F *hist_resslope = dir->make<TH1F>(name_resslope.c_str(), "", bins_resslope, min_resslope, max_resslope);
0362 TH1F *hist_residual_raw = dir->make<TH1F>(name_residual_raw.c_str(), "", bins_residual, min_residual, max_residual);
0363 TH1F *hist_resslope_raw = dir->make<TH1F>(name_resslope_raw.c_str(), "", bins_resslope, min_resslope, max_resslope);
0364 TH1F *hist_residual_cut = dir->make<TH1F>(name_residual_cut.c_str(), "", bins_residual, min_residual, max_residual);
0365 TH2F *hist_alpha = dir->make<TH2F>(name_alpha.c_str(), "", 50, min_resslope, max_resslope, 50, -50., 50.);
0366
0367 TProfile *hist_residual_trackx =
0368 dir->make<TProfile>(name_residual_trackx.c_str(), "", 50, min_trackx, max_trackx, min_residual, max_residual);
0369 TProfile *hist_resslope_trackx =
0370 dir->make<TProfile>(name_resslope_trackx.c_str(), "", 50, min_trackx, max_trackx, min_resslope, max_resslope);
0371 TProfile *hist_residual_tracky =
0372 dir->make<TProfile>(name_residual_tracky.c_str(), "", 50, min_tracky, max_tracky, min_residual, max_residual);
0373 TProfile *hist_resslope_tracky =
0374 dir->make<TProfile>(name_resslope_tracky.c_str(), "", 50, min_tracky, max_tracky, min_resslope, max_resslope);
0375 TProfile *hist_residual_trackdxdz = dir->make<TProfile>(
0376 name_residual_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_residual, max_residual);
0377 TProfile *hist_resslope_trackdxdz = dir->make<TProfile>(
0378 name_resslope_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz, min_resslope, max_resslope);
0379 TProfile *hist_residual_trackdydz = dir->make<TProfile>(
0380 name_residual_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz, min_residual, max_residual);
0381 TProfile *hist_resslope_trackdydz = dir->make<TProfile>(
0382 name_resslope_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz, min_resslope, max_resslope);
0383
0384 hist_residual_trackx->SetAxisRange(-10., 10., "Y");
0385 hist_resslope_trackx->SetAxisRange(-10., 10., "Y");
0386 hist_residual_tracky->SetAxisRange(-10., 10., "Y");
0387 hist_resslope_tracky->SetAxisRange(-10., 10., "Y");
0388 hist_residual_trackdxdz->SetAxisRange(-10., 10., "Y");
0389 hist_resslope_trackdxdz->SetAxisRange(-10., 10., "Y");
0390 hist_residual_trackdydz->SetAxisRange(-10., 10., "Y");
0391 hist_resslope_trackdydz->SetAxisRange(-10., 10., "Y");
0392
0393 name_residual += "_fit";
0394 name_resslope += "_fit";
0395 name_alpha += "_fit";
0396 name_residual_trackx += "_fit";
0397 name_resslope_trackx += "_fit";
0398 name_residual_tracky += "_fit";
0399 name_resslope_tracky += "_fit";
0400 name_residual_trackdxdz += "_fit";
0401 name_resslope_trackdxdz += "_fit";
0402 name_residual_trackdydz += "_fit";
0403 name_resslope_trackdydz += "_fit";
0404
0405 TF1 *fit_residual = nullptr;
0406 TF1 *fit_resslope = nullptr;
0407 if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) {
0408 fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_residual, max_residual, 3);
0409 fit_residual->SetParameters(
0410 sum_of_weights * (max_residual - min_residual) / bins_residual, 10. * value(kAlignX), 10. * value(kResidSigma));
0411 const double er_res[3] = {0., 10. * errorerror(kAlignX), 10. * errorerror(kResidSigma)};
0412 fit_residual->SetParErrors(er_res);
0413 fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_resslope, max_resslope, 3);
0414 fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope) / bins_resslope,
0415 1000. * value(kAlignPhiY),
0416 1000. * value(kResSlopeSigma));
0417 const double er_resslope[3] = {0., 1000. * errorerror(kAlignPhiY), 1000. * errorerror(kResSlopeSigma)};
0418 fit_resslope->SetParErrors(er_resslope);
0419 } else if (residualsModel() == kPowerLawTails) {
0420 fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_residual, max_residual, 4);
0421 fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual) / bins_residual,
0422 10. * value(kAlignX),
0423 10. * value(kResidSigma),
0424 10. * value(kResidGamma));
0425 fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_resslope, max_resslope, 4);
0426 fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope) / bins_resslope,
0427 1000. * value(kAlignPhiY),
0428 1000. * value(kResSlopeSigma),
0429 1000. * value(kResSlopeGamma));
0430 } else if (residualsModel() == kROOTVoigt) {
0431 fit_residual = new TF1(name_residual.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_residual, max_residual, 4);
0432 fit_residual->SetParameters(sum_of_weights * (max_residual - min_residual) / bins_residual,
0433 10. * value(kAlignX),
0434 10. * value(kResidSigma),
0435 10. * value(kResidGamma));
0436 fit_resslope = new TF1(name_resslope.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_resslope, max_resslope, 4);
0437 fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope) / bins_resslope,
0438 1000. * value(kAlignPhiY),
0439 1000. * value(kResSlopeSigma),
0440 1000. * value(kResSlopeGamma));
0441 } else if (residualsModel() == kGaussPowerTails) {
0442 fit_residual =
0443 new TF1(name_residual.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_residual, max_residual, 3);
0444 fit_residual->SetParameters(
0445 sum_of_weights * (max_residual - min_residual) / bins_residual, 10. * value(kAlignX), 10. * value(kResidSigma));
0446 fit_resslope =
0447 new TF1(name_resslope.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_resslope, max_resslope, 3);
0448 fit_resslope->SetParameters(sum_of_weights * (max_resslope - min_resslope) / bins_resslope,
0449 1000. * value(kAlignPhiY),
0450 1000. * value(kResSlopeSigma));
0451 } else {
0452 assert(false);
0453 }
0454
0455 fit_residual->SetLineColor(2);
0456 fit_residual->SetLineWidth(2);
0457 fit_residual->Write();
0458 fit_resslope->SetLineColor(2);
0459 fit_resslope->SetLineWidth(2);
0460 fit_resslope->Write();
0461
0462 TF1 *fit_alpha = new TF1(name_alpha.c_str(), "[0] + x*[1]", min_resslope, max_resslope);
0463 double a = 10. * value(kAlignX), b = 10. * value(kAlpha) / 1000.;
0464 if (residualsModel() == kPureGaussian2D) {
0465 double sx = 10. * value(kResidSigma), sy = 1000. * value(kResSlopeSigma), r = value(kAlpha);
0466 a = mean_residual;
0467 b = 0.;
0468 if (sx != 0.) {
0469 b = 1. / (sy / sx * r);
0470 a = -b * mean_resslope;
0471 }
0472 }
0473 fit_alpha->SetParameters(a, b);
0474 fit_alpha->SetLineColor(2);
0475 fit_alpha->SetLineWidth(2);
0476 fit_alpha->Write();
0477
0478 TProfile *fit_residual_trackx = dir->make<TProfile>(name_residual_trackx.c_str(), "", 100, min_trackx, max_trackx);
0479 TProfile *fit_resslope_trackx = dir->make<TProfile>(name_resslope_trackx.c_str(), "", 100, min_trackx, max_trackx);
0480 TProfile *fit_residual_tracky = dir->make<TProfile>(name_residual_tracky.c_str(), "", 100, min_tracky, max_tracky);
0481 TProfile *fit_resslope_tracky = dir->make<TProfile>(name_resslope_tracky.c_str(), "", 100, min_tracky, max_tracky);
0482 TProfile *fit_residual_trackdxdz =
0483 dir->make<TProfile>(name_residual_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
0484 TProfile *fit_resslope_trackdxdz =
0485 dir->make<TProfile>(name_resslope_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
0486 TProfile *fit_residual_trackdydz =
0487 dir->make<TProfile>(name_residual_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
0488 TProfile *fit_resslope_trackdydz =
0489 dir->make<TProfile>(name_resslope_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
0490
0491 fit_residual_trackx->SetLineColor(2);
0492 fit_residual_trackx->SetLineWidth(2);
0493 fit_resslope_trackx->SetLineColor(2);
0494 fit_resslope_trackx->SetLineWidth(2);
0495 fit_residual_tracky->SetLineColor(2);
0496 fit_residual_tracky->SetLineWidth(2);
0497 fit_resslope_tracky->SetLineColor(2);
0498 fit_resslope_tracky->SetLineWidth(2);
0499 fit_residual_trackdxdz->SetLineColor(2);
0500 fit_residual_trackdxdz->SetLineWidth(2);
0501 fit_resslope_trackdxdz->SetLineColor(2);
0502 fit_resslope_trackdxdz->SetLineWidth(2);
0503 fit_residual_trackdydz->SetLineColor(2);
0504 fit_residual_trackdydz->SetLineWidth(2);
0505 fit_resslope_trackdydz->SetLineColor(2);
0506 fit_resslope_trackdydz->SetLineWidth(2);
0507
0508 name_residual_trackx += "line";
0509 name_resslope_trackx += "line";
0510 name_residual_tracky += "line";
0511 name_resslope_tracky += "line";
0512 name_residual_trackdxdz += "line";
0513 name_resslope_trackdxdz += "line";
0514 name_residual_trackdydz += "line";
0515 name_resslope_trackdydz += "line";
0516
0517 TF1 *fitline_residual_trackx = new TF1(name_residual_trackx.c_str(), residual_trackx_TF1, min_trackx, max_trackx, 12);
0518 TF1 *fitline_resslope_trackx = new TF1(name_resslope_trackx.c_str(), resslope_trackx_TF1, min_trackx, max_trackx, 12);
0519 TF1 *fitline_residual_tracky = new TF1(name_residual_tracky.c_str(), residual_tracky_TF1, min_tracky, max_tracky, 12);
0520 TF1 *fitline_resslope_tracky = new TF1(name_resslope_tracky.c_str(), resslope_tracky_TF1, min_tracky, max_tracky, 12);
0521 TF1 *fitline_residual_trackdxdz =
0522 new TF1(name_residual_trackdxdz.c_str(), residual_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
0523 TF1 *fitline_resslope_trackdxdz =
0524 new TF1(name_resslope_trackdxdz.c_str(), resslope_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 12);
0525 TF1 *fitline_residual_trackdydz =
0526 new TF1(name_residual_trackdydz.c_str(), residual_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);
0527 TF1 *fitline_resslope_trackdydz =
0528 new TF1(name_resslope_trackdydz.c_str(), resslope_trackdydz_TF1, min_trackdydz, max_trackdydz, 12);
0529
0530 std::vector<TF1 *> fitlines;
0531 fitlines.push_back(fitline_residual_trackx);
0532 fitlines.push_back(fitline_resslope_trackx);
0533 fitlines.push_back(fitline_residual_tracky);
0534 fitlines.push_back(fitline_resslope_tracky);
0535 fitlines.push_back(fitline_residual_trackdxdz);
0536 fitlines.push_back(fitline_resslope_trackdxdz);
0537 fitlines.push_back(fitline_residual_trackdydz);
0538 fitlines.push_back(fitline_resslope_trackdydz);
0539
0540 double fitparameters[12] = {value(kAlignX),
0541 value(kAlignY),
0542 value(kAlignZ),
0543 value(kAlignPhiX),
0544 value(kAlignPhiY),
0545 value(kAlignPhiZ),
0546 mean_trackx,
0547 mean_tracky,
0548 mean_trackdxdz,
0549 mean_trackdydz,
0550 value(kAlpha),
0551 mean_resslope};
0552 if (residualsModel() == kPureGaussian2D)
0553 fitparameters[10] = 0.;
0554
0555 for (std::vector<TF1 *>::const_iterator itr = fitlines.begin(); itr != fitlines.end(); itr++) {
0556 (*itr)->SetParameters(fitparameters);
0557 (*itr)->SetLineColor(2);
0558 (*itr)->SetLineWidth(2);
0559 (*itr)->Write();
0560 }
0561
0562 for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0563 const double resid = (*resiter)[kResid];
0564 const double resslope = (*resiter)[kResSlope];
0565 const double positionX = (*resiter)[kPositionX];
0566 const double positionY = (*resiter)[kPositionY];
0567 const double angleX = (*resiter)[kAngleX];
0568 const double angleY = (*resiter)[kAngleY];
0569 const double redchi2 = (*resiter)[kRedChi2];
0570 double weight = 1. / redchi2;
0571 if (!m_weightAlignment)
0572 weight = 1.;
0573
0574 if (!m_weightAlignment || TMath::Prob(redchi2 * 6, 6) < 0.99) {
0575 hist_alpha->Fill(1000. * resslope, 10. * resid);
0576
0577 double coeff = value(kAlpha);
0578 if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D)
0579 coeff = 0.;
0580 double geom_resid = getResidual(value(kAlignX),
0581 value(kAlignY),
0582 value(kAlignZ),
0583 value(kAlignPhiX),
0584 value(kAlignPhiY),
0585 value(kAlignPhiZ),
0586 positionX,
0587 positionY,
0588 angleX,
0589 angleY,
0590 effectiveR(positionX, positionY),
0591 coeff,
0592 resslope);
0593 hist_residual->Fill(10. * (resid - geom_resid + value(kAlignX)), weight);
0594 hist_residual_trackx->Fill(positionX, 10. * resid, weight);
0595 hist_residual_tracky->Fill(positionY, 10. * resid, weight);
0596 hist_residual_trackdxdz->Fill(angleX, 10. * resid, weight);
0597 hist_residual_trackdydz->Fill(angleY, 10. * resid, weight);
0598 fit_residual_trackx->Fill(positionX, 10. * geom_resid, weight);
0599 fit_residual_tracky->Fill(positionY, 10. * geom_resid, weight);
0600 fit_residual_trackdxdz->Fill(angleX, 10. * geom_resid, weight);
0601 fit_residual_trackdydz->Fill(angleY, 10. * geom_resid, weight);
0602
0603 double geom_resslope = getResSlope(value(kAlignX),
0604 value(kAlignY),
0605 value(kAlignZ),
0606 value(kAlignPhiX),
0607 value(kAlignPhiY),
0608 value(kAlignPhiZ),
0609 positionX,
0610 positionY,
0611 angleX,
0612 angleY,
0613 effectiveR(positionX, positionY));
0614 hist_resslope->Fill(1000. * (resslope - geom_resslope + value(kAlignPhiY)), weight);
0615 hist_resslope_trackx->Fill(positionX, 1000. * resslope, weight);
0616 hist_resslope_tracky->Fill(positionY, 1000. * resslope, weight);
0617 hist_resslope_trackdxdz->Fill(angleX, 1000. * resslope, weight);
0618 hist_resslope_trackdydz->Fill(angleY, 1000. * resslope, weight);
0619 fit_resslope_trackx->Fill(positionX, 1000. * geom_resslope, weight);
0620 fit_resslope_tracky->Fill(positionY, 1000. * geom_resslope, weight);
0621 fit_resslope_trackdxdz->Fill(angleX, 1000. * geom_resslope, weight);
0622 fit_resslope_trackdydz->Fill(angleY, 1000. * geom_resslope, weight);
0623 }
0624
0625 hist_residual_raw->Fill(10. * resid);
0626 hist_resslope_raw->Fill(1000. * resslope);
0627 if (fabs(resslope) < 0.005)
0628 hist_residual_cut->Fill(10. * resid);
0629 }
0630
0631 double chi2 = 0.;
0632 double ndof = 0.;
0633 for (int i = 1; i <= hist_residual->GetNbinsX(); i++) {
0634 double xi = hist_residual->GetBinCenter(i);
0635 double yi = hist_residual->GetBinContent(i);
0636 double yerri = hist_residual->GetBinError(i);
0637 double yth = fit_residual->Eval(xi);
0638 if (yerri > 0.) {
0639 chi2 += pow((yth - yi) / yerri, 2);
0640 ndof += 1.;
0641 }
0642 }
0643 for (int i = 1; i <= hist_resslope->GetNbinsX(); i++) {
0644 double xi = hist_resslope->GetBinCenter(i);
0645 double yi = hist_resslope->GetBinContent(i);
0646 double yerri = hist_resslope->GetBinError(i);
0647 double yth = fit_resslope->Eval(xi);
0648 if (yerri > 0.) {
0649 chi2 += pow((yth - yi) / yerri, 2);
0650 ndof += 1.;
0651 }
0652 }
0653 ndof -= npar();
0654
0655 return (ndof > 0. ? chi2 / ndof : -1.);
0656 }
0657
0658 TTree *MuonResiduals6DOFrphiFitter::readNtuple(std::string fname,
0659 unsigned int endcap,
0660 unsigned int station,
0661 unsigned int ring,
0662 unsigned int chamber,
0663 unsigned int preselected) {
0664 TFile *f = new TFile(fname.c_str());
0665 TTree *t = (TTree *)f->Get("mual_ttree");
0666
0667
0668 TFile *tmpf = new TFile("small_tree_fit.root", "recreate");
0669 assert(tmpf != nullptr);
0670
0671
0672 TTree *tt = t->CopyTree(Form("!is_dt && is_plus==%d && ring_wheel==%d && station==%d && sector==%d && select==%d",
0673 2 - endcap,
0674 station,
0675 ring,
0676 chamber,
0677 (bool)preselected));
0678
0679 MuonAlignmentTreeRow r;
0680 tt->SetBranchAddress("res_x", &r.res_x);
0681 tt->SetBranchAddress("res_slope_x", &r.res_slope_x);
0682 tt->SetBranchAddress("pos_x", &r.pos_x);
0683 tt->SetBranchAddress("pos_y", &r.pos_y);
0684 tt->SetBranchAddress("angle_x", &r.angle_x);
0685 tt->SetBranchAddress("angle_y", &r.angle_y);
0686 tt->SetBranchAddress("pz", &r.pz);
0687 tt->SetBranchAddress("pt", &r.pt);
0688 tt->SetBranchAddress("q", &r.q);
0689
0690 Long64_t nentries = tt->GetEntries();
0691 for (Long64_t i = 0; i < nentries; i++) {
0692 tt->GetEntry(i);
0693 double *rdata = new double[MuonResiduals6DOFrphiFitter::kNData];
0694 rdata[kResid] = r.res_x;
0695 rdata[kResSlope] = r.res_slope_x;
0696 rdata[kPositionX] = r.pos_x;
0697 rdata[kPositionY] = r.pos_y;
0698 rdata[kAngleX] = r.angle_x;
0699 rdata[kAngleY] = r.angle_y;
0700 rdata[kRedChi2] = 0.1;
0701 rdata[kPz] = r.pz;
0702 rdata[kPt] = r.pt;
0703 rdata[kCharge] = r.q;
0704 fill(rdata);
0705 }
0706 delete f;
0707
0708 return tt;
0709 }