TRTTransitionRadiation Class Reference

#include <TRTTransitionRadiation.h>

Inheritance diagram for TRTTransitionRadiation:

Collaboration diagram for TRTTransitionRadiation:

Public Member Functions
	TRTTransitionRadiation (const G4String &processName="TransitionRadiation", const std::string &xmlfilename="TRgeomodelgeometry.xml")
virtual	~TRTTransitionRadiation ()
void	Initialize ()
G4bool	IsApplicable (const G4ParticleDefinition &)
G4double	GetMeanFreePath (const G4Track &aTrack, G4double, G4ForceCondition *)
G4VParticleChange *	PostStepDoIt (const G4Track &aTrack, const G4Step &aStep)
void	AddRadiatorParameters (const TRTRadiatorParameters &p)
bool	msgLvl (const MSG::Level lvl) const
	Test the output level.
MsgStream &	msg () const
	The standard message stream.
MsgStream &	msg (const MSG::Level lvl) const
	The standard message stream.
void	setLevel (MSG::Level lvl)
	Change the current logging level.

Private Member Functions
	TRTTransitionRadiation (const TRTTransitionRadiation &)
TRTTransitionRadiation &	operator= (const TRTTransitionRadiation &)
G4double	ComputePhotoAbsCof (const G4Material *Material, const G4double &GammaEnergy) const
G4double	XEmitanNevski (const G4double &PhotonEnergy, const G4double &Gamma, const G4double &GasThickness, const G4double &FoilThickness) const
G4double	XEmitanArtru (const G4double &PhotonEnergy, const G4double &Gamma, const G4double &GasThickness, const G4double &FoilThickness) const
G4double	NeffNevski (const G4double &sigGas, const G4double &sigFoil, const G4double &GasThickness, const G4double &FoilThickness, const G4int &FoilsTraversed) const
G4double	NeffArtru (const G4double &sigGas, const G4double &sigFoil, const G4double &GasThickness, const G4double &FoilThickness, const G4int &FoilsTraversed) const
G4double	XFinter (const G4double &X, const G4double *A, const G4double *F) const
G4double	XInteg (const G4double *yy, G4double *ss) const
void	initMessaging () const
	Initialize our message level and MessageSvc.

Private Attributes
TRRegionXMLHandler *	m_XMLhandler
const std::string	m_xmlfilename
std::vector< TRTRadiatorParameters >	m_radiators
G4double	m_MinEnergyTR
G4double	m_MaxEnergyTR
G4int	m_NumBins
G4double	m_WplasmaGas
G4double	m_WplasmaFoil
G4double	m_GammaMin
G4double	m_EkinMin
G4double *	m_Ey
G4double *	m_Sr
G4double *	m_om
G4double *	m_Omg
G4double *	m_sigmaGas
G4double *	m_sigmaFoil
std::string	m_nm
	Message source name.
boost::thread_specific_ptr< MsgStream >	m_msg_tls
	MsgStream instance (a std::cout like with print-out levels)
std::atomic< IMessageSvc * >	m_imsg { nullptr }
	MessageSvc pointer.
std::atomic< MSG::Level >	m_lvl { MSG::NIL }
	Current logging level.
std::atomic_flag m_initialized	ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
	Messaging initialized (initMessaging)

Detailed Description

Definition at line 20 of file TRTTransitionRadiation.h.

Constructor & Destructor Documentation

TRTTransitionRadiation() [1/2]

TRTTransitionRadiation::TRTTransitionRadiation	(	const G4String &	processName = "TransitionRadiation",
		const std::string &	xmlfilename = "TRgeomodelgeometry.xml" )

Definition at line 48 of file TRTTransitionRadiation.cxx.

                                                                                                          :
  G4VDiscreteProcess(processName,fElectromagnetic),
  AthMessaging("TRTTransitionRadiation"),
  m_XMLhandler(nullptr),m_xmlfilename(xmlfilename),
  m_MinEnergyTR(0.0),m_MaxEnergyTR(0.0),m_NumBins(0),m_WplasmaGas(0.0),
  m_WplasmaFoil(0.0),m_GammaMin(0.0),m_EkinMin(0.0),m_Ey(nullptr),m_Sr(nullptr),
  m_om(nullptr),m_Omg(nullptr),m_sigmaGas(nullptr),m_sigmaFoil(nullptr)
{
  m_radiators.clear();
  m_XMLhandler = new TRRegionXMLHandler( this );
 
  SetProcessSubType(  fTransitionRadiation );
 
  ATH_MSG_DEBUG ( "##### Constructor TRTTransitionRadiation" );
 
  Initialize();
 
  ATH_MSG_DEBUG ( "##### Constructor TRTTransitionRadiation done" );
 
}

