CMS 3D CMS Logo

VertexKinematicConstraintT.cc
Go to the documentation of this file.
4 
6 
8 
9 void VertexKinematicConstraintT::init(const std::vector<KinematicState>& states,
10  const GlobalPoint& ipoint,
11  const GlobalVector& fieldValue) {
12  int num = states.size();
13  if (num != 2)
14  throw VertexException("VertexKinematicConstraintT !=2 states passed");
15  double mfz = fieldValue.z();
16 
17  int j = 0;
18  for (std::vector<KinematicState>::const_iterator i = states.begin(); i != states.end(); i++) {
19  mom[j] = i->globalMomentum();
20  dpos[j] = ipoint - i->globalPosition();
21  a_i[j] = -i->particleCharge() * mfz;
22 
23  double pvx = mom[j].x() - a_i[j] * dpos[j].y();
24  double pvy = mom[j].y() + a_i[j] * dpos[j].x();
25  double pvt2 = pvx * pvx + pvy * pvy;
26  novera[j] = (dpos[j].x() * mom[j].x() + dpos[j].y() * mom[j].y());
27  n[j] = a_i[j] * novera[j];
28  m[j] = (pvx * mom[j].x() + pvy * mom[j].y());
29  k[j] = -mom[j].z() / (mom[j].perp2() * pvt2);
30  delta[j] = std::atan2(n[j], m[j]);
31 
32  ++j;
33  }
34 }
35 
37  //it is 2 equations per track
38  for (int j = 0; j != 2; ++j) {
39  if (a_i[j] != 0) {
40  //vector of values
41  super::vl(j * 2) = dpos[j].y() * mom[j].x() - dpos[j].x() * mom[j].y() -
42  a_i[j] * (dpos[j].x() * dpos[j].x() + dpos[j].y() * dpos[j].y()) * 0.5;
43  super::vl(j * 2 + 1) = dpos[j].z() - mom[j].z() * delta[j] / a_i[j];
44  } else {
45  //neutral particle
46  double pt2Inverse = 1. / mom[j].perp2();
47  super::vl(j * 2) = dpos[j].y() * mom[j].x() - dpos[j].x() * mom[j].y();
48  super::vl(j * 2 + 1) =
49  dpos[j].z() - mom[j].z() * (dpos[j].x() * mom[j].x() + dpos[j].y() * mom[j].y()) * pt2Inverse;
50  }
51  }
52 }
53 
55  ROOT::Math::SMatrix<double, 2, 7> el_part_d;
56  for (int j = 0; j != 2; ++j) {
57  if (a_i[j] != 0) {
58  //charged particle
59 
60  //D Jacobian matrix
61  el_part_d(0, 0) = mom[j].y() + a_i[j] * dpos[j].x();
62  el_part_d(0, 1) = -mom[j].x() + a_i[j] * dpos[j].y();
63  el_part_d(1, 0) = -k[j] * (m[j] * mom[j].x() - n[j] * mom[j].y());
64  el_part_d(1, 1) = -k[j] * (m[j] * mom[j].y() + n[j] * mom[j].x());
65  el_part_d(1, 2) = -1.;
66  el_part_d(0, 3) = dpos[j].y();
67  el_part_d(0, 4) = -dpos[j].x();
68  el_part_d(1, 3) = k[j] * (m[j] * dpos[j].x() - novera[j] * (2 * mom[j].x() - a_i[j] * dpos[j].y()));
69  el_part_d(1, 4) = k[j] * (m[j] * dpos[j].y() - novera[j] * (2 * mom[j].y() + a_i[j] * dpos[j].x()));
70  el_part_d(1, 5) = -delta[j] / a_i[j];
71  super::jac_d().Place_at(el_part_d, j * 2, j * 7);
72  } else {
73  //neutral particle
74  double pt2Inverse = 1. / mom[j].perp2();
75  el_part_d(0, 0) = mom[j].y();
76  el_part_d(0, 1) = -mom[j].x();
77  el_part_d(1, 0) = mom[j].x() * (mom[j].z() * pt2Inverse);
78  el_part_d(1, 1) = mom[j].y() * (mom[j].z() * pt2Inverse);
79  el_part_d(1, 2) = -1.;
80  el_part_d(0, 3) = dpos[j].y();
81  el_part_d(0, 4) = -dpos[j].x();
82  el_part_d(1, 3) = 2 * (dpos[j].x() * mom[j].x() + dpos[j].y() * mom[j].y()) * pt2Inverse * mom[j].x() *
83  (mom[j].z() * pt2Inverse) -
84  dpos[j].x() * (mom[j].z() * pt2Inverse);
85  el_part_d(1, 4) = 2 * (dpos[j].x() * mom[j].x() + dpos[j].y() * mom[j].y()) * pt2Inverse * mom[j].y() *
86  (mom[j].z() * pt2Inverse) -
87  dpos[j].x() * (mom[j].z() * pt2Inverse);
88  el_part_d(1, 5) = -(dpos[j].x() * mom[j].x() + dpos[j].y() * mom[j].y()) * pt2Inverse;
89  super::jac_d().Place_at(el_part_d, j * 2, j * 7);
90  }
91  }
92 }
93 
95  ROOT::Math::SMatrix<double, 2, 3> el_part_e;
96  for (int j = 0; j != 2; ++j) {
97  if (a_i[j] != 0) {
98  //charged particle
99 
100  //E jacobian matrix
101  el_part_e(0, 0) = -(mom[j].y() + a_i[j] * dpos[j].x());
102  el_part_e(0, 1) = mom[j].x() - a_i[j] * dpos[j].y();
103  el_part_e(1, 0) = k[j] * (m[j] * mom[j].x() - n[j] * mom[j].y());
104  el_part_e(1, 1) = k[j] * (m[j] * mom[j].y() + n[j] * mom[j].x());
105  el_part_e(1, 2) = 1;
106  super::jac_e().Place_at(el_part_e, 2 * j, 0);
107  } else {
108  //neutral particle
109  double pt2Inverse = 1. / mom[j].perp2();
110  el_part_e(0, 0) = -mom[j].y();
111  el_part_e(0, 1) = mom[j].x();
112  el_part_e(1, 0) = -mom[j].x() * mom[j].z() * pt2Inverse;
113  el_part_e(1, 1) = -mom[j].y() * mom[j].z() * pt2Inverse;
114  el_part_e(1, 2) = 1;
115  super::jac_e().Place_at(el_part_e, 2 * j, 0);
116  }
117  }
118 }
119 
Vector3DBase
Definition: Vector3DBase.h:8
VertexKinematicConstraintT::VertexKinematicConstraintT
VertexKinematicConstraintT()
Definition: VertexKinematicConstraintT.cc:5
DDAxes::y
mps_fire.i
i
Definition: mps_fire.py:428
MessageLogger.h
VertexKinematicConstraintT::numberOfEquations
int numberOfEquations() const override
Definition: VertexKinematicConstraintT.cc:120
VertexException
Common base class.
Definition: VertexException.h:12
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
VertexKinematicConstraintT::novera
double novera[2]
Definition: VertexKinematicConstraintT.h:21
VertexKinematicConstraintT::delta
double delta[2]
Definition: VertexKinematicConstraintT.h:21
DDAxes::x
VertexKinematicConstraintT::n
double n[2]
Definition: VertexKinematicConstraintT.h:21
MultiTrackKinematicConstraintT< 2, 4 >::vl
valueType & vl() const
Definition: MultiTrackKinematicConstraintT.h:92
VertexKinematicConstraintT::fillPositionDerivative
void fillPositionDerivative() const override
Definition: VertexKinematicConstraintT.cc:94
VertexKinematicConstraintT::dpos
GlobalVector dpos[2]
Definition: VertexKinematicConstraintT.h:23
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
MultiTrackKinematicConstraintT< 2, 4 >::jac_d
parametersDerivativeType & jac_d() const
Definition: MultiTrackKinematicConstraintT.h:93
VertexKinematicConstraintT.h
plotBeamSpotDB.ipoint
ipoint
Definition: plotBeamSpotDB.py:340
Point3DBase< float, GlobalTag >
VertexKinematicConstraintT::~VertexKinematicConstraintT
~VertexKinematicConstraintT() override
Definition: VertexKinematicConstraintT.cc:7
MultiTrackKinematicConstraintT< 2, 4 >::jac_e
positionDerivativeType & jac_e() const
Definition: MultiTrackKinematicConstraintT.h:94
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
VertexKinematicConstraintT::init
void init(const std::vector< KinematicState > &states, const GlobalPoint &point, const GlobalVector &mf) override
Definition: VertexKinematicConstraintT.cc:9
VertexKinematicConstraintT::mom
GlobalVector mom[2]
Definition: VertexKinematicConstraintT.h:22
EgammaValidation_cff.num
num
Definition: EgammaValidation_cff.py:34
VertexKinematicConstraintT::fillValue
void fillValue() const override
Definition: VertexKinematicConstraintT.cc:36
VertexKinematicConstraintT::a_i
double a_i[2]
Definition: VertexKinematicConstraintT.h:20
VertexException.h
VertexKinematicConstraintT::k
double k[2]
Definition: VertexKinematicConstraintT.h:21
VertexKinematicConstraintT::m
double m[2]
Definition: VertexKinematicConstraintT.h:21
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
VertexKinematicConstraintT::fillParametersDerivative
void fillParametersDerivative() const override
Definition: VertexKinematicConstraintT.cc:54
PV3DBase::perp2
T perp2() const
Definition: PV3DBase.h:68