EvtEvalHelAmp.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) 2000      Caltech, UCSB
//
// Module: EvtHelAmp.cc
//
// Description: Decay model for implementation of generic 2 body
//              decay specified by the partial wave amplitudes
//
//
// Modification history:
//
//    fkw        February 2, 2001     changes to satisfy KCC
//    RYD       September 7, 2000       Module created
//    Ping RG   Apr.  15,2007     change to helicity angle(@Line=210) for sequential decays
//------------------------------------------------------------------------
//
#include "EvtEvalHelAmp.hh"
#include "EvtAmp.hh"
#include "EvtConst.hh"
#include "EvtHelSys.hh"
#include "EvtId.hh"
#include "EvtPDL.hh"
#include "EvtParticle.hh"
#include "EvtPatches.hh"
#include "EvtReport.hh"
#include "EvtVector4C.hh"
#include "EvtdFunction.hh"
#include <stdlib.h>
 
using std::endl;
 
EvtEvalHelAmp::~EvtEvalHelAmp() {
 
  // allocate memory
  delete[] _lambdaA2;
  delete[] _lambdaB2;
  delete[] _lambdaC2;
 
  int ia, ib, ic;
  for ( ib = 0; ib < _nB; ib++ ) { delete[] _HBC[ib]; }
 
  delete[] _HBC;
 
  for ( ia = 0; ia < _nA; ia++ ) { delete[] _RA[ia]; }
  delete[] _RA;
 
  for ( ib = 0; ib < _nB; ib++ ) { delete[] _RB[ib]; }
  delete[] _RB;
 
  for ( ic = 0; ic < _nC; ic++ ) { delete[] _RC[ic]; }
  delete[] _RC;
 
  for ( ia = 0; ia < _nA; ia++ )
  {
    for ( ib = 0; ib < _nB; ib++ )
    {
      delete[] _amp[ia][ib];
      delete[] _amp1[ia][ib];
      delete[] _amp3[ia][ib];
    }
    delete[] _amp[ia];
    delete[] _amp1[ia];
    delete[] _amp3[ia];
  }
 
  delete[] _amp;
  delete[] _amp1;
  delete[] _amp3;
}
 
EvtEvalHelAmp::EvtEvalHelAmp( EvtSpinType::spintype typeA, EvtSpinType::spintype typeB,
                              EvtSpinType::spintype typeC, EvtComplexPtrPtr HBC ) {
 
  // find out how many states each particle have
  _nA = EvtSpinType::getSpinStates( typeA );
  _nB = EvtSpinType::getSpinStates( typeB );
  _nC = EvtSpinType::getSpinStates( typeC );
 
  // find out what 2 times the spin is
  _JA2 = EvtSpinType::getSpin2( typeA );
  _JB2 = EvtSpinType::getSpin2( typeB );
  _JC2 = EvtSpinType::getSpin2( typeC );
 
  // allocate memory
  _lambdaA2 = new int[_nA];
  _lambdaB2 = new int[_nB];
  _lambdaC2 = new int[_nC];
 
  _HBC = new EvtComplexPtr[_nB];
  int ia, ib, ic;
  for ( ib = 0; ib < _nB; ib++ ) { _HBC[ib] = new EvtComplex[_nC]; }
 
  _RA = new EvtComplexPtr[_nA];
  for ( ia = 0; ia < _nA; ia++ ) { _RA[ia] = new EvtComplex[_nA]; }
  _RB = new EvtComplexPtr[_nB];
  for ( ib = 0; ib < _nB; ib++ ) { _RB[ib] = new EvtComplex[_nB]; }
  _RC = new EvtComplexPtr[_nC];
  for ( ic = 0; ic < _nC; ic++ ) { _RC[ic] = new EvtComplex[_nC]; }
 
  _amp  = new EvtComplexPtrPtr[_nA];
  _amp1 = new EvtComplexPtrPtr[_nA];
  _amp3 = new EvtComplexPtrPtr[_nA];
  for ( ia = 0; ia < _nA; ia++ )
  {
    _amp[ia]  = new EvtComplexPtr[_nB];
    _amp1[ia] = new EvtComplexPtr[_nB];
    _amp3[ia] = new EvtComplexPtr[_nB];
    for ( ib = 0; ib < _nB; ib++ )
    {
      _amp[ia][ib]  = new EvtComplex[_nC];
      _amp1[ia][ib] = new EvtComplex[_nC];
      _amp3[ia][ib] = new EvtComplex[_nC];
    }
  }
 
  // find the allowed helicities (actually 2*times the helicity!)
 
  fillHelicity( _lambdaA2, _nA, _JA2 );
  fillHelicity( _lambdaB2, _nB, _JB2 );
  fillHelicity( _lambdaC2, _nC, _JC2 );
 
  for ( ib = 0; ib < _nB; ib++ )
  {
    for ( ic = 0; ic < _nC; ic++ ) { _HBC[ib][ic] = HBC[ib][ic]; }
  }
}
 