~TRTTransitionRadiation()

TRTTransitionRadiation::~TRTTransitionRadiation ( )

virtual

Definition at line 71 of file TRTTransitionRadiation.cxx.

                                                {
  delete [] m_Ey;
  delete [] m_Sr;
  delete [] m_Omg;
  delete [] m_om;
  delete [] m_sigmaGas;
  delete [] m_sigmaFoil;
}

TRTTransitionRadiation() [2/2]

TRTTransitionRadiation::TRTTransitionRadiation ( const TRTTransitionRadiation & )

private

Member Function Documentation

AddRadiatorParameters()

void TRTTransitionRadiation::AddRadiatorParameters

(

const TRTRadiatorParameters &

p

)

Definition at line 81 of file TRTTransitionRadiation.cxx.

                                                                                 {
 
  ATH_MSG_DEBUG(" New Radiator parameters being defined for TR process");
  ATH_MSG_DEBUG(" Volume " << p.GetLogicalVolume()->GetName());
  ATH_MSG_DEBUG(" FoilThickness " << p.GetFoilThickness());
  ATH_MSG_DEBUG(" GasThickness  " << p.GetGasThickness());
 
  m_radiators.push_back(p);
 
  return;
}

ComputePhotoAbsCof()

G4double TRTTransitionRadiation::ComputePhotoAbsCof ( const G4Material * Material, const G4double & GammaEnergy ) const

private

Definition at line 658 of file TRTTransitionRadiation.cxx.

                                                                                           {
 
  const G4double *SandiaCof(Material->GetSandiaTable()->GetSandiaCofForMaterial(PhotonEnergy));
 
  const G4double Energy2(PhotonEnergy*PhotonEnergy);
  const G4double Energy3(PhotonEnergy*Energy2);
  const G4double Energy4(Energy2*Energy2);
 
  return ( SandiaCof[0]/PhotonEnergy +
           SandiaCof[1]/Energy2 +
           SandiaCof[2]/Energy3 +
           SandiaCof[3]/Energy4 );
}

GetMeanFreePath()

G4double TRTTransitionRadiation::GetMeanFreePath	(	const G4Track &	aTrack,
		G4double	,
		G4ForceCondition *	condition )

Definition at line 246 of file TRTTransitionRadiation.cxx.

                                                                              {
  const G4DynamicParticle *aParticle(aTrack.GetDynamicParticle());
  if(std::abs(aParticle->GetDefinition()->GetPDGCharge())< std::numeric_limits<double>::epsilon() ||
     std::abs(aParticle->GetDefinition()->GetPDGMass())<std::numeric_limits<double>::epsilon()  ) {
    *condition = NotForced;
  } else {
    *condition = Forced;
  }
  return DBL_MAX;
}

Initialize()

void TRTTransitionRadiation::Initialize

(

)

