10 thePrimaryMass( mPrimary ), theSecondaryMass( mSecondary ) {}
16 const std::pair< AlgebraicMatrix, AlgebraicMatrix >
22 const std::pair< AlgebraicMatrix, AlgebraicMatrix >
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 );
30 std::pair< AlgebraicMatrix, AlgebraicMatrix >
dqsdphi = this->
dqsdphi( param );
31 std::pair< AlgebraicMatrix, AlgebraicMatrix >
dqsdm = this->
dqsdm( param );
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 );
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 );
49 return std::make_pair( dqplusdz, dqminusdz );
53 const std::pair< AlgebraicMatrix, AlgebraicMatrix >
60 const std::pair< AlgebraicMatrix, AlgebraicMatrix >
65 throw cms::Exception(
"BadConfig" ) <<
"@SUB=TwoBodyDecayDerivatives::selectedDerivatives"
66 <<
"selector has bad dimension (size=" << selector.size() <<
").";
69 int nSelected =
std::count( selector.begin(), selector.end(),
true );
74 std::pair< AlgebraicMatrix, AlgebraicMatrix >
dqsdzi;
81 dqplusdz.sub( 1, iSelected, dqsdzi.first );
82 dqminusdz.sub( 1, iSelected, dqsdzi.second );
87 return std::make_pair( dqplusdz, dqminusdz );
100 double pT2 = px*px + py*
py;
101 double p2 = pT2 + pz*
pz;
102 double pT =
sqrt( pT2 );
103 double p =
sqrt( p2 );
105 double sphi =
sin( phi );
106 double cphi =
cos( phi );
107 double stheta =
sin( theta );
108 double ctheta =
cos( theta );
112 double c2 =
sqrt( c1*c1 - 1. );
113 double c3 = 0.5*c2*ctheta/
c1;
120 pplus[2][0] = 0.5*p + c3*c4;
124 pminus[0][0] = -pplus[0][0];
125 pminus[1][0] = -pplus[1][0];
126 pminus[2][0] = 0.5*p - c3*c4;
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;
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);
139 dRotMatdpx[2][0] = -(1./pT-pT/
p2)*px/p;
140 dRotMatdpx[2][1] = 0.;
141 dRotMatdpx[2][2] = -px*pz/(p*
p2);
144 double dpplusdpx = px*( 0.5/p + c3/c4 );
147 dqplusdpx[0][0] += px*dpplusdpx/
p;
148 dqplusdpx[1][0] += py*dpplusdpx/
p;
149 dqplusdpx[2][0] += pz*dpplusdpx/
p;
152 double dpminusdpx = px*( 0.5/p - c3/c4 );
155 dqminusdpx[0][0] += px*dpminusdpx/
p;
156 dqminusdpx[1][0] += py*dpminusdpx/
p;
157 dqminusdpx[2][0] += pz*dpminusdpx/
p;
160 return std::make_pair( dqplusdpx, dqminusdpx );
173 double pT2 = px*px + py*
py;
174 double p2 = pT2 + pz*
pz;
175 double pT =
sqrt( pT2 );
176 double p =
sqrt( p2 );
178 double sphi =
sin( phi );
179 double cphi =
cos( phi );
180 double stheta =
sin( theta );
181 double ctheta =
cos( theta );
185 double c2 =
sqrt( c1*c1 - 1. );
186 double c3 = 0.5*c2*ctheta/
c1;
193 pplus[2][0] = 0.5*p + c3*c4;
197 pminus[0][0] = -pplus[0][0];
198 pminus[1][0] = -pplus[1][0];
199 pminus[2][0] = 0.5*p - c3*c4;
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);
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;
212 dRotMatdpy[2][0] = -(1./pT-pT/
p2)*py/p;
213 dRotMatdpy[2][1] = 0.;
214 dRotMatdpy[2][2] = -py*pz/(p*
p2);
217 double dpplusdpy = py*( 0.5/p + c3/c4 );
220 dqplusdpy[0][0] += px*dpplusdpy/
p;
221 dqplusdpy[1][0] += py*dpplusdpy/
p;
222 dqplusdpy[2][0] += pz*dpplusdpy/
p;
225 double dpminusdpy = py*( 0.5/p - c3/c4 );
228 dqminusdpy[0][0] += px*dpminusdpy/
p;
229 dqminusdpy[1][0] += py*dpminusdpy/
p;
230 dqminusdpy[2][0] += pz*dpminusdpy/
p;
233 return std::make_pair( dqplusdpy, dqminusdpy );
246 double pT2 = px*px + py*
py;
247 double p2 = pT2 + pz*
pz;
248 double pT =
sqrt( pT2 );
249 double p =
sqrt( p2 );
251 double sphi =
sin( phi );
252 double cphi =
cos( phi );
253 double stheta =
sin( theta );
254 double ctheta =
cos( theta );
258 double c2 =
sqrt( c1*c1 - 1. );
259 double c3 = 0.5*c2*ctheta/
c1;
266 pplus[2][0] = 0.5*p + c3*c4;
270 pminus[0][0] = -pplus[0][0];
271 pminus[1][0] = -pplus[1][0];
272 pminus[2][0] = 0.5*p - c3*c4;
277 dRotMatdpz[0][0] = px/(pT*
p)*(1.-pz*pz/p2);
278 dRotMatdpz[0][1] = 0.;
279 dRotMatdpz[0][2] = -px*pz/(p*
p2);
281 dRotMatdpz[1][0] = py/(p*pT)*(1.-pz*pz/p2);
282 dRotMatdpz[1][1] = 0.;
283 dRotMatdpz[1][2] = -py*pz/(p*
p2);
285 dRotMatdpz[2][0] = pT*pz/(p*
p2);
286 dRotMatdpz[2][1] = 0.;
287 dRotMatdpz[2][2] = (1.-pz*pz/
p2)/p;
290 double dpplusdpz = pz*( 0.5/p + c3/c4 );
293 dqplusdpz[0][0] += px*dpplusdpz/
p;
294 dqplusdpz[1][0] += py*dpplusdpz/
p;
295 dqplusdpz[2][0] += pz*dpplusdpz/
p;
298 double dpminusdpz = pz*( 0.5/p - c3/c4 );
301 dqminusdpz[0][0] += px*dpminusdpz/
p;
302 dqminusdpz[1][0] += py*dpminusdpz/
p;
303 dqminusdpz[2][0] += pz*dpminusdpz/
p;
306 return std::make_pair( dqplusdpz, dqminusdpz );
319 double pT2 = px*px + py*
py;
320 double p2 = pT2 + pz*
pz;
322 double sphi =
sin( phi );
323 double cphi =
cos( phi );
324 double stheta =
sin( theta );
325 double ctheta =
cos( theta );
329 double c2 =
sqrt( c1*c1 - 1. );
330 double c3 = -0.5*c2*stheta/
c1;
337 dpplusdtheta[2][0] = c3*c4;
343 dpminusdtheta[2][0] = -c3*c4;
351 return std::make_pair( dqplusdtheta, dqminusdtheta );
363 double sphi =
sin( phi );
364 double cphi =
cos( phi );
365 double stheta =
sin( theta );
369 double c2 =
sqrt( c1*c1 - 1. );
375 dpplusdphi[2][0] = 0.;
381 dpminusdphi[2][0] = 0.;
389 return std::make_pair( dqplusdphi, dqminusdphi );
401 double pT2 = px*px + py*
py;
402 double p2 = pT2 + pz*
pz;
404 double sphi =
sin( phi );
405 double cphi =
cos( phi );
406 double ctheta =
cos( theta );
407 double stheta =
sin( theta );
411 double c2 = 1./
sqrt( c1*c1 - 1. );
416 dpplusdm[0][0] = c2*0.5*c1*stheta*cphi;
417 dpplusdm[1][0] = c2*0.5*c1*stheta*sphi;
422 dpminusdm[0][0] = -dpplusdm[0][0];
423 dpminusdm[1][0] = -dpplusdm[1][0];
424 dpminusdm[2][0] = -dpplusdm[2][0];
432 return std::make_pair( dqplusdm, dqminusdm );
436 const std::pair< AlgebraicMatrix, AlgebraicMatrix >
457 return dqsdm( param );
460 throw cms::Exception(
"BadConfig" ) <<
"@SUB=TwoBodyDecayDerivatives::dqsdzi"
461 <<
"no decay parameter related to selector (" << i <<
").";
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpx(const TwoBodyDecayParameters ¶m) const
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdphi(const TwoBodyDecayParameters ¶m) const
Sin< T >::type sin(const T &t)
TwoBodyDecayDerivatives(double mPrimary=91.1876, double mSecondary=0.105658)
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdzi(const TwoBodyDecayParameters ¶m, 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
Cos< T >::type cos(const T &t)
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdm(const TwoBodyDecayParameters ¶m) const
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpz(const TwoBodyDecayParameters ¶m) const
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdtheta(const TwoBodyDecayParameters ¶m) const
const std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpy(const TwoBodyDecayParameters ¶m) const
const TwoBodyDecayParameters & decayParameters(void) const
const std::pair< AlgebraicMatrix, AlgebraicMatrix > selectedDerivatives(const TwoBodyDecay &tbd, const std::vector< bool > &selector) const
~TwoBodyDecayDerivatives()