#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 (TRTRadiatorParameters p)

bool	msgLvl (const MSG::Level lvl) const
	Test the output level. More...

MsgStream &	msg () const
	The standard message stream. More...

MsgStream &	msg (const MSG::Level lvl) const
	The standard message stream. More...

void	setLevel (MSG::Level lvl)
	Change the current logging level. More...

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. More...

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. More...

boost::thread_specific_ptr< MsgStream >	m_msg_tls
	MsgStream instance (a std::cout like with print-out levels) More...

std::atomic< IMessageSvc * >	m_imsg { nullptr }
	MessageSvc pointer. More...

std::atomic< MSG::Level >	m_lvl { MSG::NIL }
	Current logging level. More...

std::atomic_flag m_initialized	ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
	Messaging initialized (initMessaging) More...

Detailed Description

Definition at line 24 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 ( 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 657 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 245 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);
         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();
   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 235 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 164 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 179 of file AthMessaging.h.

180 { 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_initialized.test_and_set()) initMessaging();
   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 467 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 450 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 261 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 542 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 486 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 629 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 613 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 95 of file TRTTransitionRadiation.h.

◆ m_Ey

G4double* TRTTransitionRadiation::m_Ey

private

Definition at line 97 of file TRTTransitionRadiation.h.

◆ m_GammaMin

G4double TRTTransitionRadiation::m_GammaMin

private

Definition at line 94 of file TRTTransitionRadiation.h.

◆ m_imsg

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

mutableprivateinherited

MessageSvc pointer.

Definition at line 135 of file AthMessaging.h.

◆ m_lvl

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

mutableprivateinherited

Current logging level.

Definition at line 138 of file AthMessaging.h.

◆ m_MaxEnergyTR

G4double TRTTransitionRadiation::m_MaxEnergyTR

private

Definition at line 88 of file TRTTransitionRadiation.h.

◆ m_MinEnergyTR

G4double TRTTransitionRadiation::m_MinEnergyTR

private

Definition at line 87 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 89 of file TRTTransitionRadiation.h.

◆ m_om

G4double* TRTTransitionRadiation::m_om

private

Definition at line 99 of file TRTTransitionRadiation.h.

◆ m_Omg

G4double* TRTTransitionRadiation::m_Omg

private

Definition at line 100 of file TRTTransitionRadiation.h.

◆ m_radiators

std::vector<TRTRadiatorParameters> TRTTransitionRadiation::m_radiators

private

Definition at line 85 of file TRTTransitionRadiation.h.

◆ m_sigmaFoil

G4double* TRTTransitionRadiation::m_sigmaFoil

private

Definition at line 102 of file TRTTransitionRadiation.h.

◆ m_sigmaGas

G4double* TRTTransitionRadiation::m_sigmaGas

private

Definition at line 101 of file TRTTransitionRadiation.h.

◆ m_Sr

G4double* TRTTransitionRadiation::m_Sr

private

Definition at line 98 of file TRTTransitionRadiation.h.

◆ m_WplasmaFoil

G4double TRTTransitionRadiation::m_WplasmaFoil

private

Definition at line 92 of file TRTTransitionRadiation.h.

◆ m_WplasmaGas

G4double TRTTransitionRadiation::m_WplasmaGas

private

Definition at line 91 of file TRTTransitionRadiation.h.

◆ m_xmlfilename

const std::string TRTTransitionRadiation::m_xmlfilename

private

Definition at line 84 of file TRTTransitionRadiation.h.

◆ m_XMLhandler

TRRegionXMLHandler* TRTTransitionRadiation::m_XMLhandler

private

Definition at line 83 of file TRTTransitionRadiation.h.

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

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ TRTTransitionRadiation() [1/2]

◆ ~TRTTransitionRadiation()

◆ TRTTransitionRadiation() [2/2]

Member Function Documentation

◆ AddRadiatorParameters()

◆ ComputePhotoAbsCof()

◆ GetMeanFreePath()

◆ Initialize()

◆ initMessaging()

◆ IsApplicable()

◆ msg() [1/2]

◆ msg() [2/2]

◆ msgLvl()

◆ NeffArtru()

◆ NeffNevski()

◆ operator=()

◆ PostStepDoIt()

◆ setLevel()

◆ XEmitanArtru()

◆ XEmitanNevski()

◆ XFinter()

◆ XInteg()

Member Data Documentation

◆ ATLAS_THREAD_SAFE

◆ m_EkinMin

◆ m_Ey

◆ m_GammaMin

◆ m_imsg

◆ m_lvl

◆ m_MaxEnergyTR

◆ m_MinEnergyTR

◆ m_msg_tls

◆ m_nm

◆ m_NumBins

◆ m_om

◆ m_Omg

◆ m_radiators

◆ m_sigmaFoil

◆ m_sigmaGas

◆ m_Sr

◆ m_WplasmaFoil

◆ m_WplasmaGas

◆ m_xmlfilename

◆ m_XMLhandler