---

Material.cxx

Go to the documentation of this file.

// $Id: Material.cxx,v 1.2 2001/08/08 05:51:11 chehtman Exp $
//
// This file is part of Gismo 2
//
// Contents ----------------------------------------------------------------
//
//  Material::hatch
//  Material::Material(std::istream& fin, int iray)
//  void Material::addElement(char* name, float Z,
//     float Material::dedx(float beta, float eta)
//
//  AtomicElement::AtomicElement(char* pname, float pZ, float pA, float f,
//      float den, AtomicElement* pnext)

// Description
//
//   implementation of class Material and struct Material::AtomicElement
//
// End --------------------------------------------------------------------

#include "gismo/Material.h"

#include "gismo/Units.h"
#include "facilities/error.h"
#include <math.h>
#include <cstdio>


#include <strstream>
#include <fstream>
#include <iomanip>


// ----------------allocate static variables ----------------------------
Material* Material::firstMaterial=0;
Material* Material::vacuum = 0;


//--------------------------- path manipulation code, data ---------------------
#ifndef DEFPEGS4PATH
#  define DEFPEGS4PATH "PEGS4"
#endif

static std::string* defaultPath()
{
    std::string defPath(DEFPEGS4PATH);

    char *penv = getenv("PEGS4PATH");
    std::string *path;
    if (penv != NULL && strlen(penv) > 0)  {
          path = new std::string(penv);
          if (defPath.length() > 0) *path += ':';
    }
    else
        path = new std::string;

    *path += defPath;
    return path;
}

std::string *Material::s_dataPath = defaultPath();


void Material::dataPath(const char* path)
{
    if ( s_dataPath ) delete s_dataPath;
    s_dataPath = new std::string(path);
}

void Material::addPath(const char* path)
{
    *s_dataPath = std::string(path) + ':' + *s_dataPath;
}

const char *Material::getDataPath()
{
    return s_dataPath->c_str();
}
//---------------------- function OpenReadFile --------------------------
static std::istream * OpenReadFile(const std::string & fileName)
{
    std::ifstream * is;
  // check for absolute pathname

    if ( fileName[0] == '/' ) {
         is = new std::ifstream(fileName.c_str(), std::ios::in);
         if ( ! is->is_open() )
         {
                FATAL((std::string("Material description file '")+fileName+std::string("' does not exist.")).c_str());
                delete is;
         }
         return is;
    }

    // next, see if we can find it in the search path

    std::string fullFileName;
    std::istrstream path(Material::getDataPath());
    char directory[256];                // should be FILENAME_MAX, but gcc stdio.h
                                // for a Sun does not have it.

    while ( path ) {
         path.getline(directory, 256, ':');     //THB: changed get to getline
         fullFileName = std::string(directory) + '/' + fileName;
         is = new std::ifstream(fullFileName.c_str(), std::ios::in);
         if ( is->is_open() )  //note standard library form
            return is;
         else
            delete is;
    }

    // maybe it is in the current directory

    is = new std::ifstream(fileName.c_str(), std::ios::in);
    if ( is->is_open() )  //note standard library form
     return is;
    else         delete is;


    FATAL((std::string("PEGS4 file '") + fileName + std::string("' not found.\nPath was ")
        + std::string(Material::getDataPath())).c_str() );
    return is;
}

