ATLAS Offline Software
MuonSegmentMomentum.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 #include "MuonSegmentMomentum.h"
6 #include "MuonSegment/MuonSegment.h"
7 #include <sstream>
8 #include <iostream>
9 #include <vector>
10 
11 MuonSegmentMomentum::MuonSegmentMomentum(const std::string& type,const std::string& name,const IInterface* parent)
12  :
13  AthAlgTool(type,name,parent)
14 {
15  declareInterface<IMuonSegmentMomentumEstimator>(this);
16 }
17 
18 void MuonSegmentMomentum::fitMomentumVectorSegments( const EventContext& ctx, const std::vector <const Muon::MuonSegment*> & segments, double & signedMomentum ) const
19 {
20 
24 // MsgStream log(messageService(),name());
25  ATH_MSG_VERBOSE(" Executing MuonSegmentMomentumTool fitMomentumVectorSegments ");
26  ATH_MSG_DEBUG(" fitMomentumVectorSegments " << segments.size() << " segments ");
27 
28  std::vector<const Muon::MuonSegment*>::const_iterator it = segments.begin();
29  std::vector<const Muon::MuonSegment*>::const_iterator it2 = segments.begin();
30  std::vector<const Muon::MuonSegment*>::const_iterator it_end = segments.end();
31 
32  double pmin = 10000000000.;
33  for(; it < it_end ; ++it ) {
34  for(it2 = it; it2 < it_end ; ++it2 ) {
35  if (it == it2) continue;
36  double smom;
37  fitMomentum2Segments(ctx, *it, *it2, smom);
38  ATH_MSG_DEBUG(" Fit pair of segments signed momentum " << smom);
39  if (fabs(smom) < fabs(pmin)) pmin = smom;
40  }
41  }
42 
43  signedMomentum = pmin;
44  if (signedMomentum < 100. && signedMomentum > 0 ) signedMomentum = 100.;
45  if (signedMomentum > -100. && signedMomentum <= 0 ) signedMomentum = -100.;
46 
47  ATH_MSG_DEBUG( " Estimated signed momentum " << signedMomentum);
48 }
49 
50 
51 
52 void MuonSegmentMomentum::fitMomentum2Segments( const EventContext&, const Muon::MuonSegment* segment1, const Muon::MuonSegment* segment2, double & signedMomentum ) const
53 {
54 
58 // MsgStream log(messageService(),name());
59  ATH_MSG_VERBOSE(" Executing MuonSegmentMomentumTool fitMomentum2Segments ");
60 
61  ATH_MSG_DEBUG(" fitMomentum2Segments");
62 
63 
64 
65  // Tracking matrices
66  Amg::MatrixX model(5,3);// H projection matrix PDG
67  model.setIdentity();
68  Amg::MatrixX covT(3,3); // covariance or Track parameters here 3: z beam; r/z segment; tan/cot theta segment;
69  covT.setIdentity();
70  Amg::MatrixX v(3,1);// V = H^T cov^-1 r cov-1 = covariance of measurements
71  v.setIdentity();
72  Amg::MatrixX t(3,1);// Track parameters
73  t.setIdentity();
74  Amg::MatrixX ym(5,1); // Vector of measurememns
75  ym.setIdentity();
76  // Extrapolation Matrices
77  Amg::MatrixX resi(5,1); // residuals
78  resi.setIdentity();
79 
80  // Geometry conventions
81  double radius_cylinder = 4000.;
82  double z_cylinder = 6000.;
83  double z_end = 15000.;
84  double cos_barrel = std::cos(std::atan2(radius_cylinder,z_cylinder));
85  int imeth;
86 
87  // First Segment
88  double thetas = segment1->globalDirection().theta();
89  double rs = segment1->globalPosition().perp();
90  double zs = segment1->globalPosition().z();
91 // double phi = segment1->globalDirection().phi();
92 // double phisp = atan2( segment1->globalPosition().y(), segment1->globalPosition().x());
93  double thesp = atan2( rs, zs);
94  // angular resolution
95  double era1 = 0.002;
96  // Second Segment
97  double thetase = segment2->globalDirection().theta();
98  double rse = segment2->globalPosition().perp();
99  double zse = segment2->globalPosition().z();
100 // double phie = segment2->globalDirection().phi();
101 // double phispe = atan2( segment2->globalPosition().y(), segment2->globalPosition().x());
102  double thespe = atan2( rse,zse);
103  double era2 = 0.002;
104  ATH_MSG_DEBUG(" radius1 " << rs << " radius2 " << rse << " z1 " << zs << " z2 " << zse << " theta1 " << thetas << " theta2 " << thetase << " thetasp1 " << thesp << " thetasp2 " << thespe);
105 
106  bool barrel = false;
107  if (fabs(cos(thesp)) < cos_barrel || fabs(cos(thespe)) < cos_barrel ) barrel = true;
108  bool forward = false;
109  if (fabs(cos(thesp)) > cos_barrel || fabs(cos(thespe)) > cos_barrel ) forward = true;
110 
111  double scf = 20.;
112  ATH_MSG_DEBUG(" error scaling in Curved fit " << scf);
113  if (forward) scf = 2*scf;
114  double scfn = 20.0;
115  if (forward) scfn = 2*scfn;
116 
117 
118 // Error definitions
119 
120  // double scf = 1; // pt = 100 GeV scf = 1; pt 5 GeV scf = 20*20
121  double ers2 = 0.1*0.1+scf*scf+scfn*scfn; //error squared position
122  double ebs2 = 50*50*fabs(sin(thetas)*sin(thetase))+ers2; //error squared beam position
123  double ets21 = era1*era1 + 0.002*0.002*scf*scf; // error squared angle
124  double ets22 = era2*era2 + 0.002*0.002*scf*scf; // error squared angle
125 
126 // No beam constraint for cosmics
127 
128  if(m_cosmics) ebs2 = 100000.*100000.;
129 
130  double sign = 1;
131  if ( zs < 0 ) sign = -1.;
132 
133  if (barrel) {
134  imeth = 0;
135  // Barrel Track Model Matrix
136  model(0,0) = 1.;
137  model(1,0) = 1.;
138  model(1,1) = rs;
139  model(1,2) = (rs-radius_cylinder)*(rs-radius_cylinder);
140  model(2,1) = 1.;
141  model(2,2) = 2*(rs-radius_cylinder);
142  model(3,0) = 1.;
143  model(3,1) = rse;
144  model(3,2) = (rse-radius_cylinder)*(rse-radius_cylinder);
145  model(4,1) = 1.;
146  model(4,2) = 2*(rse-radius_cylinder);
147  // Measurements ym
148  // correspondig squared errors: zs -> ers2 cot(thetas) ->ets2
149  ym(0,0) = 0.;
150  ym(1,0) = zs;
151  ym(2,0)= cos(thetas)/sin(thetas);
152  ym(3,0) = zse;
153  ym(4,0)= cos(thetase)/sin(thetase);
154  } else {
155  imeth = 1;
156  ATH_MSG_DEBUG(" forward fit " << cos(thetas));
157  // Forward Track Model Matrix
158  model(0,0) = 1.;
159  model(1,0) = 1.;
160  model(1,1) = zs;
161  model(1,2) = sign*(zs-sign*z_cylinder)*(zs-sign*z_cylinder);
162  model(2,1) = 1;
163  model(2,2) = sign*2*(zs-sign*z_cylinder);
164  model(3,0) = 1.;
165  model(3,1) = zse;
166  model(3,2) = sign*(zse-sign*z_cylinder)*(zse-sign*z_cylinder);
167  model(4,1) = 1;
168  model(4,2) = sign*2*(zse-sign*z_cylinder);
169 
170  if (fabs(zs) > z_end+2000) {
171  model(1,2) = sign*(-z_end*z_end + z_cylinder*z_cylinder + 2*zs*sign*(z_end-z_cylinder)
172  + (zs-sign*z_end)*(zs-sign*z_end)/5.);
173  model(2,2) = sign*(2*(sign*z_end-sign*z_cylinder)
174  + (zs-sign*z_end)/5.);
175  }
176  if (fabs(zse) > z_end+2000) {
177  model(3,2) = sign*(-z_end*z_end + z_cylinder*z_cylinder + 2*zse*sign*(z_end-z_cylinder)
178  + (zse-sign*z_end)*(zse-sign*z_end)/5.);
179  model(4,2) = sign*(2*(sign*z_end-sign*z_cylinder)
180  + (zse-sign*z_end)/5.);
181  }
182  // Measurements ym
183  // correspondig squared errors: rs -> ers2 tan(thetas) ->ets2
184  ym(0,0) = 0.;
185  ym(1,0) = rs;
186  ym(2,0)= tan(thetas);
187  ym(3,0) = rse;
188  ym(4,0)= tan(thetase);
189  // std::cout << " Forward zs " << zs << std::endl;
190  }
191  ATH_MSG_DEBUG(" distance segments " << sqrt((zs-zse)*(zs-zse)+(rs-rse)*(rs-rse)) );
192 
193  for(int i = 0; i <3 ; ++i ) {
194  v(i,0)= model(0,i)*ym(0,0)/ebs2 + model(1,i)*ym(1,0)/ers2 + model(2,i)*ym(2,0)/ets21
195  + model(3,i)*ym(3,0)/ers2 + model(4,i)*ym(4,0)/ets22;
196  }
197 
198  Amg::MatrixX covTI(3,3); // covariance Inverse of Track parameters
199  covTI.setIdentity();
200  for(int i = 0; i <3 ; ++i ) {
201  for(int j = 0; j <3 ; ++j ) {
202  covTI(i,j) += model(0,i)*model(0,j)/ebs2;
203  covTI(i,j) += model(1,i)*model(1,j)/ers2;
204  covTI(i,j) += model(2,i)*model(2,j)/ets21;
205  covTI(i,j) += model(3,i)*model(3,j)/ers2;
206  covTI(i,j) += model(4,i)*model(4,j)/ets22;
207  }
208  }
209  // Solution for Track parameters
210  covT = covTI.inverse();
211  t = covT*v;
212  // std::cout << " covariance at track 00 " << covT(0,0) << " 11 " << covT(1,1) << " 12 " << covT(1,2) << " 22 " << covT(2,2) << std::endl;
213 
214  double theta=0.; double invcurvature=0.;
215  if (imeth == 0) {
216  theta = atan2(1.,t(1,0));
217  invcurvature = t(2,0)*sin(theta);
218  } else if (imeth == 1) {
219  theta = atan2(1.,1./t(1,0));
220  invcurvature = -t(2,0)*fabs(cos(theta));
221  }
222 
223  signedMomentum = (1./invcurvature)/10.;
224 
225  ATH_MSG_DEBUG(" MuonSegmentMomentum in MeV " << (1./invcurvature)/10. << " theta fit " << theta << " cos theta " << cos(theta) );
226 
227  // calculate residuals and chi2
228  for(int i = 0; i <5 ; ++i ) {
229  resi(i,0)= model(i,0)*t(0,0) + model(i,1)*t(1,0) + model(i,2)*t(2,0) - ym(i,0);
230  }
231  double chi2 = resi(0,0)*resi(0,0)/ebs2 + resi(1,0)*resi(1,0)/ers2 + resi(2,0)*resi(2,0)/ets21
232  + resi(3,0)*resi(3,0)/ers2 + resi(4,0)*resi(4,0)/ets22;
233  ATH_MSG_DEBUG("resi 00 " << resi(0,0) << " chi2 " << chi2);
234  ATH_MSG_DEBUG(" Track parameters Matrix T00 " << t(0,0) << " T10 " << t(1,0) << " T20 " << t(2,0) );
235 
236  // Reshuffle residuals res(0) -> segment position in model (1,x)
237  // res(1) -> segment angle in model (2,x)
238 
239  Amg::VectorX res(4);
240  for(int i = 0; i<4 ; ++i ) {
241  res[i]= resi(i,0);
242  }
243 
244  std::vector <double> pull(4);
245 
246  bool toobig = false;
247  for(int i = 0; i <4 ; ++i ) {
248  if ( i==0 ) pull[i] = res[i]/sqrt(ers2);
249  if ( i==1 ) pull[i] = res[i]/sqrt(ets21);
250  if ( i==2 ) pull[i] = res[i]/sqrt(ers2);
251  if ( i==3 ) pull[i] = res[i]/sqrt(ets22);
252  if ( fabs( pull[i] ) > 5 ) toobig = true;
253  }
254 
255  if (toobig) {
256  ATH_MSG_DEBUG(" Pull too BIGFIT " << " rad pos1 " << rs << " rad pos2 " << rse << " ang1 " << thetas << " ang2 " << thetase );
257  double phisp = atan2( segment1->globalPosition().y(), segment1->globalPosition().x());
258  double phispe = atan2( segment2->globalPosition().y(), segment2->globalPosition().x());
259  ATH_MSG_DEBUG(" z pos1 " << zs << " z pos2 " << zse << " phi pos 1 " << phisp << " phi pos 2 " << phispe );
260 
261  Amg::Vector3D d1 = segment1->globalDirection();
262  Amg::Vector3D d2 = segment2->globalDirection();
263  ATH_MSG_DEBUG(" diff phi " << d1.x()*d2.y() - d1.y()*d2.x() );
264  }
265  ATH_MSG_DEBUG( " Fit2MomentumSegments: residual 0 " << res[0] << " pull 0 " << pull[0] << " residual 1 " << res[1] << " pull 1 " << pull[1] );
266  ATH_MSG_DEBUG(" Fit2MomentumSegments: residual 2 " << res[2] << " pull 2 " << pull[2] << " residual 3 " << res[3] << " pull 3 " << pull[3] );
267  ATH_MSG_DEBUG(" radius 1 " << ym(1,0) << " cottan theta 1 " << ym(2,0) << " radius 2 " << ym(3,0)
268  << " cottan theta 2 " << ym(4,0) );
269  ATH_MSG_DEBUG(" radius fit 1 " << -(res[0]+ym(1,0)) << " cottan theta fit 1 " << -(res[1]+ym(2,0)) << " radius 2 fit " << -(res[2]+ym(3,0))
270  << " cottan theta 2 fit " << -(res[3]+ym(4,0)) );
271 }
272 
273 
274 
275 
Amg::VectorX
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorX
Dynamic Vector - dynamic allocation.
Definition: EventPrimitives.h:30
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:27
theta
Scalar theta() const
theta method
Definition: AmgMatrixBasePlugin.h:75
skel.it
it
Definition: skel.GENtoEVGEN.py:407
dq_defect_virtual_defect_validation.d1
d1
Definition: dq_defect_virtual_defect_validation.py:79
MuonHough::z_cylinder
constexpr double z_cylinder
length of cylinder
Definition: MuonHoughMathUtils.h:26
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
MuonSegmentMomentum::m_cosmics
Gaudi::Property< bool > m_cosmics
flag for use of cosmics, straight line model will be used, no interaction point constraint
Definition: MuonSegmentMomentum.h:32
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
MuonSegmentMomentum::fitMomentumVectorSegments
virtual void fitMomentumVectorSegments(const EventContext &, const std::vector< const Muon::MuonSegment * > &, double &signedMomentum) const override
fits a momentum to a vector of segments
Definition: MuonSegmentMomentum.cxx:18
lumiFormat.i
int i
Definition: lumiFormat.py:85
MuonSegmentMomentum::MuonSegmentMomentum
MuonSegmentMomentum(const std::string &, const std::string &, const IInterface *)
constructor
Definition: MuonSegmentMomentum.cxx:11
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
MuonSegmentMomentum.h
res
std::pair< std::vector< unsigned int >, bool > res
Definition: JetGroupProductTest.cxx:11
chi2
double chi2(TH1 *h0, TH1 *h1)
Definition: comparitor.cxx:525
test_pyathena.parent
parent
Definition: test_pyathena.py:15
sign
int sign(int a)
Definition: TRT_StrawNeighbourSvc.h:108
python.StandardJetMods.pull
pull
Definition: StandardJetMods.py:309
MuonHough::z_end
constexpr double z_end
z value whereafter no magnetic field / curvature
Definition: MuonHoughMathUtils.h:28
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
MuonSegmentMomentum::fitMomentum2Segments
virtual void fitMomentum2Segments(const EventContext &, const Muon::MuonSegment *segment1, const Muon::MuonSegment *segment2, double &signedMomentum) const override
fits a momentum to 2 segments
Definition: MuonSegmentMomentum.cxx:52
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
python.PyAthena.v
v
Definition: PyAthena.py:154
correlationModel::model
model
Definition: AsgElectronEfficiencyCorrectionTool.cxx:46
dq_defect_virtual_defect_validation.d2
d2
Definition: dq_defect_virtual_defect_validation.py:81
DetectorZone::barrel
@ barrel
MuonSegment.h
Muon::MuonSegment::globalPosition
virtual const Amg::Vector3D & globalPosition() const override final
global position
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:157
MuonHough::radius_cylinder
constexpr double radius_cylinder
radius of cylinder
Definition: MuonHoughMathUtils.h:24
Muon::MuonSegment
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:45
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
AthAlgTool
Definition: AthAlgTool.h:26
Muon::MuonSegment::globalDirection
const Amg::Vector3D & globalDirection() const
global direction
Definition: MuonSpectrometer/MuonReconstruction/MuonRecEvent/MuonSegment/MuonSegment/MuonSegment.h:163
Generated on Mon Nov 3 2025 21:15:34 for ATLAS Offline Software by   doxygen 1.8.18