Skip to content

Commit 103bb7e

Browse files
MarkusFrankATcernchMarkusFrankATcernch
authored
Fix computation of the radiation length and nuclear interaction lengt… (#11060)
* Fix computation of the radiation length and nuclear interaction length when ROOT uses G4 units * Fix computation of the radiation length and nuclear interaction length when ROOT uses G4 units
1 parent 7ba7de8 commit 103bb7e

File tree

1 file changed

+35
-53
lines changed

1 file changed

+35
-53
lines changed

‎geom/geom/src/TGeoMaterial.cxx‎

Lines changed: 35 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -488,25 +488,18 @@ void TGeoMaterial::SetRadLen(Double_t radlen, Double_t intlen)
488488
}
489489
TGeoManager::EDefaultUnits typ = TGeoManager::GetDefaultUnits();
490490
// compute radlen systematically with G3 formula for a valid material
491-
if ( typ == TGeoManager::kRootUnits && radlen>=0 ) {
491+
if (radlen >= 0) {
492492
//taken grom Geant3 routine GSMATE
493-
constexpr Double_t alr2av = 1.39621E-03*TGeoUnit::cm2;
493+
constexpr Double_t alr2av = 1.39621E-03;
494494
constexpr Double_t al183 = 5.20948;
495495
fRadLen = fA/(alr2av*fDensity*fZ*(fZ +TGeoMaterial::ScreenFactor(fZ))*
496496
(al183-TMath::Log(fZ)/3-TGeoMaterial::Coulomb(fZ)));
497-
fRadLen *= TGeoUnit::cm;
498-
}
499-
else if ( typ == TGeoManager::kG4Units && radlen>=0 ) {
500-
//taken grom Geant3 routine GSMATE
501-
constexpr Double_t alr2av = 1.39621E-03*TGeant4Unit::cm2;
502-
constexpr Double_t al183 = 5.20948;
503-
fRadLen = fA/(alr2av*fDensity*fZ*(fZ +TGeoMaterial::ScreenFactor(fZ))*
504-
(al183-TMath::Log(fZ)/3-TGeoMaterial::Coulomb(fZ)));
505-
fRadLen *= TGeant4Unit::cm;
497+
// fRadLen is in TGeo units. Apply conversion factor in requested length-units
498+
fRadLen *= (typ == TGeoManager::kRootUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
506499
}
507500
// Compute interaction length using the same formula as in GEANT4
508-
if ( typ == TGeoManager::kRootUnits && intlen>=0 ) {
509-
constexpr Double_t lambda0 = 35.*TGeoUnit::g/TGeoUnit::cm2; // [g/cm^2]
501+
if (intlen >= 0) {
502+
constexpr Double_t lambda0 = 35. * TGeoUnit::g / TGeoUnit::cm2; // [g/cm^2]
510503
Double_t nilinv = 0.0;
511504
TGeoElement *elem = GetElement();
512505
if (!elem) {
@@ -515,21 +508,10 @@ void TGeoMaterial::SetRadLen(Double_t radlen, Double_t intlen)
515508
}
516509
Double_t nbAtomsPerVolume = TGeoUnit::Avogadro*fDensity/elem->A();
517510
nilinv += nbAtomsPerVolume*TMath::Power(elem->Neff(), 0.6666667);
518-
nilinv *= TGeoUnit::amu/lambda0;
519-
fIntLen = (nilinv<=0) ? TGeoShape::Big() : (TGeoUnit::cm/nilinv);
520-
}
521-
else if ( typ == TGeoManager::kG4Units && intlen>=0 ) {
522-
constexpr Double_t lambda0 = 35.*TGeant4Unit::g/TGeant4Unit::cm2; // [g/cm^2]
523-
Double_t nilinv = 0.0;
524-
TGeoElement *elem = GetElement();
525-
if (!elem) {
526-
Fatal("SetRadLen", "Element not found for material %s", GetName());
527-
return;
528-
}
529-
Double_t nbAtomsPerVolume = TGeant4Unit::Avogadro*fDensity/elem->A();
530-
nilinv += nbAtomsPerVolume*TMath::Power(elem->Neff(), 0.6666667);
531-
nilinv *= TGeant4Unit::amu/lambda0;
532-
fIntLen = (nilinv<=0) ? TGeoShape::Big() : (TGeant4Unit::cm/nilinv);
511+
nilinv *= TGeoUnit::amu / lambda0;
512+
fIntLen = (nilinv <= 0) ? TGeoShape::Big() : (1.0 / nilinv);
513+
// fIntLen is in TGeo units. Apply conversion factor in requested length-units
514+
fIntLen *= (typ == TGeoManager::kRootUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
533515
}
534516
}
535517

@@ -800,15 +782,10 @@ TGeoMixture::~TGeoMixture()
800782

801783
void TGeoMixture::AverageProperties()
802784
{
803-
TGeoManager::EDefaultUnits typ = TGeoManager::GetDefaultUnits();
804-
const Double_t cm = (typ==TGeoManager::kRootUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
805-
const Double_t cm2 = (typ==TGeoManager::kRootUnits) ? TGeoUnit::cm2 : TGeant4Unit::cm2;
806-
const Double_t amu = (typ==TGeoManager::kRootUnits) ? TGeoUnit::amu : TGeant4Unit::amu; // [MeV/c^2]
807-
const Double_t gram = (typ==TGeoManager::kRootUnits) ? TGeoUnit::gram : TGeant4Unit::gram;
808-
const Double_t na = (typ==TGeoManager::kRootUnits) ? TGeoUnit::Avogadro : TGeant4Unit::Avogadro;
809-
const Double_t alr2av = 1.39621E-03 * cm2;
810-
const Double_t al183 = 5.20948;
811-
const Double_t lambda0 = 35.*gram/cm2; // [g/cm^2]
785+
constexpr const Double_t na = TGeoUnit::Avogadro;
786+
constexpr const Double_t alr2av = 1.39621E-03;
787+
constexpr const Double_t al183 = 5.20948;
788+
constexpr const Double_t lambda0 = 35. * TGeoUnit::g / TGeoUnit::cm2; // [g/cm^2]
812789
Double_t radinv = 0.0;
813790
Double_t nilinv = 0.0;
814791
Double_t nbAtomsPerVolume;
@@ -827,10 +804,15 @@ void TGeoMixture::AverageProperties()
827804
radinv += xinv*fWeights[j];
828805
}
829806
radinv *= alr2av*fDensity;
830-
if (radinv > 0) fRadLen = cm/radinv;
807+
fRadLen = (radinv <= 0) ? TGeoShape::Big() : 1.0 / radinv;
808+
// fRadLen is in TGeo units. Apply conversion factor in requested length-units
809+
fRadLen *= (TGeoManager::GetDefaultUnits() == TGeoManager::kRootUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
810+
831811
// Compute interaction length
832-
nilinv *= amu/lambda0;
833-
fIntLen = (nilinv<=0) ? TGeoShape::Big() : (cm/nilinv);
812+
nilinv *= TGeoUnit::amu / lambda0;
813+
fIntLen = (nilinv <= 0) ? TGeoShape::Big() : 1.0 / nilinv;
814+
// fIntLen is in TGeo units. Apply conversion factor in requested length-units
815+
fIntLen *= (TGeoManager::GetDefaultUnits() == TGeoManager::kRootUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
834816
}
835817

836818
////////////////////////////////////////////////////////////////////////////////
@@ -1312,35 +1294,35 @@ void TGeoMixture::ComputeDerivedQuantities()
13121294
void TGeoMixture::ComputeRadiationLength()
13131295
{
13141296
// Formula taken from G4Material.cxx L556
1315-
const Double_t cm = (TGeoManager::GetDefaultUnits()==TGeoManager::kRootUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
13161297
Double_t radinv = 0.0 ;
13171298
for (Int_t i=0;i<fNelements;++i) {
1318-
radinv += fVecNbOfAtomsPerVolume[i]*((TGeoElement*)fElements->At(i))->GetfRadTsai();
1299+
// GetfRadTsai is in units of cm2
1300+
radinv += fVecNbOfAtomsPerVolume[i] * ((TGeoElement *)fElements->At(i))->GetfRadTsai() / TGeoUnit::cm2;
13191301
}
1320-
fRadLen = (radinv <= 0.0 ? DBL_MAX : cm/radinv);
1302+
fRadLen = (radinv <= 0.0 ? DBL_MAX : 1.0 / radinv);
1303+
// fRadLen is in TGeo units. Apply conversion factor in requested length-units
1304+
fRadLen *= (TGeoManager::GetDefaultUnits() == TGeoManager::kRootUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
13211305
}
13221306

13231307
////////////////////////////////////////////////////////////////////////////////
13241308
/// Compute Nuclear Interaction Length based on Geant4 formula
13251309
void TGeoMixture::ComputeNuclearInterLength()
13261310
{
13271311
// Formula taken from G4Material.cxx L567
1328-
TGeoManager::EDefaultUnits typ = TGeoManager::GetDefaultUnits();
1329-
const Double_t g = (typ==TGeoManager::kRootUnits) ? TGeoUnit::g : TGeant4Unit::g;
1330-
const Double_t cm = (typ==TGeoManager::kRootUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
1331-
const Double_t amu = (typ==TGeoManager::kRootUnits) ? TGeoUnit::amu : TGeant4Unit::amu;
1332-
const Double_t lambda0 = 35*g/(cm*cm);
1312+
constexpr Double_t lambda0 = 35. * TGeoUnit::g / TGeoUnit::cm2; // [g/cm^2]
13331313
const Double_t twothird = 2.0/3.0;
13341314
Double_t NILinv = 0.0;
13351315
for (Int_t i=0; i<fNelements; ++i) {
1336-
Int_t Z = static_cast<Int_t>(((TGeoElement*)fElements->At(i))->Z()+0.5);
1316+
Int_t Z = static_cast<Int_t>(((TGeoElement *)fElements->At(i))->Z() + 0.5);
13371317
Double_t A = ((TGeoElement*)fElements->At(i))->Neff();
13381318
if(1 == Z) {
1339-
NILinv += fVecNbOfAtomsPerVolume[i]*A;
1319+
NILinv += fVecNbOfAtomsPerVolume[i] * A;
13401320
} else {
1341-
NILinv += fVecNbOfAtomsPerVolume[i]*TMath::Exp(twothird*TMath::Log(A));
1321+
NILinv += fVecNbOfAtomsPerVolume[i] * TMath::Exp(twothird * TMath::Log(A));
13421322
}
13431323
}
1344-
NILinv *= amu/lambda0;
1345-
fIntLen = (NILinv <= 0.0 ? DBL_MAX : cm/NILinv);
1324+
NILinv *= TGeoUnit::amu / lambda0;
1325+
fIntLen = (NILinv <= 0.0 ? DBL_MAX : 1.0 / NILinv);
1326+
// fIntLen is in TGeo units. Apply conversion factor in requested length-units
1327+
fIntLen *= (TGeoManager::GetDefaultUnits() == TGeoManager::kRootUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
13461328
}

0 commit comments

Comments
 (0)