// -----------------member function Material::hatch ----------------------
Material*
Material::hatch(const char* name)
{
    if( vacuum==0) vacuum = new Material();
    std::string sname(name);
    if( sname.length()==0
          || sname == "vacuum" )
        return vacuum;

    // check to see if the material has been hatched already

    Material * m = firstMaterial;
    for(; m; m = m->nextMaterial)
        if( std::string(name) ==  m->m_fileName)
            return m;

    std::string nameExt = sname + ".mat";
    std::istream*  fin = OpenReadFile(nameExt);

    // create a new Material and link it into the list of Materials
    // (should now use a map: THB)
    m = new Material(*fin, name);
    delete fin;

    if( firstMaterial ){
        Material* mat;
        for( mat=firstMaterial;  mat->nextMaterial;   mat = mat->nextMaterial);
            mat->nextMaterial = m;
    }
    else {
        firstMaterial = m;
    }
    m->nextMaterial = 0;
        
    return m;
}
// ------------constructor Material::Material, from stream --------------
Material::Material(std::istream& fin, const char* fname)
: m_fileName(fname)
, energyLoss1(0)
{
        // start reading file, in PEGS4 format, example follows
/*
 MEDIUM=Liquid Argon            ,STERNCID=LIQUID ARGON
 ELEM,RHO= 1.4000E+00,NE= 1
 ASYM=Ar,Z=18.,A=   39.948,PZ= 1.00000E+00,RHOZ= 3.99480E+01
    1.39648E+01   1.50000E+00   1.00000E-01   1.00000E+05   1.00000E+05
     0   63    0  150    0    0    0    0    0
*/
    char linebuf[90], type[7];

    fin.getline(linebuf,90);
//    sscanf(linebuf, " %6s=%15s", type, materialName);
    sscanf(linebuf, " %6s=%15[^,]", type, materialName); //materialName up to comma
    if( strcmp("MEDIUM",type) != 0) {
        std::ostrstream message;
        message << "Material file format for " << m_fileName << " is unrecognizable: should start MEDIUM=\n";
        message << "But first line was: '" << linebuf << "'\0";
        FATAL( message.str() );
    }
    int nne=0;

    fin.getline(linebuf,90);
    sscanf(linebuf, " %4s,RHO=%f,NE=%i", type, &rho, &nne);
    if( nne<0 || nne > 10 )  {
        std::ostrstream message;
        message << "Material file format problem: bad number of elements\n";
        message << "data line: '" << linebuf << "'\0";
        FATAL( message.str() );
    }


#if 1
    std::cout << "Material name: " << materialName
       << ", Type: "<< type
       << ", rho=" << rho
       << ", # elements: " << nne << '\n';
#endif

  // here loop over nne

    float pznorm = 0;
    zeff = aeff = 0;

    for (int ie = 0; ie<nne; ie++) {
        fin.getline(linebuf,90);
        AtomicElement* ae = addElement(linebuf);
        pznorm += ae->pz;
        zeff   += ae->pz * ae->z;
        aeff   += ae->pz * ae->a;

#if 1
    std::cout  << "\tElement#" << (ie+1) << ": " << ae->name
              << ", Z="     << ae->z
            << ", A="     << ae->a
            << ", pz="    << ae->pz
            << ", rhoz="  << ae->rhoz << '\n';
#endif
    }

    zeff /= pznorm;
    aeff /= pznorm;


  // the first quantitiy after the atomic element data is the radiation length. This is used
  // to calculate multiple scattering and must be done if egs is not loaded

    fin >> rlc;

  // finally, read all the cross section data, using copies of
  // mat data objects in Vacuum

    MatData_list * vaclist = &vacuum->matDataList;

    for( MatData_list::iterator it=vaclist->begin(); it<vaclist->end(); ++it){
        MatData* mdv = *it;
        MatData* md = mdv->copy();
        md->read( fin, *this );
        matDataList.push_back( md );
    }

}

// -------------member functions Material::addElement ---------------

Material::AtomicElement*
Material::addElementData(char* name, float Z,
                        float A, float fraction, float density)
{
    AtomicElement* ae = new AtomicElement(name,Z,A,fraction,density);
    elementList.push_back(ae);
    return ae;

}
// ASYM=Ar,Z=18.,A=   39.948,PZ= 1.00000E+00,RHOZ= 3.99480E+01

Material::AtomicElement*
Material::addElement(const char* buffer)
{
    char name[3];
    float Z,A,f,d;
    sscanf(buffer, " ASYM=%2c, Z=%f, A=%f, PZ=%f, RHOZ=%f",
          name, &Z, &A, &f, &d);
    name[2] = '\0';
    return addElementData(name,Z,A,f,d);
}

// -----------------member function Material::dedx ---------------------

// This formula,and the calculations of the constants,
// are adapted from the PDG section on Passage of Particle through Matter

