---

pair.cxx

Go to the documentation of this file.

// -*- C++ -*-  $Id: pair.cxx,v 1.1 2000/09/01 20:20:21 burnett Exp $
//
// This file is part of Gismo 2

#include "EGS.h"
#include "PEGSData.h"
#include "Random.h"

// used in angle sampling below
static double set_pair_rejection_function(double input,double TTEIG,double ZTARG,double ESEDEI);

// ----------------- member function pair ---------------------------//

void EGS::pairProd(const PEGSData & mat)
{
    float eig = this->e;   // energy of gamma
    float br;   // fraction to electron #2
    
    //******************************************************************
    //   for a photon energy less than 2.1 MeV, the approximation is
    //   made that one pair electron (or positron) has only its rest
    //   mass energy.   for photon energy between 2.1 mev and 50 mev the
    //   Bethe-Heitler cross section is employed.  above 50 mev the
    //   coulomb corrected Bethe-Heitler cross section is used.
    //   (Butcher and Messel, op. cit., p. 17-19, 22).
    //******************************************************************
    if (eig<=2.1)
    {
        //   below 2.1,use approximation
        br=rm/eig; //energy of secondary 'electron'
    }
    else
    {   //above 2.1, must sample
        
        //   decide whether to use bethe-heitler(lvx=0,lvl=0,2) or
        //   coulomb corrected(lvx=1,lvl=3,5) cross sections.
        //   see related comments in brems.
        int lvx = 1, lvl0=3;
        if (eig< 50.) {lvx=0;lvl0=0;}
        
        for (;;)
        {
            //retry if rejected because del out of range, or by screening
            float rnno30=Random::flat(); //we'll need at least one random number
            int lvl;
            
            //   now decide which of the two subdistributions to use.
            float rnno31=Random::flat();
            if (rnno31 >= mat.bpar(lvx) )
            {   //use the subdistribution that is
                //   proportional to 12*(br-0.5)**2.  it uses a(delta) for
                //   screening function.
                lvl=lvl0;
                //   from symmetry, only need to sample br in interval (0,.5)
                float rnno32=Random::flat();
                float rnno33=Random::flat();
                br=0.5*(1.0-std::max(rnno32,std::max(rnno33,rnno30) ) );
            }
            else
            {   //use the subdistribution that is proportional to 1,
                //   i.e., uniform.  it uses c(delta) for a screening
                //   rejection function.
                
                lvl=lvl0+2;
                br = rnno30*0.5;
            }
            
            //   the screening functions are functions of delta=delcm*del,
            //   where delcm= 136.0*exp(zg)*rm (same as for brems)
            //   and where del=1./(eg0*br*(1.0-br))
            //   with eg0 = incident photon energy and br=energy fraction.
            
            if (br==0.0) continue; //to avoid division by zero
            
            float del=1.0/(eig*br*(1.0-br));
            
            float rejf = mat.screening(del,lvx, lvl);
            if (rejf == 0.) continue;
            
            float rnscrn=Random::flat(); //random number for screening rejection
            if (rnscrn <= rejf) break ; //retry until accepted
        } // end of rejection loop
        
        
    } //end of eig.gt.2.1 else
    
    
    EGS& es1 = *this;   // use as an alias for incoming photon,
    // converted to a secondary
    EGS& es2 = *new EGS(this);   // create other secondary
    
    //   energy going to lower secondary has now been determined
    double pese2=br*eig;     //precise energy of secondary 'electron' 2
    double pese1=(double)eig-pese2; //precise energy of secondary 'electron' 1
    
    //negative kinetic energy fix
    //EGS4 DEFAULT code did not include a fix for this, it cropped up once
    //Dr. Bielajew's code was implemented, not suprising since
    //the EGS4 bug was not found until this distibution was
    //implemented (as a result of)--harrison
    
    //using Dr. Hartman's (NASA GSFC LHEA) fix
    //C    R.C. HARTMAN 94/11/19 - PROTECT AGAINST ENERGY BELOW REST MASS:            
    //C           GIVE LOWER ENERGY ELECTRON 1% OF AVAILABLE KINETIC ENERGY.          
    //      IF ((ESE2.LE.RM)) THEN                                                    
    //       ESE2 = RM + 0.01*(PEIG - 2*RM)                                          
    //      END IF 
    
    if ((pese2 <= rm2) || (pese1 < rm2)){
        pese2 = (rm2 + 0.01*((double)eig - 2*rm2));
        pese1 = (double) eig - pese2;
    }
    //end negative kinetic energy fix

    //assign energies
    es1.e = pese1;
    es2.e = pese2;
    
    //begin angular distribution sampling
    
    double theta;
    double costhe1 = 100000.0;  // set to crazy value to test
    double costhe2 = 100000.0; 
    double sinthe1 = 100000.0;
    double sinthe2 = 100000.0;
    //set IPRDST to choose which distribution to use
    //false for default EGS4, true for improved distributions
    //Should be done as a #define/#ifdef fix later....harrison
    static bool IPRDST = true;
    
    if (! IPRDST){
        //Default EGS4 routine
        
        //   this average angle of emission for both pair production and
        //   bremsstrahlung is much smaller than the average angle of
        //   multiple scattering for delta t transport=0.01 r.l.
        //   the initial and final directions are coplanar
        
        theta = (double)rm/eig;
        costhe1 = cos(theta);
        sinthe1 = sin(theta);
        sinthe2 = sinthe1;
        costhe2 = costhe1;
    }
    
    else if ((IPRDST) && (eig < 4.14)){ 
        //"lowest order approx" distribution from paper by 
        //Alex F. Bielajew:
        //"Improved angular sampling for pair production
        //in the EGS4 code system", 10/20/94, PIRS-0287R1 
        
        double E1,E2,p1,p2;
        
        E1 = (double) pese1; //Energy in MeV 
        E2 = (double) pese2;
        p1 = sqrt((E1 - rm2)*(E1+rm2));
        p2 = sqrt((E2 - rm2)*(E2+rm2));
        costhe1 = (double) Random::flat();
        costhe2 = (double) Random::flat();
        costhe1 = 1.0 - 2.0*costhe1;
        costhe2 = 1.0 - 2.0*costhe2;
        sinthe1 = rm2*sqrt((1.0-costhe1)*(1.0+costhe1))/(p1*costhe1+E1);
        sinthe2 = rm2*sqrt((1.0-costhe2)*(1.0+costhe2))/(p2*costhe2+E2);
        costhe1 = (E1 * costhe1+p1)/(p1*costhe1+E1);
        costhe2 = (E2 * costhe2+p2)/(p2*costhe2+E2);
        
    }
    
    else if ((IPRDST) && (eig >= 4.14)){
        //Schiff (Koch, Motz, & Olsen) Distribution from paper by Alex F. Bielajew:
        //"Improved angular sampling for pair production
        //in the EGS4 code system", 10/20/94, PIRS-0287R1 
        //function prototypes from other modifications:
        //from EGS.h/cxx: 
        //double set_pair_rejection_function(double input,double TTEIG,double ZTARG,double ESEDEI);
        //from PEGSData.h/cxx:
        //double hzbrang()const;
        
        int n;
        double ZTARG,TTEIG,ESEDEI,ESEDER,YA,test,XIMIN;
        double TTESE = 1000000.0;
        double GALPHA,GBETA,XITRY,XIMID,REJTST,REJMIN;
        double REJTOP = 1000000.0; // set to crazy value
        double RTEST,XITST,REJMID;
        
        ZTARG = (double) mat.hzbrang(); //((1/111)*zeff**(1/3))**2
        TTEIG = (double) eig/rm2; //energy of gamma in units of e- rest mass
        
        
        n = 0;
        while (n < 2){  
            if (n == 0) TTESE = (double)pese1/rm2; //energy of particle 1 in units of rest mass
            if (n == 1) TTESE = (double)pese2/rm2; //energy of particle 2 in units as above
            ESEDEI = TTESE/(TTEIG-TTESE); //r
            ESEDER = 1.0/ESEDEI; //1 over r
            XIMIN = 1.0/(1.0+pow((3.141593*TTESE),2));
            REJMIN = set_pair_rejection_function(XIMIN,TTEIG,ZTARG,ESEDEI);
            YA = pow((2.0/TTEIG),2);
            test = sqrt(YA/ZTARG);
            if (test >= 0.5) test = 0.5;
            if (test <= XIMIN) test = XIMIN;
            if (test <= 0.01) test = 0.01;
            XITRY = test;
            GALPHA = 1.0 + 0.25*log(YA+ZTARG*pow(XITRY,2));
            GBETA = 0.5*ZTARG*XITRY/(YA+ZTARG*pow(XITRY,2));
            GALPHA = GALPHA-GBETA*(XITRY-0.5);
            XIMID = GALPHA/(3.0*GBETA);
            if (GALPHA >= 0.0){
                XIMID = 0.5-XIMID+sqrt(pow(XIMID,2)+0.25);
            }
            if (GALPHA < 0.0){
                XIMID = 0.5-XIMID-sqrt(pow(XIMID,2)+0.25);
            }
            test = XIMID;
            if (test >= 0.5) test = 0.5;
            if (test <= XIMIN) test = XIMIN;
            if (test <= 0.01) test = 0.01;
            XIMID = test;
            REJMID = set_pair_rejection_function(XIMID,TTEIG,ZTARG,ESEDEI);
            if (REJMID <= REJMIN) REJTOP = 1.02*REJMIN;
            if (REJMIN < REJMID)  REJTOP = 1.02*REJMID;
            
            do {
                XITST = (double) Random::flat();
                RTEST = (double) Random::flat();
                REJTST = set_pair_rejection_function(XITST,TTEIG,ZTARG,ESEDEI);
                theta = sqrt((1.0/XITST)-1.0)/TTESE;
            } while ((RTEST > (REJTST/REJTOP)) || (theta >= 3.1415926535));
            if (n == 0){ 
                costhe1 = cos(theta); 
                sinthe1 = sin(theta);
                //if (sinthe1 < 0) sinthe1 = -1.0 * sinthe1; //just in case     
            }
            
            if (n == 1){ 
                costhe2 = cos(theta); 
                sinthe2 = sin(theta);
                //if (sinthe2 < 0) sinthe2 = sinthe2 * -1.0; //again, just in case
            }
            
            n++;
        }
    }
    
    
    
    //generate uniform phi for both secondaries
    double phi = Random::flat(2.*M_PI); 
    
    
    //set'm up, angles w/respect to gamma direction
    es1.rotate(costhe1, phi, sinthe1);
    es2.rotate(costhe2, phi+M_PI,sinthe2);
    
    //   now randomly decide which is the positron
    
    es2.iq = -( es1.iq = (Random::flat() <= 0.5)? 1 : -1 );
    
} // end of EGS::pair

static double set_pair_rejection_function(double input,double TTEIG,double ZTARG,double ESEDEI)
{ //harrison
    double term,ESEDER,result;
    ESEDER = 1.0/ESEDEI;
    term = pow(((1.0 + ESEDEI) * (1.0 + ESEDER)/(2.0*TTEIG)),2);
    result = 2.0 + 3.0*(ESEDEI+ESEDER) - 4.00*(ESEDEI+ESEDER + 1.0 - 4.0*pow((input - 0.5),2))*(1.0 + 0.25*log(term+ZTARG*pow(input,2)));
    return result;
}

---