EvtDecayBase.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: EvtDecayBase.cc
//
// Description: Store decay parameters for one decay.
//
// Modification history:
//
//    RYD     September 30, 1997         Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <ctype.h>
#include "EvtGenBase/EvtStatus.hh"
#include "EvtGenBase/EvtDecayBase.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtSpinType.hh"
#include <vector>
using std::endl;
using std::fstream;
void EvtDecayBase::checkQ() {
  int i;
  int q=0;
  int qpar;
 
  //If there are no daughters (jetset etc) then we do not
  //want to do this test.  Why?  Because sometimes the parent
  //will have a nonzero charge.
 
  if ( _ndaug != 0) {
    for(i=0; i<_ndaug; i++ ) {
      q += EvtPDL::chg3(_daug[i]);
    }
    qpar = EvtPDL::chg3(_parent);
 
    if ( q != qpar ) {
      report(ERROR,"EvtGen") <<_modelname.c_str()<< " generator expected "
                 << " charge to be conserved, found:"<<endl;
      report(ERROR,"EvtGen") << "Parent charge of "<<(qpar/3)<<endl;
      report(ERROR,"EvtGen") << "Sum of daughter charge of "<<(q/3)<<endl;
      report(ERROR,"EvtGen") << "The parent is "<< EvtPDL::name(_parent).c_str()<<endl;
      for(i=0; i<_ndaug; i++ ) {
      report(ERROR,"EvtGen") << "Daughter "<< EvtPDL::name(_daug[i]).c_str()<<endl;
      }
      report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
      
      ::abort();
    }
  }
}
    
 
double EvtDecayBase::getProbMax( double prob ) {
 
  int i;
 
  //diagnostics
  sum_prob+=prob;
  if (prob>max_prob) max_prob=prob;
 
 
  if ( defaultprobmax && ntimes_prob<=500 ) { 
    //We are building up probmax with this iteration
     ntimes_prob += 1;
     if ( prob > probmax ) { probmax = prob;}
     if (ntimes_prob==500) { 
       probmax*=1.2;
     }
     return 1000000.0*prob;
  }
 
  if ( prob> probmax*1.0001) {
 
    report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")";
    report(INFO,"") << "("<<_modelname.c_str()<<") ";
    report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
    for(i=0;i<_ndaug;i++){
       report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
    }
    report(INFO,"") << endl;
 
    if (defaultprobmax) probmax = prob;
 
  }
 
  ntimes_prob += 1;
 
 
  return probmax;
 
} //getProbMax
 
 
double EvtDecayBase::resetProbMax(double prob) {
  
  report(INFO,"EvtGen") << "Reseting prob max\n"; 
  report(INFO,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")";
  report(INFO,"") << "("<<_modelname.c_str()<<")";
  report(INFO,"") << EvtPDL::getStdHep(_parent)<<"->";
  
  for( int i=0;i<_ndaug;i++){
    report(INFO,"") << EvtPDL::getStdHep(_daug[i]) << " ";
  }
  report(INFO,"") << endl;
  
  probmax = 0.0;
  defaultprobmax = 0;
  ntimes_prob = 0;
  
  return prob;
 
}
 
 
std::string EvtDecayBase::commandName(){
  return std::string("");
}
void EvtDecayBase::command(std::string cmd){
  report(ERROR,"EvtGen") << "Should never call EvtDecayBase::command"<<endl;
  ::abort();
}
 
 
 
void EvtDecayBase::init() {
 
  //This default version of init does nothing;
  //A specialized version of this function can be
  //supplied for each decay model to do initialization.
 
  return;
 
}
 
void EvtDecayBase::initProbMax() {
 
  //This function is called if the decay does not have a
  //specialized initialization.  
  //The default is to set the maximum
  //probability to 0 and the number of times called to 0
  //and defaultprobmax to 1 such that the decay will be 
  //generated many many times
  //in order to generate a reasonable maximum probability
  //for the decay.
 
  defaultprobmax=1;
  ntimes_prob = 0;
  probmax = 0.0;
 
} //initProbMax
 
 
void EvtDecayBase::saveDecayInfo(EvtId ipar, int ndaug, EvtId *daug, 
                 int narg,std::vector<std::string>& args,
                 std::string name,
                 double brfr) {
 
  int i;
 
  _brfr=brfr;
  _ndaug=ndaug;
  _narg=narg;
  _parent=ipar; 
 
  _dsum=0;
 
  if (_ndaug>0) {
    _daug=new EvtId [_ndaug];
    for(i=0;i<_ndaug;i++){
      _daug[i]=daug[i];
      _dsum+=daug[i].getAlias();
    }
  }
  else{
    _daug=0;
  }
 
  if (_narg>0) {
    _args=new std::string[_narg+1];
    for(i=0;i<_narg;i++){
      _args[i]=args[i];
    }
  }
  else{
     _args = 0;
  }
 
  _modelname=name;
 
  this->init();
  this->initProbMax();
 
  if (_chkCharge){
    this->checkQ();
  }
 
 
  if (defaultprobmax){
    report(INFO,"EvtGen") << "No default probmax for ";
    report(INFO,"") << "("<<_modelname.c_str()<<") ";
    report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
    for(i=0;i<_ndaug;i++){
      report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
    }
    report(INFO,"") << endl;
    report(INFO,"") << "This is fine for development, but must be provided for production."<<endl;
    report(INFO,"EvtGen") << "Never fear though - the decay will use the \n";
    report(INFO,"EvtGen") << "500 iterations to build up a good probmax \n";
    report(INFO,"EvtGen") << "before accepting a decay. "<<endl;
  }
 
}
 
EvtDecayBase::EvtDecayBase() {
 
  //the default is that the user module does _not_ set
  // any probmax.
  defaultprobmax=1;
  ntimes_prob = 0;
  probmax = 0.0;
  
  _photos=0;
  _verbose=0;
  _summary=0;
  _parent=EvtId(-1,-1);
  _ndaug=0;
  _narg=0;
  _daug=0;
  _args=0;
  _argsD=0;
  _modelname="**********";
 
  //Default is to check that charge is conserved
  _chkCharge=1;
 
  //statistics collection!
 
  max_prob=0.0;
  sum_prob=0.0;
 
}
 
 
 
void EvtDecayBase::printSummary() {
 
  int i;
 
  if (ntimes_prob>0) {
 
    report(INFO,"EvtGen") << "Calls="<<ntimes_prob<<" eff:"<<
      sum_prob/(probmax*ntimes_prob)<<" frac. max:"<<max_prob/probmax;
    report(INFO,"") <<" probmax:"<<probmax<<" max:"<<max_prob<<" : ";
  }
 
  report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
  for(i=0;i<_ndaug;i++){
    report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
  }
  report(INFO,"") << " ("<<_modelname.c_str()<<"):"<< endl;
  
  
  
}
 
 
EvtDecayBase::~EvtDecayBase() {
 
  if (_daug!=0){
     delete [] _daug;
  }
 
  if (_args!=0){
     delete [] _args;
  }
 
  if (_argsD!=0){
     delete [] _argsD;
  }
 
}
 
void EvtDecayBase::setProbMax(double prbmx){
 
  defaultprobmax=0;
  probmax=prbmx;
 
}
 
void EvtDecayBase::noProbMax(){
 
  defaultprobmax=0;
 
}
 
 
double EvtDecayBase::findMaxMass(EvtParticle *p) {
 
  
  double maxOkMass=EvtPDL::getMaxMass(p->getId());
 
  //protect against vphotons
  if ( maxOkMass < 0.0000000001 ) return 10000000.;
  //and against already determined masses
  if ( p->hasValidP4() ) maxOkMass=p->mass();
 
  EvtParticle *par=p->getParent();
  if ( par ) {
    double maxParMass=findMaxMass(par);
    int i;
    double minDaugMass=0.;
    for(i=0;i<par->getNDaug();i++){
      EvtParticle *dau=par->getDaug(i);
      if ( dau!=p) {
    // it might already have a mass
    if ( dau->isInitialized() || dau->hasValidP4() )
      minDaugMass+=dau->mass();
    else
    //give it a bit of phase space 
      minDaugMass+=1.000001*EvtPDL::getMinMass(dau->getId());
      }
    }
    if ( maxOkMass>(maxParMass-minDaugMass)) maxOkMass=maxParMass-minDaugMass;
  }
  return maxOkMass;
}
 
 
// given list of daughters ( by number ) returns a
// list of viable masses. 
 
void EvtDecayBase::findMass(EvtParticle *p) {
 
  //Need to also check that this mass does not screw
  //up the parent
  //This code assumes that for the ith daughter, 0..i-1
  //already have a mass
  double maxOkMass=findMaxMass(p);
 
  int count=0;
  double mass;
  bool massOk=false;
  int i;
  while (!massOk) { 
    count++;
    if ( count > 10000 ) {
      report(INFO,"EvtGen") << "Can not find a valid mass for: " << EvtPDL::name(p->getId()).c_str() <<endl;
      report(INFO,"EvtGen") << "Now printing parent and/or grandparent tree\n";
      if ( p->getParent() ) {
    if ( p->getParent()->getParent() ) {
      p->getParent()->getParent()->printTree();
      report(INFO,"EvtGen") << p->getParent()->getParent()->mass() <<endl;
      report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
    }
    else{
      p->getParent()->printTree();
      report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
    }
      }
      else  p->printTree();
      report(INFO,"EvtGen") << "maxokmass=" << maxOkMass << " " << EvtPDL::getMinMass(p->getId()) << " " << EvtPDL::getMaxMass(p->getId())<<endl;
      if ( p->getNDaug() ) { 
    for (i=0; i<p->getNDaug(); i++) {
      report(INFO,"EvtGen") << p->getDaug(i)->mass()<<" ";
    }
    report(INFO,"EvtGen") << endl;
      }
      if ( maxOkMass >= EvtPDL::getMinMass(p->getId()) ) {
    report(INFO,"EvtGen") << "taking a default value\n";
    p->setMass(maxOkMass);
    return;
      } 
      assert(0);
    }
    mass = EvtPDL::getMass(p->getId());
    //Just need to check that this mass is > than
    //the mass of all daughters
    double massSum=0.;
    if ( p->getNDaug() ) { 
      for (i=0; i<p->getNDaug(); i++) {
    massSum+= p->getDaug(i)->mass();
      }
    }
    //some special cases are handled with 0 (stable) or 1 (k0->ks/kl) daughters
    if (p->getNDaug()<2)  massOk=true;
    if ( p->getParent() ) {
      if ( p->getParent()->getNDaug()==1 ) massOk=true;
    }
    if ( !massOk ) { 
      if (massSum < mass) massOk=true;
      if ( mass> maxOkMass) massOk=false;
    }
  }
 
  p->setMass(mass);
  
}
 
 
void EvtDecayBase::findMasses(EvtParticle *p, int ndaugs, 
                 EvtId daugs[10], double masses[10]) {
 
  int i;
  double mass_sum;
 
  int count=0;
 
  if (!( p->firstornot() )) {
    for (i = 0; i < ndaugs; i++ ) {
      masses[i] = p->getDaug(i)->mass();
    } //for
  } //if
  else {
    p->setFirstOrNot();
    // if only one daughter do it
 
    if (ndaugs==1) {
      masses[0]=p->mass();
      return;
    }
    
    //until we get a combo whose masses are less than _parent mass.
    do {
      mass_sum = 0.0;
 
      for (i = 0; i < ndaugs; i++ ) {
        masses[i] = EvtPDL::getMass(daugs[i]);
        mass_sum = mass_sum + masses[i];
      } 
 
      count++;
 
     
      if(count==10000) {
        report(ERROR,"EvtGen") <<"Decaying particle:"<<
      EvtPDL::name(p->getId()).c_str()<<" (m="<<p->mass()<<")"<<endl;
        report(ERROR,"EvtGen") <<"To the following daugthers"<<endl;
        for (i = 0; i < ndaugs; i++ ) {
          report(ERROR,"EvtGen") <<  
        EvtPDL::name(daugs[i]).c_str() << endl;
        } 
    report(ERROR,"EvtGen") << "Has been rejected "<<count
                   << " times, will now take minimal masses "
                   << " of daugthers"<<endl;
        
    mass_sum=0.;
    for (i = 0; i < ndaugs; i++ ) {
      masses[i] = EvtPDL::getMinMass(daugs[i]);
      mass_sum = mass_sum + masses[i];
    } 
    if (mass_sum > p->mass()){
      report(ERROR,"EvtGen") << "Parent mass="<<p->mass()
                 << "to light for daugthers."<<endl
                 << "Will throw the event away."<<endl;
      //dont terminate - start over on the event.
      EvtStatus::setRejectFlag();
      mass_sum=0.;
      //      ::abort();
    }
 
      }
    } while ( mass_sum > p->mass());
  } //else
  
  return;
}       
 
void EvtDecayBase::checkNArg(int a1, int a2, int a3, int a4) {
 
  if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) {
    report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "<<endl;
    report(ERROR,"EvtGen") << a1<<endl;; 
    if ( a2>-1) {
      report(ERROR,"EvtGen") << " or " << a2<<endl; 
    }
    if ( a3>-1) {
      report(ERROR,"EvtGen") << " or " << a3<<endl; 
    }
    if ( a4>-1) {
      report(ERROR,"EvtGen") << " or " << a4<<endl; 
    }
    report(ERROR,"EvtGen") << " arguments but found:"<< _narg << endl;
    printSummary();
    report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
    ::abort();
 
  } 
 
}
void EvtDecayBase::checkNDaug(int d1, int d2){
 
  if ( _ndaug != d1 && _ndaug != d2 ) {
    report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected ";
    report(ERROR,"EvtGen") << d1; 
    if ( d2>-1) {
      report(ERROR,"EvtGen") << " or " << d2; 
    }
    report(ERROR,"EvtGen") << " daughters but found:"<< _ndaug << endl;
    printSummary();
    report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
    ::abort();
  } 
 
}
 
