EvtHelAmp.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: EvtHelAmp.cc
//
// Description: Decay model for implementation of generic 2 body
//              decay specified by the helicity amplitudes
//
//
// Modification history:
//
//    RYD       March 14, 1999       Module created
//
//------------------------------------------------------------------------
//
#include "EvtHelAmp.hh"
#include "../EvtGenBase/EvtConst.hh"
#include "../EvtGenBase/EvtEvalHelAmp.hh"
#include "../EvtGenBase/EvtGenKine.hh"
#include "../EvtGenBase/EvtId.hh"
#include "../EvtGenBase/EvtPDL.hh"
#include "../EvtGenBase/EvtParticle.hh"
#include "../EvtGenBase/EvtPatches.hh"
#include "../EvtGenBase/EvtReport.hh"
#include <stdlib.h>
#include <string>
using std::endl;
 
EvtHelAmp::~EvtHelAmp() { delete _evalHelAmp; }
 
void EvtHelAmp::getName( std::string& model_name ) { model_name = "HELAMP"; }
 
EvtDecayBase* EvtHelAmp::clone() { return new EvtHelAmp; }
 
void EvtHelAmp::init() {
 
  checkNDaug( 2 );
 
  // find out how many states each particle have
  int _nA = EvtSpinType::getSpinStates( EvtPDL::getSpinType( getParentId() ) );
  int _nB = EvtSpinType::getSpinStates( EvtPDL::getSpinType( getDaug( 0 ) ) );
  int _nC = EvtSpinType::getSpinStates( EvtPDL::getSpinType( getDaug( 1 ) ) );
 
  if ( verbose() )
  { report( INFO, "EvtGen" ) << "_nA,_nB,_nC:" << _nA << "," << _nB << "," << _nC << endl; }
 
  // find out what 2 times the spin is
  int _JA2 = EvtSpinType::getSpin2( EvtPDL::getSpinType( getParentId() ) );
  int _JB2 = EvtSpinType::getSpin2( EvtPDL::getSpinType( getDaug( 0 ) ) );
  int _JC2 = EvtSpinType::getSpin2( EvtPDL::getSpinType( getDaug( 1 ) ) );
 
  if ( verbose() )
  {
    report( INFO, "EvtGen" ) << "_JA2,_JB2,_JC2:" << _JA2 << "," << _JB2 << "," << _JC2
                             << endl;
  }
 
  // allocate memory
  int* _lambdaA2 = new int[_nA];
  int* _lambdaB2 = new int[_nB];
  int* _lambdaC2 = new int[_nC];
 
  EvtComplexPtr* _HBC = new EvtComplexPtr[_nB];
  int            ib, ic;
  for ( ib = 0; ib < _nB; ib++ ) { _HBC[ib] = new EvtComplex[_nC]; }
 
  int i;
  // find the allowed helicities (actually 2*times the helicity!)
  fillHelicity( _lambdaA2, _nA, _JA2, getParentId() );
  fillHelicity( _lambdaB2, _nB, _JB2, getDaug( 0 ) );
  fillHelicity( _lambdaC2, _nC, _JC2, getDaug( 1 ) );
 
  if ( verbose() )
  {
    report( INFO, "EvtGen" ) << "Helicity states of particle A:" << endl;
    for ( i = 0; i < _nA; i++ ) { report( INFO, "EvtGen" ) << _lambdaA2[i] << endl; }
 
    report( INFO, "EvtGen" ) << "Helicity states of particle B:" << endl;
    for ( i = 0; i < _nB; i++ ) { report( INFO, "EvtGen" ) << _lambdaB2[i] << endl; }
 
    report( INFO, "EvtGen" ) << "Helicity states of particle C:" << endl;
    for ( i = 0; i < _nC; i++ ) { report( INFO, "EvtGen" ) << _lambdaC2[i] << endl; }
  }
 
  // now read in the helicity amplitudes
 
  int argcounter = 0;
 
  for ( ib = 0; ib < _nB; ib++ )
  {
    for ( ic = 0; ic < _nC; ic++ )
    {
      _HBC[ib][ic] = 0.0;
      if ( abs( _lambdaB2[ib] - _lambdaC2[ic] ) <= _JA2 ) argcounter += 2;
    }
  }
 
  checkNArg( argcounter );
 
  argcounter = 0;
 
  for ( ib = 0; ib < _nB; ib++ )
  {
    for ( ic = 0; ic < _nC; ic++ )
    {
      if ( abs( _lambdaB2[ib] - _lambdaC2[ic] ) <= _JA2 )
      {
        _HBC[ib][ic] =
            getArg( argcounter ) * exp( EvtComplex( 0.0, getArg( argcounter + 1 ) ) );
        ;
        argcounter += 2;
        if ( verbose() )
        {
          report( INFO, "EvtGen" )
              << "_HBC[" << ib << "][" << ic << "]=" << _HBC[ib][ic] << endl;
        }
      }
    }
  }
 
  _evalHelAmp = new EvtEvalHelAmp( EvtPDL::getSpinType( getParentId() ),
                                   EvtPDL::getSpinType( getDaug( 0 ) ),
                                   EvtPDL::getSpinType( getDaug( 1 ) ), _HBC );
 
  // Note: these are not class data members but local variables.
  delete[] _lambdaA2;
  delete[] _lambdaB2;
  delete[] _lambdaC2;
  for ( int ib = 0; ib < _nB; ib++ ) { delete[] _HBC[ib]; }
  delete[] _HBC; // _HBC is copied in ctor of EvtEvalHelAmp above.
}
 
void EvtHelAmp::initProbMax() {
 
  double maxprob = _evalHelAmp->probMax();
 
  if ( verbose() ) { report( INFO, "EvtGen" ) << "Calculated probmax" << maxprob << endl; }
 
  setProbMax( maxprob );
}
 
void EvtHelAmp::decay( EvtParticle* p ) {
 
  // first generate simple phase space
  p->initializePhaseSpace( getNDaug(), getDaugs() );
 
  _evalHelAmp->evalAmp( p, _amp2 );
 
  return;
}
 
void EvtHelAmp::fillHelicity( int* lambda2, int n, int J2, EvtId id ) {
 
  int i;
 
  // photon is special case!
  if ( n == 2 && J2 == 2 )
  {
    lambda2[0] = 2;
    lambda2[1] = -2;
    return;
  }
 
  // and so is the neutrino!
  if ( n == 1 && J2 == 1 )
  {
    if ( EvtPDL::getStdHep( id ) > 0 )
    {
      // particle i.e. lefthanded
      lambda2[0] = -1;
    }
    else
    {
      // anti particle i.e. righthanded
      lambda2[0] = 1;
    }
    return;
  }
 
  assert( n == J2 + 1 );
 
  for ( i = 0; i < n; i++ ) { lambda2[i] = n - i * 2 - 1; }
 
  return;
}