00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include <iostream>
00024 #include <math.h>
00025 #include <fstream>
00026 #include <stdio.h>
00027 #include <stdlib.h>
00028 #include <sys/stat.h>
00029 #include <strstream>
00030 #include "EvtGenBase/EvtParticle.hh"
00031 #include "EvtGenBase/EvtVector4R.hh"
00032 #include "EvtGenBase/EvtVector3R.hh"
00033 #include "EvtGenBase/EvtReport.hh"
00034 #include "EvtGenBase/EvtCPUtil.hh"
00035 #include "EvtGenBase/EvtHelSys.hh"
00036 #include "EvtGenBase/EvtKine.hh"
00037 #include "EvtGenBase/EvtdFunction.hh"
00038 using namespace std;
00039
00040
00041
00042
00043 EvtHelSys::~EvtHelSys() {
00044 }
00045
00046 EvtHelSys::EvtHelSys( const EvtVector4R& p4p, const EvtVector4R& p4d) {
00047 _p4parent=p4p;
00048 _p4daug =p4d;
00049 }
00050
00051 EvtHelSys::EvtHelSys( ) {
00052 }
00053
00054 double EvtHelSys::getHelAng(int i) {
00055 EvtVector4R b_p4p,rp4p, rp4d,boostdaug;
00056 EvtVector3R GetHelAng;
00057 while (_p4parent.d3mag()!=0) {
00058
00059
00060
00061
00062
00063
00064 double theta=Angles(_p4parent,1);
00065 double phi=Angles(_p4parent,2);
00066
00067
00068 rp4p=Helrotate(_p4parent,phi,theta);
00069 rp4d=Helrotate(_p4daug, phi,theta);
00070
00071
00072
00073
00074
00075
00076 _rotatep4p=rp4p;
00077 _rotatep4d=rp4d;
00078 _bst=rp4d;
00079
00080 return Angles(_bst,i);
00081 }
00082 return Angles(_p4daug,i);
00083 }
00084
00085 EvtVector4R EvtHelSys::checkparent(){ return _p4parent;}
00086 EvtVector4R EvtHelSys::checkdaug(){ return _p4daug;}
00087 EvtVector4R EvtHelSys::checkst(int i){
00088 getHelAng(1);
00089 if(i==0) return _bp4p;
00090 if(i==1) return _rotatep4p;
00091 if(i==2) return _rotatep4d;
00092 if(i==3) return _bst;
00093 }
00094
00095
00096 EvtVector4R EvtHelSys::Helrotate(EvtVector4R p1, double phi, double theta){
00097 EvtVector4R Rp;
00098 double cp=cos(phi);
00099 double sp=sin(phi);
00100 double ct=cos(theta);
00101 double st=sin(theta);
00102 double t=p1.get(0),x=p1.get(1), y=p1.get(2), z=p1.get(3);
00103 double xp=x*cp*ct+y*sp*ct-z*st,yp= -x*sp+y*cp,zp=x*cp*st+y*sp*st+z*ct;
00104 Rp.set(t,xp,yp,zp);
00105 return Rp;
00106 }
00107
00108 double EvtHelSys::Angles(EvtVector4R bbst,int i){
00109
00110 double rxy=sqrt(pow(bbst.get(1),2.)+pow(bbst.get(2),2.));
00111 if(bbst.d3mag()<=1e-10) {
00112 if(i==0) return 0.;
00113 else if (i==1) return 0.;
00114 else if (i==2) return 0.;
00115 else {cout<<"Angles(i): i<=2"<<endl;abort();}
00116 }
00117 else if(rxy<=1e-10){
00118 if(i==0) return bbst.d3mag();
00119 if(i==1) return 0.;
00120 if(i==2) return 0.;
00121 else {cout<<"Angles(i): i<=2"<<endl;abort();}
00122 }
00123 else {
00124 double theta=acos(bbst.get(3)/bbst.d3mag());
00125 double csphi=bbst.get(1)/bbst.d3mag()/sin(theta);
00126 if(fabs(csphi)>1.0) csphi=csphi/fabs(csphi);
00127 double phi=acos(csphi);
00128 if(bbst.get(2)<0.0) phi=2*3.1415926-phi;
00129 if(i==0) return bbst.d3mag();
00130 else if (i==1) return theta;
00131 else if (i==2) return phi;
00132 else {cout<<"Angles(i): i<=2"<<endl;abort();}
00133 }
00134
00135 }
00136
00137
00138 double djmn(int j, int m, int n, double theta){
00139 int j2=j*2,m2=m*2,n2=n*2;
00140 double temp=EvtdFunction::d(j2,m2,n2,theta);
00141 return temp;
00142 }
00143
00144 double djmn(double j, double m, double n, double theta){
00145 int j2=(int)(j*2*1.1),m2=(int)(m*2*1.1),n2=(int)(n*2*1.1);
00146 double temp=EvtdFunction::d(j2,m2,n2,theta);
00147 return temp;
00148 }
00149
00150
00151 EvtComplex Djmn(int j, int m, int n, double phi,double theta, double gamma){
00152 int j2=j*2,m2=m*2,n2=n*2;
00153 EvtComplex gp(cos(-phi*m ), -sin(phi*m));
00154 EvtComplex gm(cos(-gamma*n), -sin(gamma*n));
00155 double tp3=EvtdFunction::d(j2,m2,n2,theta);
00156
00157
00158
00159 EvtComplex temp=gp * tp3 * gm;
00160
00161 return temp;
00162 }
00163
00164
00165 EvtComplex Djmn(double j, double m, double n, double phi,double theta, double gamma){
00166 int j2=(int)(j*2*1.1),m2=(int)(m*2*1.1),n2=(int)(n*2*1.1);
00167 EvtComplex gp(cos(-phi*m ), -sin(phi*m));
00168 EvtComplex gm(cos(-gamma*n), -sin(gamma*n));
00169 double tp3=EvtdFunction::d(j2,m2,n2,theta);
00170 EvtComplex temp=gp * tp3 * gm;
00171 return temp;
00172 }
00173