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