void EvtDecayBase::checkSpinParent(EvtSpinType::spintype sp) {
 
  EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
  if ( parenttype != sp ) {
    report(ERROR,"EvtGen") << _modelname.c_str() 
               << " did not get the correct parent spin\n";
    printSummary();
    report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
    ::abort();
  } 
 
}
 
void EvtDecayBase::checkSpinDaughter(int d1, EvtSpinType::spintype sp) {
 
  EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getDaug(d1));
  if ( parenttype != sp ) {
    report(ERROR,"EvtGen") << _modelname.c_str() 
               << " did not get the correct daughter spin d=" 
               << d1 << endl;
    printSummary();
    report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
    ::abort();
  } 
 
}
 
double* EvtDecayBase::getArgs() {
 
  if ( _argsD ) return _argsD;
  //The user has asked for a list of doubles - the arguments 
  //better all be doubles...
  if ( _narg==0 ) return _argsD;
 
  _argsD = new double[_narg];
 
  int i;
  char * tc;
  for(i=0;i<_narg;i++) { 
    _argsD[i] =  strtod(_args[i].c_str(),&tc);
  }
  return _argsD;
}
 
double EvtDecayBase::getArg(int j) {
 
  // Verify string
 
  const char* str = _args[j].c_str();
  int i = 0;
  while(str[i]!=0){
    if (isalpha(str[i]) && str[i]!='e') {
 
      report(INFO,"EvtGen") << "String " << str << " is not a number" << endl;
      assert(0);
    }
    i++;
  }
  
  char** tc=0; 
  return strtod(_args[j].c_str(),tc);
}
 
 
 
 
 
bool EvtDecayBase::matchingDecay(const EvtDecayBase &other) const {
 
  if ( _ndaug != other._ndaug) return false;
  if ( _parent != other._parent) return false;
  
  std::vector<int> useDs;
  for ( unsigned int i=0; i<_ndaug; i++) useDs.push_back(0);
 
  for ( unsigned int i=0; i<_ndaug; i++) {
    bool foundIt=false;
    for ( unsigned int j=0; j<_ndaug; j++) {
      if ( useDs[j] == 1 ) continue;
      if ( _daug[i] == other._daug[j] && _daug[i].getAlias() == other._daug[j].getAlias()) {
    foundIt=true;
    useDs[j]=1;
    break;
      }
    }
    if ( foundIt==false) return false;
  }
  for ( unsigned int i=0; i<_ndaug; i++) if ( useDs[i]==0) return false;
 
  return true;
 
}