FairGeoMedium.cxx

//*-- ö : Ilse Koenig
//*-- Modified : 10/11/2003 by Ilse Koenig

////////////////////////////////////////////////////////////////////////////////
//
//  Class for tracking medium ( includes also material )
//
////////////////////////////////////////////////////////////////////////////////
#include "FairGeoMedium.h"

#include <iostream>
#include <cmath>
#include "stdlib.h"
using std::cout;
using std::log;
using std::pow;

ClassImp(FairGeoMedium)

FairGeoMedium::FairGeoMedium(const char* name)
  : TNamed(name,""),
    medId(0),
    autoflag(1),
    nComponents(0),
    weightFac(0),
    ca(NULL),
    cz(NULL),
    cw(NULL),
    density(0),
    radLen(0),
    sensFlag(0),
    fldFlag(0),
    fld(0),
    epsil(0),
    madfld(-1),
    maxstep(-1),
    maxde(-1),
    minstep(-1),
    npckov(0),
    ppckov(NULL),
    absco(NULL),
    effic(NULL),
    rindex(NULL)
{
  // Constructor for a medium with name and index id
//  SetName(name);
};

FairGeoMedium::~FairGeoMedium()
{
  // Destructor
  if (nComponents>0) {
    delete [] ca;
    ca=0;
    delete [] cz;
    cz=0;
    delete [] cw;
    cw=0;
    nComponents=0;
  }
  if (npckov>0) {
    delete [] ppckov;
    ppckov=0;
    delete [] absco;
    absco=0;
    delete [] effic;
    effic=0;
    delete [] rindex;
    rindex=0;
    npckov=0;
  }
}

void FairGeoMedium::setNComponents(Int_t n)
{
  // Sets the number of components in the material
  if (n==0) { return; }
  Int_t k=abs(n);
  if (nComponents!=0 && k!=nComponents) {
    delete [] ca;
    delete [] cz;
    delete [] cw;
    nComponents=0;
  }
  if (nComponents==0) {
    nComponents=k;
    ca=new Double_t[k];
    cz=new Double_t[k];
    cw=new Double_t[k];
  }
  weightFac=(Int_t)(n/nComponents);
}

Bool_t FairGeoMedium::setComponent (Int_t i,Double_t a,Double_t z,Double_t weight)
{
  // Defines the ith material component
  if (i<0||i>=nComponents) {
    Error("setNComponents","Wrong index");
    return kFALSE;
  }
  ca[i]=a;
  cz[i]=z;
  cw[i]=weight;
  return kTRUE;
}

void FairGeoMedium::getComponent(Int_t i,Double_t* p)
{
  // Returns the ith material component
  if (i>=0&&i<nComponents) {
    p[0]=ca[i];
    p[1]=cz[i];
    p[2]=cw[i];
//   cout << " -I p: " << p[0] << p[1] << p[2] << endl;
  } else { p[0]=p[1]=p[2]=0.; }
}

void FairGeoMedium::setNpckov (Int_t n)
{
  // Sets the number of optical parameters for the tracking of Cerenkov light
  if (n!=npckov && npckov>0) {
    delete [] ppckov;
    delete [] absco;
    delete [] effic;
    delete [] rindex;
  }
  npckov=n;
  if (n>0) {
    ppckov=new Double_t[npckov];
    absco=new Double_t[npckov];
    effic=new Double_t[npckov];
    rindex=new Double_t[npckov];
  }
}

Bool_t FairGeoMedium::setCerenkovPar(Int_t i,Double_t p,Double_t a,Double_t e ,Double_t r)
{
  // Defines the ith parameter set of the optical parameters
  if (i<0 || i>=npckov) {
    Error("setNpckov","Wrong index");
    return kFALSE;
  }
  ppckov[i]=p;
  absco[i]=a;
  effic[i]=e;
  rindex[i]=r;
  return kTRUE;
}

void FairGeoMedium::getCerenkovPar(Int_t i,Double_t* p)
{
  // returns the ith parameter set of the optical parameters
  if (i>=0&&i<npckov) {
    p[0]=ppckov[i];
    p[1]=absco[i];
    p[2]=effic[i];
    p[3]=rindex[i];
  } else { p[0]=p[1]=p[2]=p[3]=0.; }
}

void FairGeoMedium::setMediumPar(Int_t sensitivityFlag,Int_t fieldFlag,
                                 Double_t maxField,Double_t precision,Double_t maxDeviation,
                                 Double_t maxStep,Double_t maxDE,Double_t minStep )
{
  // Sets the medium parameters
  sensFlag=sensitivityFlag;
  fldFlag=fieldFlag;
  fld=maxField;
  epsil=precision;
  madfld=maxDeviation;
  maxstep=maxStep;
  maxde=maxDE;
  minstep=minStep;
}

void FairGeoMedium::getMediumPar(Double_t* params)
{
  // Returns the medium parameters
  params[0]=sensFlag;
  params[1]=fldFlag;
  params[2]=fld;
  params[3]=madfld;
  params[4]=maxstep;
  params[5]=maxde;
  params[6]=epsil;
  params[7]=minstep;
  params[8]=0.;
  params[9]=0.;
}