double EvtEvalHelAmp::probMax() {
 
  double c = 1.0 / sqrt( 4 * EvtConst::pi / ( _JA2 + 1 ) );
 
  int ia, ib, ic;
 
  double theta;
  int    itheta;
 
  double maxprob = 0.0;
 
  for ( itheta = -10; itheta <= 10; itheta++ )
  {
    theta = acos( 0.099999 * itheta );
    for ( ia = 0; ia < _nA; ia++ )
    {
      double prob = 0.0;
      for ( ib = 0; ib < _nB; ib++ )
      {
        for ( ic = 0; ic < _nC; ic++ )
        {
          _amp[ia][ib][ic] = 0.0;
          if ( abs( _lambdaB2[ib] - _lambdaC2[ic] ) <= _JA2 )
          {
            _amp[ia][ib][ic] =
                c * _HBC[ib][ic] *
                EvtdFunction::d( _JA2, _lambdaA2[ia], _lambdaB2[ib] - _lambdaC2[ic], theta );
            prob += real( _amp[ia][ib][ic] * conj( _amp[ia][ib][ic] ) );
          }
        }
      }
 
      prob *= sqrt( 1.0 * _nA );
 
      if ( prob > maxprob ) maxprob = prob;
    }
  }
 
  return maxprob;
}
 
void EvtEvalHelAmp::evalAmp( EvtParticle* p, EvtAmp& amp ) {
 
  // find theta and phi of the first daughter
 
  EvtVector4R pB = p->getDaug( 0 )->getP4();
  EvtVector4R pC = p->getDaug( 1 )->getP4();
  EvtVector4R pP = pB + pC;
 
  EvtHelSys angles( pP, pB );
  double    theta = angles.getHelAng( 1 );
  double    phi   = angles.getHelAng( 2 );
  //  double theta=acos(pB.get(3)/pB.d3mag());  //pingrg ,2007,4,15
  //  double phi=atan2(pB.get(2),pB.get(1));
 
  double c = sqrt( ( _JA2 + 1 ) / ( 4 * EvtConst::pi ) );
 
  int ia, ib, ic;
 
  double prob1 = 0.0;
 
  for ( ia = 0; ia < _nA; ia++ )
  {
    for ( ib = 0; ib < _nB; ib++ )
    {
      for ( ic = 0; ic < _nC; ic++ )
      {
        _amp[ia][ib][ic] = 0.0;
        if ( abs( _lambdaB2[ib] - _lambdaC2[ic] ) <= _JA2 )
        {
          double dfun =
              EvtdFunction::d( _JA2, _lambdaA2[ia], _lambdaB2[ib] - _lambdaC2[ic], theta );
 
          _amp[ia][ib][ic] =
              c * _HBC[ib][ic] *
              exp( EvtComplex( 0.0, phi * 0.5 *
                                        ( _lambdaA2[ia] - _lambdaB2[ib] + _lambdaC2[ic] ) ) ) *
              dfun;
        }
        prob1 += real( _amp[ia][ib][ic] * conj( _amp[ia][ib][ic] ) );
      }
    }
  }
 
  setUpRotationMatrices( p, theta, phi );
 
  applyRotationMatrices();
 
  double prob2 = 0.0;
 
  for ( ia = 0; ia < _nA; ia++ )
  {
    for ( ib = 0; ib < _nB; ib++ )
    {
      for ( ic = 0; ic < _nC; ic++ )
      {
        prob2 += real( _amp[ia][ib][ic] * conj( _amp[ia][ib][ic] ) );
        if ( _nA == 1 )
        {
          if ( _nB == 1 )
          {
            if ( _nC == 1 ) { amp.vertex( _amp[ia][ib][ic] ); }
            else { amp.vertex( ic, _amp[ia][ib][ic] ); }
          }
          else
          {
            if ( _nC == 1 ) { amp.vertex( ib, _amp[ia][ib][ic] ); }
            else
            {
              amp.vertex( ib, ic, _amp[ia][ib][ic] );
              // report(INFO,"EvtGen") << "ib,ic:"<<ib<<" "<<ic<<" "<<_amp[ia][ib][ic]<<endl;
            }
          }
        }
        else
        {
          if ( _nB == 1 )
          {
            if ( _nC == 1 ) { amp.vertex( ia, _amp[ia][ib][ic] ); }
            else { amp.vertex( ia, ic, _amp[ia][ib][ic] ); }
          }
          else
          {
            if ( _nC == 1 ) { amp.vertex( ia, ib, _amp[ia][ib][ic] ); }
            else { amp.vertex( ia, ib, ic, _amp[ia][ib][ic] ); }
          }
        }
      }
    }
  }
 
  if ( fabs( prob1 - prob2 ) > 0.000001 * prob1 )
  {
    report( INFO, "EvtGen" ) << "prob1,prob2:" << prob1 << " " << prob2 << endl;
    ::abort();
  }
 
  return;
}
 
