Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:45:33

0001 // $Id:$
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 }  // namespace
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)  // no spikes allowed
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();  // if not already initialized
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)  // no spikes allowed
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) {  // no spikes allowed
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   // Create  new temporary file
0616   TFile *tmpf = new TFile("small_tree_fit.root", "recreate");
0617   assert(tmpf != nullptr);
0618 
0619   // filter the tree (temporarily lives in the new file)
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   //delete tmpf;
0652   return tt;
0653 }