---

Elecin.cxx

Go to the documentation of this file.

// -*- C++ -*-  $Id: Elecin.cxx,v 1.1 2000/09/01 20:20:14 burnett Exp $
//
// This file is part of Gismo 2
//
#include "Elecin.h"

#include <cmath>
#include <stdlib.h>
#include <stdio.h>

void Elecin::readPEGS(std::istream& fin, int m)
{
   meke = m;
   fin >> xr0  >> teff0 >> blcc >> xcc;
   fin >> eke0 >> eke1;
   esig0  = new float[meke];
   esig1  = new float[meke];
   psig0  = new float[meke];
   psig1  = new float[meke];
   ededx0 = new float[meke];
   ededx1 = new float[meke];
   pdedx0 = new float[meke];
   pdedx1 = new float[meke];
   ebr10  = new float[meke];
   ebr11  = new float[meke];
   pbr10  = new float[meke];
   pbr11  = new float[meke];
   pbr20  = new float[meke];
   pbr21  = new float[meke];
   tmxs0  = new float[meke];
   tmxs1  = new float[meke];

   for (int i=0; i< meke; i++)
      fin >> esig0[i]  >> esig1[i]
          >> psig0[i]  >> psig1[i]
          >> ededx0[i] >> ededx1[i]
          >> pdedx0[i] >> pdedx1[i]
          >> ebr10[i]  >> ebr11[i]
          >> pbr10[i]  >> pbr11[i]
          >> pbr20[i]  >> pbr21[i]
          >> tmxs0[i]  >> tmxs1[i];
}

void Elecin::rescale(float dfact)
{

   // Modify distance-scale-dependent quantities


    float dfacti=1.0/dfact;  //convert rl**-1 to dunits**-1//

    for (int i=0; i< meke; i++)
    {  //$scale $lgn(esig,psig,ededx,pdedx(i,im)/0,1/) by dfacti;
       //$scale $lgn(tmxs(i,im)/0,1/) by dfact;

       esig0[i]  *= dfacti;
       esig1[i]  *= dfacti;
       psig0[i]  *= dfacti;
       psig1[i]  *= dfacti;
       ededx0[i] *= dfacti;
       ededx1[i] *= dfacti;
       pdedx0[i] *= dfacti;
       pdedx1[i] *= dfacti;
       tmxs0[i]  *= dfact;
       tmxs1[i]  *= dfact;
    }
    teff0 *= dfact;          //$scale teff0(im) by dfact;

    blcc  *= dfacti;         //$scale blcc(im) by dfacti
    xcc   *= sqrt(dfacti);   //$scale xcc(im) by sqrt(dfacti)
}

void Elecin::fixtmx(float estepe){
/**********************************************************************
"  this routine changes the step size algorithm used in egs so that   "
"  the step size arrays for tmxs correspond to an arbitrary,but       "
"  fixed fractional energy loss estepe.                               "
"  it is only necessary for low energy electron problems since        "
"  typically the 200*teff0 restriction on tmxs is more stringent      "
"   for electrons with energies above a few mev                       "
"                                                                     "
"  note that the $tmxs-over-ride macro is still in force in egs.      "
"                                                                     "
"  the routine changes the values only for the medium 'medium'        "
"  and it should probably be used for all media in a problem.         "
"                                                                     "
"  the routine must be called after hatch has been called and before  "
"  the simulation is begun.                                           "
"                                                                     "
"  the routine is independent of what units are being used, as long   "
"  as they are consistent( e.g. cm, rl or g/cm**2 )                   "
"                                                                     "
"  if called with estepe=0, the default algorithm is used             "
"                                                                     "
"  for a detailed discussion of the use of this routine, see          "
"  'low energy electron transport with egs' in nuclear instr. and     "
"  methods a227 (1984)535-548.  d.w.o. rogers                         "
"                                                                     "
"  v01        dec 10,1981 dave rogers nrcc                            "
"  v02        dec 1984  egs4 version                                  "
"  Adapted for Gismo's EGS  jan 1993  T. Burnett
*********************************************************************/
if(estepe == 0) return; //i.e. use the current algorithm

// set up some variables for first pass through loop"
float ei = exp( (1.-eke0)/eke1);  //energy of first table entry"
updateIndex(ei);
float eil = loge;
// this is equivalent to $setinterval eil,eke; but avoids roundoff"


// $evaluate ededx using ededx(eil);"get the electron stoppping at ei"
//"now calculate step required to cause an estepe reduction in energy"
float si=estepe*ei/ededx();
//"tabulated energies are in a fixed ratio - calc log of the ratio"
float eratio=-1./eke1;

for(int i=1; i<=meke-1; i++)
{
   float eip1=exp(((i+1)-eke0)/eke1);  //"energy at i+1"
   updateIndex(eip1);
   float eip1l=loge; index=i-1; ;//"designed this way=$setinterval"
   //$evaluate ededx using ededx(eip1l);
   float sip1=estepe*eip1/ededx();

   //"now solve these equations                                     "
   //"  si   = tmxs1 * eil   + tmxs0                                "
   //"  sip1 = tmxs1 * eip1l + tmxs0                                "

   tmxs1[i-1] = (si-sip1)/eratio;
   tmxs0[i-1] = si-tmxs1[i-1]*eil;
   // "transfer values for next loop"
   eil=eip1l; si=sip1;
}
// "now pick up last table entry which applies only to last energy"
tmxs0[meke-1] = tmxs0[meke-2];
tmxs1[meke-1] = tmxs1[meke-2];
return;
}

#define inline

inline void
Elecin::updateIndex(float enew) const
{
   if (enew == eold) return;
   Elecin * self = (Elecin *)this;
   self->loge = log(enew);
   self->index= (int)(eke1*loge + eke0 -1); // note -1 for C arrays
   if (index>=meke)
   {
      char strl[70];
      sprintf(strl, "Elecin: energy %9.3f too large for table: index %u , max %u", 
                     enew, self->index, meke);

     WARNING(strl);
     self->index = meke-1;
   }
   self->eold = enew;
}
inline float Elecin::esig()const{return  esig1[index]*loge + esig0[index];}
inline float Elecin::psig()const{return  psig1[index]*loge + psig0[index];}
inline float Elecin::ededx()const{return ededx1[index]*loge +ededx0[index];}
inline float Elecin::pdedx()const{return pdedx1[index]*loge +pdedx0[index];}
inline float Elecin::ebr1()const{return  ebr11[index]*loge + ebr10[index];}
inline float Elecin::pbr1()const{return  pbr11[index]*loge + pbr10[index];}
inline float Elecin::pbr2()const{return  pbr21[index]*loge + pbr20[index];}
inline float Elecin::tmxs()const{return  tmxs1[index]*loge + tmxs0[index];}

Elecin::~Elecin()
{
    if( meke==0 ) return;
    delete [] esig0 ;
    delete [] esig1 ;
    delete [] psig0 ;
    delete [] psig1 ;
    delete [] ededx0 ;
    delete [] ededx1;
    delete [] pdedx0;
    delete [] pdedx1;
    delete [] ebr10 ;
    delete [] ebr11 ;
    delete [] pbr10 ;
    delete [] pbr11 ;
    delete [] pbr20 ;
    delete [] pbr21 ;
    delete [] tmxs0 ;
    delete [] tmxs1 ;

}

---