The Photon Absorption Ionisation model used for the ionisation in the TRT. More...

#include <TRT_PAI_Process.h>

Inheritance diagram for TRT_PAI_Process:

Collaboration diagram for TRT_PAI_Process:

Public Member Functions
	TRT_PAI_Process (const std::string &type, const std::string &name, const IInterface *parent)
	Not much action.
virtual StatusCode	initialize () override final
	Initialization of the PAI model:
virtual StatusCode	finalize () override final
virtual double	GetMeanFreePath (double scaledKineticEnergy, double squaredCharge) const override final
	Get the mean free path in gas (CLHEP units)
virtual double	GetEnergyTransfer (double scaledKineticEnergy, CLHEP::HepRandomEngine *rndmEngine) const override final
	Get the energy transferred from the charged particle to the gas (CLHEP units).

Private Member Functions
double	ScaledEkin2GamVarTab (double scaledKineticEnergy) const
	Converting Lorentz gamma to table index (well, double)

Private Attributes
const unsigned int	m_nTabulatedGammaValues {56}
const double	m_gamExpMin {-2.}
const double	m_gamExpMax {5.}
const double	m_deltaGamExp {}
std::vector< float >	m_en_array
std::vector< std::vector< float > >	m_fn_array
std::vector< float >	m_dndx
TRT_PAI_gasMixture *	m_trtgas {}
std::string	m_gasType {"Auto"}

Detailed Description

The Photon Absorption Ionisation model used for the ionisation in the TRT.

The main reference for this implementation is

V.M.Grishin, V.K.Ermilova, and S.K.Kotelnikov, NIM A307 (1991) 273-278.

Other references to the PAI model are

V.C.Ermilova, L.P.Kotenko, and G.I.Merzon, NIM145 (1977) 555-563.
J.H.Cobb, W.W.M.Allison, and J.N.Bunch, NIM133 (1976) 315-323.
W.W.M.Allison and J.H.Cobb, Ann.Rev.Nucl.Part.Sci (1980), 253-98.

Note: Internally this class uses the original units of the fortran version. To external world sees CLHEP units.

Author: Originally implemented in age FORTRAN by Pavel Nevski. Adapted to C++ by T.Kittelmann/Mogens Dam.

Definition at line 42 of file TRT_PAI_Process.h.

Constructor & Destructor Documentation

◆ TRT_PAI_Process()

TRT_PAI_Process::TRT_PAI_Process	(	const std::string &	type,
		const std::string &	name,
		const IInterface *	parent )

Not much action.

Initializes a few variables

Definition at line 33 of file TRT_PAI_Process.cxx.

  : base_class( type, name, parent )
  , m_deltaGamExp( (m_gamExpMax-m_gamExpMin)/m_nTabulatedGammaValues )
{
  //Properties:
  declareProperty( "GasType", m_gasType, "Gas Type" );
}

Member Function Documentation

◆ finalize()

StatusCode TRT_PAI_Process::finalize ( )

finaloverridevirtual

Definition at line 234 of file TRT_PAI_Process.cxx.

                                     {
  delete m_trtgas;
  return StatusCode::SUCCESS;
}

◆ GetEnergyTransfer()

double TRT_PAI_Process::GetEnergyTransfer	(	double	scaledKineticEnergy,
		CLHEP::HepRandomEngine *	rndmEngine ) const

finaloverridevirtual

Get the energy transferred from the charged particle to the gas (CLHEP units).

The energy transfer is randomly generated according to the distributions calculated based on the input data in gasdata.h

Parameters

scaledKineticEnergy The kinetic energy a proton would have had, if it had same Lorentz gamma factor as the particle in question.

Returns: Random energy transfer consistent with PAI model calculation.

Definition at line 278 of file TRT_PAI_Process.cxx.

                                                                                                            {
 
  double gv = ScaledEkin2GamVarTab( scaledKineticEnergy );
  unsigned int tabIndx = 0;
 
  if ( gv>0 ) tabIndx = static_cast<unsigned int>( gv );
 
  tabIndx = std::min( tabIndx, m_nTabulatedGammaValues-1 );
 
  // Generate a random number uniformly in [0,1].
  //  CLHEP::HepRandomEngine* pHRengine = m_pAtRndmGenSvc->GetEngine("TRT_PAI");
  double random = CLHEP::RandFlat::shoot(rndmEngine, 0., 1.);
 
  // What we are doing next is actually to select a value of E from
  // dN^2/dXdE considered as a function of E using the standard Monte
  // Carlo method. The total area under the function is equal to
  // dN/dX, and since we can not calculate the inverse antiderivative,
  // we have tabulated the antiderivative in (m_en_array, m_fn_array)
  // and taking the inverse is as simple as switching the roles of the
  // two arrays.
 
  double Etransf = TRT_PAI_utils::Interpolate (random * m_dndx[tabIndx], m_fn_array[tabIndx], m_en_array ) ;
 
  Etransf = exp(Etransf);
 
#ifndef NDEBUG
  ATH_MSG_VERBOSE ( "Energy transfer for scaledKineticEnergy = "
                    << scaledKineticEnergy << " is " << Etransf << " eV"
                    );
#endif
 
  return Etransf * CLHEP::eV;
}

