Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:41 

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  *   Laurent Forthomme
0005  ****************************************************************************/
0006 
0007 #include "RecoPPS/ProtonReconstruction/interface/ProtonReconstructionAlgorithm.h"
0008 
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 
0011 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0012 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0013 #include "DataFormats/ProtonReco/interface/ForwardProton.h"
0014 
0015 #include "TMinuitMinimizer.h"
0016 
0017 using namespace std;
0018 using namespace edm;
0019 
0020 //----------------------------------------------------------------------------------------------------
0021 
0022 ProtonReconstructionAlgorithm::ProtonReconstructionAlgorithm(bool fit_vtx_y,
0023                                                              bool improved_estimate,
0024                                                              const std::string &multiRPAlgorithm,
0025                                                              unsigned int verbosity)
0026     : verbosity_(verbosity),
0027       fitVtxY_(fit_vtx_y),
0028       useImprovedInitialEstimate_(improved_estimate),
0029       multi_rp_algorithm_(mrpaUndefined),
0030       initialized_(false),
0031       fitter_(new ROOT::Fit::Fitter),
0032       chiSquareCalculator_(new ChiSquareCalculator) {
0033   // needed for thread safety
0034   TMinuitMinimizer::UseStaticMinuit(false);
0035 
0036   // check and set multi-RP algorithm
0037   if (multiRPAlgorithm == "chi2")
0038     multi_rp_algorithm_ = mrpaChi2;
0039   else if (multiRPAlgorithm == "newton")
0040     multi_rp_algorithm_ = mrpaNewton;
0041   else if (multiRPAlgorithm == "anal-iter")
0042     multi_rp_algorithm_ = mrpaAnalIter;
0043 
0044   if (multi_rp_algorithm_ == mrpaUndefined)
0045     throw cms::Exception("ProtonReconstructionAlgorithm") << "Algorithm '" << multiRPAlgorithm << "' not understood.";
0046 
0047   // initialise fitter
0048   double pStart[] = {0, 0, 0, 0};
0049   fitter_->SetFCN(4, *chiSquareCalculator_, pStart, 0, true);
0050 }
0051 
0052 //----------------------------------------------------------------------------------------------------
0053 
0054 void ProtonReconstructionAlgorithm::init(const LHCInterpolatedOpticalFunctionsSetCollection &opticalFunctions) {
0055   // reset cache
0056   release();
0057 
0058   // build optics data for each object
0059   for (const auto &p : opticalFunctions) {
0060     const LHCInterpolatedOpticalFunctionsSet &ofs = p.second;
0061 
0062     // make record
0063     RPOpticsData rpod;
0064     rpod.optics = &p.second;
0065     rpod.s_x_d_vs_xi = ofs.splines()[LHCOpticalFunctionsSet::exd];
0066     rpod.s_L_x_vs_xi = ofs.splines()[LHCOpticalFunctionsSet::eLx];
0067     rpod.s_y_d_vs_xi = ofs.splines()[LHCOpticalFunctionsSet::eyd];
0068     rpod.s_v_y_vs_xi = ofs.splines()[LHCOpticalFunctionsSet::evy];
0069     rpod.s_L_y_vs_xi = ofs.splines()[LHCOpticalFunctionsSet::eLy];
0070 
0071     vector<double> xiValues =
0072         ofs.getXiValues();  // local copy made since the TSpline constructor needs non-const parameters
0073     vector<double> xDValues = ofs.getFcnValues()[LHCOpticalFunctionsSet::exd];
0074     rpod.s_xi_vs_x_d = make_shared<TSpline3>("", xDValues.data(), xiValues.data(), xiValues.size());
0075 
0076     // calculate auxiliary data
0077     LHCInterpolatedOpticalFunctionsSet::Kinematics k_in = {0., 0., 0., 0., 0.};
0078     LHCInterpolatedOpticalFunctionsSet::Kinematics k_out;
0079     rpod.optics->transport(k_in, k_out);
0080     rpod.x0 = k_out.x;
0081     rpod.y0 = k_out.y;
0082 
0083     doLinearFit(ofs.getXiValues(), ofs.getFcnValues()[LHCOpticalFunctionsSet::exd], rpod.ch0, rpod.ch1);
0084     rpod.ch0 -= rpod.x0;
0085 
0086     doLinearFit(ofs.getXiValues(), ofs.getFcnValues()[LHCOpticalFunctionsSet::eLx], rpod.la0, rpod.la1);
0087 
0088     // insert record
0089     const CTPPSDetId rpId(p.first);
0090     m_rp_optics_.emplace(rpId, std::move(rpod));
0091   }
0092 
0093   // update settings
0094   initialized_ = true;
0095 }
0096 
0097 //----------------------------------------------------------------------------------------------------
0098 
0099 void ProtonReconstructionAlgorithm::doLinearFit(const std::vector<double> &vx,
0100                                                 const std::vector<double> &vy,
0101                                                 double &b,
0102                                                 double &a) {
0103   double s_1 = 0., s_x = 0., s_xx = 0., s_y = 0., s_xy = 0.;
0104   for (unsigned int i = 0; i < vx.size(); ++i) {
0105     s_1 += 1.;
0106     s_x += vx[i];
0107     s_xx += vx[i] * vx[i];
0108     s_y += vy[i];
0109     s_xy += vx[i] * vy[i];
0110   }
0111 
0112   const double d = s_xx * s_1 - s_x * s_x;
0113   a = (s_1 * s_xy - s_x * s_y) / d;
0114   b = (-s_x * s_xy + s_xx * s_y) / d;
0115 }
0116 
0117 //----------------------------------------------------------------------------------------------------
0118 
0119 void ProtonReconstructionAlgorithm::release() {
0120   initialized_ = false;
0121 
0122   m_rp_optics_.clear();
0123 }
0124 
0125 //----------------------------------------------------------------------------------------------------
0126 
0127 double ProtonReconstructionAlgorithm::ChiSquareCalculator::operator()(const double *parameters) const {
0128   // extract proton parameters
0129   const LHCInterpolatedOpticalFunctionsSet::Kinematics k_in = {
0130       0., parameters[1], parameters[3], parameters[2], parameters[0]};
0131 
0132   // calculate chi^2 by looping over hits
0133   double s2 = 0.;
0134 
0135   for (const auto &track : *tracks) {
0136     const CTPPSDetId rpId(track->rpId());
0137 
0138     // transport proton to the RP
0139     auto oit = m_rp_optics->find(rpId);
0140     LHCInterpolatedOpticalFunctionsSet::Kinematics k_out;
0141     oit->second.optics->transport(k_in, k_out);
0142 
0143     // proton position wrt. beam
0144     const double x = k_out.x - oit->second.x0;
0145     const double y = k_out.y - oit->second.y0;
0146 
0147     // calculate chi^2 contributions, convert track data mm --> cm
0148     const double x_diff_norm = (x - track->x() * 1E-1) / (track->xUnc() * 1E-1);
0149     const double y_diff_norm = (y - track->y() * 1E-1) / (track->yUnc() * 1E-1);
0150 
0151     // increase chi^2
0152     s2 += x_diff_norm * x_diff_norm + y_diff_norm * y_diff_norm;
0153   }
0154 
0155   return s2;
0156 }
0157 
0158 //----------------------------------------------------------------------------------------------------
0159 
0160 double ProtonReconstructionAlgorithm::newtonGoalFcn(double xi,
0161                                                     double x_N,
0162                                                     double x_F,
0163                                                     const ProtonReconstructionAlgorithm::RPOpticsData &i_N,
0164                                                     const ProtonReconstructionAlgorithm::RPOpticsData &i_F) {
0165   return (x_N - i_N.s_x_d_vs_xi->Eval(xi)) * i_F.s_L_x_vs_xi->Eval(xi) -
0166          (x_F - i_F.s_x_d_vs_xi->Eval(xi)) * i_N.s_L_x_vs_xi->Eval(xi);
0167 }
0168 
0169 //----------------------------------------------------------------------------------------------------
0170 
0171 reco::ForwardProton ProtonReconstructionAlgorithm::reconstructFromMultiRP(const CTPPSLocalTrackLiteRefVector &tracks,
0172                                                                           const float energy,
0173                                                                           std::ostream &os) const {
0174   // make sure optics is available for all tracks
0175   for (const auto &it : tracks) {
0176     auto oit = m_rp_optics_.find(it->rpId());
0177     if (oit == m_rp_optics_.end())
0178       throw cms::Exception("ProtonReconstructionAlgorithm")
0179           << "Optics data not available for RP " << it->rpId() << ", i.e. " << CTPPSDetId(it->rpId()) << ".";
0180   }
0181 
0182   // initial estimate of xi and th_x
0183   double xi_init = 0., th_x_init = 0.;
0184 
0185   if (useImprovedInitialEstimate_) {
0186     double x_N = tracks[0]->x() * 1E-1,  // conversion: mm --> cm
0187         x_F = tracks[1]->x() * 1E-1;
0188 
0189     const RPOpticsData &i_N = m_rp_optics_.find(tracks[0]->rpId())->second,
0190                        &i_F = m_rp_optics_.find(tracks[1]->rpId())->second;
0191 
0192     const double a = i_F.ch1 * i_N.la1 - i_N.ch1 * i_F.la1;
0193     const double b =
0194         i_F.ch0 * i_N.la1 - i_N.ch0 * i_F.la1 + i_F.ch1 * i_N.la0 - i_N.ch1 * i_F.la0 + x_N * i_F.la1 - x_F * i_N.la1;
0195     const double c = x_N * i_F.la0 - x_F * i_N.la0 + i_F.ch0 * i_N.la0 - i_N.ch0 * i_F.la0;
0196     const double D = b * b - 4. * a * c;
0197     const double sqrt_D = (D >= 0.) ? sqrt(D) : 0.;
0198 
0199     xi_init = (-b + sqrt_D) / 2. / a;
0200     th_x_init = (x_N - i_N.ch0 - i_N.ch1 * xi_init) / (i_N.la0 + i_N.la1 * xi_init);
0201   } else {
0202     double s_xi0 = 0., s_1 = 0.;
0203     for (const auto &track : tracks) {
0204       auto oit = m_rp_optics_.find(track->rpId());
0205       double xi = oit->second.s_xi_vs_x_d->Eval(track->x() * 1E-1 + oit->second.x0);  // conversion: mm --> cm
0206 
0207       s_1 += 1.;
0208       s_xi0 += xi;
0209     }
0210 
0211     xi_init = s_xi0 / s_1;
0212   }
0213 
0214   if (!std::isfinite(xi_init))
0215     xi_init = 0.;
0216   if (!std::isfinite(th_x_init))
0217     th_x_init = 0.;
0218 
0219   // initial estimate of th_y and vtx_y
0220   double y[2] = {0}, v_y[2] = {0}, L_y[2] = {0};
0221   unsigned int y_idx = 0;
0222   for (const auto &track : tracks) {
0223     if (y_idx >= 2)
0224       break;
0225 
0226     auto oit = m_rp_optics_.find(track->rpId());
0227 
0228     y[y_idx] = track->y() * 1E-1 - oit->second.s_y_d_vs_xi->Eval(xi_init);  // track y: mm --> cm
0229     v_y[y_idx] = oit->second.s_v_y_vs_xi->Eval(xi_init);
0230     L_y[y_idx] = oit->second.s_L_y_vs_xi->Eval(xi_init);
0231 
0232     y_idx++;
0233   }
0234 
0235   double vtx_y_init = 0.;
0236   double th_y_init = 0.;
0237 
0238   if (fitVtxY_) {
0239     const double det_y = v_y[0] * L_y[1] - L_y[0] * v_y[1];
0240     vtx_y_init = (det_y != 0.) ? (L_y[1] * y[0] - L_y[0] * y[1]) / det_y : 0.;
0241     th_y_init = (det_y != 0.) ? (v_y[0] * y[1] - v_y[1] * y[0]) / det_y : 0.;
0242   } else {
0243     vtx_y_init = 0.;
0244     th_y_init = (y[1] / L_y[1] + y[0] / L_y[0]) / 2.;
0245   }
0246 
0247   if (!std::isfinite(vtx_y_init))
0248     vtx_y_init = 0.;
0249   if (!std::isfinite(th_y_init))
0250     th_y_init = 0.;
0251 
0252   unsigned int armId = CTPPSDetId((*tracks.begin())->rpId()).arm();
0253 
0254   if (verbosity_)
0255     os << "ProtonReconstructionAlgorithm::reconstructFromMultiRP(" << armId << ")" << std::endl
0256        << "    initial estimate: xi_init = " << xi_init << ", th_x_init = " << th_x_init
0257        << ", th_y_init = " << th_y_init << ", vtx_y_init = " << vtx_y_init << "." << std::endl;
0258 
0259   // prepare result containers
0260   bool valid = false;
0261   double chi2 = 0.;
0262   double xi = 0., th_x = 0., th_y = 0., vtx_y = 0.;
0263 
0264   using FP = reco::ForwardProton;
0265   FP::CovarianceMatrix cm;
0266 
0267   if (multi_rp_algorithm_ == mrpaChi2) {
0268     // minimisation
0269     fitter_->Config().ParSettings(0).Set("xi", xi_init, 0.005);
0270     fitter_->Config().ParSettings(1).Set("th_x", th_x_init, 2E-6);
0271     fitter_->Config().ParSettings(2).Set("th_y", th_y_init, 2E-6);
0272     fitter_->Config().ParSettings(3).Set("vtx_y", vtx_y_init, 10E-6);
0273 
0274     if (!fitVtxY_)
0275       fitter_->Config().ParSettings(3).Fix();
0276 
0277     chiSquareCalculator_->tracks = &tracks;
0278     chiSquareCalculator_->m_rp_optics = &m_rp_optics_;
0279 
0280     fitter_->FitFCN();
0281     fitter_->FitFCN();  // second minimisation in case the first one had troubles
0282 
0283     // extract proton parameters
0284     const ROOT::Fit::FitResult &result = fitter_->Result();
0285     const double *params = result.GetParams();
0286 
0287     valid = result.IsValid();
0288     chi2 = result.Chi2();
0289     xi = params[0];
0290     th_x = params[1];
0291     th_y = params[2];
0292     vtx_y = params[3];
0293 
0294     map<unsigned int, signed int> index_map = {
0295         {(unsigned int)FP::Index::xi, 0},
0296         {(unsigned int)FP::Index::th_x, 1},
0297         {(unsigned int)FP::Index::th_y, 2},
0298         {(unsigned int)FP::Index::vtx_y, ((fitVtxY_) ? 3 : -1)},
0299         {(unsigned int)FP::Index::vtx_x, -1},
0300     };
0301 
0302     for (unsigned int i = 0; i < (unsigned int)FP::Index::num_indices; ++i) {
0303       signed int fit_i = index_map[i];
0304 
0305       for (unsigned int j = 0; j < (unsigned int)FP::Index::num_indices; ++j) {
0306         signed int fit_j = index_map[j];
0307 
0308         cm(i, j) = (fit_i >= 0 && fit_j >= 0) ? result.CovMatrix(fit_i, fit_j) : 0.;
0309       }
0310     }
0311   }
0312 
0313   else if (multi_rp_algorithm_ == mrpaNewton || multi_rp_algorithm_ == mrpaAnalIter) {
0314     // settings
0315     const unsigned int maxIterations = 100;
0316     const double maxXiDiff = 1E-6;
0317     const double xi_ep = 1E-5;
0318 
0319     // collect input
0320     double x_N = tracks[0]->x() * 1E-1,  // conversion: mm --> cm
0321         x_F = tracks[1]->x() * 1E-1;
0322 
0323     double y_N = tracks[0]->y() * 1E-1,  // conversion: mm --> cm
0324         y_F = tracks[1]->y() * 1E-1;
0325 
0326     const RPOpticsData &i_N = m_rp_optics_.find(tracks[0]->rpId())->second,
0327                        &i_F = m_rp_optics_.find(tracks[1]->rpId())->second;
0328 
0329     // horizontal reconstruction - run iterations
0330     valid = false;
0331     double xi_prev = xi_init;
0332     for (unsigned int it = 0; it < maxIterations; ++it) {
0333       if (multi_rp_algorithm_ == mrpaNewton) {
0334         const double g = newtonGoalFcn(xi_prev, x_N, x_F, i_N, i_F);
0335         const double gp = (newtonGoalFcn(xi_prev + xi_ep, x_N, x_F, i_N, i_F) - g) / xi_ep;
0336 
0337         xi = xi_prev - g / gp;
0338       }
0339 
0340       if (multi_rp_algorithm_ == mrpaAnalIter) {
0341         const double d_x_eff_N = i_N.s_x_d_vs_xi->Eval(xi_prev) / xi_prev;
0342         const double d_x_eff_F = i_F.s_x_d_vs_xi->Eval(xi_prev) / xi_prev;
0343 
0344         const double l_x_N = i_N.s_L_x_vs_xi->Eval(xi_prev);
0345         const double l_x_F = i_F.s_L_x_vs_xi->Eval(xi_prev);
0346 
0347         xi = (l_x_N * x_F - l_x_F * x_N) / (l_x_N * d_x_eff_F - l_x_F * d_x_eff_N);
0348       }
0349 
0350       if (abs(xi - xi_prev) < maxXiDiff) {
0351         valid = true;
0352         break;
0353       }
0354 
0355       xi_prev = xi;
0356     }
0357 
0358     th_x = (x_F - i_F.s_x_d_vs_xi->Eval(xi)) / i_F.s_L_x_vs_xi->Eval(xi);
0359 
0360     // vertical reconstruction
0361     const double y_eff_N = y_N - i_N.s_y_d_vs_xi->Eval(xi);
0362     const double y_eff_F = y_F - i_F.s_y_d_vs_xi->Eval(xi);
0363 
0364     const double l_y_N = i_N.s_L_y_vs_xi->Eval(xi);
0365     const double l_y_F = i_F.s_L_y_vs_xi->Eval(xi);
0366 
0367     const double v_y_N = i_N.s_v_y_vs_xi->Eval(xi);
0368     const double v_y_F = i_F.s_v_y_vs_xi->Eval(xi);
0369 
0370     const double det = l_y_N * v_y_F - l_y_F * v_y_N;
0371 
0372     if (fitVtxY_) {
0373       th_y = (y_eff_N * v_y_F - y_eff_F * v_y_N) / det;
0374       vtx_y = (l_y_N * y_eff_F - l_y_F * y_eff_N) / det;
0375     } else {
0376       th_y = (y_eff_N / l_y_N + y_eff_F / l_y_F) / 2.;
0377       vtx_y = 0.;
0378     }
0379 
0380     // covariance matrix kept empty - not to compromise the speed of these special algorithms
0381   }
0382 
0383   // print results
0384   if (verbosity_)
0385     os << "    fit valid=" << valid << std::endl
0386        << "    xi=" << xi << ", th_x=" << th_x << ", th_y=" << th_y << ", vtx_y=" << vtx_y << ", chiSq = " << chi2
0387        << std::endl;
0388 
0389   // save reco candidate
0390   const double sign_z = (armId == 0) ? +1. : -1.;  // CMS convention
0391   const FP::Point vertex(0., vtx_y, 0.);
0392   const double cos_th_sq = 1. - th_x * th_x - th_y * th_y;
0393   const double cos_th = (cos_th_sq > 0.) ? sqrt(cos_th_sq) : 1.;
0394   const double p = energy * (1. - xi);
0395   const FP::Vector momentum(-p * th_x,  // the signs reflect change LHC --> CMS convention
0396                             +p * th_y,
0397                             sign_z * p * cos_th);
0398   signed int ndf = 2. * tracks.size() - ((fitVtxY_) ? 4. : 3.);
0399 
0400   return reco::ForwardProton(chi2, ndf, vertex, momentum, xi, cm, FP::ReconstructionMethod::multiRP, tracks, valid);
0401 }
0402 
0403 //----------------------------------------------------------------------------------------------------
0404 
0405 reco::ForwardProton ProtonReconstructionAlgorithm::reconstructFromSingleRP(const CTPPSLocalTrackLiteRef &track,
0406                                                                            const float energy,
0407                                                                            std::ostream &os) const {
0408   CTPPSDetId rpId(track->rpId());
0409 
0410   if (verbosity_)
0411     os << "reconstructFromSingleRP(" << rpId.arm() * 100 + rpId.station() * 10 + rpId.rp() << ")" << std::endl;
0412 
0413   // make sure optics is available for the track
0414   auto oit = m_rp_optics_.find(track->rpId());
0415   if (oit == m_rp_optics_.end())
0416     throw cms::Exception("ProtonReconstructionAlgorithm")
0417         << "Optics data not available for RP " << track->rpId() << ", i.e. " << rpId << ".";
0418 
0419   // rough estimate of xi and th_y from each track
0420   const double x_full = track->x() * 1E-1 + oit->second.x0;  // conversion mm --> cm
0421   const double xi = oit->second.s_xi_vs_x_d->Eval(x_full);
0422   const double L_y = oit->second.s_L_y_vs_xi->Eval(xi);
0423   const double th_y = track->y() * 1E-1 / L_y;  // conversion mm --> cm
0424 
0425   const double ep_x = 1E-6;
0426   const double dxi_dx = (oit->second.s_xi_vs_x_d->Eval(x_full + ep_x) - xi) / ep_x;
0427   const double xi_unc = abs(dxi_dx) * track->xUnc() * 1E-1;  // conversion mm --> cm
0428 
0429   const double ep_xi = 1E-4;
0430   const double dL_y_dxi = (oit->second.s_L_y_vs_xi->Eval(xi + ep_xi) - L_y) / ep_xi;
0431   const double th_y_unc_sq = th_y * th_y * (pow(track->yUnc() / track->y(), 2.) + pow(dL_y_dxi * xi_unc / L_y, 2.));
0432 
0433   if (verbosity_)
0434     os << "    xi = " << xi << " +- " << xi_unc << ", th_y = " << th_y << " +- " << sqrt(th_y_unc_sq) << "."
0435        << std::endl;
0436 
0437   using FP = reco::ForwardProton;
0438 
0439   // save proton candidate
0440   const double sign_z = (CTPPSDetId(track->rpId()).arm() == 0) ? +1. : -1.;  // CMS convention
0441   const FP::Point vertex(0., 0., 0.);
0442   const double cos_th_sq = 1. - th_y * th_y;
0443   const double cos_th = (cos_th_sq > 0.) ? sqrt(cos_th_sq) : 1.;
0444   const double p = energy * (1. - xi);
0445   const FP::Vector momentum(0., p * th_y, sign_z * p * cos_th);
0446 
0447   FP::CovarianceMatrix cm;
0448   cm((int)FP::Index::xi, (int)FP::Index::xi) = xi_unc * xi_unc;
0449   cm((int)FP::Index::th_y, (int)FP::Index::th_y) = th_y_unc_sq;
0450 
0451   CTPPSLocalTrackLiteRefVector trk;
0452   trk.push_back(track);
0453 
0454   return reco::ForwardProton(0., 0, vertex, momentum, xi, cm, FP::ReconstructionMethod::singleRP, trk, true);
0455 }
[ 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!