---

annih.cxx

Go to the documentation of this file.

// -*- C++ -*-  $Id: annih.cxx,v 1.1 2000/09/01 20:20:17 burnett Exp $
//
// This file is part of Gismo 2
//
// file: annih.cpp
// created 1-Dec-1991  T. Burnett

// Contents -------------------------------------------------------
//
//  EGS:annih
//
// End -----------------------------------------------------------


#include "EGS.h"
#include "Random.h"
#include <cmath>
#include <algorithm>


void EGS::annih()
{

    //                                version 4.00  --  24 jul 1986/1400
    //******************************************************************
    //   gamma spectrum for two gamma in-flight positron annihilation.
    //   using scheme based on heitler's p269-27o formulae
    //   this routine should give the correct distribution, but more
    //   thought could be put into devising a faster scheme.  however,
    //   since positron annihilation in flight is relatively infrequent
    //   this may not be worthwhile.
    //******************************************************************

    // Modified to handle annihilation at rest (THB 4-jan-1991)
    double avip = this->e + (double)rm; //available energy of incident positron
    double a  = avip/rm,
          ai = 1.0/a,
          g  = a-1.0,
          t  = std::max(0.,g-1.0),
          p  = sqrt(a*t),
          ep0= 1.0/(a+p),
          ep;

    //   sample 1/ep from ep=ep0 to 1.0-ep0
    if( ::fabs(ep0-0.5) > 1e-5)
     for(;;)
     {   float rnno01= Random::flat();
        ep=ep0*exp(rnno01*log((1.0-ep0)/ep0));
        //   now decide whether to accept
        float rnno02= Random::flat();
        float rejf=1.0-ep+ai*ai*(2.0*g-1.0/ep);
        if (rnno02 <=rejf) break ;
     }
     else
       ep = ep0;
    //   this completes sampling of a distribution which is asymmetric
    //   about ep=1/2, but which when symmetrized is the symmetric
    //   annihilation distribution. pick ep in (1/2,1-ep0).

    if (ep>0.5) ep=1.0-ep;

    //   set up energies
    double pesg1 = avip*ep; //energy of secondary gamma 1
    double pesg2 = avip - pesg1;

    EGS& g1 = *this;        // alias for input particle
    EGS& g2 = *new EGS(this);

    g1.e = pesg1; g1.iq = 0;
    g2.e = pesg2; g2.iq = 0;

    //   set up gamma lab angles
    float phi = Random::flat(2.*M_PI);   // random phi to use for both
    float costh1, costh2;

    if (t<=0.)  // at rest
    {  costh1 = Random::flat(-1.,1.);
       costh2 = -costh1;
    }
    else
    {
       costh1 = std::min(1.0,(pesg1-rm)*p/t/pesg1);
       costh2 = std::min(1.0,(pesg2-rm)*p/t/pesg2);
    }
    g1.rotate(costh1, phi);

    g2.rotate(costh2, phi+M_PI);

} //end of subroutine annih

---