CMS 3D CMS Logo

List of all members | Classes | Public Member Functions | Private Types | Static Private Member Functions | Private Attributes
ProtonReconstructionAlgorithm Class Reference

#include <ProtonReconstructionAlgorithm.h>

Classes

class  ChiSquareCalculator
 class for calculation of chi^2 More...
 
struct  RPOpticsData
 optics data associated with 1 RP More...
 

Public Member Functions

void init (const LHCInterpolatedOpticalFunctionsSetCollection &opticalFunctions)
 
 ProtonReconstructionAlgorithm (bool fit_vtx_y, bool improved_estimate, const std::string &multiRPAlgorithm, unsigned int verbosity)
 
reco::ForwardProton reconstructFromMultiRP (const CTPPSLocalTrackLiteRefVector &tracks, const float energy, std::ostream &os) const
 run proton reconstruction using multiple-RP strategy More...
 
reco::ForwardProton reconstructFromSingleRP (const CTPPSLocalTrackLiteRef &track, const float energy, std::ostream &os) const
 run proton reconstruction using single-RP strategy More...
 
void release ()
 
 ~ProtonReconstructionAlgorithm ()=default
 

Private Types

enum  { mrpaUndefined, mrpaChi2, mrpaNewton, mrpaAnalIter }
 

Static Private Member Functions

static void doLinearFit (const std::vector< double > &vx, const std::vector< double > &vy, double &b, double &a)
 
static double newtonGoalFcn (double xi, double x_N, double x_F, const RPOpticsData &i_N, const RPOpticsData &i_F)
 

Private Attributes

std::unique_ptr< ChiSquareCalculatorchiSquareCalculator_
 object to calculate chi^2 More...
 
std::unique_ptr< ROOT::Fit::Fitter > fitter_
 fitter object More...
 
bool fitVtxY_
 
bool initialized_
 
std::map< unsigned int, RPOpticsDatam_rp_optics_
 map: RP id –> optics data More...
 
enum ProtonReconstructionAlgorithm:: { ... }  multi_rp_algorithm_
 
bool useImprovedInitialEstimate_
 
unsigned int verbosity_
 

Detailed Description

Definition at line 25 of file ProtonReconstructionAlgorithm.h.

Member Enumeration Documentation

◆ anonymous enum

anonymous enum
private

Constructor & Destructor Documentation

◆ ProtonReconstructionAlgorithm()

ProtonReconstructionAlgorithm::ProtonReconstructionAlgorithm ( bool  fit_vtx_y,
bool  improved_estimate,
const std::string &  multiRPAlgorithm,
unsigned int  verbosity 
)

Definition at line 22 of file ProtonReconstructionAlgorithm.cc.

References chiSquareCalculator_, Exception, fitter_, mrpaAnalIter, mrpaChi2, mrpaNewton, mrpaUndefined, and multi_rp_algorithm_.

27  fitVtxY_(fit_vtx_y),
28  useImprovedInitialEstimate_(improved_estimate),
30  initialized_(false),
32  chiSquareCalculator_(new ChiSquareCalculator) {
33  // needed for thread safety
34  TMinuitMinimizer::UseStaticMinuit(false);
35 
36  // check and set multi-RP algorithm
37  if (multiRPAlgorithm == "chi2")
39  else if (multiRPAlgorithm == "newton")
41  else if (multiRPAlgorithm == "anal-iter")
43 
45  throw cms::Exception("ProtonReconstructionAlgorithm") << "Algorithm '" << multiRPAlgorithm << "' not understood.";
46 
47  // initialise fitter
48  double pStart[] = {0, 0, 0, 0};
49  fitter_->SetFCN(4, *chiSquareCalculator_, pStart, 0, true);
50 }
std::unique_ptr< ROOT::Fit::Fitter > fitter_
fitter object
enum ProtonReconstructionAlgorithm::@992 multi_rp_algorithm_
std::unique_ptr< ChiSquareCalculator > chiSquareCalculator_
object to calculate chi^2
const int verbosity

