EvtAbsLineShape.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: EvtLineShape.cc
//
// Description: Store particle properties for one particle.
//
// Modification history:
//
//    Lange     March 10, 2001        Module created
//
//------------------------------------------------------------------------
//
#include "EvtAbsLineShape.hh"
#include "EvtComplex.hh"
#include "EvtPDL.hh"
#include "EvtPatches.hh"
#include "EvtPropBreitWigner.hh"
#include "EvtRandom.hh"
#include "EvtReport.hh"
#include "EvtTwoBodyVertex.hh"
#include <algorithm>
#include <ctype.h>
#include <fstream>
#include <iostream>
#include <math.h>
#include <stdlib.h>
using std::fstream;
using namespace std; ///////////
 
EvtAbsLineShape::EvtAbsLineShape() {}
 
EvtAbsLineShape::~EvtAbsLineShape() {}
 
EvtAbsLineShape::EvtAbsLineShape( double mass, double width, double maxRange,
                                  EvtSpinType::spintype sp ) {
 
  _includeDecayFact = false;
  _includeBirthFact = false;
  _applyFixForSP8   = false;
  _mass             = mass;
  _width            = width;
  _spin             = sp;
  _maxRange         = maxRange;
  double maxdelta   = 15.0 * width;
  // if ( width>0.001 ) {
  //   if ( 5.0*width < 0.6 ) maxdelta = 0.6;
  // }
  if ( maxRange > 0.00001 )
  {
    _massMax = mass + maxdelta;
    _massMin = mass - maxRange;
  }
  else
  {
    _massMax = mass + maxdelta;
    _massMin = mass - 15.0 * width;
  }
  if ( _massMin < 0. ) _massMin = 0.;
  _massMax = mass + maxdelta;
}
 
EvtAbsLineShape::EvtAbsLineShape( const EvtAbsLineShape& x ) {
 
  _includeDecayFact = x._includeDecayFact;
  _includeBirthFact = x._includeBirthFact;
  _mass             = x._mass;
  _massMax          = x._massMax;
  _massMin          = x._massMin;
  _width            = x._width;
  _spin             = x._spin;
  _maxRange         = x._maxRange;
  _applyFixForSP8   = x._applyFixForSP8;
}
 
EvtAbsLineShape& EvtAbsLineShape::operator=( const EvtAbsLineShape& x ) {
 
  _includeDecayFact = x._includeDecayFact;
  _includeBirthFact = x._includeBirthFact;
  _mass             = x._mass;
  _massMax          = x._massMax;
  _massMin          = x._massMin;
  _width            = x._width;
  _spin             = x._spin;
  _maxRange         = x._maxRange;
  _applyFixForSP8   = x._applyFixForSP8;
  return *this;
}
 
EvtAbsLineShape* EvtAbsLineShape::clone() { return new EvtAbsLineShape( *this ); }
 
double EvtAbsLineShape::rollMass() {
 
  double ymin, ymax;
  double temp;
 
  if ( _width < 0.0001 ) { return _mass; }
  else
  {
    ymin = atan( 2.0 * ( _massMin - _mass ) / _width );
    ymax = atan( 2.0 * ( _massMax - _mass ) / _width );
 
    temp = ( _mass + ( ( _width / 2.0 ) * tan( EvtRandom::Flat( ymin, ymax ) ) ) );
 
    return temp;
  }
}
double EvtAbsLineShape::getRandMass( EvtId* parId, int nDaug, EvtId* dauId, EvtId* othDaugId,
                                     double maxMass, double* dauMasses ) {
 
  if ( _width < 0.0001 ) return _mass;
  // its not flat - but generated according to a BW
 
  if ( maxMass > 0 && maxMass < _massMin )
  {
    report( ERROR, "EvtGen" ) << "In EvtAbsLineShape::getRandMass" << endl;
    report( ERROR, "EvtGen" ) << "Decaying particle " << EvtPDL::name( *parId )
                              << " with mass " << maxMass << endl;
    report( ERROR, "EvtGen" ) << " to particle"
                              << " with a minimal mass of " << _massMin << endl;
  }
 
  double mMin = _massMin;
  double mMax = _massMax;
  if ( maxMass > -0.5 && maxMass < mMax ) mMax = maxMass;
  double ymin = atan( 2.0 * ( mMin - _mass ) / _width );
  double ymax = atan( 2.0 * ( mMax - _mass ) / _width );
 
loop:
  double themass = ( _mass + ( ( _width / 2.0 ) * tan( EvtRandom::Flat( ymin, ymax ) ) ) );
 
  if ( fabs( _addFactorPn ) > 0.00000001 )
  { // addFactorPn, pingrg-2010-1-10
    double massOthD   = EvtPDL::getMeanMass( *othDaugId );
    double massParent = EvtPDL::getMeanMass( *parId );
    double phsp, maxp, maxp1, maxp2;
    if ( themass + massOthD < massParent )
    {
      EvtTwoBodyVertex vb( themass, massOthD, massParent, 1.0 );
      phsp = vb.pD();
    }
    else { return themass; }
 
    if ( ( massOthD + mMax ) < massParent )
    {
      EvtTwoBodyVertex vb1( massOthD, mMax, massParent, 1 );
      EvtTwoBodyVertex vb2( massOthD, mMin, massParent, 1 );
      maxp  = vb1.pD();
      maxp1 = pow( maxp, _addFactorPn * 2.0 );
      maxp  = vb2.pD();
      maxp2 = pow( maxp, _addFactorPn * 2.0 );
      maxp  = max( maxp1, maxp2 );
    }
    else { return themass; }
 
    double wt  = pow( phsp, _addFactorPn * 2.0 ) / maxp;
    double rdm = EvtRandom::Flat( 0.0, 1.0 );
 
    if ( rdm > wt ) goto loop;
  }
 
  return themass;
  //  return EvtRandom::Flat(_massMin,_massMax);
}
 
double EvtAbsLineShape::getMassProb( double mass, double massPar, int nDaug,
                                     double* massDau ) {
 
  double dTotMass = 0.;
  if ( nDaug > 1 )
  {
    int i;
    for ( i = 0; i < nDaug; i++ ) { dTotMass += massDau[i]; }
    // report(INFO,"EvtGen") << mass << " " << massPar << " " << dTotMass << " "<< endl;
    //     if ( (mass-dTotMass)<0.0001 ) return 0.;
    if ( ( mass < dTotMass ) ) return 0.;
  }
  if ( _width < 0.0001 ) return 1.;
 
  // no parent - lets not panic
  if ( massPar > 0.0000000001 )
  {
    if ( mass > massPar ) return 0.;
  }
  // Otherwise return the right value.
  // Fortunately we have generated events according to a non-rel BW, so
  // just return..
  // EvtPropBreitWigner bw(_mass,_width);
  // EvtPropFactor<EvtTwoBodyVertex> f(bw);
  // EvtComplex fm=f.eval(mass);
  // EvtComplex fm0=f.eval(_mass);
  // return (abs(fm)*abs(fm))/(abs(fm0)*abs(fm0));
  return 1.0;
}