Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:39:51 

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