CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
TwoBodyDecayDerivatives.cc
Go to the documentation of this file.
1 
4 
6 
7 #include <algorithm>
8 
9 TwoBodyDecayDerivatives::TwoBodyDecayDerivatives(double mPrimary, double mSecondary)
10  : thePrimaryMass(mPrimary), theSecondaryMass(mSecondary) {}
11 
13 
14 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::derivatives(const TwoBodyDecay &tbd) const {
15  return derivatives(tbd.decayParameters());
16 }
17 
18 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::derivatives(
19  const TwoBodyDecayParameters &param) const {
20  // get the derivatives with respect to all parameters
21  std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdpx = this->dqsdpx(param);
22  std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdpy = this->dqsdpy(param);
23  std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdpz = this->dqsdpz(param);
24  std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdtheta = this->dqsdtheta(param);
25  std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdphi = this->dqsdphi(param);
26  std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdm = this->dqsdm(param);
27 
28  AlgebraicMatrix dqplusdz(3, dimension);
29  dqplusdz.sub(1, px, dqsdpx.first);
30  dqplusdz.sub(1, py, dqsdpy.first);
31  dqplusdz.sub(1, pz, dqsdpz.first);
32  dqplusdz.sub(1, theta, dqsdtheta.first);
33  dqplusdz.sub(1, phi, dqsdphi.first);
34  dqplusdz.sub(1, mass, dqsdm.first);
35 
36  AlgebraicMatrix dqminusdz(3, dimension);
37  dqminusdz.sub(1, px, dqsdpx.second);
38  dqminusdz.sub(1, py, dqsdpy.second);
39  dqminusdz.sub(1, pz, dqsdpz.second);
40  dqminusdz.sub(1, theta, dqsdtheta.second);
41  dqminusdz.sub(1, phi, dqsdphi.second);
42  dqminusdz.sub(1, mass, dqsdm.second);
43 
44  return std::make_pair(dqplusdz, dqminusdz);
45 }
46 
47 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::selectedDerivatives(
48  const TwoBodyDecay &tbd, const std::vector<bool> &selector) const {
49  return selectedDerivatives(tbd.decayParameters(), selector);
50 }
51 
52 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::selectedDerivatives(
53  const TwoBodyDecayParameters &param, const std::vector<bool> &selector) const {
54  if (selector.size() != dimension) {
55  throw cms::Exception("BadConfig") << "@SUB=TwoBodyDecayDerivatives::selectedDerivatives"
56  << "selector has bad dimension (size=" << selector.size() << ").";
57  }
58 
59  int nSelected = std::count(selector.begin(), selector.end(), true);
60  int iSelected = 1;
61 
62  AlgebraicMatrix dqplusdz(3, nSelected);
63  AlgebraicMatrix dqminusdz(3, nSelected);
64  std::pair<AlgebraicMatrix, AlgebraicMatrix> dqsdzi;
65 
66  for (unsigned int i = 1; i <= dimension; i++) {
67  if (selector[i]) {
68  dqsdzi = this->dqsdzi(param, DerivativeParameterName(i));
69  dqplusdz.sub(1, iSelected, dqsdzi.first);
70  dqminusdz.sub(1, iSelected, dqsdzi.second);
71  iSelected++;
72  }
73  }
74 
75  return std::make_pair(dqplusdz, dqminusdz);
76 }
77 
78 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdpx(
79  const TwoBodyDecayParameters &param) const {
80  double px = param[TwoBodyDecayParameters::px];
81  double py = param[TwoBodyDecayParameters::py];
82  double pz = param[TwoBodyDecayParameters::pz];
83  double theta = param[TwoBodyDecayParameters::theta];
84  double phi = param[TwoBodyDecayParameters::phi];
85 
86  // compute transverse and absolute momentum
87  double pT2 = px * px + py * py;
88  double p2 = pT2 + pz * pz;
89  double pT = sqrt(pT2);
90  double p = sqrt(p2);
91 
92  double sphi = sin(phi);
93  double cphi = cos(phi);
94  double stheta = sin(theta);
95  double ctheta = cos(theta);
96 
97  // some constants from kinematics
98  double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
99  double c2 = sqrt(c1 * c1 - 1.);
100  double c3 = 0.5 * c2 * ctheta / c1;
101  double c4 = sqrt(p2 + thePrimaryMass * thePrimaryMass);
102 
103  // momentum of decay particle 1 in the primary's boosted frame
104  AlgebraicMatrix pplus(3, 1);
105  pplus[0][0] = theSecondaryMass * c2 * stheta * cphi;
106  pplus[1][0] = theSecondaryMass * c2 * stheta * sphi;
107  pplus[2][0] = 0.5 * p + c3 * c4;
108 
109  // momentum of decay particle 2 in the primary's boosted frame
110  AlgebraicMatrix pminus(3, 1);
111  pminus[0][0] = -pplus[0][0];
112  pminus[1][0] = -pplus[1][0];
113  pminus[2][0] = 0.5 * p - c3 * c4;
114 
115  // derivative of rotation matrix w.r.t. px
116  AlgebraicMatrix dRotMatdpx(3, 3);
117 
118  dRotMatdpx[0][0] = pz / (pT * p) * (1. - px * px * (1. / pT2 + 1. / p2));
119  dRotMatdpx[0][1] = px * py / (pT * pT2);
120  dRotMatdpx[0][2] = (1. - px * px / p2) / p;
121 
122  dRotMatdpx[1][0] = -px * py * pz / (pT * p) * (1. / pT2 + 1. / p2);
123  dRotMatdpx[1][1] = (1. - px * px / pT2) / pT;
124  dRotMatdpx[1][2] = -px * py / (p * p2);
125 
126  dRotMatdpx[2][0] = -(1. / pT - pT / p2) * px / p;
127  dRotMatdpx[2][1] = 0.;
128  dRotMatdpx[2][2] = -px * pz / (p * p2);
129 
130  // derivative of the momentum of particle 1 in the lab frame w.r.t. px
131  double dpplusdpx = px * (0.5 / p + c3 / c4);
132 
133  AlgebraicMatrix dqplusdpx = dRotMatdpx * pplus;
134  dqplusdpx[0][0] += px * dpplusdpx / p;
135  dqplusdpx[1][0] += py * dpplusdpx / p;
136  dqplusdpx[2][0] += pz * dpplusdpx / p;
137 
138  // derivative of the momentum of particle 2 in the lab frame w.r.t. px
139  double dpminusdpx = px * (0.5 / p - c3 / c4);
140 
141  AlgebraicMatrix dqminusdpx = dRotMatdpx * pminus;
142  dqminusdpx[0][0] += px * dpminusdpx / p;
143  dqminusdpx[1][0] += py * dpminusdpx / p;
144  dqminusdpx[2][0] += pz * dpminusdpx / p;
145 
146  // return result
147  return std::make_pair(dqplusdpx, dqminusdpx);
148 }
149 
150 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdpy(
151  const TwoBodyDecayParameters &param) const {
152  double px = param[TwoBodyDecayParameters::px];
153  double py = param[TwoBodyDecayParameters::py];
154  double pz = param[TwoBodyDecayParameters::pz];
155  double theta = param[TwoBodyDecayParameters::theta];
156  double phi = param[TwoBodyDecayParameters::phi];
157 
158  // compute transverse and absolute momentum
159  double pT2 = px * px + py * py;
160  double p2 = pT2 + pz * pz;
161  double pT = sqrt(pT2);
162  double p = sqrt(p2);
163 
164  double sphi = sin(phi);
165  double cphi = cos(phi);
166  double stheta = sin(theta);
167  double ctheta = cos(theta);
168 
169  // some constants from kinematics
170  double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
171  double c2 = sqrt(c1 * c1 - 1.);
172  double c3 = 0.5 * c2 * ctheta / c1;
173  double c4 = sqrt(p2 + thePrimaryMass * thePrimaryMass);
174 
175  // momentum of decay particle 1 in the rest frame of the primary
176  AlgebraicMatrix pplus(3, 1);
177  pplus[0][0] = theSecondaryMass * c2 * stheta * cphi;
178  pplus[1][0] = theSecondaryMass * c2 * stheta * sphi;
179  pplus[2][0] = 0.5 * p + c3 * c4;
180 
181  // momentum of decay particle 2 in the rest frame of the primary
182  AlgebraicMatrix pminus(3, 1);
183  pminus[0][0] = -pplus[0][0];
184  pminus[1][0] = -pplus[1][0];
185  pminus[2][0] = 0.5 * p - c3 * c4;
186 
187  // derivative of rotation matrix w.r.t. py
188  AlgebraicMatrix dRotMatdpy(3, 3);
189 
190  dRotMatdpy[0][0] = -px * py * pz / (pT * p) * (1. / pT2 + 1. / p2);
191  dRotMatdpy[0][1] = (py * py / pT2 - 1.) / pT;
192  dRotMatdpy[0][2] = -px * py / (p * p2);
193 
194  dRotMatdpy[1][0] = pz / (pT * p) * (1. - py * py * (1. / pT2 + 1. / p2));
195  dRotMatdpy[1][1] = -px * py / (pT * pT2);
196  dRotMatdpy[1][2] = (1. - py * py / p2) / p;
197 
198  dRotMatdpy[2][0] = -(1. / pT - pT / p2) * py / p;
199  dRotMatdpy[2][1] = 0.;
200  dRotMatdpy[2][2] = -py * pz / (p * p2);
201 
202  // derivative of the momentum of particle 1 in the lab frame w.r.t. py
203  double dpplusdpy = py * (0.5 / p + c3 / c4);
204 
205  AlgebraicMatrix dqplusdpy = dRotMatdpy * pplus;
206  dqplusdpy[0][0] += px * dpplusdpy / p;
207  dqplusdpy[1][0] += py * dpplusdpy / p;
208  dqplusdpy[2][0] += pz * dpplusdpy / p;
209 
210  // derivative of the momentum of particle 2 in the lab frame w.r.t. py
211  double dpminusdpy = py * (0.5 / p - c3 / c4);
212 
213  AlgebraicMatrix dqminusdpy = dRotMatdpy * pminus;
214  dqminusdpy[0][0] += px * dpminusdpy / p;
215  dqminusdpy[1][0] += py * dpminusdpy / p;
216  dqminusdpy[2][0] += pz * dpminusdpy / p;
217 
218  // return result
219  return std::make_pair(dqplusdpy, dqminusdpy);
220 }
221 
222 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdpz(
223  const TwoBodyDecayParameters &param) const {
224  double px = param[TwoBodyDecayParameters::px];
225  double py = param[TwoBodyDecayParameters::py];
226  double pz = param[TwoBodyDecayParameters::pz];
227  double theta = param[TwoBodyDecayParameters::theta];
228  double phi = param[TwoBodyDecayParameters::phi];
229 
230  // compute transverse and absolute momentum
231  double pT2 = px * px + py * py;
232  double p2 = pT2 + pz * pz;
233  double pT = sqrt(pT2);
234  double p = sqrt(p2);
235 
236  double sphi = sin(phi);
237  double cphi = cos(phi);
238  double stheta = sin(theta);
239  double ctheta = cos(theta);
240 
241  // some constants from kinematics
242  double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
243  double c2 = sqrt(c1 * c1 - 1.);
244  double c3 = 0.5 * c2 * ctheta / c1;
245  double c4 = sqrt(p2 + thePrimaryMass * thePrimaryMass);
246 
247  // momentum of decay particle 1 in the rest frame of the primary
248  AlgebraicMatrix pplus(3, 1);
249  pplus[0][0] = theSecondaryMass * c2 * stheta * cphi;
250  pplus[1][0] = theSecondaryMass * c2 * stheta * sphi;
251  pplus[2][0] = 0.5 * p + c3 * c4;
252 
253  // momentum of decay particle 2 in the rest frame of the primary
254  AlgebraicMatrix pminus(3, 1);
255  pminus[0][0] = -pplus[0][0];
256  pminus[1][0] = -pplus[1][0];
257  pminus[2][0] = 0.5 * p - c3 * c4;
258 
259  // derivative of rotation matrix w.r.t. py
260  AlgebraicMatrix dRotMatdpz(3, 3);
261 
262  dRotMatdpz[0][0] = px / (pT * p) * (1. - pz * pz / p2);
263  dRotMatdpz[0][1] = 0.;
264  dRotMatdpz[0][2] = -px * pz / (p * p2);
265 
266  dRotMatdpz[1][0] = py / (p * pT) * (1. - pz * pz / p2);
267  dRotMatdpz[1][1] = 0.;
268  dRotMatdpz[1][2] = -py * pz / (p * p2);
269 
270  dRotMatdpz[2][0] = pT * pz / (p * p2);
271  dRotMatdpz[2][1] = 0.;
272  dRotMatdpz[2][2] = (1. - pz * pz / p2) / p;
273 
274  // derivative of the momentum of particle 1 in the lab frame w.r.t. pz
275  double dpplusdpz = pz * (0.5 / p + c3 / c4);
276 
277  AlgebraicMatrix dqplusdpz = dRotMatdpz * pplus;
278  dqplusdpz[0][0] += px * dpplusdpz / p;
279  dqplusdpz[1][0] += py * dpplusdpz / p;
280  dqplusdpz[2][0] += pz * dpplusdpz / p;
281 
282  // derivative of the momentum of particle 2 in the lab frame w.r.t. pz
283  double dpminusdpz = pz * (0.5 / p - c3 / c4);
284 
285  AlgebraicMatrix dqminusdpz = dRotMatdpz * pminus;
286  dqminusdpz[0][0] += px * dpminusdpz / p;
287  dqminusdpz[1][0] += py * dpminusdpz / p;
288  dqminusdpz[2][0] += pz * dpminusdpz / p;
289 
290  // return result
291  return std::make_pair(dqplusdpz, dqminusdpz);
292 }
293 
294 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdtheta(
295  const TwoBodyDecayParameters &param) const {
296  double px = param[TwoBodyDecayParameters::px];
297  double py = param[TwoBodyDecayParameters::py];
298  double pz = param[TwoBodyDecayParameters::pz];
299  double theta = param[TwoBodyDecayParameters::theta];
300  double phi = param[TwoBodyDecayParameters::phi];
301 
302  // compute transverse and absolute momentum
303  double pT2 = px * px + py * py;
304  double p2 = pT2 + pz * pz;
305 
306  double sphi = sin(phi);
307  double cphi = cos(phi);
308  double stheta = sin(theta);
309  double ctheta = cos(theta);
310 
311  // some constants from kinematics
312  double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
313  double c2 = sqrt(c1 * c1 - 1.);
314  double c3 = -0.5 * c2 * stheta / c1;
315  double c4 = sqrt(p2 + thePrimaryMass * thePrimaryMass);
316 
317  // derivative of the momentum of particle 1 in the primary's rest frame w.r.t.
318  // theta
319  AlgebraicMatrix dpplusdtheta(3, 1);
320  dpplusdtheta[0][0] = theSecondaryMass * c2 * ctheta * cphi;
321  dpplusdtheta[1][0] = theSecondaryMass * c2 * ctheta * sphi;
322  dpplusdtheta[2][0] = c3 * c4;
323 
324  // derivative of the momentum of particle 2 in the primary's rest frame w.r.t.
325  // theta
326  AlgebraicMatrix dpminusdtheta(3, 1);
327  dpminusdtheta[0][0] = -theSecondaryMass * c2 * ctheta * cphi;
328  dpminusdtheta[1][0] = -theSecondaryMass * c2 * ctheta * sphi;
329  dpminusdtheta[2][0] = -c3 * c4;
330 
331  TwoBodyDecayModel decayModel;
332  AlgebraicMatrix rotMat = decayModel.rotationMatrix(px, py, pz);
333 
334  AlgebraicMatrix dqplusdtheta = rotMat * dpplusdtheta;
335  AlgebraicMatrix dqminusdtheta = rotMat * dpminusdtheta;
336 
337  return std::make_pair(dqplusdtheta, dqminusdtheta);
338 }
339 
340 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdphi(
341  const TwoBodyDecayParameters &param) const {
342  double px = param[TwoBodyDecayParameters::px];
343  double py = param[TwoBodyDecayParameters::py];
344  double pz = param[TwoBodyDecayParameters::pz];
345  double theta = param[TwoBodyDecayParameters::theta];
346  double phi = param[TwoBodyDecayParameters::phi];
347 
348  double sphi = sin(phi);
349  double cphi = cos(phi);
350  double stheta = sin(theta);
351 
352  // some constants from kinematics
353  double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
354  double c2 = sqrt(c1 * c1 - 1.);
355 
356  // derivative of the momentum of particle 1 in the primary's rest frame w.r.t.
357  // phi
358  AlgebraicMatrix dpplusdphi(3, 1);
359  dpplusdphi[0][0] = -theSecondaryMass * c2 * stheta * sphi;
360  dpplusdphi[1][0] = theSecondaryMass * c2 * stheta * cphi;
361  dpplusdphi[2][0] = 0.;
362 
363  // derivative of the momentum of particle 2 in the primary's rest frame w.r.t.
364  // phi
365  AlgebraicMatrix dpminusdphi(3, 1);
366  dpminusdphi[0][0] = theSecondaryMass * c2 * stheta * sphi;
367  dpminusdphi[1][0] = -theSecondaryMass * c2 * stheta * cphi;
368  dpminusdphi[2][0] = 0.;
369 
370  TwoBodyDecayModel decayModel;
371  AlgebraicMatrix rotMat = decayModel.rotationMatrix(px, py, pz);
372 
373  AlgebraicMatrix dqplusdphi = rotMat * dpplusdphi;
374  AlgebraicMatrix dqminusdphi = rotMat * dpminusdphi;
375 
376  return std::make_pair(dqplusdphi, dqminusdphi);
377 }
378 
379 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdm(
380  const TwoBodyDecayParameters &param) const {
381  double px = param[TwoBodyDecayParameters::px];
382  double py = param[TwoBodyDecayParameters::py];
383  double pz = param[TwoBodyDecayParameters::pz];
384  double theta = param[TwoBodyDecayParameters::theta];
385  double phi = param[TwoBodyDecayParameters::phi];
386 
387  double pT2 = px * px + py * py;
388  double p2 = pT2 + pz * pz;
389 
390  double sphi = sin(phi);
391  double cphi = cos(phi);
392  double ctheta = cos(theta);
393  double stheta = sin(theta);
394 
395  // some constants from kinematics
396  double c1 = 0.5 * thePrimaryMass / theSecondaryMass;
397  double c2 = 1. / sqrt(c1 * c1 - 1.);
398  double m2 = thePrimaryMass * thePrimaryMass;
399 
400  // derivative of the momentum of particle 1 in the primary's rest frame w.r.t.
401  // the primary's mass
402  AlgebraicMatrix dpplusdm(3, 1);
403  dpplusdm[0][0] = c2 * 0.5 * c1 * stheta * cphi;
404  dpplusdm[1][0] = c2 * 0.5 * c1 * stheta * sphi;
405  dpplusdm[2][0] = c2 * theSecondaryMass * (c1 * c1 + p2 / m2) / sqrt(p2 + m2) * ctheta;
406 
407  // derivative of the momentum of particle 2 in the primary's rest frame w.r.t.
408  // the primary's mass
409  AlgebraicMatrix dpminusdm(3, 1);
410  dpminusdm[0][0] = -dpplusdm[0][0];
411  dpminusdm[1][0] = -dpplusdm[1][0];
412  dpminusdm[2][0] = -dpplusdm[2][0];
413 
414  TwoBodyDecayModel decayModel;
415  AlgebraicMatrix rotMat = decayModel.rotationMatrix(px, py, pz);
416 
417  AlgebraicMatrix dqplusdm = rotMat * dpplusdm;
418  AlgebraicMatrix dqminusdm = rotMat * dpminusdm;
419 
420  return std::make_pair(dqplusdm, dqminusdm);
421 }
422 
423 const std::pair<AlgebraicMatrix, AlgebraicMatrix> TwoBodyDecayDerivatives::dqsdzi(
424  const TwoBodyDecayParameters &param, const DerivativeParameterName &i) const {
425  switch (i) {
427  return dqsdpx(param);
428  break;
430  return dqsdpy(param);
431  break;
433  return dqsdpz(param);
434  break;
436  return dqsdtheta(param);
437  break;
439  return dqsdphi(param);
440  break;
442  return dqsdm(param);
443  break;
444  default:
445  throw cms::Exception("BadConfig") << "@SUB=TwoBodyDecayDerivatives::dqsdzi"
446  << "no decay parameter related to selector (" << i << ").";
447  };
448 
449  return std::make_pair(AlgebraicMatrix(3, 1, 0), AlgebraicMatrix(3, 1, 0));
450 }
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpx(const TwoBodyDecayParameters &param) const
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdphi(const TwoBodyDecayParameters &param) const
const TString p2
Definition: fwPaths.cc:13
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
TwoBodyDecayDerivatives(double mPrimary=91.1876, double mSecondary=0.105658)
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdzi(const TwoBodyDecayParameters &param, const DerivativeParameterName &i) const
tuple m2
Definition: callgraph.py:57
const std::pair< AlgebraicMatrix, AlgebraicMatrix > derivatives(const TwoBodyDecay &tbd) const
AlgebraicMatrix rotationMatrix(double px, double py, double pz)
CLHEP::HepMatrix AlgebraicMatrix
T sqrt(T t)
Definition: SSEVec.h:19
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdm(const TwoBodyDecayParameters &param) const
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpz(const TwoBodyDecayParameters &param) const
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdtheta(const TwoBodyDecayParameters &param) const
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpy(const TwoBodyDecayParameters &param) const
const TwoBodyDecayParameters & decayParameters(void) const
Definition: TwoBodyDecay.h:31
const std::pair< AlgebraicMatrix, AlgebraicMatrix > selectedDerivatives(const TwoBodyDecay &tbd, const std::vector< bool > &selector) const