---

GheishaData.cxx

Go to the documentation of this file.

// $Id: GheishaData.cxx,v 1.1 2000/09/01 20:20:23 burnett Exp $
// Author: T. Burnett
// Contents ----------------------------------------------------------------
//
//
// Description
//
//      GheishaData implementations
//
// End --------------------------------------------------------------------


#include "GheishaData.h"

#include "gismo/Material.h"

#include <cmath>
#include <cstring> // for gcc
#include <cstdlib> // for gcc

#define gsmate gsmate_
#define gsmixt gsmixt_

extern "C" {
   void gsmate(const int* imat, /* const char * namate ,*/ 
                const float* A,
              const float* Z, const float* dens, const float* radl, const float * absl);
//THB: following change required by unexplained change in GHEISHA fortran code
   
   void gsmixt(const int* imat,/* const char * namate, */
          const int* nlmat, const float * wmat,
              const int * what, const float* dens);          
}


void
GheishaData::read(std::istream&, Material& mat)
{

   float  wtNorm = 0.,
          numNorm = 0.;
   int numElem=0;

    Material::Element_list::iterator it=mat.elementList.begin();
    for(; it != mat.elementList.end(); ++ it) {
        Material::AtomicElement* el = *it;
        wtNorm += el->rhoz;
        numNorm += el->pz;
        numElem++;
    }
   float* weights = new float[numElem];
   int*   indices = new int[numElem];
   char hadMatName[20];

   int i=0;
   for(it = mat.elementList.begin(); it != mat.elementList.end(); ++it) {
        Material::AtomicElement* el = *it;
       strncpy(hadMatName, el->name, 18);
      strncat(hadMatName,"$",1);
      weights[i] = el->rhoz/wtNorm;
      float density;
      if(numElem == 1) density = mat.density();
      else             density = mat.density()* weights[i] * numNorm/el->pz;

      float lamda;
      switch( (int)el->z)
      {
       case 1:
           if(el->a < 1.5)     lamda = 790.;      // Hydrogen
            else                lamda = 342.;      // Deuterium
          break;
       case 2:  lamda = 478.f; break;     // Helium
       case 3:  lamda = 120.6f; break;    // Lithium
       case 4:  lamda = 36.7f;  break;    // Berillium
       case 6:  lamda = 49.9f;  break;    // Carbon
       case 7:  lamda = 99.4f;  break;    // Nitrogen
       case 10: lamda = 74.9f;  break;    // Neon
       case 13: lamda = 37.2f;  break;    // Aluminium
       case 26: lamda = 17.1f;  break;    // Iron
       case 29: lamda = 14.8f;  break;    // Copper
       case 74: lamda = 10.3f;  break;    // Tungsten
       case 82: lamda = 18.5f;  break;    // Lead
       case 92: lamda = 12.0f;  break;    // Uranium
       default:                                     // Everything else....
         lamda = 31.62 * pow(el->a, .335f) * density;
      }
      hadMatCounter++;
      indices[i] = hadMatCounter;
      float rlc =mat.radiationLength();

      gsmate(&hadMatCounter,/* hadMatName,*/ &el->a, &el->z, &density, &rlc,
             &lamda);

      i++;
   }

// Take care of mixtures here....
   if(numElem > 1) {
      strncpy(hadMatName, mat.name(), 18);
      strncat(hadMatName,"$",1);
      hadMatCounter++;
      float density = mat.density();
// cout << "Number of elements: " << numElem << '\n';
      gsmixt(&hadMatCounter, /* mat.name(),  */
               &numElem, weights, indices, &density);

   }

   hadMatIndex = hadMatCounter;
   delete [] weights;
   delete [] indices;

}

int GheishaData::hadMatCounter=0;

---