Definition at line 94 of file TRTTransitionRadiation.cxx.

                                        {
  ATH_MSG_DEBUG ( "Constructor TRTTransitionRadiation::Initialize()" );
  m_NumBins     = 100;         // Bins for TR generation (Nevski has 25 bins)
  m_MinEnergyTR =   2.0*CLHEP::keV;   // Minimum TR energy
  m_MaxEnergyTR =  50.0*CLHEP::keV;   // Maximum TR energy
 
  m_GammaMin = 500.;          // Min. value of Lorentz gamma for TR generation
  m_EkinMin  = (m_GammaMin-1.)*CLHEP::electron_mass_c2;
 
  //We get the materials from GeoMaterial manager and convert them to G4Materials
  //
  //NB: Maybe we should control this from TRT_GeoModel instead?
  //
  //    On a related note we should get rid of the xml files too -
  //    they are a loophole in the versioning.
 
  G4Material* g4mat_Gas{nullptr};
  G4Material* g4mat_FoilMaterial{nullptr};
  std::string errorMessage;
 
  // Get material information from storegate.
  ISvcLocator *svcLocator = Gaudi::svcLocator(); // from Bootstrap
  SmartIF<StoreGateSvc> detStore{svcLocator->service( "DetectorStore")};
  if (!detStore) {
    errorMessage = "Can not access Detector Store";
    ATH_MSG_FATAL (errorMessage);
    throw std::runtime_error(errorMessage);
  }
 
  StoredMaterialManager* materialManager = detStore->tryRetrieve<StoredMaterialManager>("MATERIALS");
  if(materialManager) {
    Geo2G4MaterialFactory geo2g4_material_fact;//Note - this is a very lightweight class!
    g4mat_Gas = geo2g4_material_fact.Build(materialManager->getMaterial("trt::CO2"));
    g4mat_FoilMaterial = geo2g4_material_fact.Build(materialManager->getMaterial("trt::CH2"));
 
    //Note: One might wonder if we now own g4mat_Gas and
    //g4mat_FoilMaterial, but since the elementfactory used inside
    //Geo2G4MaterialFactory is static, we must NOT delete them!!
  }
  else {
    ATH_MSG_INFO("GeoModel Material is not available. Construct TR materials using IRDBAccessSvc interface");
    SmartIF<IGeoDbTagSvc> geoDbTagSvc{svcLocator->service("GeoDbTagSvc")};
    if (!geoDbTagSvc) {
      errorMessage = "Can not access GeoDbTagSvc";
      ATH_MSG_FATAL (errorMessage);
      throw std::runtime_error(errorMessage);
 
      return;
    };
    SmartIF<IRDBAccessSvc> pAccessSvc{svcLocator->service(geoDbTagSvc->getParamSvcName())};
    if ( !pAccessSvc) {
      errorMessage = "Can not access " + geoDbTagSvc->getParamSvcName();
      ATH_MSG_FATAL (errorMessage);
      throw std::runtime_error(errorMessage);
    }
 
    std::map<std::string,int> materialComponentsMap; // Number of elements for each material
    std::map<std::string,G4Material*> materialMap;
    IRDBRecordset_ptr materialsRec = pAccessSvc->getRecordsetPtr("TRMaterials","","");
    // Step #1. Count elements per material
    for (const IRDBRecord_ptr& material : *materialsRec) {
      std::string key = material->getString("NAME");
      auto mapIt = materialComponentsMap.find(key);
      if (mapIt == materialComponentsMap.end()) {
        materialComponentsMap.emplace(key,1);
      }
      else {
        materialComponentsMap.at(key) = mapIt->second + 1;
      }
    }
    // Step #2. Build materials
    for (const IRDBRecord_ptr& material :*materialsRec) {
      std::string key = material->getString("NAME");
      G4Material* g4Material{nullptr};
      if(materialMap.find(key)==materialMap.end()) {
        auto itNElements = materialComponentsMap.find(key);
        if (itNElements==materialComponentsMap.end()) continue;
        g4Material = new G4Material(key
                                    ,material->getDouble("DENSITY")*(CLHEP::gram / CLHEP::cm3)
                                    ,itNElements->second);
        materialMap.emplace(key,g4Material);
      }
      else {
        g4Material = materialMap.at(key);
      }
 
      G4Element* g4Element = G4Element::GetElement(material->getString("ELEMENT_NAME"));
      if(!g4Element) {
        errorMessage = "Wrong name for element " + material->getString("ELEMENT_NAME")
          + " found in the description of material " + key;
        ATH_MSG_FATAL (errorMessage);
        throw std::runtime_error(errorMessage);
      }
      g4Material->AddElement(g4Element,(G4int)material->getInt("N"));
    }
    g4mat_Gas = materialMap.at("CO2");
    g4mat_FoilMaterial = materialMap.at("CH2");
  }
 
  G4double PlasmaCof(4.0*M_PI*CLHEP::fine_structure_const*CLHEP::hbarc*CLHEP::hbarc*CLHEP::hbarc/CLHEP::electron_mass_c2);
 
  // Plate and gas indicies
  m_WplasmaFoil = sqrt( PlasmaCof * g4mat_FoilMaterial->GetElectronDensity() );
  m_WplasmaGas  = sqrt( PlasmaCof * g4mat_Gas->GetElectronDensity() );
 
  ATH_MSG_DEBUG("Foil material : " << g4mat_FoilMaterial->GetName()
                << "; plasma energy : " << m_WplasmaFoil/CLHEP::eV << " eV");
  ATH_MSG_DEBUG("Gas : " << g4mat_Gas->GetName()
                << "; plasma energy : " << m_WplasmaGas/CLHEP::eV << " eV");
 
  m_Omg       = new G4double[m_NumBins];
  m_om        = new G4double[m_NumBins];
  m_Ey        = new G4double[m_NumBins];
  m_Sr        = new G4double[m_NumBins];
  m_sigmaGas  = new G4double[m_NumBins];
  m_sigmaFoil = new G4double[m_NumBins];
 
  for(G4int j(0); j<m_NumBins; ++j) {
    m_Omg[j] = log( m_MinEnergyTR ) +
      j*log( m_MaxEnergyTR/m_MinEnergyTR )/(m_NumBins-1);
    m_om [j] = exp( m_Omg[j] );
    m_sigmaGas[j]  = ComputePhotoAbsCof( g4mat_Gas, m_om[j] );
    m_sigmaFoil[j] = ComputePhotoAbsCof( g4mat_FoilMaterial, m_om[j] );
  }
 
  const std::string file=PathResolver::find_file(m_xmlfilename,"DATAPATH");
  if (!file.empty())
    {
      m_XMLhandler->Process(file);
    }
  else ATH_MSG_WARNING("File "<<m_xmlfilename<<" not found! Expect problems.");
 
  ATH_MSG_DEBUG ( "Constructor TRTTransitionRadiation::Initialize() DONE!" );
  return;
 
}