◆ ~ProtonReconstructionAlgorithm()

ProtonReconstructionAlgorithm::~ProtonReconstructionAlgorithm ( )
default

Member Function Documentation

◆ doLinearFit()

void ProtonReconstructionAlgorithm::doLinearFit ( const std::vector< double > &  vx,
const std::vector< double > &  vy,
double &  b,
double &  a 
)
staticprivate

Definition at line 99 of file ProtonReconstructionAlgorithm.cc.

References a, b, ztail::d, mps_fire::i, run3scouting_cff::vx, and run3scouting_cff::vy.

Referenced by init().

102  {
103  double s_1 = 0., s_x = 0., s_xx = 0., s_y = 0., s_xy = 0.;
104  for (unsigned int i = 0; i < vx.size(); ++i) {
105  s_1 += 1.;
106  s_x += vx[i];
107  s_xx += vx[i] * vx[i];
108  s_y += vy[i];
109  s_xy += vx[i] * vy[i];
110  }
111 
112  const double d = s_xx * s_1 - s_x * s_x;
113  a = (s_1 * s_xy - s_x * s_y) / d;
114  b = (-s_x * s_xy + s_xx * s_y) / d;
115 }
d
Definition: ztail.py:151
double b
Definition: hdecay.h:120
double a
Definition: hdecay.h:121

◆ init()

void ProtonReconstructionAlgorithm::init ( const LHCInterpolatedOpticalFunctionsSetCollection opticalFunctions)

Definition at line 54 of file ProtonReconstructionAlgorithm.cc.

References ProtonReconstructionAlgorithm::RPOpticsData::ch0, ProtonReconstructionAlgorithm::RPOpticsData::ch1, doLinearFit(), LHCOpticalFunctionsSet::eLx, LHCOpticalFunctionsSet::eLy, LHCOpticalFunctionsSet::evy, LHCOpticalFunctionsSet::exd, LHCOpticalFunctionsSet::eyd, LHCOpticalFunctionsSet::getFcnValues(), LHCOpticalFunctionsSet::getXiValues(), initialized_, ProtonReconstructionAlgorithm::RPOpticsData::la0, ProtonReconstructionAlgorithm::RPOpticsData::la1, m_rp_optics_, eostools::move(), profile_base_cff::opticalFunctions, ProtonReconstructionAlgorithm::RPOpticsData::optics, AlCaHLTBitMon_ParallelJobs::p, release(), ProtonReconstructionAlgorithm::RPOpticsData::s_L_x_vs_xi, ProtonReconstructionAlgorithm::RPOpticsData::s_L_y_vs_xi, ProtonReconstructionAlgorithm::RPOpticsData::s_v_y_vs_xi, ProtonReconstructionAlgorithm::RPOpticsData::s_x_d_vs_xi, ProtonReconstructionAlgorithm::RPOpticsData::s_xi_vs_x_d, ProtonReconstructionAlgorithm::RPOpticsData::s_y_d_vs_xi, LHCInterpolatedOpticalFunctionsSet::splines(), LHCInterpolatedOpticalFunctionsSet::transport(), LHCInterpolatedOpticalFunctionsSet::Kinematics::x, ProtonReconstructionAlgorithm::RPOpticsData::x0, LHCInterpolatedOpticalFunctionsSet::Kinematics::y, and ProtonReconstructionAlgorithm::RPOpticsData::y0.

Referenced by CTPPSProtonProducer::produce().