void EvtEvalHelAmp::fillHelicity( int* lambda2, int n, int J2 ) {
 
  int i;
 
  // photon is special case!
  if ( n == 2 && J2 == 2 )
  {
    lambda2[0] = 2;
    lambda2[1] = -2;
    return;
  }
 
  assert( n == J2 + 1 );
 
  for ( i = 0; i < n; i++ ) { lambda2[i] = n - i * 2 - 1; }
 
  return;
}
 
void EvtEvalHelAmp::setUpRotationMatrices( EvtParticle* p, double theta, double phi ) {
 
  switch ( _JA2 )
  {
 
  case 0:
  case 1:
  case 2:
  case 3:
  case 4:
  case 5:
  case 6:
  case 7:
  case 8:
 
  {
 
    EvtSpinDensity R = p->rotateToHelicityBasis();
 
    int i, j, n;
 
    n = R.GetDim();
 
    assert( n == _nA );
 
    for ( i = 0; i < n; i++ )
    {
      for ( j = 0; j < n; j++ ) { _RA[i][j] = R.Get( i, j ); }
    }
  }
 
  break;
 
  default:
    report( ERROR, "EvtGen" ) << "Spin2(_JA2)=" << _JA2 << " not supported!" << endl;
    ::abort();
  }
 
  switch ( _JB2 )
  {
 
  case 0:
  case 1:
  case 2:
  case 3:
  case 4:
  case 5:
  case 6:
  case 7:
  case 8:
 
  {
 
    int i, j, n;
 
    EvtSpinDensity R = p->getDaug( 0 )->rotateToHelicityBasis( phi, theta, -phi );
 
    n = R.GetDim();
 
    assert( n == _nB );
 
    // report(INFO,"EvtGen") << "RB:"<< R <<endl;
 
    for ( i = 0; i < n; i++ )
    {
      for ( j = 0; j < n; j++ ) { _RB[i][j] = conj( R.Get( i, j ) ); }
    }
  }
 
  break;
 
  default:
    report( ERROR, "EvtGen" ) << "Spin2(_JB2)=" << _JB2 << " not supported!" << endl;
    ::abort();
  }
 
  switch ( _JC2 )
  {
 
  case 0:
  case 1:
  case 2:
  case 3:
  case 4:
  case 5:
  case 6:
  case 7:
  case 8:
 
  {
 
    int i, j, n;
 
    EvtSpinDensity R = p->getDaug( 1 )->rotateToHelicityBasis( phi, EvtConst::pi + theta,
                                                               phi - EvtConst::pi );
 
    n = R.GetDim();
 
    // report(INFO,"EvtGen") << "RC:"<< R <<endl;
 
    assert( n == _nC );
 
    for ( i = 0; i < n; i++ )
    {
      for ( j = 0; j < n; j++ ) { _RC[i][j] = conj( R.Get( i, j ) ); }
    }
  }
 
  break;
 
  default:
    report( ERROR, "EvtGen" ) << "Spin2(_JC2)=" << _JC2 << " not supported!" << endl;
    ::abort();
  }
}
 
void EvtEvalHelAmp::applyRotationMatrices() {
 
  int ia, ib, ic, i;
 
  EvtComplex temp;
 
  for ( ia = 0; ia < _nA; ia++ )
  {
    for ( ib = 0; ib < _nB; ib++ )
    {
      for ( ic = 0; ic < _nC; ic++ )
      {
        temp = 0;
        for ( i = 0; i < _nC; i++ ) { temp += _RC[i][ic] * _amp[ia][ib][i]; }
        _amp1[ia][ib][ic] = temp;
      }
    }
  }
 
  for ( ia = 0; ia < _nA; ia++ )
  {
    for ( ic = 0; ic < _nC; ic++ )
    {
      for ( ib = 0; ib < _nB; ib++ )
      {
        temp = 0;
        for ( i = 0; i < _nB; i++ ) { temp += _RB[i][ib] * _amp1[ia][i][ic]; }
        _amp3[ia][ib][ic] = temp;
      }
    }
  }
 
  for ( ib = 0; ib < _nB; ib++ )
  {
    for ( ic = 0; ic < _nC; ic++ )
    {
      for ( ia = 0; ia < _nA; ia++ )
      {
        temp = 0;
        for ( i = 0; i < _nA; i++ ) { temp += _RA[i][ia] * _amp3[i][ib][ic]; }
        _amp[ia][ib][ic] = temp;
      }
    }
  }
}