initMessaging()

void AthMessaging::initMessaging ( ) const

privateinherited

Initialize our message level and MessageSvc.

This method should only be called once.

Definition at line 39 of file AthMessaging.cxx.

{
  m_imsg = Athena::getMessageSvc();
  // If user did not set an explicit level, set a default
  if (m_lvl == MSG::NIL) {
    m_lvl = m_imsg ?
      static_cast<MSG::Level>( m_imsg.load()->outputLevel(m_nm) ) :
      MSG::INFO;
  }
}

IsApplicable()

G4bool TRTTransitionRadiation::IsApplicable

(

const G4ParticleDefinition &

particle

)

Definition at line 236 of file TRTTransitionRadiation.cxx.

                                                                                {
  //return true if PDG Charge and PDG Mass are non-zero.
  return  ( std::abs(particle.GetPDGCharge())>std::numeric_limits<double>::epsilon() && std::abs(particle.GetPDGMass())>std::numeric_limits<double>::epsilon() ) ;
}

msg() [1/2]

MsgStream & AthMessaging::msg ( ) const

inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 163 of file AthMessaging.h.

{
  MsgStream* ms = m_msg_tls.get();
  if (!ms) {
    if (!m_initialized.test_and_set()) initMessaging();
    ms = new MsgStream(m_imsg,m_nm);
    m_msg_tls.reset( ms );
  }
 
  ms->setLevel (m_lvl);
  return *ms;
}

msg() [2/2]

MsgStream & AthMessaging::msg ( const MSG::Level lvl ) const

inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 178 of file AthMessaging.h.

{ return msg() << lvl; }

msgLvl()

bool AthMessaging::msgLvl ( const MSG::Level lvl ) const

inlineinherited

Test the output level.

Parameters

lvl	The message level to test against

Returns

boolean Indicating if messages at given level will be printed

Return values

true	Messages at level "lvl" will be printed

Definition at line 151 of file AthMessaging.h.

{
  if (m_lvl <= lvl) {
    msg() << lvl;
    return true;
  } else {
    return false;
  }
}

NeffArtru()

G4double TRTTransitionRadiation::NeffArtru ( const G4double & sigGas, const G4double & sigFoil, const G4double & GasThickness, const G4double & FoilThickness, const G4int & FoilsTraversed ) const

inlineprivate

Definition at line 468 of file TRTTransitionRadiation.cxx.

                                                                                           {
  if ( FoilsTraversed <= 0 )
    return 0.;
 
  G4double sigma = FoilThickness*sigFoil + GasThickness*sigGas;
  G4double xNF = FoilsTraversed;
  return (1.-exp(-xNF*sigma))/(1.-exp(-sigma));
}

NeffNevski()

G4double TRTTransitionRadiation::NeffNevski ( const G4double & sigGas, const G4double & sigFoil, const G4double & GasThickness, const G4double & FoilThickness, const G4int & FoilsTraversed ) const

inlineprivate

Definition at line 451 of file TRTTransitionRadiation.cxx.

                                                                                            {
  if ( FoilsTraversed <= 0 )
    return 0.;
 
  G4double C1 = sigGas * GasThickness*FoilsTraversed;
  G4double C2 = sigFoil*FoilThickness*FoilsTraversed;
 
  G4double W1 = exp(-C1);
  G4double W2 = exp(-C2);
 
  return FoilsTraversed * (1-W1*W2)/(C1+C2);
}

operator=()

TRTTransitionRadiation & TRTTransitionRadiation::operator= ( const TRTTransitionRadiation & )

private

PostStepDoIt()

G4VParticleChange * TRTTransitionRadiation::PostStepDoIt	(	const G4Track &	aTrack,
		const G4Step &	aStep )