◆ GetMeanFreePath()

double TRT_PAI_Process::GetMeanFreePath	(	double	scaledKineticEnergy,
		double	squaredCharge ) const

finaloverridevirtual

Get the mean free path in gas (CLHEP units)

Parameters

scaledKineticEnergy	The kinetic energy a proton would have had, if it had same Lorentz gamma factor as the particle in question
squaredCharge	Charge squared

Returns: Mean free path

Definition at line 247 of file TRT_PAI_Process.cxx.

                                                                    {
 
  if ( squaredCharge == 0. ) return 99999.0 * CLHEP::km;
 
  double gv = ScaledEkin2GamVarTab( scaledKineticEnergy );
  unsigned int index = 0;
  if ( gv > 0 ) index = static_cast<unsigned int >(gv);
 
  double d;
  if ( gv < 0 ) d = m_dndx[0];
  else if ( index >= m_nTabulatedGammaValues-1 ) d = m_dndx[m_nTabulatedGammaValues-1];
  else {
    double remain = gv-index;
    d = (1.-remain)*m_dndx[index] + remain*m_dndx[index+1];
  }
 
  double freepath = 1./(d*squaredCharge);
 
#ifndef NDEBUG
  ATH_MSG_VERBOSE ( "Mean free path for scaledKineticEnergy = " << scaledKineticEnergy
                    << " is " << freepath/CLHEP::mm << " mm " );
#endif
 
  // All internal calculations are performed without initialization with
  //   CLHEP units, so we put it here instead
 
  return freepath * CLHEP::cm;
}

◆ initialize()

StatusCode TRT_PAI_Process::initialize ( )

finaloverridevirtual

Initialization of the PAI model:

Construct chemical elements;
Construct gas components (molecules) from elements;
Construct gas mixture from gas components;
Construct "effective gas" from the gas mixture, i.e. calculate everything necessary for the current gas mixture to do the PAI simulation

