---

PEGSData.cxx

Go to the documentation of this file.

// $Header: /nfs/slac/g/glast/ground/cvs/gismo/src/egs/PEGSData.cxx,v 1.1 2000/09/01 20:20:15 burnett Exp $
//
// This file is part of Gismo 2


#include "PEGSData.h"

#include "facilities/error.h"
#include <string.h>
#include <stdlib.h>
#include <iomanip>
#include <strstream>

void
PEGSData::read(std::istream& fin, Material& mat)
{
   char linebuf[101];
   iraylm = irayl;      // set local Rayleigh flag from current static

   // initialize data for bremsstrahlung production
   // the quantity zbrang is ( (1/111)*zeff**(1/3) )**2
   // where zeff is defined in equation (7) of pirs0203

    zbrang = 0.0;
    double pznorm=0.0;

    for( Material::Element_list::iterator pel= mat.elementList.begin();
            pel != mat.elementList.end(); ++pel) {
        const Material::AtomicElement* tel = *pel;
        zbrang += tel->pz*tel->z*(tel->z+1.0);
        pznorm += tel->pz;
    }

   zbrang = 8.116224E-05 * pow(zbrang/pznorm , 1./3.);

   // read media and thresh
//   float rlc;
//   fin >> rlc ;
//   mat.setRadiationLength(rlc);

   fin>> ae >> ap >> ue >> up;
#if 0
   cout << "Ae="<<ae<<", Ap="<<ap<<", ue="<<ue<<", up="<<up<<'\n';
#endif
   te = ae-rm;
   thmoll = 2.*te+rm;

   // array sizes from PEGS
   int iray1=0, mseke, meke, mleke, mcmfp, mrange;
   int msge, mge;
   fin >> msge >> mge
       >> mseke>> meke >> mleke >> mcmfp >> mrange  >> iray1;

   if( iray1==0 && iraylm ) {
       iraylm = 0;
#if 0  // this message is generated in the original
       std::ostrstream message;
       message << "PEGSData::read -- no Rayleigh data for material \'"
           << mat.name() << "\'" << '\0';
       WARNING(message.str());
#endif
   }

   fin.getline(linebuf,100);   // skip rest of line

   bp.readPEGS(fin);         // read brempr data

   el.readPEGS(fin, meke);   //read electin data

   ph.readPEGS(fin, mge, iray1);  // read photin data

   // change to given distance units

   rldu   = mat.radiationLength()/dunit;
   el.rescale(rldu);        //converts rl to dunits
   ph.rescale(rldu);

   el.fixtmx(estepe);         // apply low E correction to tmxs
}

// --------------- member function PEGSData::meanFreePath ----------

float PEGSData::meanFreePath(float e, int iq) const
{
    switch (iq)    {
        case -1: return (e>ae)? 1.0 / esig(e-rm) : 0;
        case  1: return (e>ae)? 1.0 / psig(e-rm) : 0;
        default:  {
            if( e<=ap ) return 0;
            float mfp = gmfp(e);
            if( iraylm ) mfp *= cohe(e);
            return mfp;
        }
    }
}

// --------------- member function PEGSData::dedx -----------------------

float PEGSData::dedx(float e, int iq) const
{   switch (iq)

    { case -1: return ededx(e-rm);
      case  1: return pdedx(e-rm);
      default: return 0.;
    }
}

// --------------- member function PEGSData::chooseChannel -----------------------

PEGSData::channel
PEGSData::chooseChannel(float rnno, float e, int iq) const
{

    switch (iq)
    {
      case -1:    //electron
        if (e<=ae)                 return STOP;             //stop
        if     (rnno <=ebr1(e-rm)) return BREMS;
        else                       return MOLLER;  //e- -> e- + e-

      case 1:    // positron
        if     (e<=ae || ae==0)    return ANNIH;   //e+ => gamma+gamma at rest
        else if(rnno <=pbr1(e-rm)) return BREMS;   //e+ -> e+ + gamma
        else if(rnno < pbr2(e-rm)) return BHABHA;  //e+ -> e+ + e-
        else                       return ANNIH;   //e+ -> gamma+gamma

      default:    // photon


        if( iraylm ) {
            // rayleigh-scattering option: cohe(e) is the probability to not scatter
            float cohfac = cohe(e);
            if( rnno <= (1.0-cohfac) )  return RAYLEIGH;  // gamma -> gamma
            // didn't scatter: rescale rnno to the range (0,1]
            rnno = (rnno - (1-cohfac))/cohfac;
        }
        if     (rnno <= gbr1(e))   return PAIR;   //gamma-> e+ + e-
        else if(rnno <  gbr2(e))   return COMPT;  //gamma-> gamma + e-
        else                       return PHOTO;  //gamma-> e- (+gamma)
    }
}
//--------------------------- printing stuff ------------------------------
void
PEGSData::printHead1(std::ostream& out)
{
    out << std::setw(58)
        << "---------------------------PEGS4 data ----------------------";
}

void
PEGSData::printHead2(std::ostream& out)
{
    out
        << std::setw(10) << "rl(cm)"
        << std::setw(10) << "ecut(MeV)"
        << std::setw(10) << "pcut(Mev)"
        << std::setw(10) << "emax(Mev)"
        << std::setw(10) << "pmax(Mev)"
        << std::setw(8)  << "iraylm"
        ;
}
void
PEGSData::printOn(std::ostream& out)
{
    out
        << std::setw(10) << rl()
        << std::setw(10) << ecut()
        << std::setw(10) << pcut()
        << std::setw(10) << emax()
        << std::setw(10) << pmax()
        << std::setw(8)  << iraylm
        ;
}

PEGSData::~PEGSData()
{
}

//  -------------------- allocate static constants ----------------------

const float PEGSData::dunit = 1.0;  // distance unit in cm
const float PEGSData::rm = 0.511;   // electron Mass in MeV
float PEGSData::estepe=0.50;  // fractional eloss allowed
int PEGSData::irayl=1;        // Rayleigh flag

---