Definition at line 262 of file TRTTransitionRadiation.cxx.

                                                                               {
  aParticleChange.Initialize(aTrack);
 
  // Obtain kinetic energy
  const G4DynamicParticle *aParticle(aTrack.GetDynamicParticle());
  G4double KineticEnergy(aParticle->GetKineticEnergy());
 
  // Kinetic energy too low for TR generation for any kind of particle
  if ( KineticEnergy < m_EkinMin )
    return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
 
  // Particle properties ...
  G4double ParticleMass(aParticle->GetDefinition()->GetPDGMass());
  G4double Gamma(1.0 + KineticEnergy/ParticleMass);
 
  if( Gamma < m_GammaMin ) {
    ATH_MSG_VERBOSE ( "Leaving PostStepDoIt(): Gamma < " << m_GammaMin );
    return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
  }
 
  G4double FoilThickness(0.);
  G4double GasThickness(0.);
  BEflag BEflag(BARREL); //enum?
 
  // Check whether particle is inside radiator
  std::vector<TRTRadiatorParameters>::const_iterator currentRadiator(m_radiators.begin());
  const std::vector<TRTRadiatorParameters>::const_iterator endOfRadiators(m_radiators.end());
  while(currentRadiator!=endOfRadiators) {
    if ( currentRadiator->GetLogicalVolume() ==
         aTrack.GetVolume()->GetLogicalVolume() ) {
      FoilThickness = currentRadiator->GetFoilThickness();
      GasThickness  = currentRadiator->GetGasThickness();
      BEflag        = currentRadiator->GetBEflag();
      break;
    }
    ++currentRadiator;
  }
  // Particle not inside radiator - return
  if ( currentRadiator == endOfRadiators ) {
    return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
  }
 
  ATH_MSG_VERBOSE ( "PostStepDoIt(): Inside radiator "
                    << currentRadiator->GetLogicalVolume()->GetName() );
 
  G4ThreeVector ParticleDirection(aParticle->GetMomentumDirection());
 
  if ( BEflag == ENDCAP ) {
    G4double costh(std::abs(ParticleDirection[2]));
    FoilThickness /= costh;
    GasThickness  /= costh;
  }
 
  G4double StepLength(aStep.GetStepLength());
  G4int FoilsTraversed(static_cast<G4int>(StepLength/(FoilThickness+GasThickness)+0.5));
 
  // Steplength so short, that no foils are crossed.
  // Forget it...
  // Actually the treatment we do is only valid for "many" foils.
  // Should we put this requirement higher?
  if ( FoilsTraversed == 0 ) {
    ATH_MSG_VERBOSE ( "Leaving PostStepDoIt(): No foils traversed." );
    return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
  }
 
  G4double Gy(0.), Gmany(0.), Neff(0.);
 
  // Generate, in the language of Artru:
  //   1) G_y
  //   2) N_{eff}
  //   3) G_{many}
  // Divide by photon energy to get photon spectrum (m_Ey)
 
  for( G4int i(0); i<m_NumBins; ++i ) {
 
#ifdef ARTRU                                    // Formulas from Artru
    Gy = XEmitanArtru( m_om[i],
                       Gamma,
                       GasThickness,
                       FoilThickness );
    Neff = NeffArtru( m_sigmaGas[i],
                      m_sigmaFoil[i],
                      GasThickness,
                      FoilThickness,
                      FoilsTraversed );
#else
    Gy = XEmitanNevski( m_om[i],                   // Formulas from Cherry via
                        Gamma,                     // Pavel Nevski in G3
                        GasThickness,
                        FoilThickness );
    Neff = NeffNevski( m_sigmaGas[i],
                       m_sigmaFoil[i],
                       GasThickness,
                       FoilThickness,
                       FoilsTraversed );
#endif
 
    Gmany = Neff*Gy;                              // Energy spectrum
    m_Ey[i] = Gmany / m_om[i];                    // Photon spectrum
 
  }
 
  // Integrate spectrum
  G4double Sph(XInteg( m_Ey, m_Sr ));
  G4double Sph_org(Sph);
 
  if ( ( BEflag == BARREL ) && (Gamma > 0.0) ) {
    G4double tempLog10Gamma(log10(Gamma)); //Hardcoded numbers from TestBeam study (E.Klinkby). Better to move elsewhere.
    G4double pHT_dEdxData = -0.001 + 0.0005*tempLog10Gamma;
    G4double pHT_TRData   = 0.15/(1.0 + exp(-(tempLog10Gamma - 3.3)/0.27));
    G4double pHT_dEdxMC = -0.0031 + 0.00016*tempLog10Gamma;
    G4double pHT_TRMC   = 0.1289/(1.0 + exp(-(tempLog10Gamma - 3.01)/0.1253));
 
    G4double fudge = std::min(2.0,(pHT_dEdxData+pHT_TRData)/( pHT_dEdxMC+pHT_TRMC));
 
    Sph *= fudge;
  }
 
  // Now how many photons are generated?
  G4long NumberOfPhotonsGenerated(CLHEP::RandPoisson::shoot(Sph));
 
  //reset Sph to original value.
  Sph = Sph_org;
 
  ATH_MSG_VERBOSE ( ">>>> FoilsTraversed = " << FoilsTraversed
                    << ", Gamma = " << Gamma
                    << ", Sph = " << Sph
                    << ", Nph = " << NumberOfPhotonsGenerated );
  // No photons, return
  if ( NumberOfPhotonsGenerated <= 0 ) {
    ATH_MSG_VERBOSE ( "Leaving TRTTransitionRadiation::PostStepDoIt: "
                      << "No photons generated." );
 
    return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
  }
 
  // Now, generate the photons
  aParticleChange.SetNumberOfSecondaries( NumberOfPhotonsGenerated );
  G4StepPoint *pPostStepPoint(aStep.GetPostStepPoint());
  G4ThreeVector x0(pPostStepPoint->GetPosition());
 
  G4double Esum(0.), XRannu(0.), Ephot(0.), Theta(0.), Phi(0.),
    X(0.), Y(0.), Z(0.);
 
  G4ThreeVector pos(0),TRPhotonDirection(0);
 
  for( G4int i(0); i<NumberOfPhotonsGenerated; ++i ) {
    XRannu = CLHEP::RandFlat::shoot();
    Ephot = exp( XFinter( XRannu*Sph, m_Sr, m_Omg ) );
 
    Esum += Ephot;
 
    pos  = pPostStepPoint->GetPosition();
 
    ATH_MSG_VERBOSE ( " Photon generated with energy : " << Ephot/CLHEP::keV
                      << " keV; position: ( " << pos.x()/CLHEP::cm << " cm, "
                      << pos.y()/CLHEP::cm << " cm, " << pos.z()/CLHEP::cm << " cm )" );
 
    // Angle w.r.t. electron direction (this is not correct but anyway...)
    Theta = std::abs( CLHEP::RandGaussZiggurat::shoot( 0.0, M_PI/Gamma ) );
    if( Theta >= 0.1 ) Theta = 0.1;
 
    Phi = (2.*M_PI)*CLHEP::RandFlat::shoot();
 
    X = sin(Theta)*cos(Phi);
    Y = sin(Theta)*sin(Phi);
    Z = cos(Theta);
 
    TRPhotonDirection.set(X,Y,Z);
    TRPhotonDirection.rotateUz( ParticleDirection );
    TRPhotonDirection.unit();
 
    G4DynamicParticle*
      TRPhoton = new G4DynamicParticle( G4Gamma::Gamma(),
                                        TRPhotonDirection,
                                        Ephot );
    aParticleChange.AddSecondary( TRPhoton );
 
  }
 
  aParticleChange.ProposeEnergy( KineticEnergy-Esum );
 
  ATH_MSG_VERBOSE ( "Leaving TRTTransitionRadiation::PostStepDoIt; "
                    << "Number of TR photons generated: " << NumberOfPhotonsGenerated );
 
  return &aParticleChange;
}

