Page MenuHomeHEPForge
Files F19251953

EvtEvalHelAmp.cpp
No OneTemporary
Actions

Size
12 KB
Referenced Files
None
Subscribers
None

EvtEvalHelAmp.cpp
View Options

/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
* *
* This file is part of EvtGen. *
* *
* EvtGen is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* EvtGen is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/
#include "EvtGenBase/EvtEvalHelAmp.hh"
#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtdFunction.hh"
#include <stdlib.h>
using std::endl;
EvtEvalHelAmp::~EvtEvalHelAmp()
{
//deallocate 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( EvtId idA, EvtId idB, EvtId idC,
EvtComplexPtrPtr HBC )
{
EvtSpinType::spintype typeA = EvtPDL::getSpinType( idA );
EvtSpinType::spintype typeB = EvtPDL::getSpinType( idB );
EvtSpinType::spintype typeC = EvtPDL::getSpinType( idC );
//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, idA );
fillHelicity( _lambdaB2, _nB, _JB2, idB );
fillHelicity( _lambdaC2, _nC, _JC2, idC );
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();
double theta = acos( pB.get( 3 ) / pB.d3mag() );
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] );
}
}
} 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 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "prob1,prob2:" << prob1 << " " << prob2 << endl;
::abort();
}
return;
}
void EvtEvalHelAmp::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;
}
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:
EvtGenReport( EVTGEN_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 );
for ( i = 0; i < n; i++ ) {
for ( j = 0; j < n; j++ ) {
_RB[i][j] = conj( R.get( i, j ) );
}
}
}
break;
default:
EvtGenReport( EVTGEN_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();
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:
EvtGenReport( EVTGEN_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;
}
}
}
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 6:12 AM (15 h, 23 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6550544
Default Alt Text
EvtEvalHelAmp.cpp (12 KB)

Event Timeline

HEPForge.org · .