54  {
55  // reset cache
56  release();
57 
58  // build optics data for each object
59  for (const auto &p : opticalFunctions) {
60  const LHCInterpolatedOpticalFunctionsSet &ofs = p.second;
61 
62  // make record
63  RPOpticsData rpod;
64  rpod.optics = &p.second;
65  rpod.s_x_d_vs_xi = ofs.splines()[LHCOpticalFunctionsSet::exd];
66  rpod.s_L_x_vs_xi = ofs.splines()[LHCOpticalFunctionsSet::eLx];
67  rpod.s_y_d_vs_xi = ofs.splines()[LHCOpticalFunctionsSet::eyd];
68  rpod.s_v_y_vs_xi = ofs.splines()[LHCOpticalFunctionsSet::evy];
69  rpod.s_L_y_vs_xi = ofs.splines()[LHCOpticalFunctionsSet::eLy];
70 
71  vector<double> xiValues =
72  ofs.getXiValues(); // local copy made since the TSpline constructor needs non-const parameters
73  vector<double> xDValues = ofs.getFcnValues()[LHCOpticalFunctionsSet::exd];
74  rpod.s_xi_vs_x_d = make_shared<TSpline3>("", xDValues.data(), xiValues.data(), xiValues.size());
75 
76  // calculate auxiliary data
77  LHCInterpolatedOpticalFunctionsSet::Kinematics k_in = {0., 0., 0., 0., 0.};
79  rpod.optics->transport(k_in, k_out);
80  rpod.x0 = k_out.x;
81  rpod.y0 = k_out.y;
82 
83  doLinearFit(ofs.getXiValues(), ofs.getFcnValues()[LHCOpticalFunctionsSet::exd], rpod.ch0, rpod.ch1);
84  rpod.ch0 -= rpod.x0;
85 
86  doLinearFit(ofs.getXiValues(), ofs.getFcnValues()[LHCOpticalFunctionsSet::eLx], rpod.la0, rpod.la1);
87 
88  // insert record
89  const CTPPSDetId rpId(p.first);
90  m_rp_optics_.emplace(rpId, std::move(rpod));
91  }
92 
93  // update settings
94  initialized_ = true;
95 }
std::map< unsigned int, RPOpticsData > m_rp_optics_
map: RP id –> optics data
static void doLinearFit(const std::vector< double > &vx, const std::vector< double > &vy, double &b, double &a)
const std::vector< std::shared_ptr< const TSpline3 > > & splines() const
const std::vector< std::vector< double > > & getFcnValues() const
Set of optical functions corresponding to one scoring plane along LHC, including splines for interpol...
const std::vector< double > & getXiValues() const
Base class for CTPPS detector IDs.
Definition: CTPPSDetId.h:32
def move(src, dest)
Definition: eostools.py:511

◆ newtonGoalFcn()

double ProtonReconstructionAlgorithm::newtonGoalFcn ( double  xi,
double  x_N,
double  x_F,
const RPOpticsData i_N,
const RPOpticsData i_F 
)
staticprivate

Definition at line 160 of file ProtonReconstructionAlgorithm.cc.

References ProtonReconstructionAlgorithm::RPOpticsData::s_L_x_vs_xi, ProtonReconstructionAlgorithm::RPOpticsData::s_x_d_vs_xi, and protons_cff::xi.

Referenced by reconstructFromMultiRP().

164  {
165  return (x_N - i_N.s_x_d_vs_xi->Eval(xi)) * i_F.s_L_x_vs_xi->Eval(xi) -
166  (x_F - i_F.s_x_d_vs_xi->Eval(xi)) * i_N.s_L_x_vs_xi->Eval(xi);
167 }

◆ reconstructFromMultiRP()

reco::ForwardProton ProtonReconstructionAlgorithm::reconstructFromMultiRP ( const CTPPSLocalTrackLiteRefVector tracks,
const float  energy,
std::ostream &  os 
) const

run proton reconstruction using multiple-RP strategy

Definition at line 171 of file ProtonReconstructionAlgorithm.cc.

