Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:56:49 

0001 // $Id:$
0002 
0003 #ifdef STANDALONE_FITTER
0004 #include "MuonResiduals6DOFFitter.h"
0005 #else
0006 #include "Alignment/MuonAlignmentAlgorithms/interface/MuonResiduals6DOFFitter.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_y,
0023                     double delta_z,
0024                     double delta_phix,
0025                     double delta_phiy,
0026                     double delta_phiz,
0027                     double track_x,
0028                     double track_y,
0029                     double track_dxdz,
0030                     double track_dydz,
0031                     double alphax,
0032                     double residual_dxdz) {
0033     return delta_x - track_dxdz * delta_z - track_y * track_dxdz * delta_phix + track_x * track_dxdz * delta_phiy -
0034            track_y * delta_phiz + residual_dxdz * alphax;
0035   }
0036 
0037   double residual_y(double delta_x,
0038                     double delta_y,
0039                     double delta_z,
0040                     double delta_phix,
0041                     double delta_phiy,
0042                     double delta_phiz,
0043                     double track_x,
0044                     double track_y,
0045                     double track_dxdz,
0046                     double track_dydz,
0047                     double alphay,
0048                     double residual_dydz) {
0049     return delta_y - track_dydz * delta_z - track_y * track_dydz * delta_phix + track_x * track_dydz * delta_phiy +
0050            track_x * delta_phiz + residual_dydz * alphay;
0051   }
0052 
0053   double residual_dxdz(double delta_x,
0054                        double delta_y,
0055                        double delta_z,
0056                        double delta_phix,
0057                        double delta_phiy,
0058                        double delta_phiz,
0059                        double track_x,
0060                        double track_y,
0061                        double track_dxdz,
0062                        double track_dydz) {
0063     return -track_dxdz * track_dydz * delta_phix + (1. + track_dxdz * track_dxdz) * delta_phiy -
0064            track_dydz * delta_phiz;
0065   }
0066 
0067   double residual_dydz(double delta_x,
0068                        double delta_y,
0069                        double delta_z,
0070                        double delta_phix,
0071                        double delta_phiy,
0072                        double delta_phiz,
0073                        double track_x,
0074                        double track_y,
0075                        double track_dxdz,
0076                        double track_dydz) {
0077     return -(1. + track_dydz * track_dydz) * delta_phix + track_dxdz * track_dydz * delta_phiy +
0078            track_dxdz * delta_phiz;
0079   }
0080 
0081   Double_t residual_x_trackx_TF1(Double_t *xx, Double_t *p) {
0082     return 10. * residual_x(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9], p[10], p[11]);
0083   }
0084   Double_t residual_x_tracky_TF1(Double_t *xx, Double_t *p) {
0085     return 10. * residual_x(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9], p[10], p[11]);
0086   }
0087   Double_t residual_x_trackdxdz_TF1(Double_t *xx, Double_t *p) {
0088     return 10. * residual_x(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9], p[10], p[11]);
0089   }
0090   Double_t residual_x_trackdydz_TF1(Double_t *xx, Double_t *p) {
0091     return 10. * residual_x(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0], p[10], p[11]);
0092   }
0093 
0094   Double_t residual_y_trackx_TF1(Double_t *xx, Double_t *p) {
0095     return 10. * residual_y(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9], p[12], p[13]);
0096   }
0097   Double_t residual_y_tracky_TF1(Double_t *xx, Double_t *p) {
0098     return 10. * residual_y(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9], p[12], p[13]);
0099   }
0100   Double_t residual_y_trackdxdz_TF1(Double_t *xx, Double_t *p) {
0101     return 10. * residual_y(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9], p[12], p[13]);
0102   }
0103   Double_t residual_y_trackdydz_TF1(Double_t *xx, Double_t *p) {
0104     return 10. * residual_y(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0], p[12], p[13]);
0105   }
0106 
0107   Double_t residual_dxdz_trackx_TF1(Double_t *xx, Double_t *p) {
0108     return 1000. * residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9]);
0109   }
0110   Double_t residual_dxdz_tracky_TF1(Double_t *xx, Double_t *p) {
0111     return 1000. * residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9]);
0112   }
0113   Double_t residual_dxdz_trackdxdz_TF1(Double_t *xx, Double_t *p) {
0114     return 1000. * residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9]);
0115   }
0116   Double_t residual_dxdz_trackdydz_TF1(Double_t *xx, Double_t *p) {
0117     return 1000. * residual_dxdz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0]);
0118   }
0119 
0120   Double_t residual_dydz_trackx_TF1(Double_t *xx, Double_t *p) {
0121     return 1000. * residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], xx[0], p[7], p[8], p[9]);
0122   }
0123   Double_t residual_dydz_tracky_TF1(Double_t *xx, Double_t *p) {
0124     return 1000. * residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], xx[0], p[8], p[9]);
0125   }
0126   Double_t residual_dydz_trackdxdz_TF1(Double_t *xx, Double_t *p) {
0127     return 1000. * residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], xx[0], p[9]);
0128   }
0129   Double_t residual_dydz_trackdydz_TF1(Double_t *xx, Double_t *p) {
0130     return 1000. * residual_dydz(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], xx[0]);
0131   }
0132 }  // namespace
0133 
0134 void MuonResiduals6DOFFitter_FCN(int &npar, double *gin, double &fval, double *par, int iflag) {
0135   MuonResidualsFitterFitInfo *fitinfo = (MuonResidualsFitterFitInfo *)(minuit->GetObjectFit());
0136   MuonResidualsFitter *fitter = fitinfo->fitter();
0137 
0138   fval = 0.;
0139   for (std::vector<double *>::const_iterator resiter = fitter->residuals_begin(); resiter != fitter->residuals_end();
0140        ++resiter) {
0141     const double residX = (*resiter)[MuonResiduals6DOFFitter::kResidX];
0142     const double residY = (*resiter)[MuonResiduals6DOFFitter::kResidY];
0143     const double resslopeX = (*resiter)[MuonResiduals6DOFFitter::kResSlopeX];
0144     const double resslopeY = (*resiter)[MuonResiduals6DOFFitter::kResSlopeY];
0145     const double positionX = (*resiter)[MuonResiduals6DOFFitter::kPositionX];
0146     const double positionY = (*resiter)[MuonResiduals6DOFFitter::kPositionY];
0147     const double angleX = (*resiter)[MuonResiduals6DOFFitter::kAngleX];
0148     const double angleY = (*resiter)[MuonResiduals6DOFFitter::kAngleY];
0149     const double redchi2 = (*resiter)[MuonResiduals6DOFFitter::kRedChi2];
0150 
0151     const double alignx = par[MuonResiduals6DOFFitter::kAlignX];
0152     const double aligny = par[MuonResiduals6DOFFitter::kAlignY];
0153     const double alignz = par[MuonResiduals6DOFFitter::kAlignZ];
0154     const double alignphix = par[MuonResiduals6DOFFitter::kAlignPhiX];
0155     const double alignphiy = par[MuonResiduals6DOFFitter::kAlignPhiY];
0156     const double alignphiz = par[MuonResiduals6DOFFitter::kAlignPhiZ];
0157     const double resXsigma = par[MuonResiduals6DOFFitter::kResidXSigma];
0158     const double resYsigma = par[MuonResiduals6DOFFitter::kResidYSigma];
0159     const double slopeXsigma = par[MuonResiduals6DOFFitter::kResSlopeXSigma];
0160     const double slopeYsigma = par[MuonResiduals6DOFFitter::kResSlopeYSigma];
0161     const double alphax = par[MuonResiduals6DOFFitter::kAlphaX];
0162     const double alphay = par[MuonResiduals6DOFFitter::kAlphaY];
0163     const double resXgamma = par[MuonResiduals6DOFFitter::kResidXGamma];
0164     const double resYgamma = par[MuonResiduals6DOFFitter::kResidYGamma];
0165     const double slopeXgamma = par[MuonResiduals6DOFFitter::kResSlopeXGamma];
0166     const double slopeYgamma = par[MuonResiduals6DOFFitter::kResSlopeYGamma];
0167 
0168     double coefX = alphax, coefY = alphay;
0169     if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian ||
0170         fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D)
0171       coefX = coefY = 0.;
0172     double residXpeak = residual_x(
0173         alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coefX, resslopeX);
0174     double residYpeak = residual_y(
0175         alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY, coefY, resslopeY);
0176     double slopeXpeak =
0177         residual_dxdz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
0178     double slopeYpeak =
0179         residual_dydz(alignx, aligny, alignz, alignphix, alignphiy, alignphiz, positionX, positionY, angleX, angleY);
0180 
0181     double weight = (1. / redchi2) * number_of_hits / sum_of_weights;
0182     if (!weight_alignment)
0183       weight = 1.;
0184 
0185     if (!weight_alignment || TMath::Prob(redchi2 * 12, 12) < 0.99)  // no spikes allowed
0186     {
0187       if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian) {
0188         if (fitter->useRes() == MuonResidualsFitter::k1111) {
0189           fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
0190           fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
0191           fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
0192           fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeY, slopeYpeak, slopeYsigma);
0193         } else if (fitter->useRes() == MuonResidualsFitter::k1110) {
0194           fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
0195           fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
0196           fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
0197         } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
0198           fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
0199           fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
0200         } else if (fitter->useRes() == MuonResidualsFitter::k1010) {
0201           fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
0202           fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
0203         } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
0204           fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
0205         }
0206         //std::cout<<"FCNx("<<residX<<","<<residXpeak<<","<<resXsigma<<") = "<<MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma)<<std::endl;
0207         //std::cout<<"FCNy("<<residY<<","<<residYpeak<<","<<resYsigma<<") = "<<MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma)<<std::endl;
0208         //std::cout<<"FCNsx("<<resslopeX<<","<<slopeXpeak<<","<<slopeXsigma<<") = "<<MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma)<<std::endl;
0209         //std::cout<<"FCNsy("<<resslopeY<<","<<slopeYpeak<<","<<slopeYsigma<<") = "<<MuonResidualsFitter_logPureGaussian(resslopeY, slopeYpeak, slopeYsigma)<<std::endl;
0210       } else if (fitter->residualsModel() == MuonResidualsFitter::kPureGaussian2D) {
0211         if (fitter->useRes() == MuonResidualsFitter::k1111) {
0212           fval += -weight * MuonResidualsFitter_logPureGaussian2D(
0213                                 residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
0214           fval += -weight * MuonResidualsFitter_logPureGaussian2D(
0215                                 residY, resslopeY, residYpeak, slopeYpeak, resYsigma, slopeYsigma, alphay);
0216         } else if (fitter->useRes() == MuonResidualsFitter::k1110) {
0217           fval += -weight * MuonResidualsFitter_logPureGaussian2D(
0218                                 residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
0219           fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
0220         } else if (fitter->useRes() == MuonResidualsFitter::k1100) {
0221           fval += -weight * MuonResidualsFitter_logPureGaussian(residX, residXpeak, resXsigma);
0222           fval += -weight * MuonResidualsFitter_logPureGaussian(residY, residYpeak, resYsigma);
0223         } else if (fitter->useRes() == MuonResidualsFitter::k1010) {
0224           fval += -weight * MuonResidualsFitter_logPureGaussian2D(
0225                                 residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax);
0226         } else if (fitter->useRes() == MuonResidualsFitter::k0010) {
0227           fval += -weight * MuonResidualsFitter_logPureGaussian(resslopeX, slopeXpeak, slopeXsigma);
0228         }
0229         //std::cout<<"FCNx("<<residX<<","<<resslopeX<<","<<residXpeak<<","<<slopeXpeak<<","<<resXsigma<<","<<slopeXsigma<<","<<alphax<<") = "<<MuonResidualsFitter_logPureGaussian2D(residX, resslopeX, residXpeak, slopeXpeak, resXsigma, slopeXsigma, alphax)<<std::endl;
0230         //std::cout<<"FCNy("<<residY<<","<<resslopeY<<","<<residYpeak<<","<<slopeYpeak<<","<<resYsigma<<","<<slopeYsigma<<","<<alphay<<") = "<<MuonResidualsFitter_logPureGaussian2D(residY, resslopeY, residYpeak, slopeYpeak, resYsigma, slopeYsigma, alphay)<<std::endl;
0231       } else if (fitter->residualsModel() == MuonResidualsFitter::kPowerLawTails) {
0232         fval += -weight * MuonResidualsFitter_logPowerLawTails(residX, residXpeak, resXsigma, resXgamma);
0233         fval += -weight * MuonResidualsFitter_logPowerLawTails(residY, residYpeak, resYsigma, resYgamma);
0234         fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
0235         fval += -weight * MuonResidualsFitter_logPowerLawTails(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
0236       } else if (fitter->residualsModel() == MuonResidualsFitter::kROOTVoigt) {
0237         fval += -weight * MuonResidualsFitter_logROOTVoigt(residX, residXpeak, resXsigma, resXgamma);
0238         fval += -weight * MuonResidualsFitter_logROOTVoigt(residY, residYpeak, resYsigma, resYgamma);
0239         fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeX, slopeXpeak, slopeXsigma, slopeXgamma);
0240         fval += -weight * MuonResidualsFitter_logROOTVoigt(resslopeY, slopeYpeak, slopeYsigma, slopeYgamma);
0241       } else if (fitter->residualsModel() == MuonResidualsFitter::kGaussPowerTails) {
0242         fval += -weight * MuonResidualsFitter_logGaussPowerTails(residX, residXpeak, resXsigma);
0243         fval += -weight * MuonResidualsFitter_logGaussPowerTails(residY, residYpeak, resYsigma);
0244         fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeX, slopeXpeak, slopeXsigma);
0245         fval += -weight * MuonResidualsFitter_logGaussPowerTails(resslopeY, slopeYpeak, slopeYsigma);
0246       } else {
0247         assert(false);
0248       }
0249     }
0250   }
0251   /*
0252   const double resXsigma = par[MuonResiduals6DOFFitter::kResidXSigma];
0253   const double resYsigma = par[MuonResiduals6DOFFitter::kResidYSigma];
0254   const double slopeXsigma = par[MuonResiduals6DOFFitter::kResSlopeXSigma];
0255   const double slopeYsigma = par[MuonResiduals6DOFFitter::kResSlopeYSigma];
0256   const double alphax = par[MuonResiduals6DOFFitter::kAlphaX];
0257   const double alphay = par[MuonResiduals6DOFFitter::kAlphaY];
0258   std::cout<<"fval="<<fval<<","<<resXsigma<<","<<slopeXsigma<<","<<alphax<<","<<resYsigma<<","<<slopeYsigma<<","<<alphay<<std::endl;
0259 */
0260 }
0261 
0262 void MuonResiduals6DOFFitter::inform(TMinuit *tMinuit) { minuit = tMinuit; }
0263 
0264 void MuonResiduals6DOFFitter::correctBField() { MuonResidualsFitter::correctBField(kPt, kCharge); }
0265 
0266 double MuonResiduals6DOFFitter::sumofweights() {
0267   sum_of_weights = 0.;
0268   number_of_hits = 0.;
0269   weight_alignment = m_weightAlignment;
0270   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0271     if (m_weightAlignment) {
0272       double redchi2 = (*resiter)[kRedChi2];
0273       if (TMath::Prob(redchi2 * 12, 12) < 0.99) {  // no spikes allowed
0274         sum_of_weights += 1. / redchi2;
0275         number_of_hits += 1.;
0276       }
0277     } else {
0278       sum_of_weights += 1.;
0279       number_of_hits += 1.;
0280     }
0281   }
0282   return sum_of_weights;
0283 }
0284 
0285 bool MuonResiduals6DOFFitter::fit(Alignable *ali) {
0286   initialize_table();  // if not already initialized
0287   sumofweights();
0288 
0289   double resx_std = 0.5;
0290   double resy_std = 3.0;
0291   double resslopex_std = 0.002;
0292   double resslopey_std = 0.005;
0293 
0294   int nums[16] = {kAlignX,
0295                   kAlignY,
0296                   kAlignZ,
0297                   kAlignPhiX,
0298                   kAlignPhiY,
0299                   kAlignPhiZ,
0300                   kResidXSigma,
0301                   kResidYSigma,
0302                   kResSlopeXSigma,
0303                   kResSlopeYSigma,
0304                   kAlphaX,
0305                   kAlphaY,
0306                   kResidXGamma,
0307                   kResidYGamma,
0308                   kResSlopeXGamma,
0309                   kResSlopeYGamma};
0310   std::string names[16] = {"AlignX",
0311                            "AlignY",
0312                            "AlignZ",
0313                            "AlignPhiX",
0314                            "AlignPhiY",
0315                            "AlignPhiZ",
0316                            "ResidXSigma",
0317                            "ResidYSigma",
0318                            "ResSlopeXSigma",
0319                            "ResSlopeYSigma",
0320                            "AlphaX",
0321                            "AlphaY",
0322                            "ResidXGamma",
0323                            "ResidYGamma",
0324                            "ResSlopeXGamma",
0325                            "ResSlopeYGamma"};
0326   double starts[16] = {0.,
0327                        0.,
0328                        0.,
0329                        0.,
0330                        0.,
0331                        0.,
0332                        resx_std,
0333                        resy_std,
0334                        resslopex_std,
0335                        resslopey_std,
0336                        0.,
0337                        0.,
0338                        0.1 * resx_std,
0339                        0.1 * resy_std,
0340                        0.1 * resslopex_std,
0341                        0.1 * resslopey_std};
0342   double steps[16] = {0.1,
0343                       0.1,
0344                       0.1,
0345                       0.001,
0346                       0.001,
0347                       0.001,
0348                       0.001 * resx_std,
0349                       0.001 * resy_std,
0350                       0.001 * resslopex_std,
0351                       0.001 * resslopey_std,
0352                       0.001,
0353                       0.001,
0354                       0.01 * resx_std,
0355                       0.01 * resy_std,
0356                       0.01 * resslopex_std,
0357                       0.01 * resslopey_std};
0358   double lows[16] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1., -1., 0., 0., 0., 0.};
0359   double highs[16] = {0., 0., 0., 0., 0., 0., 10., 10., 0.1, 0.1, 1., 1., 0., 0., 0., 0.};
0360 
0361   std::vector<int> num(nums, nums + 6);
0362   std::vector<std::string> name(names, names + 6);
0363   std::vector<double> start(starts, starts + 6);
0364   std::vector<double> step(steps, steps + 6);
0365   std::vector<double> low(lows, lows + 6);
0366   std::vector<double> high(highs, highs + 6);
0367 
0368   bool add_alpha = (residualsModel() == kPureGaussian2D);
0369   bool add_gamma = (residualsModel() == kROOTVoigt || residualsModel() == kPowerLawTails);
0370 
0371   int idx[8], ni = 0;
0372   if (useRes() == k1111) {
0373     for (ni = 0; ni < 4; ni++)
0374       idx[ni] = ni + 6;
0375     if (add_alpha)
0376       for (; ni < 6; ni++)
0377         idx[ni] = ni + 6;
0378     else if (add_gamma)
0379       for (; ni < 8; ni++)
0380         idx[ni] = ni + 8;
0381     if (!add_alpha)
0382       fix(kAlphaX);
0383     if (!add_alpha)
0384       fix(kAlphaY);
0385   } else if (useRes() == k1110) {
0386     for (ni = 0; ni < 3; ni++)
0387       idx[ni] = ni + 6;
0388     if (add_alpha)
0389       idx[ni++] = 10;
0390     else if (add_gamma)
0391       for (; ni < 6; ni++)
0392         idx[ni] = ni + 9;
0393     fix(kResSlopeYSigma);
0394     fix(kAlphaY);
0395     if (!add_alpha)
0396       fix(kAlphaX);
0397   } else if (useRes() == k1100) {
0398     for (ni = 0; ni < 2; ni++)
0399       idx[ni] = ni + 6;
0400     if (add_gamma)
0401       for (; ni < 4; ni++)
0402         idx[ni] = ni + 10;
0403     fix(kResSlopeXSigma);
0404     fix(kResSlopeYSigma);
0405     fix(kAlphaX);
0406     fix(kAlphaY);
0407   } else if (useRes() == k1010) {
0408     idx[ni++] = 6;
0409     idx[ni++] = 8;
0410     if (add_alpha)
0411       idx[ni++] = 10;
0412     if (add_gamma) {
0413       idx[ni++] = 12;
0414       idx[ni++] = 14;
0415     }
0416     fix(kResidYSigma);
0417     fix(kResSlopeYSigma);
0418     fix(kAlphaY);
0419     if (!add_alpha)
0420       fix(kAlphaX);
0421   } else if (useRes() == k0010) {
0422     idx[ni++] = 8;
0423     if (add_gamma)
0424       idx[ni++] = 14;
0425     fix(kResidXSigma);
0426     fix(kResidYSigma);
0427     fix(kResSlopeYSigma);
0428     fix(kAlphaX);
0429     fix(kAlphaY);
0430   }
0431   for (int i = 0; i < ni; i++) {
0432     num.push_back(nums[idx[i]]);
0433     name.push_back(names[idx[i]]);
0434     start.push_back(starts[idx[i]]);
0435     step.push_back(steps[idx[i]]);
0436     low.push_back(lows[idx[i]]);
0437     high.push_back(highs[idx[i]]);
0438   }
0439 
0440   return dofit(&MuonResiduals6DOFFitter_FCN, num, name, start, step, low, high);
0441 }
0442 
0443 double MuonResiduals6DOFFitter::plot(std::string name, TFileDirectory *dir, Alignable *ali) {
0444   sumofweights();
0445 
0446   double mean_residualx = 0., mean_residualy = 0., mean_resslopex = 0., mean_resslopey = 0.;
0447   double mean_trackx = 0., mean_tracky = 0., mean_trackdxdz = 0., mean_trackdydz = 0.;
0448   double sum_w = 0.;
0449 
0450   for (std::vector<double *>::const_iterator rit = residuals_begin(); rit != residuals_end(); ++rit) {
0451     const double redchi2 = (*rit)[kRedChi2];
0452     double weight = 1. / redchi2;
0453     if (!m_weightAlignment)
0454       weight = 1.;
0455 
0456     if (!m_weightAlignment || TMath::Prob(redchi2 * 12, 12) < 0.99)  // no spikes allowed
0457     {
0458       double factor_w = 1. / (sum_w + weight);
0459       mean_residualx = factor_w * (sum_w * mean_residualx + weight * (*rit)[kResidX]);
0460       mean_residualy = factor_w * (sum_w * mean_residualy + weight * (*rit)[kResidY]);
0461       mean_resslopex = factor_w * (sum_w * mean_resslopex + weight * (*rit)[kResSlopeX]);
0462       mean_resslopey = factor_w * (sum_w * mean_resslopey + weight * (*rit)[kResSlopeY]);
0463       mean_trackx = factor_w * (sum_w * mean_trackx + weight * (*rit)[kPositionX]);
0464       mean_tracky = factor_w * (sum_w * mean_tracky + weight * (*rit)[kPositionY]);
0465       mean_trackdxdz = factor_w * (sum_w * mean_trackdxdz + weight * (*rit)[kAngleX]);
0466       mean_trackdydz = factor_w * (sum_w * mean_trackdydz + weight * (*rit)[kAngleY]);
0467       sum_w += weight;
0468     }
0469   }
0470 
0471   std::string name_x, name_y, name_dxdz, name_dydz, name_x_raw, name_y_raw, name_dxdz_raw, name_dydz_raw, name_x_cut,
0472       name_y_cut, name_alphax, name_alphay;
0473   std::string name_x_trackx, name_y_trackx, name_dxdz_trackx, name_dydz_trackx;
0474   std::string name_x_tracky, name_y_tracky, name_dxdz_tracky, name_dydz_tracky;
0475   std::string name_x_trackdxdz, name_y_trackdxdz, name_dxdz_trackdxdz, name_dydz_trackdxdz;
0476   std::string name_x_trackdydz, name_y_trackdydz, name_dxdz_trackdydz, name_dydz_trackdydz;
0477 
0478   name_x = name + "_x";
0479   name_y = name + "_y";
0480   name_dxdz = name + "_dxdz";
0481   name_dydz = name + "_dydz";
0482   name_x_raw = name + "_x_raw";
0483   name_y_raw = name + "_y_raw";
0484   name_dxdz_raw = name + "_dxdz_raw";
0485   name_dydz_raw = name + "_dydz_raw";
0486   name_x_cut = name + "_x_cut";
0487   name_y_cut = name + "_y_cut";
0488   name_alphax = name + "_alphax";
0489   name_alphay = name + "_alphay";
0490   name_x_trackx = name + "_x_trackx";
0491   name_y_trackx = name + "_y_trackx";
0492   name_dxdz_trackx = name + "_dxdz_trackx";
0493   name_dydz_trackx = name + "_dydz_trackx";
0494   name_x_tracky = name + "_x_tracky";
0495   name_y_tracky = name + "_y_tracky";
0496   name_dxdz_tracky = name + "_dxdz_tracky";
0497   name_dydz_tracky = name + "_dydz_tracky";
0498   name_x_trackdxdz = name + "_x_trackdxdz";
0499   name_y_trackdxdz = name + "_y_trackdxdz";
0500   name_dxdz_trackdxdz = name + "_dxdz_trackdxdz";
0501   name_dydz_trackdxdz = name + "_dydz_trackdxdz";
0502   name_x_trackdydz = name + "_x_trackdydz";
0503   name_y_trackdydz = name + "_y_trackdydz";
0504   name_dxdz_trackdydz = name + "_dxdz_trackdydz";
0505   name_dydz_trackdydz = name + "_dydz_trackdydz";
0506 
0507   double width = ali->surface().width();
0508   double length = ali->surface().length();
0509   int bins = 200;
0510   double min_x = -100., max_x = 100.;
0511   double min_y = -100., max_y = 100.;
0512   double min_dxdz = -75., max_dxdz = 75.;
0513   double min_dydz = -150., max_dydz = 150.;
0514   double min_trackx = -width / 2., max_trackx = width / 2.;
0515   double min_tracky = -length / 2., max_tracky = length / 2.;
0516   double min_trackdxdz = -1.5, max_trackdxdz = 1.5;
0517   double min_trackdydz = -1.5, max_trackdydz = 1.5;
0518 
0519   TH1F *hist_x = dir->make<TH1F>(name_x.c_str(), "", bins, min_x, max_x);
0520   TH1F *hist_y = dir->make<TH1F>(name_y.c_str(), "", bins, min_y, max_y);
0521   TH1F *hist_dxdz = dir->make<TH1F>(name_dxdz.c_str(), "", bins, min_dxdz, max_dxdz);
0522   TH1F *hist_dydz = dir->make<TH1F>(name_dydz.c_str(), "", bins, min_dydz, max_dydz);
0523   TH1F *hist_x_raw = dir->make<TH1F>(name_x_raw.c_str(), "", bins, min_x, max_x);
0524   TH1F *hist_y_raw = dir->make<TH1F>(name_y_raw.c_str(), "", bins, min_y, max_y);
0525   TH1F *hist_dxdz_raw = dir->make<TH1F>(name_dxdz_raw.c_str(), "", bins, min_dxdz, max_dxdz);
0526   TH1F *hist_dydz_raw = dir->make<TH1F>(name_dydz_raw.c_str(), "", bins, min_dydz, max_dydz);
0527   TH1F *hist_x_cut = dir->make<TH1F>(name_x_cut.c_str(), "", bins, min_x, max_x);
0528   TH1F *hist_y_cut = dir->make<TH1F>(name_y_cut.c_str(), "", bins, min_y, max_y);
0529   TH2F *hist_alphax = dir->make<TH2F>(name_alphax.c_str(), "", 50, 50, 50, 50, -50., 50.);
0530   TH2F *hist_alphay = dir->make<TH2F>(name_alphay.c_str(), "", 75, 100, 100, 75, -75., 75.);
0531 
0532   TProfile *hist_x_trackx = dir->make<TProfile>(name_x_trackx.c_str(), "", 50, min_trackx, max_trackx, min_x, max_x);
0533   TProfile *hist_y_trackx = dir->make<TProfile>(name_y_trackx.c_str(), "", 50, min_trackx, max_trackx, min_y, max_y);
0534   TProfile *hist_dxdz_trackx =
0535       dir->make<TProfile>(name_dxdz_trackx.c_str(), "", 50, min_trackx, max_trackx, min_dxdz, max_dxdz);
0536   TProfile *hist_dydz_trackx =
0537       dir->make<TProfile>(name_dydz_trackx.c_str(), "", 50, min_trackx, max_trackx, min_dydz, max_dydz);
0538   TProfile *hist_x_tracky = dir->make<TProfile>(name_x_tracky.c_str(), "", 50, min_tracky, max_tracky, min_x, max_x);
0539   TProfile *hist_y_tracky = dir->make<TProfile>(name_y_tracky.c_str(), "", 50, min_tracky, max_tracky, min_y, max_y);
0540   TProfile *hist_dxdz_tracky =
0541       dir->make<TProfile>(name_dxdz_tracky.c_str(), "", 50, min_tracky, max_tracky, min_dxdz, max_dxdz);
0542   TProfile *hist_dydz_tracky =
0543       dir->make<TProfile>(name_dydz_tracky.c_str(), "", 50, min_tracky, max_tracky, min_dydz, max_dydz);
0544   TProfile *hist_x_trackdxdz =
0545       dir->make<TProfile>(name_x_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_x, max_x);
0546   TProfile *hist_y_trackdxdz =
0547       dir->make<TProfile>(name_y_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_y, max_y);
0548   TProfile *hist_dxdz_trackdxdz =
0549       dir->make<TProfile>(name_dxdz_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_dxdz, max_dxdz);
0550   TProfile *hist_dydz_trackdxdz =
0551       dir->make<TProfile>(name_dydz_trackdxdz.c_str(), "", 250, min_trackdxdz, max_trackdxdz, min_dydz, max_dydz);
0552   TProfile *hist_x_trackdydz =
0553       dir->make<TProfile>(name_x_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_x, max_x);
0554   TProfile *hist_y_trackdydz =
0555       dir->make<TProfile>(name_y_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_y, max_y);
0556   TProfile *hist_dxdz_trackdydz =
0557       dir->make<TProfile>(name_dxdz_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_dxdz, max_dxdz);
0558   TProfile *hist_dydz_trackdydz =
0559       dir->make<TProfile>(name_dydz_trackdydz.c_str(), "", 250, min_trackdydz, max_trackdydz, min_dydz, max_dydz);
0560 
0561   hist_x_trackx->SetAxisRange(-10., 10., "Y");
0562   hist_y_trackx->SetAxisRange(-20., 20., "Y");
0563   hist_dxdz_trackx->SetAxisRange(-10., 10., "Y");
0564   hist_dydz_trackx->SetAxisRange(-20., 20., "Y");
0565   hist_x_tracky->SetAxisRange(-10., 10., "Y");
0566   hist_y_tracky->SetAxisRange(-20., 20., "Y");
0567   hist_dxdz_tracky->SetAxisRange(-10., 10., "Y");
0568   hist_dydz_tracky->SetAxisRange(-20., 20., "Y");
0569   hist_x_trackdxdz->SetAxisRange(-10., 10., "Y");
0570   hist_y_trackdxdz->SetAxisRange(-20., 20., "Y");
0571   hist_dxdz_trackdxdz->SetAxisRange(-10., 10., "Y");
0572   hist_dydz_trackdxdz->SetAxisRange(-20., 20., "Y");
0573   hist_x_trackdydz->SetAxisRange(-10., 10., "Y");
0574   hist_y_trackdydz->SetAxisRange(-20., 20., "Y");
0575   hist_dxdz_trackdydz->SetAxisRange(-10., 10., "Y");
0576   hist_dydz_trackdydz->SetAxisRange(-20., 20., "Y");
0577 
0578   name_x += "_fit";
0579   name_y += "_fit";
0580   name_dxdz += "_fit";
0581   name_dydz += "_fit";
0582   name_alphax += "_fit";
0583   name_alphay += "_fit";
0584   name_x_trackx += "_fit";
0585   name_y_trackx += "_fit";
0586   name_dxdz_trackx += "_fit";
0587   name_dydz_trackx += "_fit";
0588   name_x_tracky += "_fit";
0589   name_y_tracky += "_fit";
0590   name_dxdz_tracky += "_fit";
0591   name_dydz_tracky += "_fit";
0592   name_x_trackdxdz += "_fit";
0593   name_y_trackdxdz += "_fit";
0594   name_dxdz_trackdxdz += "_fit";
0595   name_dydz_trackdxdz += "_fit";
0596   name_x_trackdydz += "_fit";
0597   name_y_trackdydz += "_fit";
0598   name_dxdz_trackdydz += "_fit";
0599   name_dydz_trackdydz += "_fit";
0600 
0601   TF1 *fit_x = nullptr;
0602   TF1 *fit_y = nullptr;
0603   TF1 *fit_dxdz = nullptr;
0604   TF1 *fit_dydz = nullptr;
0605   if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D) {
0606     fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_x, max_x, 3);
0607     fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins, 10. * value(kAlignX), 10. * value(kResidXSigma));
0608     const double er_x[3] = {0., 10. * errorerror(kAlignX), 10. * errorerror(kResidXSigma)};
0609     fit_x->SetParErrors(er_x);
0610     fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_y, max_y, 3);
0611     fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins, 10. * value(kAlignY), 10. * value(kResidYSigma));
0612     const double er_y[3] = {0., 10. * errorerror(kAlignY), 10. * errorerror(kResidYSigma)};
0613     fit_y->SetParErrors(er_y);
0614     fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dxdz, max_dxdz, 3);
0615     fit_dxdz->SetParameters(
0616         sum_of_weights * (max_dxdz - min_dxdz) / bins, 1000. * value(kAlignPhiY), 1000. * value(kResSlopeXSigma));
0617     const double er_dxdz[3] = {0., 1000. * errorerror(kAlignPhiY), 1000. * errorerror(kResSlopeXSigma)};
0618     fit_dxdz->SetParErrors(er_dxdz);
0619     fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_pureGaussian_TF1, min_dydz, max_dydz, 3);
0620     fit_dydz->SetParameters(
0621         sum_of_weights * (max_dydz - min_dydz) / bins, -1000. * value(kAlignPhiX), 1000. * value(kResSlopeYSigma));
0622     const double er_dydz[3] = {0., 1000. * errorerror(kAlignPhiX), 1000. * errorerror(kResSlopeYSigma)};
0623     fit_dydz->SetParErrors(er_dydz);
0624   } else if (residualsModel() == kPowerLawTails) {
0625     fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_x, max_x, 4);
0626     fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins,
0627                          10. * value(kAlignX),
0628                          10. * value(kResidXSigma),
0629                          10. * value(kResidXGamma));
0630     fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_y, max_y, 4);
0631     fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins,
0632                          10. * value(kAlignY),
0633                          10. * value(kResidYSigma),
0634                          10. * value(kResidYGamma));
0635     fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dxdz, max_dxdz, 4);
0636     fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz) / bins,
0637                             1000. * value(kAlignPhiY),
0638                             1000. * value(kResSlopeXSigma),
0639                             1000. * value(kResSlopeXGamma));
0640     fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_powerLawTails_TF1, min_dydz, max_dydz, 4);
0641     fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz) / bins,
0642                             -1000. * value(kAlignPhiX),
0643                             1000. * value(kResSlopeYSigma),
0644                             1000. * value(kResSlopeYGamma));
0645   } else if (residualsModel() == kROOTVoigt) {
0646     fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_x, max_x, 4);
0647     fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins,
0648                          10. * value(kAlignX),
0649                          10. * value(kResidXSigma),
0650                          10. * value(kResidXGamma));
0651     fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_y, max_y, 4);
0652     fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins,
0653                          10. * value(kAlignY),
0654                          10. * value(kResidYSigma),
0655                          10. * value(kResidYGamma));
0656     fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dxdz, max_dxdz, 4);
0657     fit_dxdz->SetParameters(sum_of_weights * (max_dxdz - min_dxdz) / bins,
0658                             1000. * value(kAlignPhiY),
0659                             1000. * value(kResSlopeXSigma),
0660                             1000. * value(kResSlopeXGamma));
0661     fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_ROOTVoigt_TF1, min_dydz, max_dydz, 4);
0662     fit_dydz->SetParameters(sum_of_weights * (max_dydz - min_dydz) / bins,
0663                             -1000. * value(kAlignPhiX),
0664                             1000. * value(kResSlopeYSigma),
0665                             1000. * value(kResSlopeYGamma));
0666   } else if (residualsModel() == kGaussPowerTails) {
0667     fit_x = new TF1(name_x.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_x, max_x, 3);
0668     fit_x->SetParameters(sum_of_weights * (max_x - min_x) / bins, 10. * value(kAlignX), 10. * value(kResidXSigma));
0669     fit_y = new TF1(name_y.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_y, max_y, 3);
0670     fit_y->SetParameters(sum_of_weights * (max_y - min_y) / bins, 10. * value(kAlignY), 10. * value(kResidYSigma));
0671     fit_dxdz = new TF1(name_dxdz.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dxdz, max_dxdz, 3);
0672     fit_dxdz->SetParameters(
0673         sum_of_weights * (max_dxdz - min_dxdz) / bins, 1000. * value(kAlignPhiY), 1000. * value(kResSlopeXSigma));
0674     fit_dydz = new TF1(name_dydz.c_str(), MuonResidualsFitter_GaussPowerTails_TF1, min_dydz, max_dydz, 3);
0675     fit_dydz->SetParameters(
0676         sum_of_weights * (max_dydz - min_dydz) / bins, -1000. * value(kAlignPhiX), 1000. * value(kResSlopeYSigma));
0677   } else {
0678     assert(false);
0679   }
0680 
0681   fit_x->SetLineColor(2);
0682   fit_x->SetLineWidth(2);
0683   fit_x->Write();
0684   fit_y->SetLineColor(2);
0685   fit_y->SetLineWidth(2);
0686   fit_y->Write();
0687   fit_dxdz->SetLineColor(2);
0688   fit_dxdz->SetLineWidth(2);
0689   fit_dxdz->Write();
0690   fit_dydz->SetLineColor(2);
0691   fit_dydz->SetLineWidth(2);
0692   fit_dydz->Write();
0693 
0694   TF1 *fit_alphax = new TF1(name_alphax.c_str(), "[0] + x*[1]", min_dxdz, max_dxdz);
0695   TF1 *fit_alphay = new TF1(name_alphay.c_str(), "[0] + x*[1]", min_dydz, max_dydz);
0696   double aX = 10. * value(kAlignX), bX = 10. * value(kAlphaX) / 1000.;
0697   double aY = 10. * value(kAlignY), bY = 10. * value(kAlphaY) / 1000.;
0698   if (residualsModel() == kPureGaussian2D) {
0699     double sx = 10. * value(kResidXSigma), sy = 1000. * value(kResSlopeXSigma), r = value(kAlphaX);
0700     aX = mean_residualx;
0701     bX = 0.;
0702     if (sx != 0.) {
0703       bX = 1. / (sy / sx * r);
0704       aX = -bX * mean_resslopex;
0705     }
0706     sx = 10. * value(kResidYSigma);
0707     sy = 1000. * value(kResSlopeYSigma);
0708     r = value(kAlphaY);
0709     aY = mean_residualx;
0710     bY = 0.;
0711     if (sx != 0.) {
0712       bY = 1. / (sy / sx * r);
0713       aY = -bY * mean_resslopey;
0714     }
0715   }
0716   fit_alphax->SetParameters(aX, bX);
0717   fit_alphay->SetParameters(aY, bY);
0718 
0719   fit_alphax->SetLineColor(2);
0720   fit_alphax->SetLineWidth(2);
0721   fit_alphax->Write();
0722   fit_alphay->SetLineColor(2);
0723   fit_alphay->SetLineWidth(2);
0724   fit_alphay->Write();
0725 
0726   TProfile *fit_x_trackx = dir->make<TProfile>(name_x_trackx.c_str(), "", 100, min_trackx, max_trackx);
0727   TProfile *fit_y_trackx = dir->make<TProfile>(name_y_trackx.c_str(), "", 100, min_trackx, max_trackx);
0728   TProfile *fit_dxdz_trackx = dir->make<TProfile>(name_dxdz_trackx.c_str(), "", 100, min_trackx, max_trackx);
0729   TProfile *fit_dydz_trackx = dir->make<TProfile>(name_dydz_trackx.c_str(), "", 100, min_trackx, max_trackx);
0730   TProfile *fit_x_tracky = dir->make<TProfile>(name_x_tracky.c_str(), "", 100, min_tracky, max_tracky);
0731   TProfile *fit_y_tracky = dir->make<TProfile>(name_y_tracky.c_str(), "", 100, min_tracky, max_tracky);
0732   TProfile *fit_dxdz_tracky = dir->make<TProfile>(name_dxdz_tracky.c_str(), "", 100, min_tracky, max_tracky);
0733   TProfile *fit_dydz_tracky = dir->make<TProfile>(name_dydz_tracky.c_str(), "", 100, min_tracky, max_tracky);
0734   TProfile *fit_x_trackdxdz = dir->make<TProfile>(name_x_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
0735   TProfile *fit_y_trackdxdz = dir->make<TProfile>(name_y_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
0736   TProfile *fit_dxdz_trackdxdz =
0737       dir->make<TProfile>(name_dxdz_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
0738   TProfile *fit_dydz_trackdxdz =
0739       dir->make<TProfile>(name_dydz_trackdxdz.c_str(), "", 500, min_trackdxdz, max_trackdxdz);
0740   TProfile *fit_x_trackdydz = dir->make<TProfile>(name_x_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
0741   TProfile *fit_y_trackdydz = dir->make<TProfile>(name_y_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
0742   TProfile *fit_dxdz_trackdydz =
0743       dir->make<TProfile>(name_dxdz_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
0744   TProfile *fit_dydz_trackdydz =
0745       dir->make<TProfile>(name_dydz_trackdydz.c_str(), "", 500, min_trackdydz, max_trackdydz);
0746 
0747   fit_x_trackx->SetLineColor(2);
0748   fit_x_trackx->SetLineWidth(2);
0749   fit_y_trackx->SetLineColor(2);
0750   fit_y_trackx->SetLineWidth(2);
0751   fit_dxdz_trackx->SetLineColor(2);
0752   fit_dxdz_trackx->SetLineWidth(2);
0753   fit_dydz_trackx->SetLineColor(2);
0754   fit_dydz_trackx->SetLineWidth(2);
0755   fit_x_tracky->SetLineColor(2);
0756   fit_x_tracky->SetLineWidth(2);
0757   fit_y_tracky->SetLineColor(2);
0758   fit_y_tracky->SetLineWidth(2);
0759   fit_dxdz_tracky->SetLineColor(2);
0760   fit_dxdz_tracky->SetLineWidth(2);
0761   fit_dydz_tracky->SetLineColor(2);
0762   fit_dydz_tracky->SetLineWidth(2);
0763   fit_x_trackdxdz->SetLineColor(2);
0764   fit_x_trackdxdz->SetLineWidth(2);
0765   fit_y_trackdxdz->SetLineColor(2);
0766   fit_y_trackdxdz->SetLineWidth(2);
0767   fit_dxdz_trackdxdz->SetLineColor(2);
0768   fit_dxdz_trackdxdz->SetLineWidth(2);
0769   fit_dydz_trackdxdz->SetLineColor(2);
0770   fit_dydz_trackdxdz->SetLineWidth(2);
0771   fit_x_trackdydz->SetLineColor(2);
0772   fit_x_trackdydz->SetLineWidth(2);
0773   fit_y_trackdydz->SetLineColor(2);
0774   fit_y_trackdydz->SetLineWidth(2);
0775   fit_dxdz_trackdydz->SetLineColor(2);
0776   fit_dxdz_trackdydz->SetLineWidth(2);
0777   fit_dydz_trackdydz->SetLineColor(2);
0778   fit_dydz_trackdydz->SetLineWidth(2);
0779 
0780   name_x_trackx += "line";
0781   name_y_trackx += "line";
0782   name_dxdz_trackx += "line";
0783   name_dydz_trackx += "line";
0784   name_x_tracky += "line";
0785   name_y_tracky += "line";
0786   name_dxdz_tracky += "line";
0787   name_dydz_tracky += "line";
0788   name_x_trackdxdz += "line";
0789   name_y_trackdxdz += "line";
0790   name_dxdz_trackdxdz += "line";
0791   name_dydz_trackdxdz += "line";
0792   name_x_trackdydz += "line";
0793   name_y_trackdydz += "line";
0794   name_dxdz_trackdydz += "line";
0795   name_dydz_trackdydz += "line";
0796 
0797   TF1 *fitline_x_trackx = new TF1(name_x_trackx.c_str(), residual_x_trackx_TF1, min_trackx, max_trackx, 14);
0798   TF1 *fitline_y_trackx = new TF1(name_y_trackx.c_str(), residual_y_trackx_TF1, min_trackx, max_trackx, 14);
0799   TF1 *fitline_dxdz_trackx = new TF1(name_dxdz_trackx.c_str(), residual_dxdz_trackx_TF1, min_trackx, max_trackx, 14);
0800   TF1 *fitline_dydz_trackx = new TF1(name_dydz_trackx.c_str(), residual_dydz_trackx_TF1, min_trackx, max_trackx, 14);
0801   TF1 *fitline_x_tracky = new TF1(name_x_tracky.c_str(), residual_x_tracky_TF1, min_tracky, max_tracky, 14);
0802   TF1 *fitline_y_tracky = new TF1(name_y_tracky.c_str(), residual_y_tracky_TF1, min_tracky, max_tracky, 14);
0803   TF1 *fitline_dxdz_tracky = new TF1(name_dxdz_tracky.c_str(), residual_dxdz_tracky_TF1, min_tracky, max_tracky, 14);
0804   TF1 *fitline_dydz_tracky = new TF1(name_dydz_tracky.c_str(), residual_dydz_tracky_TF1, min_tracky, max_tracky, 14);
0805   TF1 *fitline_x_trackdxdz =
0806       new TF1(name_x_trackdxdz.c_str(), residual_x_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
0807   TF1 *fitline_y_trackdxdz =
0808       new TF1(name_y_trackdxdz.c_str(), residual_y_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
0809   TF1 *fitline_dxdz_trackdxdz =
0810       new TF1(name_dxdz_trackdxdz.c_str(), residual_dxdz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
0811   TF1 *fitline_dydz_trackdxdz =
0812       new TF1(name_dydz_trackdxdz.c_str(), residual_dydz_trackdxdz_TF1, min_trackdxdz, max_trackdxdz, 14);
0813   TF1 *fitline_x_trackdydz =
0814       new TF1(name_x_trackdydz.c_str(), residual_x_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
0815   TF1 *fitline_y_trackdydz =
0816       new TF1(name_y_trackdydz.c_str(), residual_y_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
0817   TF1 *fitline_dxdz_trackdydz =
0818       new TF1(name_dxdz_trackdydz.c_str(), residual_dxdz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
0819   TF1 *fitline_dydz_trackdydz =
0820       new TF1(name_dydz_trackdydz.c_str(), residual_dydz_trackdydz_TF1, min_trackdydz, max_trackdydz, 14);
0821 
0822   std::vector<TF1 *> fitlines;
0823   fitlines.push_back(fitline_x_trackx);
0824   fitlines.push_back(fitline_y_trackx);
0825   fitlines.push_back(fitline_dxdz_trackx);
0826   fitlines.push_back(fitline_dydz_trackx);
0827   fitlines.push_back(fitline_x_tracky);
0828   fitlines.push_back(fitline_y_tracky);
0829   fitlines.push_back(fitline_dxdz_tracky);
0830   fitlines.push_back(fitline_dydz_tracky);
0831   fitlines.push_back(fitline_x_trackdxdz);
0832   fitlines.push_back(fitline_y_trackdxdz);
0833   fitlines.push_back(fitline_dxdz_trackdxdz);
0834   fitlines.push_back(fitline_dydz_trackdxdz);
0835   fitlines.push_back(fitline_x_trackdydz);
0836   fitlines.push_back(fitline_y_trackdydz);
0837   fitlines.push_back(fitline_dxdz_trackdydz);
0838   fitlines.push_back(fitline_dydz_trackdydz);
0839 
0840   double fitparameters[14] = {value(kAlignX),
0841                               value(kAlignY),
0842                               value(kAlignZ),
0843                               value(kAlignPhiX),
0844                               value(kAlignPhiY),
0845                               value(kAlignPhiZ),
0846                               mean_trackx,
0847                               mean_tracky,
0848                               mean_trackdxdz,
0849                               mean_trackdydz,
0850                               value(kAlphaX),
0851                               mean_resslopex,
0852                               value(kAlphaY),
0853                               mean_resslopey};
0854   if (residualsModel() == kPureGaussian2D)
0855     fitparameters[10] = fitparameters[12] = 0.;
0856 
0857   for (std::vector<TF1 *>::const_iterator itr = fitlines.begin(); itr != fitlines.end(); itr++) {
0858     (*itr)->SetParameters(fitparameters);
0859     (*itr)->SetLineColor(2);
0860     (*itr)->SetLineWidth(2);
0861     (*itr)->Write();
0862   }
0863 
0864   for (std::vector<double *>::const_iterator resiter = residuals_begin(); resiter != residuals_end(); ++resiter) {
0865     const double residX = (*resiter)[kResidX];
0866     const double residY = (*resiter)[kResidY];
0867     const double resslopeX = (*resiter)[kResSlopeX];
0868     const double resslopeY = (*resiter)[kResSlopeY];
0869     const double positionX = (*resiter)[kPositionX];
0870     const double positionY = (*resiter)[kPositionY];
0871     const double angleX = (*resiter)[kAngleX];
0872     const double angleY = (*resiter)[kAngleY];
0873     const double redchi2 = (*resiter)[kRedChi2];
0874     double weight = 1. / redchi2;
0875     if (!m_weightAlignment)
0876       weight = 1.;
0877 
0878     if (!m_weightAlignment || TMath::Prob(redchi2 * 12, 12) < 0.99) {  // no spikes allowed
0879 
0880       hist_alphax->Fill(1000. * resslopeX, 10. * residX);
0881       hist_alphay->Fill(1000. * resslopeY, 10. * residY);
0882 
0883       double coefX = value(kAlphaX), coefY = value(kAlphaY);
0884       if (residualsModel() == kPureGaussian || residualsModel() == kPureGaussian2D)
0885         coefX = coefY = 0.;
0886       double geom_residX = residual_x(value(kAlignX),
0887                                       value(kAlignY),
0888                                       value(kAlignZ),
0889                                       value(kAlignPhiX),
0890                                       value(kAlignPhiY),
0891                                       value(kAlignPhiZ),
0892                                       positionX,
0893                                       positionY,
0894                                       angleX,
0895                                       angleY,
0896                                       coefX,
0897                                       resslopeX);
0898       hist_x->Fill(10. * (residX - geom_residX + value(kAlignX)), weight);
0899       hist_x_trackx->Fill(positionX, 10. * residX, weight);
0900       hist_x_tracky->Fill(positionY, 10. * residX, weight);
0901       hist_x_trackdxdz->Fill(angleX, 10. * residX, weight);
0902       hist_x_trackdydz->Fill(angleY, 10. * residX, weight);
0903       fit_x_trackx->Fill(positionX, 10. * geom_residX, weight);
0904       fit_x_tracky->Fill(positionY, 10. * geom_residX, weight);
0905       fit_x_trackdxdz->Fill(angleX, 10. * geom_residX, weight);
0906       fit_x_trackdydz->Fill(angleY, 10. * geom_residX, weight);
0907 
0908       double geom_residY = residual_y(value(kAlignX),
0909                                       value(kAlignY),
0910                                       value(kAlignZ),
0911                                       value(kAlignPhiX),
0912                                       value(kAlignPhiY),
0913                                       value(kAlignPhiZ),
0914                                       positionX,
0915                                       positionY,
0916                                       angleX,
0917                                       angleY,
0918                                       coefY,
0919                                       resslopeY);
0920       hist_y->Fill(10. * (residY - geom_residY + value(kAlignY)), weight);
0921       hist_y_trackx->Fill(positionX, 10. * residY, weight);
0922       hist_y_tracky->Fill(positionY, 10. * residY, weight);
0923       hist_y_trackdxdz->Fill(angleX, 10. * residY, weight);
0924       hist_y_trackdydz->Fill(angleY, 10. * residY, weight);
0925       fit_y_trackx->Fill(positionX, 10. * geom_residY, weight);
0926       fit_y_tracky->Fill(positionY, 10. * geom_residY, weight);
0927       fit_y_trackdxdz->Fill(angleX, 10. * geom_residY, weight);
0928       fit_y_trackdydz->Fill(angleY, 10. * geom_residY, weight);
0929 
0930       double geom_resslopeX = residual_dxdz(value(kAlignX),
0931                                             value(kAlignY),
0932                                             value(kAlignZ),
0933                                             value(kAlignPhiX),
0934                                             value(kAlignPhiY),
0935                                             value(kAlignPhiZ),
0936                                             positionX,
0937                                             positionY,
0938                                             angleX,
0939                                             angleY);
0940       hist_dxdz->Fill(1000. * (resslopeX - geom_resslopeX + value(kAlignPhiY)), weight);
0941       hist_dxdz_trackx->Fill(positionX, 1000. * resslopeX, weight);
0942       hist_dxdz_tracky->Fill(positionY, 1000. * resslopeX, weight);
0943       hist_dxdz_trackdxdz->Fill(angleX, 1000. * resslopeX, weight);
0944       hist_dxdz_trackdydz->Fill(angleY, 1000. * resslopeX, weight);
0945       fit_dxdz_trackx->Fill(positionX, 1000. * geom_resslopeX, weight);
0946       fit_dxdz_tracky->Fill(positionY, 1000. * geom_resslopeX, weight);
0947       fit_dxdz_trackdxdz->Fill(angleX, 1000. * geom_resslopeX, weight);
0948       fit_dxdz_trackdydz->Fill(angleY, 1000. * geom_resslopeX, weight);
0949 
0950       double geom_resslopeY = residual_dydz(value(kAlignX),
0951                                             value(kAlignY),
0952                                             value(kAlignZ),
0953                                             value(kAlignPhiX),
0954                                             value(kAlignPhiY),
0955                                             value(kAlignPhiZ),
0956                                             positionX,
0957                                             positionY,
0958                                             angleX,
0959                                             angleY);
0960       hist_dydz->Fill(1000. * (resslopeY - geom_resslopeY - value(kAlignPhiX)), weight);
0961       hist_dydz_trackx->Fill(positionX, 1000. * resslopeY, weight);
0962       hist_dydz_tracky->Fill(positionY, 1000. * resslopeY, weight);
0963       hist_dydz_trackdxdz->Fill(angleX, 1000. * resslopeY, weight);
0964       hist_dydz_trackdydz->Fill(angleY, 1000. * resslopeY, weight);
0965       fit_dydz_trackx->Fill(positionX, 1000. * geom_resslopeY, weight);
0966       fit_dydz_tracky->Fill(positionY, 1000. * geom_resslopeY, weight);
0967       fit_dydz_trackdxdz->Fill(angleX, 1000. * geom_resslopeY, weight);
0968       fit_dydz_trackdydz->Fill(angleY, 1000. * geom_resslopeY, weight);
0969     }
0970 
0971     hist_x_raw->Fill(10. * residX);
0972     hist_y_raw->Fill(10. * residY);
0973     hist_dxdz_raw->Fill(1000. * resslopeX);
0974     hist_dydz_raw->Fill(1000. * resslopeY);
0975     if (fabs(resslopeX) < 0.005)
0976       hist_x_cut->Fill(10. * residX);
0977     if (fabs(resslopeY) < 0.030)
0978       hist_y_cut->Fill(10. * residY);
0979   }
0980 
0981   double chi2 = 0.;
0982   double ndof = 0.;
0983   for (int i = 1; i <= hist_x->GetNbinsX(); i++) {
0984     double xi = hist_x->GetBinCenter(i);
0985     double yi = hist_x->GetBinContent(i);
0986     double yerri = hist_x->GetBinError(i);
0987     double yth = fit_x->Eval(xi);
0988     if (yerri > 0.) {
0989       chi2 += pow((yth - yi) / yerri, 2);
0990       ndof += 1.;
0991     }
0992   }
0993   for (int i = 1; i <= hist_y->GetNbinsX(); i++) {
0994     double xi = hist_y->GetBinCenter(i);
0995     double yi = hist_y->GetBinContent(i);
0996     double yerri = hist_y->GetBinError(i);
0997     double yth = fit_y->Eval(xi);
0998     if (yerri > 0.) {
0999       chi2 += pow((yth - yi) / yerri, 2);
1000       ndof += 1.;
1001     }
1002   }
1003   for (int i = 1; i <= hist_dxdz->GetNbinsX(); i++) {
1004     double xi = hist_dxdz->GetBinCenter(i);
1005     double yi = hist_dxdz->GetBinContent(i);
1006     double yerri = hist_dxdz->GetBinError(i);
1007     double yth = fit_dxdz->Eval(xi);
1008     if (yerri > 0.) {
1009       chi2 += pow((yth - yi) / yerri, 2);
1010       ndof += 1.;
1011     }
1012   }
1013   for (int i = 1; i <= hist_dydz->GetNbinsX(); i++) {
1014     double xi = hist_dydz->GetBinCenter(i);
1015     double yi = hist_dydz->GetBinContent(i);
1016     double yerri = hist_dydz->GetBinError(i);
1017     double yth = fit_dydz->Eval(xi);
1018     if (yerri > 0.) {
1019       chi2 += pow((yth - yi) / yerri, 2);
1020       ndof += 1.;
1021     }
1022   }
1023   ndof -= npar();
1024 
1025   return (ndof > 0. ? chi2 / ndof : -1.);
1026 }
1027 
1028 TTree *MuonResiduals6DOFFitter::readNtuple(
1029     std::string fname, unsigned int wheel, unsigned int station, unsigned int sector, unsigned int preselected) {
1030   TFile *f = new TFile(fname.c_str());
1031   TTree *t = (TTree *)f->Get("mual_ttree");
1032 
1033   // Create  new temporary file
1034   TFile *tmpf = new TFile("small_tree_fit.root", "recreate");
1035   assert(tmpf != nullptr);
1036 
1037   // filter the tree (temporarily lives in the new file)
1038   TTree *tt = t->CopyTree(Form(
1039       "is_dt && ring_wheel==%d && station==%d && sector==%d && select==%d", wheel, station, sector, (bool)preselected));
1040 
1041   MuonAlignmentTreeRow r;
1042   tt->SetBranchAddress("res_x", &r.res_x);
1043   tt->SetBranchAddress("res_slope_x", &r.res_slope_x);
1044   tt->SetBranchAddress("res_y", &r.res_y);
1045   tt->SetBranchAddress("res_slope_y", &r.res_slope_y);
1046   tt->SetBranchAddress("pos_x", &r.pos_x);
1047   tt->SetBranchAddress("pos_y", &r.pos_y);
1048   tt->SetBranchAddress("angle_x", &r.angle_x);
1049   tt->SetBranchAddress("angle_y", &r.angle_y);
1050   tt->SetBranchAddress("pz", &r.pz);
1051   tt->SetBranchAddress("pt", &r.pt);
1052   tt->SetBranchAddress("q", &r.q);
1053 
1054   Long64_t nentries = tt->GetEntries();
1055   for (Long64_t i = 0; i < nentries; i++) {
1056     tt->GetEntry(i);
1057     double *rdata = new double[MuonResiduals6DOFFitter::kNData];
1058     rdata[kResidX] = r.res_x;
1059     rdata[kResidY] = r.res_y;
1060     rdata[kResSlopeX] = r.res_slope_x;
1061     rdata[kResSlopeY] = r.res_slope_y;
1062     rdata[kPositionX] = r.pos_x;
1063     rdata[kPositionY] = r.pos_y;
1064     rdata[kAngleX] = r.angle_x;
1065     rdata[kAngleY] = r.angle_y;
1066     rdata[kRedChi2] = 0.1;
1067     rdata[kPz] = r.pz;
1068     rdata[kPt] = r.pt;
1069     rdata[kCharge] = r.q;
1070     fill(rdata);
1071   }
1072   delete f;
1073   //delete tmpf;
1074   return tt;
1075 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!