CMS 3D CMS Logo

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