EvtDMix.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtDMix.cc
//
// Description: Routine to decay a particle according th phase space
//
// Modification history:
//
//    RYD       January 8, 1997       Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtDMix.hh"
#include "EvtGenBase/EvtRandom.hh"
#include <string>
 
EvtDMix::~EvtDMix() {}
 
void EvtDMix::getName(std::string& model_name){
 
  model_name="DMIX";     
 
}
 
EvtDecayBase* EvtDMix::clone(){
 
  return new EvtDMix;
 
}
 
 
void EvtDMix::init(){
 
  // check that there are 0 arguments
  checkNArg(3);
  _rd=getArg(0);
  _xpr=getArg(1);
  _ypr=getArg(2);
}
 
void EvtDMix::initProbMax(){
 
  noProbMax();
 
}
 
void EvtDMix::decay( EvtParticle *p ){
 
  //unneeded - lange - may13-02
  //if ( p->getNDaug() != 0 ) {
    //Will end up here because maxrate multiplies by 1.2
  //  report(DEBUG,"EvtGen") << "In EvtDMix: has "
  //               <<" daugthers should not be here!"<<endl;
  //  return;
  //}
 
  p->initializePhaseSpace(getNDaug(),getDaugs());
 
  double ctau=EvtPDL::getctau(p->getId());
  if ( ctau==0. ) return;
 
 
  double pdf, random, gt, weight;
 
  double maxPdf=_rd + sqrt(_rd)*_ypr*50. + 2500.0*(_xpr*_xpr+_ypr*_ypr)/4.0;
  bool keepGoing=true;
  while ( keepGoing ) {
    random=EvtRandom::Flat();
    gt=-log(random);
    weight=random;
    pdf=_rd + sqrt(_rd)*_ypr*gt + gt*gt*(_xpr*_xpr+_ypr*_ypr)/4.0;
    pdf*=exp(-1.0*gt);
    pdf/=weight;
    if ( pdf > maxPdf )
      std::cout << pdf << " " << weight << " " << maxPdf << " " << gt << std::endl;
    if ( pdf > maxPdf*EvtRandom::Flat()) keepGoing=false;
  }
 
 
  p->setLifetime(gt*ctau);
 
  return ;
}