TEvtGen/EvtGenBase/EvtPhotonParticle.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: EvtPhotonParticle.cc
  12 //
  13 // Description: Class to describe massless vectors
  14 //
  15 // Modification history:
  16 //
  17 //    DJL/RYD   September 25, 1996           Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <stdlib.h>
  23 #include <iostream>
  24 #include <math.h>
  25 #include "EvtGenBase/EvtComplex.hh"
  26 #include "EvtGenBase/EvtPhotonParticle.hh"
  27 #include "EvtGenBase/EvtVector4C.hh"
  28 #include "EvtGenBase/EvtReport.hh"
  29 using std::endl;
  30
  31 EvtPhotonParticle::~EvtPhotonParticle(){}
  32
  33 void EvtPhotonParticle::init(EvtId part_n,const EvtVector4R& p4){
  34
  35   init(part_n,p4.get(0),p4.get(1),
  36        p4.get(2),p4.get(3));
  37 }
  38
  39 void EvtPhotonParticle::init(EvtId part_n,double e,double px,double py,double pz){
  40
  41   _validP4=true;
  42   setp(e,px,py,pz);
  43   setpart_num(part_n);
  44
  45   setLifetime();
  46
  47   //defere calculation of basis vectors untill they are needed!
  48   _evalBasis=0;
  49
  50 }
  51
  52
  53 EvtVector4C EvtPhotonParticle::epsParentPhoton(int i){
  54
  55   if (!_evalBasis){
  56
  57     _evalBasis=1;
  58     eps1.set(EvtComplex(0.0,0.0), EvtComplex(-1.0/sqrt(2.0),0.0),
  59              EvtComplex(0.0,-1.0/sqrt(2.0)), EvtComplex(0.0,0.0));
  60     eps2.set(EvtComplex(0.0,0.0), EvtComplex(1.0/sqrt(2.0),0.0),
  61              EvtComplex(0.0,-1.0/sqrt(2.0)), EvtComplex(0.0,0.0));
  62
  63     // These are for photon along z axis.  Rotate to get
  64     // correct direction...
  65
  66     double phi,theta;
  67
  68     EvtVector4R p=this->getP4();
  69
  70     double px=p.get(1);
  71     double py=p.get(2);
  72     double pz=p.get(3);
  73
  74     phi = atan2(py,px);
  75     theta = acos(pz/sqrt(px*px+py*py+pz*pz));
  76     eps1.applyRotateEuler(phi,theta,-phi);
  77     eps2.applyRotateEuler(phi,theta,-phi);
  78
  79   }
  80
  81
  82   EvtVector4C temp;
  83
  84   switch(i) {
  85
  86   case 0:
  87     temp=eps1;
  88     break;
  89   case 1:
  90     temp=eps2;
  91     break;
  92   default:
  93     report(ERROR,"EvtGen") << "EvtPhotonParticle.cc: Asked "
  94                            << "for state:"<<i<<endl;
  95     ::abort();
  96     break;
  97   }
  98
  99   return temp;
 100 }
 101
 102 EvtVector4C EvtPhotonParticle::epsPhoton(int ){
 103
 104   report(ERROR,"EvtGen") << "EvtPhotonParticle.cc: Can not get "
 105                          << "state in photons restframe."<<endl;;
 106   ::abort();
 107   return EvtVector4C();
 108
 109 }
 110
 111
 112 EvtSpinDensity EvtPhotonParticle::rotateToHelicityBasis() const {
 113
 114   EvtVector4C eplus(0.0,-1.0/sqrt(2.0),
 115                     EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
 116   EvtVector4C eminus(0.0,1.0/sqrt(2.0),
 117                      EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
 118
 119   //Really uggly have to cast away constness because the
 120   //function epsParentPhoton caches the state vectors...
 121   EvtVector4C e1=((EvtParticle*)this)->epsParentPhoton(0);
 122   EvtVector4C e2=((EvtParticle*)this)->epsParentPhoton(1);
 123
 124
 125   EvtSpinDensity R;
 126   R.setDim(2);
 127
 128   R.set(0,0,(eplus.conj())*e1);
 129   R.set(0,1,(eplus.conj())*e2);
 130
 131   R.set(1,0,(eminus.conj())*e1);
 132   R.set(1,1,(eminus.conj())*e2);
 133
 134   return R;
 135
 136 }
 137
 138
 139 EvtSpinDensity EvtPhotonParticle::rotateToHelicityBasis(double alpha,
 140                                                         double beta,
 141                                                         double gamma) const{
 142
 143   EvtVector4C eplus(0.0,-1.0/sqrt(2.0),
 144                     EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
 145   EvtVector4C eminus(0.0,1.0/sqrt(2.0),
 146                      EvtComplex(0.0,-1.0/sqrt(2.0)),0.0);
 147
 148   eplus.applyRotateEuler(alpha,beta,gamma);
 149   eminus.applyRotateEuler(alpha,beta,gamma);
 150
 151
 152   //Really uggly have to cast away constness because the
 153   //function epsParentPhoton caches the state vectors...
 154   EvtVector4C e1=((EvtParticle*)this)->epsParentPhoton(0);
 155   EvtVector4C e2=((EvtParticle*)this)->epsParentPhoton(1);
 156
 157   EvtSpinDensity R;
 158   R.setDim(2);
 159
 160   R.set(0,0,(eplus.conj())*e1);
 161   R.set(0,1,(eplus.conj())*e2);
 162
 163   R.set(1,0,(eminus.conj())*e1);
 164   R.set(1,1,(eminus.conj())*e2);
 165
 166   return R;
 167
 168 }
 169
 170