Definition at line 43 of file TRT_PAI_Process.cxx.

                                       {
 
  using namespace TRT_PAI_gasdata;
  using namespace TRT_PAI_physicsConstants;
 
  ATH_MSG_VERBOSE ( "TRT_PAI_Process::initialize()" );
 
  // Supported gasType's are 70%/27%/03%:
  // 1) "Xenon"  - Xe/CO2/O2
  // 2) "Argon"  - Ar/CO2/O2
  // 3) "Kryton" - Kr/CO2/O2
  //
  // Figure out which type of gas we are having:
 
  ATH_MSG_DEBUG ( "Gastype read as " << m_gasType );
 
  std::string gasType = "NotSet";
  // Here m_gasType="Auto" was previously used to prompt a
  // check of InDetDD::TRT_DetectorManager for whether we want
  // the old Xenon mix: Xe/CO2/CF4 or the new one: Xe/CO2/O2
  // Nowadays, it is always the new one, so we don't check this anymore.
 
  if ( m_gasType == "Auto" )
    {
      gasType = "Xenon";
    }
  else if ( m_gasType == "Xenon" || m_gasType == "Argon" || m_gasType == "Krypton" )
    {
      ATH_MSG_DEBUG ( "Gastype is overriden from joboptions to be " << m_gasType );
      gasType = m_gasType;
    }
  else
    {
      ATH_MSG_FATAL ( "GasType property set to '" << m_gasType
                     << "'. Must be one of 'Auto', 'Xenon', 'Argon' or 'Krypton'." );
      return StatusCode::FAILURE;
    };
 
  ATH_MSG_INFO ( name() << " will use gas type = " << gasType );
 
  // Done with "trivial" initialization.
  // Finally we can start to construct gas
 
  ATH_MSG_VERBOSE ( "Constructing gas mixture." );
 
  // Need the gas temperature
  const double tempK = 289.;  // At 289. degrees, we get same densities as Nevski had
 
  // Define elements
  using element_type = std::map<std::string, TRT_PAI_element, std::less<std::string>>;
 
  element_type elements;
  elements["Xe"] = TRT_PAI_element( "Xe", EXe , SXe , NXe , ZXe,  AXe);
  elements["C"]  = TRT_PAI_element( "C" , EC  , SC  , NC  , ZC ,  AC );
  elements["F"]  = TRT_PAI_element( "F" , EF  , SF  , NF  , ZF ,  AF );
  elements["O"]  = TRT_PAI_element( "O" , EO  , SO  , NO  , ZO ,  AO );
  elements["Ar"] = TRT_PAI_element( "Ar", EAr , SAr , NAr , ZAr,  AAr);
  elements["Kr"] = TRT_PAI_element( "Kr", EKr , SKr , NKr , ZKr,  AKr);
 
  // Print out elements
  {
    for (element_type::iterator ei=elements.begin(); ei!=elements.end(); ++ei) {
      ATH_MSG_DEBUG ( ". Element "  << (*ei).second.getName()
                      << ", A= "    << (*ei).second.getAtomicA()
                      << ", Z= "    << (*ei).second.getAtomicZ()
                      << ", rho="   << (*ei).second.getDensity(tempK)
                      << " (" << tempK << " Kelvin)"
                      );
    }
  }
 
  // Define gas components
 
  using component_type = std::map<std::string, TRT_PAI_gasComponent, std::less<std::string>>;
 
  component_type components;
 
  components["Xe"]  = TRT_PAI_gasComponent("Xe");
  components["Xe"].addElement(&elements["Xe"],1);
 
  components["CO2"] = TRT_PAI_gasComponent("CO2");
  components["CO2"].addElement(&elements["C"],1);
  components["CO2"].addElement(&elements["O"],2);
 
  components["CF4"] = TRT_PAI_gasComponent("CF4");
  components["CF4"].addElement(&elements["C"],1);
  components["CF4"].addElement(&elements["F"],4);
 
  components["O2"]  = TRT_PAI_gasComponent("O2");
  components["O2"].addElement(&elements["O"],2);
 
  components["Ar"]  = TRT_PAI_gasComponent("Ar");
  components["Ar"].addElement(&elements["Ar"],1);
 
  components["Kr"]  = TRT_PAI_gasComponent("Kr");
  components["Kr"].addElement(&elements["Kr"],1);
 
  // Print out gas components
 
  {
    std::vector<std::string> cnam(6);
    cnam[0] = "Xe"; cnam[1] = "CO2"; cnam[2] = "CF4"; cnam[3] = "CO2"; cnam[4] = "Ar"; cnam[5] = "Kr";
 
    int n;
    for ( int ic=0; ic<6; ++ic ) {
      n = components[cnam[ic]].getNElementTypes();
      ATH_MSG_DEBUG ( ". Gas component " << components[cnam[ic]].getName() << " contains");
      for ( int ie=0; ie<n; ++ie ) {
        ATH_MSG_DEBUG ( "   - "        <<  components[cnam[ic]].getElementMultiplicity(ie)
                        << " atoms "   << (components[cnam[ic]].getElement(ie))->getName()
                        << " with Z= " << (components[cnam[ic]].getElement(ie))->getAtomicZ()
                        << " and A= "  << (components[cnam[ic]].getElement(ie))->getAtomicA()
                        );
      }
      ATH_MSG_DEBUG ( "   > density: " << components[cnam[ic]].getDensity(tempK) << " (" << tempK << " Kelvin)" );
    }
  }
 
  // Construct TRT gas mixture
  std::string mixtureName="TRT  Gas Mixture";
  mixtureName.insert(4, gasType);
  m_trtgas = new TRT_PAI_gasMixture(mixtureName);
 
  if ( gasType == "Xenon" ) {
    ATH_MSG_DEBUG ( "Using new Xenon gas mixture (Xe/CO2/O2 - 70/27/3)." );
    m_trtgas->addComponent( &components["Xe"] , 0.70);
    m_trtgas->addComponent( &components["CO2"], 0.27);
    m_trtgas->addComponent( &components["O2"] , 0.03);
  }
 
  if (( gasType == "XenonOld" )){
    ATH_MSG_DEBUG ( "Using old Xenon gas mixture (Xe/CO2/CF4 - 70/10/20)." );
    m_trtgas->addComponent( &components["Xe"] , 0.70);
    m_trtgas->addComponent( &components["CO2"], 0.10);
    m_trtgas->addComponent( &components["CF4"], 0.20);
  }
 
  if ( gasType == "Argon" ) {
    ATH_MSG_DEBUG ( "Using Argon gas mixture (Ar/CO2/O2 - 70/27/3)." );
    m_trtgas->addComponent( &components["Ar"] , 0.70);
    m_trtgas->addComponent( &components["CO2"], 0.27);
    m_trtgas->addComponent( &components["O2"] , 0.03);
  }
 
  if ( gasType == "Krypton" ) {
    ATH_MSG_DEBUG ( "Using Krypton gas mixture (Kr/CO2/O2 - 70/27/3)." );
    m_trtgas->addComponent( &components["Kr"] , 0.70);
    m_trtgas->addComponent( &components["CO2"], 0.27);
    m_trtgas->addComponent( &components["O2"] , 0.03);
  }
 
  m_trtgas->freezeGas();
 
  m_trtgas->showStructure();
 
  // Create and initialize effective gas
  // Sasha: to check: does Emin depends from Energy levels of active gas?
  // for Xenon lowest energetic level is 12.08 eV
  // for Argon - 15.83 eV. Does Emin should be different for Argon???
 
  const double Emin = 12.0;        // [eV]
  const double Emax = 1e+07;       // [eV]
                                   // Don't change Emax from 10 MeV without
                                   // talking to Pavel Nevski!
  const double eps   = 0.01;       // Epsilon for numeric integration
 
  TRT_PAI_effectiveGas effectiveGas(m_trtgas, Emin, Emax, tempK, eps);
 
  // Now tabulate...
 
  ATH_MSG_DEBUG ( "Making tables for various gamma factors." );
 
  std::vector<float> gamvec(m_nTabulatedGammaValues);
  std::vector<float> lnE;
 
  double gamvar = m_gamExpMin - 0.5*m_deltaGamExp;
  for ( unsigned int ig = 0; ig < m_nTabulatedGammaValues; ++ig ) {
    gamvar    += m_deltaGamExp;
    gamvec[ig] = 1. + pow(10.,gamvar);
  }
 
  // Here is the real action!
 
  effectiveGas.GasTab(gamvec, m_en_array, m_fn_array, m_dndx );
 
  ATH_MSG_DEBUG ( "Initialization completed." );
 
  return StatusCode::SUCCESS;
}