void FairGeoMedium::read(std::fstream& fin, Int_t aflag )
{
  // Reads the parameters from file
  autoflag=aflag;
  Int_t n;
  fin>>n;
  setNComponents(n);
  for(Int_t ik=0; ik<nComponents; ik++) {
    fin>>ca[ik];
  }
  for(Int_t i=0; i<nComponents; i++) {
    fin>>cz[i];
  }
  fin>>density;
  if (nComponents==1) {
    cw[0]=1.;
    calcRadiationLength();
  } else {
    for(Int_t i=0; i<nComponents; i++) {
      fin>>cw[i];
    }
  }
  fin>>sensFlag>>fldFlag>>fld>>epsil ;
  if (autoflag<1) { fin>>madfld>>maxstep>>maxde>>minstep; }
  else {
    //to use this feature one has to set TGeant3::SetAUTO(0), thus if the media does not
    // defined these values one can force Geant3 to calculate them by given them a value
    // of -1
    madfld=-1;
    maxstep=-1;
    maxde=-1;
    minstep=-1;
  }
  fin>>n;
  setNpckov(n);
  if (n>0) {
    for(Int_t i=0; i<n; i++) {
      fin>>ppckov[i]>>absco[i]>>effic[i]>>rindex[i];
    }
  }

}

void FairGeoMedium::print()
{
  // Prints the medium definition
  const char* bl="  ";
  cout<<GetName()<<'\n'<<nComponents* weightFac<<bl;
  for(Int_t ii=0; ii<nComponents; ii++) { cout<<ca[ii]<<bl ;}
  for(Int_t j=0; j<nComponents; j++) { cout<<cz[j]<<bl ;}
  cout<<density<<bl;
  if (nComponents<2) { cout<<radLen; }
  else for(Int_t iik=0; iik<nComponents; iik++) { cout<<cw[iik]<<bl ;}
  cout<<'\n'<<sensFlag<<bl<<fldFlag<<bl<<fld<<bl<<epsil<<'\n';
  if (autoflag<1) {
    cout<<madfld<<bl<<maxstep<<bl<<maxde<<bl<<minstep<<'\n';
  }
  cout<<npckov<<'\n';
  if (npckov>0) {
    for(Int_t i=0; i<npckov; i++) {
      cout<<ppckov[i]<< bl<<absco[i]<<bl<<effic[i]<<bl<<rindex[i]<<'\n';
    }
  }
  cout<<'\n';
}

void FairGeoMedium::write (std::fstream& fout)
{
  // Writes the medium definition into stream
  const char* bl="  ";
  fout<<GetName()<<'\n'<<nComponents* weightFac<<bl;
  for(Int_t i=0; i<nComponents; i++) { fout<<ca[i]<<bl ;}
  for(Int_t j=0; j<nComponents; j++) { fout<<cz[j]<<bl ;}
  fout<<density<<bl;
  if (nComponents<2) { fout<<radLen; }
  else for(Int_t i=0; i<nComponents; i++) { fout<<cw[i]<<bl ;}
  fout<<'\n'<<sensFlag<<bl<<fldFlag<<bl<<fld<<bl<<epsil<<'\n';
  if (autoflag<1) {
    fout<<madfld<<bl<<maxstep<<bl<<maxde<<bl<<minstep<<'\n';
  }
  fout<<npckov<<'\n';
  if (npckov>0) {
    for(Int_t i=0; i<npckov; i++) {
      fout<<ppckov[i]<< bl<<absco[i]<<bl<<effic[i]<<bl<<rindex[i]<<'\n';
    }
  }
  fout<<'\n';
}

Bool_t FairGeoMedium::calcRadiationLength()
{
  // calculates radiation length
  // formula in GEANT manual CONS110
  if (cz[0]<1.e-6) { // VAKUUM$
    radLen=1.e+16;
    return kTRUE;
  }
  Double_t alpha=1/137.;     // fine structure constant
  Double_t fac=.1912821;     // 4*((electron radius)**2)*(avogadro's number)
  Double_t z, a, w, az2, fc, y, xi, x0i, amol=0., x0itot=0.;
  if (weightFac>0) { amol=1.; }
  else {
    for (int i=0; i<nComponents; i++) {
      amol+=cw[i]*ca[i];
      if (amol==0.) {
        Error("calcRadiationLength()","amol==0 for medium %s",fName.Data());
        return kFALSE;
      }
    }
  }
  for (int i=0; i<nComponents; i++) {
    z=cz[i];
    a=ca[i];
    if(weightFac>0) {
      w=cw[i]/amol;
    } else {
      w=a*cw[i]/amol;
    }
    az2=alpha*alpha*z*z;
    fc=az2 * (1./(1.+az2) + 0.20206 - 0.0369*az2 + 0.0083*az2*az2
              - .002F*az2*az2*az2);
    y=log(183./pow(z,1./3.)) - fc;
    xi=(float)(log(1440./pow(z,2./3.)) / y);
    x0i=fac*alpha/a*z*(z+xi)*y;
    x0itot+=(x0i*w);
  }
  if (x0itot==0. || density==0.) {
    Error("calcRadiationLength()","x0itot=0 or density=0 for medium %s",fName.Data());
    return kFALSE;
  }
  radLen=1/density/x0itot;
  return kTRUE;
}