setLevel()

void AthMessaging::setLevel ( MSG::Level lvl )

inherited

Change the current logging level.

Use this rather than msg().setLevel() for proper operation with MT.

Definition at line 28 of file AthMessaging.cxx.

{
  m_lvl = lvl;
}

XEmitanArtru()

G4double TRTTransitionRadiation::XEmitanArtru ( const G4double & PhotonEnergy, const G4double & Gamma, const G4double & GasThickness, const G4double & FoilThickness ) const

inlineprivate

Definition at line 543 of file TRTTransitionRadiation.cxx.

                                                                                  {
  G4double phimax = 4.*M_PI;
  G4double ginv2 = 1./(Gamma*Gamma);
  G4double tau = d2/d1;
 
  G4double mom = omega/CLHEP::hbarc;
 
  G4double xi1 = m_WplasmaFoil/omega;
  G4double xi2 = m_WplasmaGas/omega;
 
  G4double xi1sq = xi1*xi1;
  G4double xi2sq = xi2*xi2;
 
  G4double Z1 = 2./(mom*(ginv2+xi1sq));
  G4double Z2 = 2./(mom*(ginv2+xi2sq));
 
  G4double a = d1/Z2;
  G4double V = d1*(1./Z1-1./Z2);
 
  G4double p0   = mom*d1*(xi1sq-xi2sq)/(4.*M_PI);
  G4double pmin = p0 + a*(1.+tau)/(2.*M_PI);
  G4double ymax = tau*phimax;
  G4double pmax = p0 + ymax*(1.+tau)/(2.*M_PI);
 
  G4int ipmin = (int)pmin+1;
  G4int ipmax = (int)pmax;
 
  G4double sum=0.;
  G4double y0 = -0.5*mom*d1*(xi1sq-xi2sq)/(1+tau);
  G4double dy = (2.*M_PI)/(1+tau);
  G4double y;
 
  G4double term = 0.;
  G4double t1   = 0.;
  G4double t2   = 0.;
  G4double tx   = 0.;
 
  for ( G4int ip(ipmin); ip<=ipmax; ++ip ) {
    t2 = t1;
    t1 = term;
 
    y = y0 + ip*dy;
    term  = (1./y-1./(y+V))*sin(0.5*(y+V));
    term *= term;
    term *= y-a;
    sum  += term;
 
    //G4double theta2 = 2.*y*tau/(mom*d2) - ginv2 - xi2sq;
    //G4double theta  = sqrt(theta2);
 
    tx = (tx>term) ? tx : term;
 
    // Locate local maxima (previous term)
 
    if ( t2 > 0. && t1 > term && t1 > t2 ) {
      if ( t1 < 0.0001*sum ) break;
    }
  }
 
  return ( sum*4.*2.*CLHEP::fine_structure_const/(1.+tau));
}

