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