EvtRelBreitWignerBarrierFact.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
//    Dvoretskii June  03, 2002        Reimplemented rollMass()
//
//------------------------------------------------------------------------
#include "EvtPatches.hh"
 
#include "EvtPatches.hh"
 
#include "EvtPredGen.hh"
 
#include "EvtAmpPdf.hh"
#include "EvtIntervalFlatPdf.hh"
#include "EvtMassAmp.hh"
#include "EvtPDL.hh"
#include "EvtPropBreitWignerRel.hh"
#include "EvtRelBreitWignerBarrierFact.hh"
#include "EvtSpinType.hh"
#include "EvtTwoBodyVertex.hh"
#include <algorithm>
EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact() {}
 
EvtRelBreitWignerBarrierFact::~EvtRelBreitWignerBarrierFact() {}
 
EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact( double mass, double width,
                                                            double                maxRange,
                                                            EvtSpinType::spintype sp )
    : EvtAbsLineShape( mass, width, maxRange, sp ) { // double mDaug1, double mDaug2, int l) {
 
  _includeDecayFact = true;
  _includeBirthFact = true;
  _mass             = mass;
  _width            = width;
  _spin             = sp;
  _blatt            = 3.0;
  _maxRange         = maxRange;
  _errorCond        = false;
 
  double maxdelta = 15.0 * width;
 
  if ( maxRange > 0.00001 )
  {
    _massMax = mass + maxdelta;
    _massMin = mass - maxRange;
  }
  else
  {
    _massMax = mass + maxdelta;
    _massMin = mass - 15.0 * width;
  }
 
  _massMax = mass + maxdelta;
  if ( _massMin < 0. ) _massMin = 0.;
}
 
EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(
    const EvtRelBreitWignerBarrierFact& x )
    : EvtAbsLineShape( x ) {
  _massMax          = x._massMax;
  _massMin          = x._massMin;
  _blatt            = x._blatt;
  _maxRange         = x._maxRange;
  _includeDecayFact = x._includeDecayFact;
  _includeBirthFact = x._includeBirthFact;
  _errorCond        = x._errorCond;
}
 
EvtRelBreitWignerBarrierFact&
EvtRelBreitWignerBarrierFact::operator=( const EvtRelBreitWignerBarrierFact& x ) {
  _mass             = x._mass;
  _width            = x._width;
  _spin             = x._spin;
  _massMax          = x._massMax;
  _massMin          = x._massMin;
  _blatt            = x._blatt;
  _maxRange         = x._maxRange;
  _includeDecayFact = x._includeDecayFact;
  _includeBirthFact = x._includeBirthFact;
  _errorCond        = x._errorCond;
 
  return *this;
}
 
EvtAbsLineShape* EvtRelBreitWignerBarrierFact::clone() {
 
  return new EvtRelBreitWignerBarrierFact( *this );
}
 
double EvtRelBreitWignerBarrierFact::getMassProb( double mass, double massPar, int nDaug,
                                                  double* massDau ) {
 
  _errorCond = false;
  // return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);
  if ( nDaug != 2 ) return EvtAbsLineShape::getMassProb( mass, massPar, nDaug, massDau );
 
  double dTotMass = 0.;
 
  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.;
  // report(INFO,"EvtGen") << mass << " " << dTotMass << endl;
  if ( ( mass < dTotMass ) ) return 0.;
 
  if ( _width < 0.0001 ) return 1.;
 
  if ( massPar > 0.0000000001 )
  {
    if ( mass > massPar ) return 0.;
  }
 
  if ( _errorCond ) return 0.;
 
  // we did all the work in getRandMass
  return 1.;
}
 