XEmitanNevski()

G4double TRTTransitionRadiation::XEmitanNevski ( const G4double & PhotonEnergy, const G4double & Gamma, const G4double & GasThickness, const G4double & FoilThickness ) const

inlineprivate

Definition at line 487 of file TRTTransitionRadiation.cxx.

                                                                                    {
  //FIXME hardcoded numbers are not good.
  G4double eps = 0.05;
 
  G4double Al = Al1 + Al2;
 
  G4double CK = 1/(2*M_PI*CLHEP::hbarc);
 
  // Plasma frequencies:
  G4double om1 = m_WplasmaGas;
  G4double om2 = m_WplasmaFoil;
 
  G4double om12 = (om1*om1-om2*om2)/PhotonEnergy;
  G4double oms = (Al1*om1*om1+Al2*om2*om2)/Al;
  G4double Ao = 0.5*Al1*om12*CK;
  G4double Bo = 0.5*Al2*om12*CK;
  G4double ro = sqrt(oms)/Gamma;
  G4double co = 0.5*(PhotonEnergy/(Gamma*Gamma)+oms/PhotonEnergy);
  G4double cc = co*Al*CK;
  G4double cb = M_PI*Al2/Al;
  G4int i1 = static_cast<G4int>(1.0+(ro>co?ro:co)*Al*CK);//(G4int)(1.0+max(ro,co)*Al*CK);
  G4int i2 = static_cast<G4int>(Gamma*ro*Al*CK);
 
  G4double Sum = 0.0;
  G4double Smx = 0.0;
  G4double Sm1 = 0.0;
  G4double S = 0.0;
  G4double So = 0.0;
 
  G4double Sx = 0.0;
 
  for(G4int i(i1); i<=i2; ++i) {
    Sx = So;
    So = S;
    S = (static_cast<G4double>(i)-cc)*(sin(cb*(static_cast<G4double>(i)-Ao))/((static_cast<G4double>(i)-Ao)*(static_cast<G4double>(i)+Bo)))*
      (sin(cb*(static_cast<G4double>(i)-Ao))/((static_cast<G4double>(i)-Ao)*(static_cast<G4double>(i)+Bo)));
    Smx=Smx>S?Smx:S;//max(Smx,S);
    Sm1=Sm1>S?Sm1:S;//max(Sm1,S);
    Sum += S;
 
    if( (So<Sx)&&(So<=S) ) {
      if( Sm1<eps*Smx ) break;
      Sm1 = 0.0;
    }
  }
 
  return ( CLHEP::fine_structure_const*om12*om12*Al*Al*CK*CK*Sum/M_PI );
}

XFinter()

G4double TRTTransitionRadiation::XFinter ( const G4double & X, const G4double * A, const G4double * F ) const

private

