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/EvtGenKine.hh"
00025 #include "EvtGenBase/EvtPDL.hh"
00026 #include "EvtGenBase/EvtVector4C.hh"
00027 #include "EvtGenBase/EvtVector4R.hh"
00028 #include "EvtGenBase/EvtTensor4C.hh"
00029 #include "EvtGenBase/EvtDiracParticle.hh"
00030 #include "EvtGenBase/EvtScalarParticle.hh"
00031 #include "EvtGenBase/EvtVectorParticle.hh"
00032 #include "EvtGenBase/EvtTensorParticle.hh"
00033 #include "EvtGenBase/EvtPhotonParticle.hh"
00034 #include "EvtGenBase/EvtNeutrinoParticle.hh"
00035 #include "EvtGenBase/EvtStringParticle.hh"
00036 #include "EvtGenBase/EvtRaritaSchwingerParticle.hh"
00037 #include "EvtGenBase/EvtHighSpinParticle.hh"
00038 #include "EvtGenBase/EvtReport.hh"
00039 #include "EvtGenModels/EvtRhoPi.hh"
00040 #include "EvtGenBase/EvtRandom.hh"
00041 #include "EvtGenModels/EvtVVS.hh"
00042 #include <string>
00043 using namespace std;
00044
00045 EvtRhoPi::~EvtRhoPi() {}
00046
00047 void EvtRhoPi::getName(std::string& model_name){
00048
00049 model_name="RhoPi";
00050
00051 }
00052
00053 EvtDecayBase* EvtRhoPi::clone(){
00054
00055 return new EvtRhoPi;
00056
00057 }
00058
00059
00060 void EvtRhoPi::init(){
00061
00062
00063 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
00064
00065
00066
00067
00068 }
00069 void EvtRhoPi::initProbMax(){
00070 noProbMax();
00071
00072 }
00073
00074
00075 static int nrun=1;
00076 static double max_amps=0.0;
00077
00078 void EvtRhoPi::decay( EvtParticle *p ){
00079
00081 double ResonanceMass =0.7710;
00082 double ResonanceWidth=0.1492;
00083 double r1 =0.9;
00084 double phase1 =0;
00085 double r2 =1.1;
00086 double phase2 =1.5;
00087 if(getNArg()>0){
00088 ResonanceMass =getArg(0);
00089 ResonanceWidth=getArg(1);
00090 r1 =getArg(2);
00091 phase1 =getArg(3);
00092 r2 =getArg(4);
00093 phase2 =getArg(5);
00094 }
00095 double amps,SamAmps,rd1;
00096
00097 if(nrun==1){
00098 int ir,nd;
00099 for(ir=0;ir<=20000;ir++){
00100 p->initializePhaseSpace(getNDaug(),getDaugs());
00101 int nd=p->getNDaug(),i;
00102 _nd=nd;
00103 for(i=0;i<=nd-1;i++){
00104 _p4Lab[i]=p->getDaug(i)->getP4Lab();
00105 _p4CM[i]=p->getDaug(i)->getP4();
00106 }
00107 amps=AmplitudeSquare(ResonanceMass,ResonanceWidth,r1,r2,phase1,phase2);
00108 if(amps>max_amps) max_amps=amps*1.01;
00109 nrun++;
00110 }
00111 }
00112 if(max_amps==0.0) {report(ERROR,"EvtGen")<<"The decay amplitude square should be positive number"<<endl;abort();}
00113
00114
00115
00116 loop:
00117 p->initializePhaseSpace(getNDaug(),getDaugs());
00118 int i;
00119 for(i=0;i<=p->getNDaug()-1;i++){
00120 _p4Lab[i]=p->getDaug(i)->getP4Lab();
00121 _p4CM[i]=p->getDaug(i)->getP4();
00122 }
00123
00124 amps=AmplitudeSquare(ResonanceMass,ResonanceWidth,r1,r2,phase1,phase2);
00125 SamAmps=amps/max_amps;
00126 rd1=EvtRandom::Flat(0.0, 1.0);
00127 if(rd1>=SamAmps) goto loop;
00128 return ;
00129 }
00130
00131
00132 double EvtRhoPi::AmplitudeSquare(double ResonanceMass, double ResonanceWidth,double r1,double
00133 r2,double phase1,double phase2){
00134 EvtVector4R dp1=GetDaugMomLab(0);
00135 EvtVector4R dp2=GetDaugMomLab(1);
00136 EvtVector4R dp3=GetDaugMomLab(2);
00137 VVS Jpsi2rhopi(dp1,dp2,dp3,ResonanceMass,ResonanceWidth,r1,r2,phase1,phase2);
00138 return Jpsi2rhopi.amps();
00139
00140 }
00141