References a, funct::abs(), protons_cff::arm, b, HltBtagPostValidation_cff::c, ProtonReconstructionAlgorithm::RPOpticsData::ch0, ProtonReconstructionAlgorithm::RPOpticsData::ch1, nano_mu_local_reco_cff::chi2, chiSquareCalculator_, hcalRecHitTable_cff::energy, fitter_, fitVtxY_, g, runTauDisplay::gp, mps_fire::i, createfilelist::int, dqmiolumiharvest::j, ProtonReconstructionAlgorithm::RPOpticsData::la0, ProtonReconstructionAlgorithm::RPOpticsData::la1, m_rp_optics_, HLT_2023v12_cff::maxIterations, mrpaAnalIter, mrpaChi2, mrpaNewton, multi_rp_algorithm_, newtonGoalFcn(), AlCaHLTBitMon_ParallelJobs::p, submitPVValidationJobs::params, mps_fire::result, ProtonReconstructionAlgorithm::RPOpticsData::s_L_x_vs_xi, ProtonReconstructionAlgorithm::RPOpticsData::s_L_y_vs_xi, ProtonReconstructionAlgorithm::RPOpticsData::s_v_y_vs_xi, ProtonReconstructionAlgorithm::RPOpticsData::s_x_d_vs_xi, ProtonReconstructionAlgorithm::RPOpticsData::s_y_d_vs_xi, mathSSE::sqrt(), HLT_2023v12_cff::track, pwdgSkimBPark_cfi::tracks, useImprovedInitialEstimate_, validateGeometry_cfg::valid, verbosity_, bphysicsOniaDQM_cfi::vertex, protons_cff::xi, and y.

Referenced by CTPPSProtonProducer::produce().

