TEvtGen/EvtGenBase/EvtKine.cxx

   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
  58   p4_d1p=boostTo(p4_d1,p4_p);
  59   p4_d2p=boostTo(p4_d2,p4_p);
  60   p4_h1p=boostTo(p4_h1,p4_p);
  61   p4_h2p=boostTo(p4_h2,p4_p);
  62
  63
  64   EvtVector4R d1_perp,d1_prime,h1_perp;
  65   EvtVector4R D;
  66
  67   D=p4_d1p+p4_d2p;
  68
  69   d1_perp=p4_d1p-(D.dot(p4_d1p)/D.dot(D))*D;
  70   h1_perp=p4_h1p-(D.dot(p4_h1p)/D.dot(D))*D;
  71
  72   // orthogonal to both D and d1_perp
  73
  74   d1_prime=D.cross(d1_perp);
  75
  76   d1_perp= d1_perp/d1_perp.d3mag();
  77   d1_prime= d1_prime/d1_prime.d3mag();
  78
  79   double x,y;
  80
  81   x=d1_perp.dot(h1_perp);
  82   y=d1_prime.dot(h1_perp);
  83
  84   double chi=atan2(y,x);
  85
  86   if (chi<0.0) chi+=EvtConst::twoPi;
  87
  88   return chi;
  89
  90 }
  91
  92
  93
  94 double EvtDecayPlaneNormalAngle(const EvtVector4R& p,const EvtVector4R& q,
  95                           const EvtVector4R& d1,const EvtVector4R& d2){
  96
  97   EvtVector4C lc=dual(directProd(d1,d2)).cont2(q);
  98
  99   EvtVector4R l(real(lc.get(0)),real(lc.get(1)),
 100                 real(lc.get(2)),real(lc.get(3)));
 101
 102   double pq=p*q;
 103
 104   return q.mass()*(p*l)/sqrt(-(pq*pq-p.mass2()*q.mass2())*l.mass2());
 105
 106
 107 }
 108
 109
 110 // Calculate phi using the given 4 vectors (all in the same frame)
 111 double EvtDecayAnglePhi( const EvtVector4R& z, const EvtVector4R& p, const
 112         EvtVector4R& q, const EvtVector4R& d )
 113 {
 114     double eq = (p * q) / p.mass();
 115     double ed = (p * d) / p.mass();
 116     double mq = q.mass();
 117     double q2 = p.mag2r3(q);
 118     double qd = p.dotr3(q,d);
 119     double zq = p.dotr3(z,q);
 120     double zd = p.dotr3(z,d);
 121     double alpha = (eq - mq)/(q2 * mq) * qd - ed/mq;
 122
 123     double y = p.scalartripler3(z,q,d) + alpha * p.scalartripler3(z,q,q);
 124     double x = (zq * (qd + alpha * q2) - q2 * (zd + alpha * zq)) / sqrt(q2);
 125
 126     double phi = atan2(y,x);
 127
 128     return phi<0 ? (phi+EvtConst::twoPi) : phi;
 129 }
 130
 131 EvtComplex wignerD( int j, int m1, int m2, double phi,
 132                     double theta, double gamma )
 133 {
 134
 135     EvtComplex gp(0.0, -phi*m1);
 136     EvtComplex gm(0.0, -gamma*m2);
 137
 138     return exp( gp ) * EvtdFunction::d(j, m1, m2, theta) * exp( gm );
 139 }
 140
 141