00001 //-------------------------------------------------------------------------- 00002 // 00003 // Environment: 00004 // This software is part of the EvtGen package developed jointly 00005 // for the BaBar and CLEO collaborations. If you use all or part 00006 // of it, please give an appropriate acknowledgement. 00007 // 00008 // Copyright Information: See EvtGen/COPYRIGHT 00009 // Copyright (C) 1998 Caltech, UCSB 00010 // 00011 // Module: EvtDMix.cc 00012 // 00013 // Description: Routine to decay a particle according th phase space 00014 // 00015 // Modification history: 00016 // 00017 // RYD January 8, 1997 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/EvtReport.hh" 00027 #include "EvtGenModels/EvtDMix.hh" 00028 #include "EvtGenBase/EvtRandom.hh" 00029 #include <string> 00030 00031 EvtDMix::~EvtDMix() {} 00032 00033 void EvtDMix::getName(std::string& model_name){ 00034 00035 model_name="DMIX"; 00036 00037 } 00038 00039 EvtDecayBase* EvtDMix::clone(){ 00040 00041 return new EvtDMix; 00042 00043 } 00044 00045 00046 void EvtDMix::init(){ 00047 00048 // check that there are 0 arguments 00049 checkNArg(3); 00050 _rd=getArg(0); 00051 _xpr=getArg(1); 00052 _ypr=getArg(2); 00053 } 00054 00055 void EvtDMix::initProbMax(){ 00056 00057 noProbMax(); 00058 00059 } 00060 00061 void EvtDMix::decay( EvtParticle *p ){ 00062 00063 //unneeded - lange - may13-02 00064 //if ( p->getNDaug() != 0 ) { 00065 //Will end up here because maxrate multiplies by 1.2 00066 // report(DEBUG,"EvtGen") << "In EvtDMix: has " 00067 // <<" daugthers should not be here!"<<endl; 00068 // return; 00069 //} 00070 00071 p->initializePhaseSpace(getNDaug(),getDaugs()); 00072 00073 double ctau=EvtPDL::getctau(p->getId()); 00074 if ( ctau==0. ) return; 00075 00076 00077 double pdf, random, gt, weight; 00078 00079 double maxPdf=_rd + sqrt(_rd)*_ypr*50. + 2500.0*(_xpr*_xpr+_ypr*_ypr)/4.0; 00080 bool keepGoing=true; 00081 while ( keepGoing ) { 00082 random=EvtRandom::Flat(); 00083 gt=-log(random); 00084 weight=random; 00085 pdf=_rd + sqrt(_rd)*_ypr*gt + gt*gt*(_xpr*_xpr+_ypr*_ypr)/4.0; 00086 pdf*=exp(-1.0*gt); 00087 pdf/=weight; 00088 if ( pdf > maxPdf ) 00089 std::cout << pdf << " " << weight << " " << maxPdf << " " << gt << std::endl; 00090 if ( pdf > maxPdf*EvtRandom::Flat()) keepGoing=false; 00091 } 00092 00093 00094 p->setLifetime(gt*ctau); 00095 00096 return ; 00097 } 00098 00099