173  {
174  // make sure optics is available for all tracks
175  for (const auto &it : tracks) {
176  auto oit = m_rp_optics_.find(it->rpId());
177  if (oit == m_rp_optics_.end())
178  throw cms::Exception("ProtonReconstructionAlgorithm")
179  << "Optics data not available for RP " << it->rpId() << ", i.e. " << CTPPSDetId(it->rpId()) << ".";
180  }
181 
182  // initial estimate of xi and th_x
183  double xi_init = 0., th_x_init = 0.;
184 
186  double x_N = tracks[0]->x() * 1E-1, // conversion: mm --> cm
187  x_F = tracks[1]->x() * 1E-1;
188 
189  const RPOpticsData &i_N = m_rp_optics_.find(tracks[0]->rpId())->second,
190  &i_F = m_rp_optics_.find(tracks[1]->rpId())->second;
191 
192  const double a = i_F.ch1 * i_N.la1 - i_N.ch1 * i_F.la1;
193  const double b =
194  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;
195  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;
196  const double D = b * b - 4. * a * c;
197  const double sqrt_D = (D >= 0.) ? sqrt(D) : 0.;
198 
199  xi_init = (-b + sqrt_D) / 2. / a;
200  th_x_init = (x_N - i_N.ch0 - i_N.ch1 * xi_init) / (i_N.la0 + i_N.la1 * xi_init);
201  } else {
202  double s_xi0 = 0., s_1 = 0.;
203  for (const auto &track : tracks) {
204  auto oit = m_rp_optics_.find(track->rpId());
205  double xi = oit->second.s_xi_vs_x_d->Eval(track->x() * 1E-1 + oit->second.x0); // conversion: mm --> cm
206 
207  s_1 += 1.;
208  s_xi0 += xi;
209  }
210 
211  xi_init = s_xi0 / s_1;
212  }
213 
214  if (!std::isfinite(xi_init))
215  xi_init = 0.;
216  if (!std::isfinite(th_x_init))
217  th_x_init = 0.;
218 
219  // initial estimate of th_y and vtx_y
220  double y[2] = {0}, v_y[2] = {0}, L_y[2] = {0};
221  unsigned int y_idx = 0;
222  for (const auto &track : tracks) {
223  if (y_idx >= 2)
224  break;
225 
226  auto oit = m_rp_optics_.find(track->rpId());
227 
228  y[y_idx] = track->y() * 1E-1 - oit->second.s_y_d_vs_xi->Eval(xi_init); // track y: mm --> cm
229  v_y[y_idx] = oit->second.s_v_y_vs_xi->Eval(xi_init);
230  L_y[y_idx] = oit->second.s_L_y_vs_xi->Eval(xi_init);
231 
232  y_idx++;
233  }
234 
235  double vtx_y_init = 0.;
236  double th_y_init = 0.;
237 
238  if (fitVtxY_) {
239  const double det_y = v_y[0] * L_y[1] - L_y[0] * v_y[1];
240  vtx_y_init = (det_y != 0.) ? (L_y[1] * y[0] - L_y[0] * y[1]) / det_y : 0.;
241  th_y_init = (det_y != 0.) ? (v_y[0] * y[1] - v_y[1] * y[0]) / det_y : 0.;
242  } else {
243  vtx_y_init = 0.;
244  th_y_init = (y[1] / L_y[1] + y[0] / L_y[0]) / 2.;
245  }
246 
247  if (!std::isfinite(vtx_y_init))
248  vtx_y_init = 0.;
249  if (!std::isfinite(th_y_init))
250  th_y_init = 0.;
251 
252  unsigned int armId = CTPPSDetId((*tracks.begin())->rpId()).arm();
253 
254  if (verbosity_)
255  os << "ProtonReconstructionAlgorithm::reconstructFromMultiRP(" << armId << ")" << std::endl
256  << " initial estimate: xi_init = " << xi_init << ", th_x_init = " << th_x_init
257  << ", th_y_init = " << th_y_init << ", vtx_y_init = " << vtx_y_init << "." << std::endl;
258 
259  // prepare result containers
260  bool valid = false;
261  double chi2 = 0.;
262  double xi = 0., th_x = 0., th_y = 0., vtx_y = 0.;
263 
264  using FP = reco::ForwardProton;
266 
267  if (multi_rp_algorithm_ == mrpaChi2) {
268  // minimisation
269  fitter_->Config().ParSettings(0).Set("xi", xi_init, 0.005);
270  fitter_->Config().ParSettings(1).Set("th_x", th_x_init, 2E-6);
271  fitter_->Config().ParSettings(2).Set("th_y", th_y_init, 2E-6);
272  fitter_->Config().ParSettings(3).Set("vtx_y", vtx_y_init, 10E-6);
273 
274  if (!fitVtxY_)
275  fitter_->Config().ParSettings(3).Fix();
276 
277  chiSquareCalculator_->tracks = &tracks;
278  chiSquareCalculator_->m_rp_optics = &m_rp_optics_;
279 
280  fitter_->FitFCN();
281  fitter_->FitFCN(); // second minimisation in case the first one had troubles
282 
283  // extract proton parameters
284  const ROOT::Fit::FitResult &result = fitter_->Result();
285  const double *params = result.GetParams();
286 
287  valid = result.IsValid();
288  chi2 = result.Chi2();
289  xi = params[0];
290  th_x = params[1];
291  th_y = params[2];
292  vtx_y = params[3];
293 
294  map<unsigned int, signed int> index_map = {
295  {(unsigned int)FP::Index::xi, 0},
296  {(unsigned int)FP::Index::th_x, 1},
297  {(unsigned int)FP::Index::th_y, 2},
298  {(unsigned int)FP::Index::vtx_y, ((fitVtxY_) ? 3 : -1)},
299  {(unsigned int)FP::Index::vtx_x, -1},
300  };
301 
302  for (unsigned int i = 0; i < (unsigned int)FP::Index::num_indices; ++i) {
303  signed int fit_i = index_map[i];
304 
305  for (unsigned int j = 0; j < (unsigned int)FP::Index::num_indices; ++j) {
306  signed int fit_j = index_map[j];
307 
308  cm(i, j) = (fit_i >= 0 && fit_j >= 0) ? result.CovMatrix(fit_i, fit_j) : 0.;
309  }
310  }
311  }
312 
314  // settings
315  const unsigned int maxIterations = 100;
316  const double maxXiDiff = 1E-6;
317  const double xi_ep = 1E-5;
318 
319  // collect input
320  double x_N = tracks[0]->x() * 1E-1, // conversion: mm --> cm
321  x_F = tracks[1]->x() * 1E-1;
322 
323  double y_N = tracks[0]->y() * 1E-1, // conversion: mm --> cm
324  y_F = tracks[1]->y() * 1E-1;
325 
326  const RPOpticsData &i_N = m_rp_optics_.find(tracks[0]->rpId())->second,
327  &i_F = m_rp_optics_.find(tracks[1]->rpId())->second;
328 
329  // horizontal reconstruction - run iterations
330  valid = false;
331  double xi_prev = xi_init;
332  for (unsigned int it = 0; it < maxIterations; ++it) {
334  const double g = newtonGoalFcn(xi_prev, x_N, x_F, i_N, i_F);
335  const double gp = (newtonGoalFcn(xi_prev + xi_ep, x_N, x_F, i_N, i_F) - g) / xi_ep;
336 
337  xi = xi_prev - g / gp;
338  }
339 
341  const double d_x_eff_N = i_N.s_x_d_vs_xi->Eval(xi_prev) / xi_prev;
342  const double d_x_eff_F = i_F.s_x_d_vs_xi->Eval(xi_prev) / xi_prev;
343 
344  const double l_x_N = i_N.s_L_x_vs_xi->Eval(xi_prev);
345  const double l_x_F = i_F.s_L_x_vs_xi->Eval(xi_prev);
346 
347  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);
348  }
349 
350  if (abs(xi - xi_prev) < maxXiDiff) {
351  valid = true;
352  break;
353  }
354 
355  xi_prev = xi;
356  }
357 
358  th_x = (x_F - i_F.s_x_d_vs_xi->Eval(xi)) / i_F.s_L_x_vs_xi->Eval(xi);
359 
360  // vertical reconstruction
361  const double y_eff_N = y_N - i_N.s_y_d_vs_xi->Eval(xi);
362  const double y_eff_F = y_F - i_F.s_y_d_vs_xi->Eval(xi);
363 
364  const double l_y_N = i_N.s_L_y_vs_xi->Eval(xi);
365  const double l_y_F = i_F.s_L_y_vs_xi->Eval(xi);
366 
367  const double v_y_N = i_N.s_v_y_vs_xi->Eval(xi);
368  const double v_y_F = i_F.s_v_y_vs_xi->Eval(xi);
369 
370  const double det = l_y_N * v_y_F - l_y_F * v_y_N;
371 
372  if (fitVtxY_) {
373  th_y = (y_eff_N * v_y_F - y_eff_F * v_y_N) / det;
374  vtx_y = (l_y_N * y_eff_F - l_y_F * y_eff_N) / det;
375  } else {
376  th_y = (y_eff_N / l_y_N + y_eff_F / l_y_F) / 2.;
377  vtx_y = 0.;
378  }
379 
380  // covariance matrix kept empty - not to compromise the speed of these special algorithms
381  }
382 
383  // print results
384  if (verbosity_)
385  os << " fit valid=" << valid << std::endl
386  << " xi=" << xi << ", th_x=" << th_x << ", th_y=" << th_y << ", vtx_y=" << vtx_y << ", chiSq = " << chi2
387  << std::endl;
388 
389  // save reco candidate
390  const double sign_z = (armId == 0) ? +1. : -1.; // CMS convention
391  const FP::Point vertex(0., vtx_y, 0.);
392  const double cos_th_sq = 1. - th_x * th_x - th_y * th_y;
393  const double cos_th = (cos_th_sq > 0.) ? sqrt(cos_th_sq) : 1.;
394  const double p = energy * (1. - xi);
395  const FP::Vector momentum(-p * th_x, // the signs reflect change LHC --> CMS convention
396  +p * th_y,
397  sign_z * p * cos_th);
398  signed int ndf = 2. * tracks.size() - ((fitVtxY_) ? 4. : 3.);
399 
400  return reco::ForwardProton(chi2, ndf, vertex, momentum, xi, cm, FP::ReconstructionMethod::multiRP, tracks, valid);
401 }
std::unique_ptr< ROOT::Fit::Fitter > fitter_
fitter object
enum ProtonReconstructionAlgorithm::@992 multi_rp_algorithm_
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
std::map< unsigned int, RPOpticsData > m_rp_optics_
map: RP id –> optics data
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
std::unique_ptr< ChiSquareCalculator > chiSquareCalculator_
object to calculate chi^2
T sqrt(T t)
Definition: SSEVec.h:19
static double newtonGoalFcn(double xi, double x_N, double x_F, const RPOpticsData &i_N, const RPOpticsData &i_F)
math::XYZPoint Point
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:141
double b
Definition: hdecay.h:120
Base class for CTPPS detector IDs.
Definition: CTPPSDetId.h:32
double a
Definition: hdecay.h:121
Eigen::Matrix< float, EcalPulseShape::TEMPLATESAMPLES, EcalPulseShape::TEMPLATESAMPLES > CovarianceMatrix

