TEvtGen/EvtGenBase/EvtKine.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtKine.cc
  12 //
  13 // Description: routines to calculate decay angles.
  14 //
  15 // Modification history:
  16 //
  17 //    DJL/RYD     September 25, 1996         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <math.h>
  23 #include "EvtGenBase/EvtKine.hh"
  24 #include "EvtGenBase/EvtConst.hh"
  25 #include "EvtGenBase/EvtVector4R.hh"
  26 #include "EvtGenBase/EvtVector4C.hh"
  27 #include "EvtGenBase/EvtTensor4C.hh"
  28 #include "EvtGenBase/EvtdFunction.hh"
  29 #include "EvtGenBase/EvtReport.hh"
  30
  31
  32
  33 double EvtDecayAngle(const EvtVector4R& p,const EvtVector4R& q,
  34                      const EvtVector4R& d) {
  35
  36   double pd=p*d;
  37   double pq=p*q;
  38   double qd=q*d;
  39   double mp2=p.mass2();
  40   double mq2=q.mass2();
  41   double md2=d.mass2();
  42
  43   double cost=(pd*mq2-pq*qd)/sqrt((pq*pq-mq2*mp2)*(qd*qd-mq2*md2));
  44
  45   return cost;
  46
  47 }
  48
  49 double EvtDecayAngleChi(const EvtVector4R& p4_p,const EvtVector4R& p4_d1,
  50                         const EvtVector4R& p4_d2,const EvtVector4R& p4_h1,
  51                         const EvtVector4R& p4_h2 ) {
  52
  53   EvtVector4R p4_d1p,p4_h1p,p4_h2p,p4_d2p;
  54
  55
  56   // boost all vectors parent restframe
  57   // This does not boost particle to parent rest frame !!!
  58   // It goes from parents rest frame to frame where parent has given momentum.
  59   p4_d1p=boostTo(p4_d1,p4_p,true);
  60   p4_d2p=boostTo(p4_d2,p4_p,true);
  61   p4_h1p=boostTo(p4_h1,p4_p,true);
  62   p4_h2p=boostTo(p4_h2,p4_p,true);
  63
  64
  65   EvtVector4R d1_perp,d1_prime,h1_perp;
  66   EvtVector4R D;
  67
  68   D=p4_d1p+p4_d2p;
  69
  70   d1_perp=p4_d1p-(D.dot(p4_d1p)/D.dot(D))*D;
  71   h1_perp=p4_h1p-(D.dot(p4_h1p)/D.dot(D))*D;
  72
  73   // orthogonal to both D and d1_perp
  74
  75   d1_prime=D.cross(d1_perp);
  76
  77   d1_perp= d1_perp/d1_perp.d3mag();
  78   d1_prime= d1_prime/d1_prime.d3mag();
  79
  80   double x,y;
  81
  82   x=d1_perp.dot(h1_perp);
  83   y=d1_prime.dot(h1_perp);
  84
  85   double chi=atan2(y,x);
  86
  87   if (chi<0.0) chi+=EvtConst::twoPi;
  88
  89   return chi;
  90
  91 }
  92
  93
  94
  95 double EvtDecayPlaneNormalAngle(const EvtVector4R& p,const EvtVector4R& q,
  96                           const EvtVector4R& d1,const EvtVector4R& d2){
  97
  98   EvtVector4C lc=dual(EvtGenFunctions::directProd(d1,d2)).cont2(q);
  99
 100   EvtVector4R l(real(lc.get(0)),real(lc.get(1)),
 101                 real(lc.get(2)),real(lc.get(3)));
 102
 103   double pq=p*q;
 104
 105   return q.mass()*(p*l)/sqrt(-(pq*pq-p.mass2()*q.mass2())*l.mass2());
 106
 107
 108 }
 109
 110
 111 // Calculate phi using the given 4 vectors (all in the same frame)
 112 double EvtDecayAnglePhi( const EvtVector4R& z, const EvtVector4R& p, const
 113         EvtVector4R& q, const EvtVector4R& d )
 114 {
 115     double eq = (p * q) / p.mass();
 116     double ed = (p * d) / p.mass();
 117     double mq = q.mass();
 118     double q2 = p.mag2r3(q);
 119     double qd = p.dotr3(q,d);
 120     double zq = p.dotr3(z,q);
 121     double zd = p.dotr3(z,d);
 122     double alpha = (eq - mq)/(q2 * mq) * qd - ed/mq;
 123
 124     double y = p.scalartripler3(z,q,d) + alpha * p.scalartripler3(z,q,q);
 125     double x = (zq * (qd + alpha * q2) - q2 * (zd + alpha * zq)) / sqrt(q2);
 126
 127     double phi = atan2(y,x);
 128
 129     return phi<0 ? (phi+EvtConst::twoPi) : phi;
 130 }
 131
 132 EvtComplex wignerD( int j, int m1, int m2, double phi,
 133                     double theta, double gamma )
 134 {
 135
 136     EvtComplex gp(0.0, -phi*m1);
 137     EvtComplex gm(0.0, -gamma*m2);
 138
 139     return exp( gp ) * EvtdFunction::d(j, m1, m2, theta) * exp( gm );
 140 }
 141
 142