15 #include "CLHEP/Random/RandGaussZiggurat.h"
28 #include "G4DynamicParticle.hh"
30 #include "G4Material.hh"
31 #include "G4Element.hh"
32 #include "G4EmProcessSubType.hh"
33 #include "G4ProcessType.hh"
35 #include "Randomize.hh"
38 #include "CLHEP/Random/RandPoisson.h"
49 G4VDiscreteProcess(processName,fElectromagnetic),
51 m_XMLhandler(nullptr),m_xmlfilename(xmlfilename),
52 m_MinEnergyTR(0.0),m_MaxEnergyTR(0.0),m_NumBins(0),m_WplasmaGas(0.0),
53 m_WplasmaFoil(0.0),m_GammaMin(0.0),m_EkinMin(0.0),m_Ey(nullptr),m_Sr(nullptr),
54 m_om(nullptr),m_Omg(nullptr),m_sigmaGas(nullptr),m_sigmaFoil(nullptr)
59 SetProcessSubType( fTransitionRadiation );
65 ATH_MSG_DEBUG (
"##### Constructor TRTTransitionRadiation done" );
83 ATH_MSG_DEBUG(
" New Radiator parameters being defined for TR process");
95 ATH_MSG_DEBUG (
"Constructor TRTTransitionRadiation::Initialize()" );
110 G4Material* g4mat_Gas{
nullptr};
111 G4Material* g4mat_FoilMaterial{
nullptr};
112 std::string errorMessage;
115 ISvcLocator *svcLocator = Gaudi::svcLocator();
117 if( StatusCode::SUCCESS != svcLocator->service(
"DetectorStore",
detStore ) ) {
118 errorMessage =
"Can not access Detector Store";
120 throw std::runtime_error(errorMessage);
123 if(materialManager) {
125 g4mat_Gas = geo2g4_material_fact.
Build(materialManager->
getMaterial(
"trt::CO2"));
126 g4mat_FoilMaterial = geo2g4_material_fact.
Build(materialManager->
getMaterial(
"trt::CH2"));
133 ATH_MSG_INFO(
"GeoModel Material is not available. Construct TR materials using IRDBAccessSvc interface");
135 if(StatusCode::SUCCESS != svcLocator->service(
"GeoDbTagSvc",geoDbTagSvc)) {
136 errorMessage =
"Can not access GeoDbTagSvc";
138 throw std::runtime_error(errorMessage);
143 if(StatusCode::SUCCESS != svcLocator->service(geoDbTagSvc->getParamSvcName(),pAccessSvc)) {
144 errorMessage =
"Can not access " + geoDbTagSvc->getParamSvcName();
146 throw std::runtime_error(errorMessage);
149 std::map<std::string,int> materialComponentsMap;
150 std::map<std::string,G4Material*> materialMap;
151 IRDBRecordset_ptr materialsRec = pAccessSvc->getRecordsetPtr(
"TRMaterials",
"",
"");
154 std::string
key = material->getString(
"NAME");
155 auto mapIt = materialComponentsMap.find(
key);
156 if(mapIt == materialComponentsMap.end()) {
157 materialComponentsMap.emplace(
key,1);
160 materialComponentsMap.at(
key) = mapIt->second + 1;
165 std::string
key = material->getString(
"NAME");
166 G4Material* g4Material{
nullptr};
167 if(materialMap.find(
key)==materialMap.end()) {
168 auto itNElements = materialComponentsMap.find(
key);
169 g4Material =
new G4Material(
key
171 ,itNElements->second);
172 materialMap.emplace(
key,g4Material);
175 g4Material = materialMap.at(
key);
178 G4Element* g4Element = G4Element::GetElement(material->getString(
"ELEMENT_NAME"));
180 errorMessage =
"Wrong name for element " + material->getString(
"ELEMENT_NAME")
181 +
" found in the description of material " +
key;
183 throw std::runtime_error(errorMessage);
185 g4Material->AddElement(g4Element,(G4int)material->getInt(
"N"));
187 g4mat_Gas = materialMap.at(
"CO2");
188 g4mat_FoilMaterial = materialMap.at(
"CH2");
194 m_WplasmaFoil = sqrt( PlasmaCof * g4mat_FoilMaterial->GetElectronDensity() );
195 m_WplasmaGas = sqrt( PlasmaCof * g4mat_Gas->GetElectronDensity() );
197 ATH_MSG_DEBUG(
"Foil material : " << g4mat_FoilMaterial->GetName()
224 ATH_MSG_DEBUG (
"Constructor TRTTransitionRadiation::Initialize() DONE!" );
236 return ( std::abs(
particle.GetPDGCharge())>std::numeric_limits<double>::epsilon() && std::abs(
particle.GetPDGMass())>std::numeric_limits<double>::epsilon() ) ;
247 const G4DynamicParticle *aParticle(aTrack.GetDynamicParticle());
248 if(std::abs(aParticle->GetDefinition()->GetPDGCharge())< std::numeric_limits<double>::epsilon() ||
249 std::abs(aParticle->GetDefinition()->GetPDGMass())<std::numeric_limits<double>::epsilon() ) {
261 const G4Step& aStep ) {
262 aParticleChange.Initialize(aTrack);
265 const G4DynamicParticle *aParticle(aTrack.GetDynamicParticle());
266 G4double KineticEnergy(aParticle->GetKineticEnergy());
270 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
273 G4double
ParticleMass(aParticle->GetDefinition()->GetPDGMass());
278 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
281 G4double FoilThickness(0.);
282 G4double GasThickness(0.);
286 std::vector<TRTRadiatorParameters>::const_iterator currentRadiator(
m_radiators.begin());
287 const std::vector<TRTRadiatorParameters>::const_iterator endOfRadiators(
m_radiators.end());
288 while(currentRadiator!=endOfRadiators) {
289 if ( currentRadiator->GetLogicalVolume() ==
290 aTrack.GetVolume()->GetLogicalVolume() ) {
291 FoilThickness = currentRadiator->GetFoilThickness();
292 GasThickness = currentRadiator->GetGasThickness();
293 BEflag = currentRadiator->GetBEflag();
299 if ( currentRadiator == endOfRadiators ) {
300 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
304 << currentRadiator->GetLogicalVolume()->GetName() );
306 G4ThreeVector ParticleDirection(aParticle->GetMomentumDirection());
309 G4double costh(std::abs(ParticleDirection[2]));
310 FoilThickness /= costh;
311 GasThickness /= costh;
315 G4int FoilsTraversed(
static_cast<G4int
>(
StepLength/(FoilThickness+GasThickness)+0.5));
321 if ( FoilsTraversed == 0 ) {
323 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
326 G4double Gy(0.), Gmany(0.), Neff(0.);
336 #ifdef ARTRU // Formulas from Artru
365 G4double Sph_org(Sph);
368 G4double tempLog10Gamma(log10(
Gamma));
369 G4double pHT_dEdxData = -0.001 + 0.0005*tempLog10Gamma;
370 G4double pHT_TRData = 0.15/(1.0 +
exp(-(tempLog10Gamma - 3.3)/0.27));
371 G4double pHT_dEdxMC = -0.0031 + 0.00016*tempLog10Gamma;
372 G4double pHT_TRMC = 0.1289/(1.0 +
exp(-(tempLog10Gamma - 3.01)/0.1253));
374 G4double fudge =
std::min(2.0,(pHT_dEdxData+pHT_TRData)/( pHT_dEdxMC+pHT_TRMC));
380 G4long NumberOfPhotonsGenerated(CLHEP::RandPoisson::shoot(Sph));
386 <<
", Gamma = " <<
Gamma
388 <<
", Nph = " << NumberOfPhotonsGenerated );
390 if ( NumberOfPhotonsGenerated <= 0 ) {
392 <<
"No photons generated." );
394 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
398 aParticleChange.SetNumberOfSecondaries( NumberOfPhotonsGenerated );
399 G4StepPoint *pPostStepPoint(aStep.GetPostStepPoint());
400 G4ThreeVector x0(pPostStepPoint->GetPosition());
402 G4double Esum(0.), XRannu(0.), Ephot(0.), Theta(0.),
Phi(0.),
405 G4ThreeVector
pos(0),TRPhotonDirection(0);
407 for( G4int
i(0);
i<NumberOfPhotonsGenerated; ++
i ) {
408 XRannu = CLHEP::RandFlat::shoot();
413 pos = pPostStepPoint->GetPosition();
420 Theta = std::abs( CLHEP::RandGaussZiggurat::shoot( 0.0,
M_PI/
Gamma ) );
421 if( Theta >= 0.1 ) Theta = 0.1;
423 Phi = (2.*
M_PI)*CLHEP::RandFlat::shoot();
429 TRPhotonDirection.set(
X,
Y,
Z);
430 TRPhotonDirection.rotateUz( ParticleDirection );
431 TRPhotonDirection.unit();
437 aParticleChange.AddSecondary( TRPhoton );
441 aParticleChange.ProposeEnergy( KineticEnergy-Esum );
444 <<
"Number of TR photons generated: " << NumberOfPhotonsGenerated );
446 return &aParticleChange;
450 const G4double & sigFoil,
451 const G4double & GasThickness,
452 const G4double & FoilThickness,
453 const G4int & FoilsTraversed )
const {
454 if ( FoilsTraversed <= 0 )
457 G4double C1 = sigGas * GasThickness*FoilsTraversed;
458 G4double C2 = sigFoil*FoilThickness*FoilsTraversed;
460 G4double W1 =
exp(-C1);
461 G4double W2 =
exp(-C2);
463 return FoilsTraversed * (1-W1*W2)/(C1+C2);
467 const G4double & sigFoil,
468 const G4double & GasThickness,
469 const G4double & FoilThickness,
470 const G4int & FoilsTraversed )
const {
471 if ( FoilsTraversed <= 0 )
474 G4double
sigma = FoilThickness*sigFoil + GasThickness*sigGas;
475 G4double xNF = FoilsTraversed;
486 const G4double &
Gamma,
487 const G4double & Al1,
488 const G4double & Al2 )
const {
492 G4double Al = Al1 + Al2;
500 G4double om12 = (om1*om1-om2*om2)/PhotonEnergy;
501 G4double oms = (Al1*om1*om1+Al2*om2*om2)/Al;
502 G4double Ao = 0.5*Al1*om12*CK;
503 G4double Bo = 0.5*Al2*om12*CK;
504 G4double ro = sqrt(oms)/
Gamma;
505 G4double co = 0.5*(PhotonEnergy/(
Gamma*
Gamma)+oms/PhotonEnergy);
506 G4double
cc = co*Al*CK;
507 G4double cb =
M_PI*Al2/Al;
508 G4int i1 =
static_cast<G4int
>(1.0+(ro>co?ro:co)*Al*CK);
509 G4int i2 =
static_cast<G4int
>(
Gamma*ro*Al*CK);
519 for(G4int
i(i1);
i<=i2; ++
i) {
522 S = (
static_cast<G4double
>(
i)-
cc)*(
sin(cb*(
static_cast<G4double
>(
i)-Ao))/((
static_cast<G4double
>(
i)-Ao)*(
static_cast<G4double
>(
i)+Bo)))*
523 (
sin(cb*(
static_cast<G4double
>(
i)-Ao))/((
static_cast<G4double
>(
i)-Ao)*(
static_cast<G4double
>(
i)+Bo)));
528 if( (So<Sx)&&(So<=
S) ) {
529 if( Sm1<eps*Smx )
break;
542 const G4double &
Gamma,
544 const G4double &
d1 )
const {
545 G4double phimax = 4.*
M_PI;
547 G4double tau =
d2/
d1;
554 G4double xi1sq = xi1*xi1;
555 G4double xi2sq = xi2*xi2;
557 G4double Z1 = 2./(
mom*(ginv2+xi1sq));
558 G4double Z2 = 2./(
mom*(ginv2+xi2sq));
561 G4double V =
d1*(1./Z1-1./Z2);
563 G4double p0 =
mom*
d1*(xi1sq-xi2sq)/(4.*
M_PI);
564 G4double pmin = p0 +
a*(1.+tau)/(2.*
M_PI);
565 G4double
ymax = tau*phimax;
566 G4double pmax = p0 +
ymax*(1.+tau)/(2.*
M_PI);
568 G4int ipmin = (
int)pmin+1;
569 G4int ipmax = (
int)pmax;
572 G4double y0 = -0.5*
mom*
d1*(xi1sq-xi2sq)/(1+tau);
573 G4double
dy = (2.*
M_PI)/(1+tau);
581 for ( G4int
ip(ipmin);
ip<=ipmax; ++
ip ) {
586 term = (1./
y-1./(
y+V))*
sin(0.5*(
y+V));
594 tx = (
tx>term) ?
tx : term;
598 if (
t2 > 0. &&
t1 > term &&
t1 >
t2 ) {
599 if (
t1 < 0.0001*
sum )
break;
637 K =
static_cast<G4int
>((K1+K2)/2);
646 G4double X2(
A[K1+1]);
648 G4double xfinter((
F[K1]*(
X-X2)+
F[K1+1]*(X1-
X))/(X1-X2));
657 const G4double & PhotonEnergy )
const {
659 const G4double *SandiaCof(
Material->GetSandiaTable()->GetSandiaCofForMaterial(PhotonEnergy));
661 const G4double Energy2(PhotonEnergy*PhotonEnergy);
662 const G4double Energy3(PhotonEnergy*Energy2);
663 const G4double Energy4(Energy2*Energy2);
665 return ( SandiaCof[0]/PhotonEnergy +
666 SandiaCof[1]/Energy2 +
667 SandiaCof[2]/Energy3 +
668 SandiaCof[3]/Energy4 );