TEvtGen/EvtGen/EvtGenModels/EvtDMix.cpp

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtDMix.cc
  12 //
  13 // Description: Routine to decay a particle according th phase space
  14 //
  15 // Modification history:
  16 //
  17 //    RYD       January 8, 1997       Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <stdlib.h>
  23 #include "EvtGenBase/EvtParticle.hh"
  24 #include "EvtGenBase/EvtGenKine.hh"
  25 #include "EvtGenBase/EvtPDL.hh"
  26 #include "EvtGenBase/EvtReport.hh"
  27 #include "EvtGenModels/EvtDMix.hh"
  28 #include "EvtGenBase/EvtRandom.hh"
  29 #include <string>
  30
  31 EvtDMix::~EvtDMix() {}
  32
  33 std::string EvtDMix::getName(){
  34
  35   return "DMIX";
  36
  37 }
  38
  39 EvtDecayBase* EvtDMix::clone(){
  40
  41   return new EvtDMix;
  42
  43 }
  44
  45
  46 void EvtDMix::init(){
  47
  48   // check that there are 0 arguments
  49   checkNArg(3);
  50   _rd=getArg(0);
  51   _xpr=getArg(1);
  52   _ypr=getArg(2);
  53 }
  54
  55 void EvtDMix::initProbMax(){
  56
  57   noProbMax();
  58
  59 }
  60
  61 void EvtDMix::decay( EvtParticle *p ){
  62
  63   //unneeded - lange - may13-02
  64   //if ( p->getNDaug() != 0 ) {
  65     //Will end up here because maxrate multiplies by 1.2
  66   //  report(DEBUG,"EvtGen") << "In EvtDMix: has "
  67   //                       <<" daugthers should not be here!"<<endl;
  68   //  return;
  69   //}
  70
  71   p->initializePhaseSpace(getNDaug(),getDaugs());
  72
  73   double ctau=EvtPDL::getctau(p->getId());
  74   if ( ctau==0. ) return;
  75
  76
  77   double pdf, random, gt, weight;
  78
  79   double maxPdf=_rd + sqrt(_rd)*_ypr*50. + 2500.0*(_xpr*_xpr+_ypr*_ypr)/4.0;
  80   bool keepGoing=true;
  81   while ( keepGoing ) {
  82     random=EvtRandom::Flat();
  83     gt=-log(random);
  84     weight=random;
  85     pdf=_rd + sqrt(_rd)*_ypr*gt + gt*gt*(_xpr*_xpr+_ypr*_ypr)/4.0;
  86     pdf*=exp(-1.0*gt);
  87     pdf/=weight;
  88     if ( pdf > maxPdf )
  89       std::cout << pdf << " " << weight << " " << maxPdf << " " << gt << std::endl;
  90     if ( pdf > maxPdf*EvtRandom::Flat()) keepGoing=false;
  91   }
  92
  93
  94   p->setLifetime(gt*ctau);
  95
  96   return ;
  97 }
  98
  99