Definition at line 630 of file TRTTransitionRadiation.cxx.

                                                                                                         {
  G4int K1(0);
  G4int K2(0);
 
  if( (K1>=m_NumBins) || ( X < A[K1] ) ) K1 = 0; //FIXME This line has no effect?
  if( (K2>=m_NumBins) || ( X > A[K2] ) ) K2 = m_NumBins-1;
 
  G4int K(0);
  while( K2-K1>1 ) {
    K = static_cast<G4int>((K1+K2)/2);
    if( A[K]<X ) {
      K1 = K;
    } else {
      K2 = K;
    }
  }
 
  G4double X1(A[K1]);
  G4double X2(A[K1+1]);
 
  G4double xfinter((F[K1]*(X-X2)+F[K1+1]*(X1-X))/(X1-X2));
  return xfinter;
}

XInteg()

G4double TRTTransitionRadiation::XInteg ( const G4double * yy, G4double * ss ) const

inlineprivate

Definition at line 614 of file TRTTransitionRadiation.cxx.

                                                                                         {
  G4double S(0.);
  ss[0] = 0.;
 
  for( G4int i(1); i<m_NumBins; ++i ) {
    S += 0.5*( m_Omg[i]-m_Omg[i-1] )*( yy[i-1]*m_om[i-1] + yy[i]*m_om[i] );
    ss[i] = S;
  }
 
  return S;
}

Member Data Documentation

ATLAS_THREAD_SAFE

std::atomic_flag m_initialized AthMessaging::ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT

mutableprivateinherited

Messaging initialized (initMessaging)

Definition at line 141 of file AthMessaging.h.

m_EkinMin

G4double TRTTransitionRadiation::m_EkinMin

private

Definition at line 91 of file TRTTransitionRadiation.h.

m_Ey

G4double* TRTTransitionRadiation::m_Ey

private

Definition at line 93 of file TRTTransitionRadiation.h.

m_GammaMin

G4double TRTTransitionRadiation::m_GammaMin

private

Definition at line 90 of file TRTTransitionRadiation.h.

m_imsg

std::atomic<IMessageSvc*> AthMessaging::m_imsg { nullptr }

mutableprivateinherited

MessageSvc pointer.

Definition at line 135 of file AthMessaging.h.

{ nullptr };

m_lvl

std::atomic<MSG::Level> AthMessaging::m_lvl { MSG::NIL }

mutableprivateinherited

Current logging level.

Definition at line 138 of file AthMessaging.h.

{ MSG::NIL };

m_MaxEnergyTR

G4double TRTTransitionRadiation::m_MaxEnergyTR

private

Definition at line 84 of file TRTTransitionRadiation.h.

m_MinEnergyTR

G4double TRTTransitionRadiation::m_MinEnergyTR

private

Definition at line 83 of file TRTTransitionRadiation.h.

m_msg_tls

boost::thread_specific_ptr<MsgStream> AthMessaging::m_msg_tls

mutableprivateinherited

MsgStream instance (a std::cout like with print-out levels)

Definition at line 132 of file AthMessaging.h.

m_nm

std::string AthMessaging::m_nm

privateinherited

Message source name.

Definition at line 129 of file AthMessaging.h.

m_NumBins

G4int TRTTransitionRadiation::m_NumBins

private

Definition at line 85 of file TRTTransitionRadiation.h.

m_om

G4double* TRTTransitionRadiation::m_om

private

Definition at line 95 of file TRTTransitionRadiation.h.

m_Omg

G4double* TRTTransitionRadiation::m_Omg

private

Definition at line 96 of file TRTTransitionRadiation.h.

m_radiators

std::vector<TRTRadiatorParameters> TRTTransitionRadiation::m_radiators

private

Definition at line 81 of file TRTTransitionRadiation.h.

m_sigmaFoil

G4double* TRTTransitionRadiation::m_sigmaFoil

private

Definition at line 98 of file TRTTransitionRadiation.h.

m_sigmaGas

G4double* TRTTransitionRadiation::m_sigmaGas

private

Definition at line 97 of file TRTTransitionRadiation.h.

m_Sr

G4double* TRTTransitionRadiation::m_Sr

private

Definition at line 94 of file TRTTransitionRadiation.h.

m_WplasmaFoil

G4double TRTTransitionRadiation::m_WplasmaFoil

private

Definition at line 88 of file TRTTransitionRadiation.h.

m_WplasmaGas

G4double TRTTransitionRadiation::m_WplasmaGas

private

Definition at line 87 of file TRTTransitionRadiation.h.

m_xmlfilename

const std::string TRTTransitionRadiation::m_xmlfilename

private

Definition at line 80 of file TRTTransitionRadiation.h.

m_XMLhandler

TRRegionXMLHandler* TRTTransitionRadiation::m_XMLhandler

private

Definition at line 79 of file TRTTransitionRadiation.h.

---

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

• TRTTransitionRadiation.h
• TRTTransitionRadiation.cxx