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/EvtAngSamLab.hh"
00042 #include "EvtGenBase/EvtRandom.hh"
00043 #include <string>
00044 using std::endl;
00045
00046 EvtAngSamLab::~EvtAngSamLab() {}
00047
00048 void EvtAngSamLab::getName(std::string& model_name){
00049
00050 model_name="AngSamLab";
00051
00052 }
00053
00054 EvtDecayBase* EvtAngSamLab::clone(){
00055
00056 return new EvtAngSamLab;
00057
00058 }
00059
00060
00061 void EvtAngSamLab::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 EvtAngSamLab::initProbMax(){
00076
00077 noProbMax();
00078
00079 }
00080
00081 void EvtAngSamLab::decay( EvtParticle *p ){
00082
00083 loop:
00084 p->initializePhaseSpace(getNDaug(),getDaugs());
00085
00086 EvtParticle *v,*s1;
00087 EvtVector4R pv;
00088
00089 v =p->getDaug(0);
00090 s1=p->getDaug(1);
00091 pv=v->getP4Lab();
00092
00093
00094 double theta =pv.theta();
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