float Material::dedx(float beta, float eta) const
{

    if( !energyLoss1 ) {
        Material* This = const_cast<Material*>(this); // cast from const
        This->energyLoss1= Units::fromMeV(0.3071f) * Z()/A() * density();
        if( Z()>1 )  This->IP = 16. * pow((Z()), 0.9F);
        else         This->IP = 19.2f;   // for H2 gas: 21.8 for liquid

        This->energyLoss2= 1.022e6 / IP ;  // const is 2 * me in eV
    This->energyLoss3= log(IP/(28.816*sqrt(density()))); //density correction
    }


    float arg2,   arg1 = eta*eta * energyLoss2*beta;
    if(arg1 > 1)  arg2 = log(arg1)- beta*beta
        -0.5*density_correction(eta);
    else          arg2 = 1.;
    if(arg2 < 0.) arg2 = 1.;

    return energyLoss1/(beta*beta) *arg2;

}

float Material::density_correction(float eta) const {
    float x0,x1;
    float cc = 1.0+ 2.0*energyLoss3;

    if(density() > 0.05)
// condensed material   
    {
        if(IP < 100.){ x0 = (cc<3.681) ? 0.2 : 0.326*cc-1.; x1=2.;}
        else         { x0 = (cc<5.215) ? 0.2 : 0.326*cc-1.5; x1=3.;}

    }else{
// gas
        if(cc<12.25){
            int i=int((cc-10)/0.5)+1;
            if(i<0)i=0;
            if(i>4)i=4;
            x0=1.6+0.1*i; x1=4.;
        }else{
            if(cc < 13.804){
                x0=2.;x1=5.;
            }else{
                x0=0.326*cc-2.5; x1=5;
            }
        }
    }

    float xa = cc/4.606; float xm = 3.;
    float aa = 4.606*(xa-x0)/pow(x1-x0,xm);
    float x = log10(eta);
    float delta = 0;
    if(x>x0){        
        delta=4.606*x-cc;
        if(x<x1)delta+=aa*pow(x1-x,xm);
    }

    return delta;

}



// -------------static member function Material::declareMatData----------
unsigned
Material::declareMatData(MatData* mdnew)
{
  if( Material::vacuum ==0 )
    Material::vacuum = new Material();

  if( firstMaterial )
     FATAL("Material: cannot add new matdata object after a call to hatch");
  // add this object to the list attached to Vacuum
  vacuum->matDataList.push_back(mdnew);
  return vacuum->matDataList.size()-1;
}

// ---------------- constructor AtomicElement::AtomicElement ----------

Material::AtomicElement::AtomicElement(char* pname, float pZ, float pA, float f,        float den)
  : z(pZ)
  , a(pA)
  , pz(f)
  , rhoz(den)
{
  strncpy(name, pname, 3);
}
// -------------static member function Material::printAll----------
void
Material::printAll(std::ostream& out)
{
    out << "\n\t===========List of defined materials======================";

    Material * mat = Material::firstMaterial;
    if( ! mat ) {
        out << "\n====================\n"
            <<   "No materials defined\n"
            <<   "====================\n";
        return;
    }
    unsigned i, n = mat->matDataList.size();

    // first header line:
    out<< '\n' << std::setw (36) << "                                    ";
    for(i=0; i< n ; i++ )
        mat->matDataList[i]->printHead1(out);

    // second header line:
    out << '\n' <<std::setw(36) <<  "            name    Z    A   density  ";
    //                          12345678901234567890123456789012345678
    for(i=0; i< n ; i++ )
        mat->matDataList[i]->printHead2(out);

    out << '\n';

    while (mat ) {
        mat->printOn(out);
        mat = mat->nextMaterial;
    }
}
void
Material::printOn(std::ostream& out)
{
    out << std::setw(16) << m_fileName.c_str() //title()
        << std::setw(6) << std::setprecision(3) << Z()
        << std::setw(6) << std::setprecision(3) << A()
        << std::setw(8) << density();
    for(unsigned i=0; i< matDataList.size() ; i++ )
        matDataList[i]->printOn(out);

    out << '\n';
}

Material::~Material()
{

    unsigned i, n = matDataList.size();

    for(i=0; i< n ; i++ )
        delete matDataList[i];

    if( this==vacuum) return;

    Element_list::iterator eit = elementList.begin();
    while (eit != elementList.end() )
        delete (*eit++);

}

void Material::deleteAll()
{
    Material * mat = Material::firstMaterial;

    while (mat ) {
        Material* next = mat->nextMaterial;
        delete mat;
        mat = next;
    }
    delete s_dataPath;
    s_dataPath=0;
    delete vacuum;
}

---