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 <stdlib.h>
00024 #include "EvtGenBase/EvtParticle.hh"
00025 #include "EvtGenBase/EvtGenKine.hh"
00026 #include "EvtGenBase/EvtPDL.hh"
00027 #include "EvtGenBase/EvtVector4C.hh"
00028 #include "EvtGenBase/EvtVector4R.hh"
00029 #include "EvtGenBase/EvtTensor4C.hh"
00030 #include "EvtGenBase/EvtDiracParticle.hh"
00031 #include "EvtGenBase/EvtScalarParticle.hh"
00032 #include "EvtGenBase/EvtVectorParticle.hh"
00033 #include "EvtGenBase/EvtTensorParticle.hh"
00034 #include "EvtGenBase/EvtPhotonParticle.hh"
00035 #include "EvtGenBase/EvtNeutrinoParticle.hh"
00036 #include "EvtGenBase/EvtStringParticle.hh"
00037 #include "EvtGenBase/EvtRaritaSchwingerParticle.hh"
00038 #include "EvtGenBase/EvtHighSpinParticle.hh"
00039 #include "EvtGenBase/EvtReport.hh"
00040 #include "EvtGenBase/EvtHelSys.hh"
00041 #include "EvtGenModels/EvtAngSamX.hh"
00042 #include "EvtGenBase/EvtRandom.hh"
00043 #include <string>
00044 using std::endl;
00045
00046 EvtAngSamX::~EvtAngSamX() {}
00047
00048 void EvtAngSamX::getName(std::string& model_name){
00049
00050 model_name="AngSamX";
00051
00052 }
00053
00054 EvtDecayBase* EvtAngSamX::clone(){
00055
00056 return new EvtAngSamX;
00057
00058 }
00059
00060
00061 void EvtAngSamX::init(){
00062
00063 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
00064
00065
00066
00067
00068 if(getNArg()>6 || getNArg()<2 ) {
00069 report(ERROR,"EvtGen")<<" 2<= expected parameters <=6, but it is "<< getNArg()<<endl;
00070 ::abort();
00071 }
00072 }
00073 void EvtAngSamX::initProbMax(){
00074
00075 noProbMax();
00076
00077 }
00078
00079 void EvtAngSamX::decay( EvtParticle *p ){
00080
00081 loop:
00082 p->initializePhaseSpace(getNDaug(),getDaugs());
00083
00084 EvtParticle *v,*s1;
00085 EvtVector4R pv,ps,ppr;
00086
00087 v =p->getDaug(0);
00088 pv=v->getP4();
00089 s1=p->getDaug(1);
00090 ps=s1->getP4();
00091 ppr=p->getP4();
00092
00093 EvtHelSys angles(ppr,pv);
00094 double theta =angles.getHelAng(1);
00095
00096
00097
00098 int narg=getNArg();
00099 double c0,c2,c4,c6,c8,c10;
00100 c0=0;c2=0;c4=0;c6=0;c8=0;c10=0;
00101 if(narg==2) {
00102 c0=getArg(0);
00103 c2=getArg(1);
00104 } else if (narg==3){
00105 c0=getArg(0);
00106 c2=getArg(1);
00107 c4=getArg(2);
00108 } else if (narg==4){
00109 c0=getArg(0);
00110 c2=getArg(1);
00111 c4=getArg(2);
00112 c6=getArg(3);
00113 } else if (narg==5){
00114 c0=getArg(0);
00115 c2=getArg(1);
00116 c4=getArg(2);
00117 c6=getArg(3);
00118 c8=getArg(4);
00119 } else if (narg==6){
00120 c0=getArg(0);
00121 c2=getArg(1);
00122 c4=getArg(2);
00123 c6=getArg(3);
00124 c8=getArg(4);
00125 c10=getArg(5);
00126 }
00127
00128 double costheta= cos(theta);
00129 double costheta2= costheta*costheta;
00130 double costheta4= costheta2*costheta2;
00131 double costheta6= costheta4*costheta2;
00132 double costheta8= costheta6*costheta2;
00133 double costheta10= costheta8*costheta2;
00134
00135 double amp1=(c0+c2*costheta2+c4*costheta4+c6*costheta6+c8*costheta8+c10*costheta10);
00136 double a0,a2,a4,a6,a8,a10;
00137 double ampflag;
00138 if(c0<0) {a0=0;}else{a0=c0;}
00139 if(c2<0) {a2=0;}else{a2=c2;}
00140 if(c4<0) {a4=0;}else{a4=c4;}
00141 if(c6<0) {a6=0;}else{a6=c6;}
00142 if(c8<0) {a8=0;}else{a8=c8;}
00143 if(c10<0) {a10=0;}else{a10=c10;}
00144 ampflag=a0+a2+a4+a6+a8+a10;
00145 if( ampflag<=0 ) {
00146 report(ERROR,"EvtGen")<<" The maxium value of amplitude square should be positive, but it is "<< ampflag<<endl;
00147 ::abort();
00148 }
00149 ampflag = amp1/ampflag;
00150 double rd1=EvtRandom::Flat(0.0, 1.0);
00151 if(rd1>=ampflag) goto loop;
00152 return ;
00153 }
00154
00155