double EvtRelBreitWignerBarrierFact::getRandMass( EvtId* parId, int nDaug, EvtId* dauId,
                                                  EvtId* othDaugId, double maxMass,
                                                  double* dauMasses ) {
  if ( nDaug != 2 )
  {
    if ( fabs( _addFactorPn ) > 0.00000001 ) EvtAbsLineShape::addFactorPn( _addFactorPn );
    double themass =
        EvtAbsLineShape::getRandMass( parId, nDaug, dauId, othDaugId, maxMass, dauMasses );
    return themass;
  }
 
  if ( _width < 0.00001 ) return _mass;
 
  // first figure out L - take the lowest allowed.
 
  EvtSpinType::spintype spinD1 = EvtPDL::getSpinType( dauId[0] );
  EvtSpinType::spintype spinD2 = EvtPDL::getSpinType( dauId[1] );
 
  int t1 = EvtSpinType::getSpin2( spinD1 );
  int t2 = EvtSpinType::getSpin2( spinD2 );
  int t3 = EvtSpinType::getSpin2( _spin );
 
  int Lmin = -10;
 
  // the user has overridden the partial wave to use.
  for ( int vC = 0; vC < _userSetPW.size(); vC++ )
  {
    if ( dauId[0] == _userSetPWD1[vC] && dauId[1] == _userSetPWD2[vC] )
      Lmin = 2 * _userSetPW[vC];
    if ( dauId[0] == _userSetPWD2[vC] && dauId[1] == _userSetPWD1[vC] )
      Lmin = 2 * _userSetPW[vC];
  }
 
  // allow for special cases.
  if ( Lmin < -1 )
  {
 
    // There are some things I don't know how to deal with
    if ( t3 > 4 )
      return EvtAbsLineShape::getRandMass( parId, nDaug, dauId, othDaugId, maxMass,
                                           dauMasses );
    if ( t1 > 4 )
      return EvtAbsLineShape::getRandMass( parId, nDaug, dauId, othDaugId, maxMass,
                                           dauMasses );
    if ( t2 > 4 )
      return EvtAbsLineShape::getRandMass( parId, nDaug, dauId, othDaugId, maxMass,
                                           dauMasses );
 
    // figure the min and max allowwed "spins" for the daughters state
    Lmin = std::max( t3 - t2 - t1, std::max( t2 - t3 - t1, t1 - t3 - t2 ) );
    if ( Lmin < 0 ) Lmin = 0;
    assert( Lmin == 0 || Lmin == 2 || Lmin == 4 );
  }
 
  // double massD1=EvtPDL::getMeanMass(dauId[0]);
  // double massD2=EvtPDL::getMeanMass(dauId[1]);
  double massD1 = dauMasses[0];
  double massD2 = dauMasses[1];
 
  // I'm not sure how to define the vertex factor here - so retreat to nonRel code.
  if ( ( massD1 + massD2 ) > _mass )
    return EvtAbsLineShape::getRandMass( parId, nDaug, dauId, othDaugId, maxMass, dauMasses );
 
  // parent vertex factor not yet implemented
  double massOthD   = -10.;
  double massParent = -10.;
  int    birthl     = -10;
  if ( othDaugId )
  {
    EvtSpinType::spintype spinOth = EvtPDL::getSpinType( *othDaugId );
    EvtSpinType::spintype spinPar = EvtPDL::getSpinType( *parId );
 
    int tt1 = EvtSpinType::getSpin2( spinOth );
    int tt2 = EvtSpinType::getSpin2( spinPar );
    int tt3 = EvtSpinType::getSpin2( _spin );
 
    // figure the min and max allowwed "spins" for the daughters state
    if ( ( tt1 <= 4 ) && ( tt2 <= 4 ) )
    {
      birthl = std::max( tt3 - tt2 - tt1, std::max( tt2 - tt3 - tt1, tt1 - tt3 - tt2 ) );
      if ( birthl < 0 ) birthl = 0;
 
      massOthD   = EvtPDL::getMeanMass( *othDaugId );
      massParent = EvtPDL::getMeanMass( *parId );
    }
 
    // allow user to override
    for ( size_t vC = 0; vC < _userSetBirthPW.size(); vC++ )
    {
      if ( *othDaugId == _userSetBirthOthD[vC] && *parId == _userSetBirthPar[vC] )
      { birthl = 2 * _userSetBirthPW[vC]; }
    }
  }
  double massM = _massMax;
  if ( ( maxMass > -0.5 ) && ( maxMass < massM ) ) massM = maxMass;
 
  // special case... if the parent mass is _fixed_ we can do a little better
  // and only for a two body decay as that seems to be where we have problems
 
  // Define relativistic propagator amplitude
 
  EvtTwoBodyVertex vd( massD1, massD2, _mass, Lmin / 2 );
  vd.set_f( _blatt );
  EvtPropBreitWignerRel bw( _mass, _width );
  EvtMassAmp            amp( bw, vd );
  if ( _includeDecayFact )
  {
    amp.addDeathFact();
    amp.addDeathFactFF();
  }
 
  if ( fabs( _addFactorPn ) > 0.00000001 )
  {
    //    std::cout<<"EvtRelBreitWignerBarrierFact "<< _addFactorPn<<std::endl;
    amp.addFactorPn( _addFactorPn );
  }
  if ( massParent > -1. )
  {
    if ( _includeBirthFact )
    {
 
      EvtTwoBodyVertex vb( _mass, massOthD, massParent, birthl / 2 );
      if ( _applyFixForSP8 ) vb.set_f( _blatt );
      amp.setBirthVtx( vb );
      amp.addBirthFact();
      amp.addBirthFactFF();
    }
  }
 
  EvtAmpPdf<EvtPoint1D> pdf( amp );
 
  // Estimate maximum and create predicate for accept reject
 
  double tempMaxLoc = _mass;
  if ( maxMass > -0.5 && maxMass < _mass ) tempMaxLoc = maxMass;
  double tempMax = _massMax;
  if ( maxMass > -0.5 && maxMass < _massMax ) tempMax = maxMass;
  double tempMinMass = _massMin;
  if ( massD1 + massD2 > _massMin ) tempMinMass = massD1 + massD2;
 
  // redo sanity check - is there a solution to our problem.
  // if not return an error condition that is caught by the
  // mass prob calculation above.
  if ( tempMinMass > tempMax )
  {
    _errorCond = true;
    return tempMinMass;
  }
 
  if ( tempMaxLoc < tempMinMass ) tempMaxLoc = tempMinMass;
 
  double safetyFactor = 1.2;
  if ( _applyFixForSP8 ) safetyFactor = 1.4;
 
  EvtPdfMax<EvtPoint1D> max( safetyFactor *
                             pdf.evaluate( EvtPoint1D( tempMinMass, tempMax, tempMaxLoc ) ) );
 
  EvtPdfPred<EvtPoint1D> pred( pdf );
  pred.setMax( max );
 
  EvtIntervalFlatPdf                                        flat( tempMinMass, tempMax );
  EvtPdfGen<EvtPoint1D>                                     gen( flat );
  EvtPredGen<EvtPdfGen<EvtPoint1D>, EvtPdfPred<EvtPoint1D>> predgen( gen, pred );
 
  EvtPoint1D point = predgen();
  return point.value();
}