◆ ScaledEkin2GamVarTab()

double TRT_PAI_Process::ScaledEkin2GamVarTab ( double scaledKineticEnergy ) const

inlineprivate

Converting Lorentz gamma to table index (well, double)

Definition at line 240 of file TRT_PAI_Process.cxx.

        {
  using namespace TRT_PAI_physicsConstants;
  return (log(scaledKineticEnergy/MProtonMeV)*invlog10 - m_gamExpMin)/m_deltaGamExp;
}

Member Data Documentation

◆ m_deltaGamExp

const double TRT_PAI_Process::m_deltaGamExp {}

private

Definition at line 92 of file TRT_PAI_Process.h.

92{};

◆ m_dndx

std::vector<float> TRT_PAI_Process::m_dndx

private

Definition at line 97 of file TRT_PAI_Process.h.

◆ m_en_array

std::vector<float> TRT_PAI_Process::m_en_array

private

Definition at line 95 of file TRT_PAI_Process.h.

◆ m_fn_array

std::vector< std::vector<float> > TRT_PAI_Process::m_fn_array

private

Definition at line 96 of file TRT_PAI_Process.h.

◆ m_gamExpMax

const double TRT_PAI_Process::m_gamExpMax {5.}

private

Definition at line 91 of file TRT_PAI_Process.h.

91{5.}; // Max Lorentz gamma: 1+10^(m_gamExpMax)

◆ m_gamExpMin

const double TRT_PAI_Process::m_gamExpMin {-2.}

private

Definition at line 90 of file TRT_PAI_Process.h.

90{-2.}; // Min Lorentz gamma: 1+10^(m_gamExpMin)

◆ m_gasType

std::string TRT_PAI_Process::m_gasType {"Auto"}

private

Definition at line 106 of file TRT_PAI_Process.h.

106{"Auto"};

◆ m_nTabulatedGammaValues

const unsigned int TRT_PAI_Process::m_nTabulatedGammaValues {56}

private

Definition at line 89 of file TRT_PAI_Process.h.

89{56};

◆ m_trtgas

TRT_PAI_gasMixture* TRT_PAI_Process::m_trtgas {}

private

Definition at line 98 of file TRT_PAI_Process.h.

98{};

The documentation for this class was generated from the following files:

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ TRT_PAI_Process()

Member Function Documentation

◆ finalize()

◆ GetEnergyTransfer()

◆ GetMeanFreePath()

◆ initialize()

◆ ScaledEkin2GamVarTab()

Member Data Documentation

◆ m_deltaGamExp

◆ m_dndx

◆ m_en_array

◆ m_fn_array

◆ m_gamExpMax

◆ m_gamExpMin

◆ m_gasType

◆ m_nTabulatedGammaValues

◆ m_trtgas