EvtNeutrinoParticle.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/EvtNeutrinoParticle.hh"
 
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4R.hh"
 
#include <iostream>
#include <math.h>
#include <stdlib.h>
using std::endl;
 
void EvtNeutrinoParticle::init( EvtId part_n, const EvtVector4R& p4 )
{
    m_validP4 = true;
    setp( p4 );
    setpart_num( part_n );
 
    double e, px, py, pz;
    e = p4.get( 0 );
    px = p4.get( 1 );
    py = p4.get( 2 );
    pz = p4.get( 3 );
 
    if ( EvtPDL::getStdHep( part_n ) == 0 ) {
        EvtGenReport( EVTGEN_ERROR, "EvtGen" )
            << "Error in EvtNeutrinoParticle::init, part_n=" << part_n.getId()
            << endl;
    }
 
    if ( EvtPDL::getStdHep( part_n ) > 0 ) {
        double beta, alpha, p2, norm;
 
        // See Sakurai p. 167-169
        // and Renton p. 126
 
        p2 = px * px + py * py + pz * pz;
 
        beta = acos( pz / sqrt( p2 ) );
        alpha = atan2( py, px );
 
        norm = sqrt( 2 * e );
 
        double cosb, sinb, cosa, sina;
 
        cosb = cos( 0.5 * beta );
        sinb = sin( 0.5 * beta );
 
        cosa = cos( 0.5 * alpha );
        sina = sin( 0.5 * alpha );
 
        m_spinor_parent.set( -norm * sinb * EvtComplex( cosa, -sina ),
                             norm * cosb * EvtComplex( cosa, sina ),
                             norm * sinb * EvtComplex( cosa, -sina ),
                             -norm * cosb * EvtComplex( cosa, sina ) );
 
    } else {
        px = -p4.get( 1 );
        py = -p4.get( 2 );
        pz = -p4.get( 3 );
 
        double pn, sqrpn;
 
        pn = e;
        sqrpn = sqrt( pn - pz );
 
        m_spinor_parent.set( ( 1.0 / sqrpn ) * EvtComplex( px, -py ),
                             EvtComplex( sqrpn, 0.0 ),
                             ( -1.0 / sqrpn ) * EvtComplex( px, -py ),
                             -EvtComplex( sqrpn, 0.0 ) );
    }
 
    setLifetime();
}
 
EvtDiracSpinor EvtNeutrinoParticle::spParentNeutrino() const
{
    return m_spinor_parent;
}
 
EvtDiracSpinor EvtNeutrinoParticle::spNeutrino() const
{
    EvtGenReport( EVTGEN_ERROR, "EvtGen" )
        << "Tried to get neutrino spinor in restframe";
    EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Will terminate execution.";
 
    ::abort();
 
    return m_spinor_rest;
}
 
EvtSpinDensity EvtNeutrinoParticle::rotateToHelicityBasis() const
{
    EvtGenReport( EVTGEN_ERROR, "EvtGen" )
        << "rotateToHelicityBasis not implemented for neutrino.";
    EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Will terminate execution.";
 
    ::abort();
 
    EvtSpinDensity rho;
    return rho;
}
 
EvtSpinDensity EvtNeutrinoParticle::rotateToHelicityBasis( double, double,
                                                           double ) const
{
    EvtGenReport( EVTGEN_ERROR, "EvtGen" )
        << "rotateToHelicityBasis(alpha,beta,gama) not implemented for neutrino.";
    EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Will terminate execution.";
 
    ::abort();
 
    EvtSpinDensity R;
    R.setDiag( 1 );
 
    return R;
}