00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <stdlib.h>
00023 #include "EvtGenBase/EvtParticle.hh"
00024 #include "EvtGenBase/EvtRandom.hh"
00025 #include "EvtGenBase/EvtPDL.hh"
00026 #include "EvtGenModels/EvtSingleParticle.hh"
00027 #include "EvtGenBase/EvtReport.hh"
00028 #include <string>
00029 #include "EvtGenBase/EvtConst.hh"
00030 using std::endl;
00031
00032 EvtSingleParticle::~EvtSingleParticle() {}
00033
00034 void EvtSingleParticle::getName(std::string& model_name){
00035
00036 model_name="SINGLE";
00037
00038 }
00039
00040 EvtDecayBase* EvtSingleParticle::clone(){
00041
00042 return new EvtSingleParticle();
00043
00044 }
00045
00046 void EvtSingleParticle::init(){
00047
00048
00049
00050 disableCheckQ();
00051
00052 if ((getNArg()==6)||(getNArg()==4)||(getNArg()==2)) {
00053
00054 if (getNArg()==6){
00055
00056
00057 pmin=getArg(0);
00058 pmax=getArg(1);
00059
00060 cthetamin=getArg(2);
00061 cthetamax=getArg(3);
00062
00063 phimin=getArg(4);
00064 phimax=getArg(5);
00065
00066 }
00067
00068 if (getNArg()==4){
00069
00070
00071 pmin=getArg(0);
00072 pmax=getArg(1);
00073
00074 cthetamin=getArg(2);
00075 cthetamax=getArg(3);
00076
00077 phimin=0.0;
00078 phimax=EvtConst::twoPi;
00079
00080 }
00081
00082 if (getNArg()==2){
00083
00084
00085 pmin=getArg(0);
00086 pmax=getArg(1);
00087
00088 cthetamin=-1.0;
00089 cthetamax=1.0;
00090
00091 phimin=0.0;
00092 phimax=EvtConst::twoPi;
00093
00094 }
00095
00096
00097 }else{
00098
00099 report(ERROR,"EvtGen") << "EvtSingleParticle generator expected "
00100 << " 6, 4, or 2 arguments but found:"<<getNArg()<<endl;
00101 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00102 ::abort();
00103
00104 }
00105
00106
00107 report(INFO,"EvtGen") << "The single particle generator has been configured:"
00108 <<endl;
00109 report(INFO,"EvtGen") << pmax << " > p > " << pmin <<endl;
00110 report(INFO,"EvtGen") << cthetamax << " > costheta > " << cthetamin <<endl;
00111 report(INFO,"EvtGen") << phimax << " > phi > " << phimin <<endl;
00112
00113 }
00114
00115 void EvtSingleParticle::decay( EvtParticle *p ){
00116
00117 EvtParticle *d;
00118 EvtVector4R p4;
00119
00120 double mass=EvtPDL::getMass(getDaug(0));
00121
00122 p->makeDaughters(getNDaug(),getDaugs());
00123 d=p->getDaug(0);
00124
00125
00126
00127
00128 double pcm;
00129 if (pmin==pmax) pcm=pmin;
00130 if (pmin< pmax) pcm=EvtRandom::Flat(pmin,pmax);
00131
00132
00133 double phi;
00134 if(phimin==phimax) phi=phimin;
00135 else phi=EvtRandom::Flat(phimin,phimax);
00136
00137 double cthetalab;
00138
00139 do{
00140
00141 double ctheta;
00142 if(cthetamin==cthetamax) ctheta=cthetamin;
00143 else ctheta=EvtRandom::Flat(cthetamin,cthetamax);
00144 double stheta=sqrt(1.0-ctheta*ctheta);
00145 p4.set(sqrt(mass*mass+pcm*pcm),pcm*cos(phi)*stheta,
00146 pcm*sin(phi)*stheta,pcm*ctheta);
00147
00148 d->init( getDaug(0),p4);
00149
00150
00151 EvtVector4R p4lab=d->getP4Lab();
00152 cthetalab=p4lab.get(3)/p4lab.d3mag();
00153 }while (cthetalab>cthetamax||cthetalab<cthetamin);
00154
00155 return ;
00156 }
00157
00158
00159