00001 //-------------------------------------------------------------------------- 00002 // 00003 // Environment: 00004 // This software is part of models developed at BES collaboration 00005 // based on the EvtGen framework. If you use all or part 00006 // of it, please give an appropriate acknowledgement. 00007 // 00008 // Copyright Information: See EvtGen/BesCopyright 00009 // Copyright (A) 2006 Ping Rong-Gang @IHEP 00010 // 00011 // Module: EvtDIY.cc 00012 // 00013 // Description: Model provided by user, see the mannual 00014 // 00015 // Modification history: 00016 // 00017 // Ping R.-G. December, 2006 Module created 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/EvtDIY.hh" 00040 #include "EvtGenBase/EvtRandom.hh" 00041 #include <string> 00042 using namespace std; //::endl; 00043 00044 EvtDIY::~EvtDIY() {} 00045 00046 void EvtDIY::getName(std::string& model_name){ 00047 00048 model_name="DIY"; 00049 00050 } 00051 00052 EvtDecayBase* EvtDIY::clone(){ 00053 00054 return new EvtDIY; 00055 00056 } 00057 00058 00059 void EvtDIY::init(){ 00060 00061 // check that there are 1 arguments:angular distribution parameter 00062 // checkNArg(0); 00063 // checkNDaug(2); 00064 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId()); 00065 // if ( parenttype == EvtSpinType::SCALAR ){ 00066 // report(ERROR,"EvtGen")<<"Scalar decays with flat distribution"<<endl; 00067 // ::abort(); 00068 // } 00069 } 00070 void EvtDIY::initProbMax(){ 00071 noProbMax(); 00072 00073 } 00074 00075 00076 static int nrun=1; 00077 static double max_amps=0.0; 00078 00079 void EvtDIY::decay( EvtParticle *p ){ 00080 00081 double amps,SamAmps,rd1; 00082 // calculated the max amplitude square in 20000 events, is it enough larger? 00083 if(nrun==1){ 00084 int ir,nd; 00085 for(ir=0;ir<=200000;ir++){ 00086 loop0: 00087 p->initializePhaseSpace(getNDaug(),getDaugs()); 00088 int nd=p->getNDaug(),i; 00089 _nd=nd; 00090 for(i=0;i<=nd-1;i++){ 00091 _p4Lab[i]=p->getDaug(i)->getP4Lab(); 00092 _p4CM[i]=p->getDaug(i)->getP4(); 00093 } 00094 amps=AmplitudeSquare(); 00095 if(amps<0) goto loop0; 00096 if(amps>max_amps) max_amps=amps*1.01; 00097 nrun++; 00098 } 00099 } 00100 if(max_amps==0.0) {report(ERROR,"EvtGen")<<"The decay amplitude square should be positive number"<<endl;abort();} 00101 //cout<<"max_amp="<<max_amps<<endl; 00102 00103 00104 loop: 00105 p->initializePhaseSpace(getNDaug(),getDaugs()); 00106 int i; 00107 for(i=0;i<=p->getNDaug()-1;i++){ 00108 _p4Lab[i]=p->getDaug(i)->getP4Lab(); 00109 _p4CM[i]=p->getDaug(i)->getP4(); 00110 } 00111 // Put phase space results into the daughters. 00112 amps=AmplitudeSquare(); 00113 if(amps < 0) goto loop; 00114 SamAmps=amps/max_amps; 00115 rd1=EvtRandom::Flat(0.0, 1.0); 00116 if(rd1>=SamAmps) goto loop; 00117 return ; 00118 } 00119 00120