◆ reconstructFromSingleRP()

reco::ForwardProton ProtonReconstructionAlgorithm::reconstructFromSingleRP ( const CTPPSLocalTrackLiteRef track,
const float  energy,
std::ostream &  os 
) const

run proton reconstruction using single-RP strategy

Definition at line 405 of file ProtonReconstructionAlgorithm.cc.

References funct::abs(), protons_cff::arm, CTPPSDetId::arm(), hcalRecHitTable_cff::energy, m_rp_optics_, AlCaHLTBitMon_ParallelJobs::p, conifer::pow(), edm::RefVector< C, T, F >::push_back(), mathSSE::sqrt(), HLT_2023v12_cff::track, verbosity_, bphysicsOniaDQM_cfi::vertex, and protons_cff::xi.

Referenced by CTPPSProtonProducer::produce().

407  {
408  CTPPSDetId rpId(track->rpId());
409 
410  if (verbosity_)
411  os << "reconstructFromSingleRP(" << rpId.arm() * 100 + rpId.station() * 10 + rpId.rp() << ")" << std::endl;
412 
413  // make sure optics is available for the track
414  auto oit = m_rp_optics_.find(track->rpId());
415  if (oit == m_rp_optics_.end())
416  throw cms::Exception("ProtonReconstructionAlgorithm")
417  << "Optics data not available for RP " << track->rpId() << ", i.e. " << rpId << ".";
418 
419  // rough estimate of xi and th_y from each track
420  const double x_full = track->x() * 1E-1 + oit->second.x0; // conversion mm --> cm
421  const double xi = oit->second.s_xi_vs_x_d->Eval(x_full);
422  const double L_y = oit->second.s_L_y_vs_xi->Eval(xi);
423  const double th_y = track->y() * 1E-1 / L_y; // conversion mm --> cm
424 
425  const double ep_x = 1E-6;
426  const double dxi_dx = (oit->second.s_xi_vs_x_d->Eval(x_full + ep_x) - xi) / ep_x;
427  const double xi_unc = abs(dxi_dx) * track->xUnc() * 1E-1; // conversion mm --> cm
428 
429  const double ep_xi = 1E-4;
430  const double dL_y_dxi = (oit->second.s_L_y_vs_xi->Eval(xi + ep_xi) - L_y) / ep_xi;
431  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.));
432 
433  if (verbosity_)
434  os << " xi = " << xi << " +- " << xi_unc << ", th_y = " << th_y << " +- " << sqrt(th_y_unc_sq) << "."
435  << std::endl;
436 
437  using FP = reco::ForwardProton;
438 
439  // save proton candidate
440  const double sign_z = (CTPPSDetId(track->rpId()).arm() == 0) ? +1. : -1.; // CMS convention
441  const FP::Point vertex(0., 0., 0.);
442  const double cos_th_sq = 1. - th_y * th_y;
443  const double cos_th = (cos_th_sq > 0.) ? sqrt(cos_th_sq) : 1.;
444  const double p = energy * (1. - xi);
445  const FP::Vector momentum(0., p * th_y, sign_z * p * cos_th);
446 
448  cm((int)FP::Index::xi, (int)FP::Index::xi) = xi_unc * xi_unc;
449  cm((int)FP::Index::th_y, (int)FP::Index::th_y) = th_y_unc_sq;
450 
452  trk.push_back(track);
453 
454  return reco::ForwardProton(0., 0, vertex, momentum, xi, cm, FP::ReconstructionMethod::singleRP, trk, true);
455 }
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
uint32_t arm() const
Definition: CTPPSDetId.h:57
constexpr int pow(int x)
Definition: conifer.h:24
std::map< unsigned int, RPOpticsData > m_rp_optics_
map: RP id –> optics data
T sqrt(T t)
Definition: SSEVec.h:19
math::XYZPoint Point
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
Base class for CTPPS detector IDs.
Definition: CTPPSDetId.h:32
Eigen::Matrix< float, EcalPulseShape::TEMPLATESAMPLES, EcalPulseShape::TEMPLATESAMPLES > CovarianceMatrix
void push_back(value_type const &ref)
Add a Ref<C, T> to the RefVector.
Definition: RefVector.h:67

