TEvtGen/EvtGenBase/EvtDecayProb.cxx

   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: EvtGen/EvtDecayProb.cc
  12 //
  13 // Description:
  14 //
  15 // Modification history:
  16 //
  17 //    DJL/RYD     August 11, 1998         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 #include "EvtGenBase/EvtPatches.hh"
  21
  22 #include "EvtGenBase/EvtDecayBase.hh"
  23 #include "EvtGenBase/EvtDecayProb.hh"
  24 #include "EvtGenBase/EvtParticle.hh"
  25 #include "EvtGenBase/EvtRadCorr.hh"
  26 #include "EvtGenBase/EvtRandom.hh"
  27 #include "EvtGenBase/EvtPDL.hh"
  28 #include "EvtGenBase/EvtReport.hh"
  29 using std::endl;
  30
  31 void EvtDecayProb::makeDecay(EvtParticle* p, bool recursive){
  32
  33   int ntimes=10000;
  34
  35   double dummy;
  36
  37   do{
  38     _weight=1.0;
  39     _daugsDecayedByParentModel=false;
  40
  41     decay(p);
  42
  43     ntimes--;
  44
  45     _prob = _prob/_weight;
  46
  47     dummy=getProbMax(_prob)*EvtRandom::Flat();
  48     p->setDecayProb(_prob/getProbMax(_prob));
  49
  50   }while(ntimes&&(_prob<dummy));
  51
  52   if (ntimes==0){
  53     report(DEBUG,"EvtGen") << "Tried accept/reject:10000"
  54                            <<" times, and rejected all the times!"<<endl;
  55     report(DEBUG,"EvtGen") << "Is therefore accepting the last event!"<<endl;
  56     report(DEBUG,"EvtGen") << "Decay of particle:"<<
  57       EvtPDL::name(p->getId()).c_str()<<"(channel:"<<
  58       p->getChannel()<<") with mass "<<p->mass()<<endl;
  59
  60     for(size_t ii=0;ii<p->getNDaug();ii++){
  61       report(DEBUG,"EvtGen") <<"Daughter "<<ii<<":"<<
  62         EvtPDL::name(p->getDaug(ii)->getId()).c_str()<<" with mass "<<
  63         p->getDaug(ii)->mass()<<endl;
  64     }
  65   }
  66
  67
  68   EvtSpinDensity rho;
  69   rho.setDiag(p->getSpinStates());
  70   p->setSpinDensityBackward(rho);
  71   if (getPHOTOS() || EvtRadCorr::alwaysRadCorr()) {
  72     EvtRadCorr::doRadCorr(p);
  73   }
  74
  75   if(!recursive) return;
  76
  77   //Now decay the daughters.
  78   if ( !daugsDecayedByParentModel()) {
  79     for(size_t i=0;i<p->getNDaug();i++){
  80       //Need to set the spin density of the daughters to be
  81       //diagonal.
  82       rho.setDiag(p->getDaug(i)->getSpinStates());
  83       p->getDaug(i)->setSpinDensityForward(rho);
  84
  85       //Now decay the daughter.  Really!
  86       p->getDaug(i)->decay();
  87     }
  88   }
  89
  90 }
  91
  92
  93
  94