EvtDalitzPlot.cpp

Go to the documentation of this file.

/***********************************************************************
* 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/EvtDalitzPlot.hh"

#include "EvtGenBase/EvtDecayMode.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtTwoBodyVertex.hh"

#include <assert.h>
#include <math.h>
#include <stdio.h>

using namespace EvtCyclic3;

EvtDalitzPlot::EvtDalitzPlot() :
    _mA( 0. ), _mB( 0. ), _mC( 0. ), _bigM( 0. ), _ldel( 0. ), _rdel( 0. )
{
}

EvtDalitzPlot::EvtDalitzPlot( double mA, double mB, double mC, double bigM,
                              double ldel, double rdel ) :
    _mA( mA ), _mB( mB ), _mC( mC ), _bigM( bigM ), _ldel( ldel ), _rdel( rdel )
{
    sanityCheck();
}

EvtDalitzPlot::EvtDalitzPlot( const EvtDecayMode& mode, double ldel, double rdel )
{
    _mA = EvtPDL::getMeanMass( EvtPDL::getId( mode.dau( A ) ) );
    _mB = EvtPDL::getMeanMass( EvtPDL::getId( mode.dau( B ) ) );
    _mC = EvtPDL::getMeanMass( EvtPDL::getId( mode.dau( C ) ) );
    _bigM = EvtPDL::getMeanMass( EvtPDL::getId( mode.mother() ) );

    _ldel = ldel;
    _rdel = rdel;

    sanityCheck();
}

bool EvtDalitzPlot::operator==( const EvtDalitzPlot& other ) const
{
    bool ret = false;
    if ( _mA == other._mA && _mB == other._mB && _mC == other._mC &&
         _bigM == other._bigM )
        ret = true;

    return ret;
}

const EvtDalitzPlot* EvtDalitzPlot::clone() const
{
    return new EvtDalitzPlot( *this );
}

void EvtDalitzPlot::sanityCheck() const
{
    if ( _mA < 0 || _mB < 0 || _mC < 0 || _bigM <= 0 ||
         _bigM - _mA - _mB - _mC < 0. ) {
        printf( "Invalid Dalitz plot %f %f %f %f\n", _mA, _mB, _mC, _bigM );
        assert( 0 );
    }
    assert( _ldel <= 0. );
    assert( _rdel >= 0. );
}

double EvtDalitzPlot::m( Index i ) const
{
    double m = _mA;
    if ( i == B )
        m = _mB;
    else if ( i == C )
        m = _mC;

    return m;
}

double EvtDalitzPlot::sum() const
{
    return _mA * _mA + _mB * _mB + _mC * _mC + _bigM * _bigM;
}

double EvtDalitzPlot::qAbsMin( Pair i ) const
{
    Index j = first( i );
    Index k = second( i );

    return ( m( j ) + m( k ) ) * ( m( j ) + m( k ) );
}

double EvtDalitzPlot::qAbsMax( Pair i ) const
{
    Index j = other( i );
    return ( _bigM - m( j ) ) * ( _bigM - m( j ) );
}

double EvtDalitzPlot::qResAbsMin( EvtCyclic3::Pair i ) const
{
    return qAbsMin( i ) - sum() / 3.;
}

double EvtDalitzPlot::qResAbsMax( EvtCyclic3::Pair i ) const
{
    return qAbsMax( i ) - sum() / 3.;
}

double EvtDalitzPlot::qHelAbsMin( EvtCyclic3::Pair i ) const
{
    Pair j = next( i );
    Pair k = prev( i );
    return ( qAbsMin( j ) - qAbsMax( k ) ) / 2.;
}

double EvtDalitzPlot::qHelAbsMax( EvtCyclic3::Pair i ) const
{
    Pair j = next( i );
    Pair k = prev( i );
    return ( qAbsMax( j ) - qAbsMin( k ) ) / 2.;
}

double EvtDalitzPlot::mAbsMin( Pair i ) const
{
    return sqrt( qAbsMin( i ) );
}

double EvtDalitzPlot::mAbsMax( Pair i ) const
{
    return sqrt( qAbsMax( i ) );
}

// parallel

double EvtDalitzPlot::qMin( Pair i, Pair j, double q ) const
{
    if ( i == j )
        return q;

    else {
        // Particle pair j defines the rest-frame
        // 0 - particle common to r.f. and angle calculations
        // 1 - particle belonging to r.f. but not angle
        // 2 - particle not belonging to r.f.

        Index k0 = common( i, j );
        Index k2 = other( j );
        Index k1 = other( k0, k2 );

        // Energy, momentum of particle common to rest-frame and angle
        EvtTwoBodyKine jpair( m( k0 ), m( k1 ), sqrt( q ) );
        double pk = jpair.p();
        double ek = jpair.e( EvtTwoBodyKine::A, EvtTwoBodyKine::AB );

        // Energy and momentum of the other particle
        EvtTwoBodyKine mother( sqrt( q ), m( k2 ), bigM() );
        double ej = mother.e( EvtTwoBodyKine::B, EvtTwoBodyKine::A );
        double pj = mother.p( EvtTwoBodyKine::A );

        // See PDG 34.4.3.1
        return ( ek + ej ) * ( ek + ej ) - ( pk + pj ) * ( pk + pj );
    }
}

// antiparallel

double EvtDalitzPlot::qMax( Pair i, Pair j, double q ) const
{
    if ( i == j )
        return q;
    else {
        // Particle pair j defines the rest-frame
        // 0 - particle common to r.f. and angle calculations
        // 1 - particle belonging to r.f. but not angle
        // 2 - particle not belonging to r.f.

        Index k0 = common( i, j );
        Index k2 = other( j );
        Index k1 = other( k0, k2 );

        // Energy, momentum of particle common to rest-frame and angle
        EvtTwoBodyKine jpair( m( k0 ), m( k1 ), sqrt( q ) );
        double ek = jpair.e( EvtTwoBodyKine::A, EvtTwoBodyKine::AB );
        double pk = jpair.p();

        // Energy and momentum of the other particle
        EvtTwoBodyKine mother( sqrt( q ), m( k2 ), bigM() );
        double ej = mother.e( EvtTwoBodyKine::B, EvtTwoBodyKine::A );
        double pj = mother.p( EvtTwoBodyKine::A );

        // See PDG 34.4.3.1
        return ( ek + ej ) * ( ek + ej ) - ( pk - pj ) * ( pk - pj );
    }
}

double EvtDalitzPlot::getArea( int N, Pair i, Pair j ) const
{
    // Trapezoidal integral over qi. qj can be calculated.
    // The first and the last point are zero, so they are not counted

    double dh = ( qAbsMax( i ) - qAbsMin( i ) ) / ( (double)N );
    double sum = 0;

    int ii;
    for ( ii = 1; ii < N; ii++ ) {
        double x = qAbsMin( i ) + ii * dh;
        double dy = qMax( j, i, x ) - qMin( j, i, x );
        sum += dy;
    }

    return sum * dh;
}

double EvtDalitzPlot::cosTh( EvtCyclic3::Pair i1, double q1,
                             EvtCyclic3::Pair i2, double q2 ) const
{
    if ( i1 == i2 )
        return 1.;

    double qmax = qMax( i1, i2, q2 );
    double qmin = qMin( i1, i2, q2 );

    double cos = ( qmax + qmin - 2 * q1 ) / ( qmax - qmin );

    return cos;
}

double EvtDalitzPlot::e( Index i, Pair j, double q ) const
{
    if ( i == other( j ) ) {
        // i does not belong to pair j

        return ( bigM() * bigM() - q - m( i ) * m( i ) ) / 2 / sqrt( q );
    } else {
        // i and k make pair j

        Index k;
        if ( first( j ) == i )
            k = second( j );
        else
            k = first( j );

        double e = ( q + m( i ) * m( i ) - m( k ) * m( k ) ) / 2 / sqrt( q );
        return e;
    }
}

double EvtDalitzPlot::p( Index i, Pair j, double q ) const
{
    double en = e( i, j, q );
    double p2 = en * en - m( i ) * m( i );

    if ( p2 < 0 ) {
        printf( "Bad value of p2 %f %d %d %f %f\n", p2, i, j, en, m( i ) );
        assert( 0 );
    }

    return sqrt( p2 );
}

double EvtDalitzPlot::q( EvtCyclic3::Pair i1, double cosTh, EvtCyclic3::Pair i2,
                         double q2 ) const
{
    if ( i1 == i2 )
        return q2;

    EvtCyclic3::Index f = first( i1 );
    EvtCyclic3::Index s = second( i1 );
    return m( f ) * m( f ) + m( s ) * m( s ) +
           2 * e( f, i2, q2 ) * e( s, i2, q2 ) -
           2 * p( f, i2, q2 ) * p( s, i2, q2 ) * cosTh;
}

double EvtDalitzPlot::jacobian( EvtCyclic3::Pair i, double q ) const
{
    return 2 * p( first( i ), i, q ) *
           p( other( i ), i, q );    // J(BC) = 2pA*pB = 2pA*pC
}

EvtTwoBodyVertex EvtDalitzPlot::vD( Pair iRes, double m0, int L ) const
{
    return EvtTwoBodyVertex( m( first( iRes ) ), m( second( iRes ) ), m0, L );
}

EvtTwoBodyVertex EvtDalitzPlot::vB( Pair iRes, double m0, int L ) const
{
    return EvtTwoBodyVertex( m0, m( other( iRes ) ), bigM(), L );
}

void EvtDalitzPlot::print() const
{
    // masses
    printf( "Mass  M    %f\n", bigM() );
    printf( "Mass mA    %f\n", _mA );
    printf( "Mass mB    %f\n", _mB );
    printf( "Mass mC    %f\n", _mC );
    // limits
    printf( "Limits qAB %f : %f\n", qAbsMin( AB ), qAbsMax( AB ) );
    printf( "Limits qBC %f : %f\n", qAbsMin( BC ), qAbsMax( BC ) );
    printf( "Limits qCA %f : %f\n", qAbsMin( CA ), qAbsMax( CA ) );
    printf( "Sum q       %f\n", sum() );
    printf( "Limit qsum  %f : %f\n", qSumMin(), qSumMax() );
}