◆ release()

void ProtonReconstructionAlgorithm::release ( )

Definition at line 119 of file ProtonReconstructionAlgorithm.cc.

References initialized_, and m_rp_optics_.

Referenced by init(), and CTPPSProtonProducer::produce().

119  {
120  initialized_ = false;
121 
122  m_rp_optics_.clear();
123 }
std::map< unsigned int, RPOpticsData > m_rp_optics_
map: RP id –> optics data

Member Data Documentation

◆ chiSquareCalculator_

std::unique_ptr<ChiSquareCalculator> ProtonReconstructionAlgorithm::chiSquareCalculator_
private

object to calculate chi^2

Definition at line 83 of file ProtonReconstructionAlgorithm.h.

Referenced by ProtonReconstructionAlgorithm(), and reconstructFromMultiRP().

◆ fitter_

std::unique_ptr<ROOT::Fit::Fitter> ProtonReconstructionAlgorithm::fitter_
private

fitter object

Definition at line 80 of file ProtonReconstructionAlgorithm.h.

Referenced by ProtonReconstructionAlgorithm(), and reconstructFromMultiRP().

◆ fitVtxY_

bool ProtonReconstructionAlgorithm::fitVtxY_
private

Definition at line 48 of file ProtonReconstructionAlgorithm.h.

Referenced by reconstructFromMultiRP().

◆ initialized_

bool ProtonReconstructionAlgorithm::initialized_
private

Definition at line 51 of file ProtonReconstructionAlgorithm.h.

Referenced by init(), and release().

◆ m_rp_optics_

std::map<unsigned int, RPOpticsData> ProtonReconstructionAlgorithm::m_rp_optics_
private

map: RP id –> optics data

Definition at line 66 of file ProtonReconstructionAlgorithm.h.

Referenced by init(), reconstructFromMultiRP(), reconstructFromSingleRP(), and release().

◆ multi_rp_algorithm_

enum { ... } ProtonReconstructionAlgorithm::multi_rp_algorithm_

◆ useImprovedInitialEstimate_

bool ProtonReconstructionAlgorithm::useImprovedInitialEstimate_
private

Definition at line 49 of file ProtonReconstructionAlgorithm.h.

Referenced by reconstructFromMultiRP().

◆ verbosity_

unsigned int ProtonReconstructionAlgorithm::verbosity_
private