Material table in ROOT, radiation and nuclear interaction lengths.

[MarkusFrankATcernch]

Since this is an issue coming up over and over again, I want to create here an entry so that the discussion is not lost.

[MarkusFrankATcernch]

For a speciality of the TGeoMaterial's interface SetRadLen, ROOT was always computing the radiation length and the nuclear interaction length. The values specified in the XML description read in DD4hep was ALWAYS ignored. see: https://root.cern.ch/root/html534/TGeoMaterial.html#TGeoMaterial:SetRadLen

Hence: If there is a discrepancy between the values in Geant4 we have a problem, which must be sorted out.

[gaede]

I just found the cause of the problem with the interaction lengths: In Compact2Objects.cpp you call:

    // Update estimated density if not provided.
    if ( has_density )   {
      mix->SetDensity(dens_val);
    }

after the Material (TGeoMixture) had benn created correctly.

However this call has a nasty side effect (in TGeoMixture.h):

SetDensity(Double_t density) {fDensity = density; SetRadLen(0);}

and SetRadLen() is only implemented for TGeoMaterial and re-conmputes radLen and IntLen for the first Element.

[MarkusFrankATcernch]

This call to SetDensity should be OK. SetRadLen re-computes in this event both values according to the new density. We want the computed values and ignore the XML entries. Hence, this is the correct behaviour.

[MarkusFrankATcernch]

UUuuhhh. I see what you mean. This means that mixtures are entirely mis-treated and effectively reduced to their first compound. Is this correct? If yes, I have to propagate this bug to Andrei.

[gaede]

Yes, this is a bug (oversight) in ROOT. I believe the correct fix would be to implement TGeoMixture::SetDensity() as :

SetDensity(Double_t density) {fDensity = density;  AverageProperties() ;}

as AverageProperties() is the equivalent for SetRadLen(0.) for individual element materials.

On the other hand, I am not sure we need the SetDensity() call in the first place, as the density should always be given in the material definition ( we could throw exception if it isn't ...)

[agheata]

Right, calling TGeoMaterial::SetDensity() from a mixture was wrong. Fixed now in the ROOT master, the call will trigger recomputing the average properties as proposed by Frank. Please let me know if any unforeseen side effect.

[petricm]

There is still stuff to check, in case the material is gas you have:

    <material name="beam" state="gas">
        <P unit="pascal" value="6.25e-06"/>
        <MEE unit="eV" value="38.5760755714278"/>
        <D unit="g/cm3" value="1.7e-14"/>
        <fraction n="0.36264" ref="H"/>
        <fraction n="0.36276" ref="N"/>
        <fraction n="0.117748421296248" ref="C"/>
        <fraction n="0.156851578703752" ref="O"/>
    </material>

So what happens with values of preasure and MEE

[petricm]

It seems that the RadL and Nucl.Int.Length from ROOT and Geant4 are not the same.

I used the dump from ROOT where:

    <element Z="14" N="28" formula="SI" name="SI">
        <atom type="A" unit="g/mole" value="2.80855000e+01" />
    </element>

and defined a material

    <material formula="Si" name="Silicon" state="solid" >
        <D type="density" unit="g/cm3" value="2.329" />
        <composite n="1" ref="Si" />
    </material>

when you call materialScan on this you get RadL=9.3536 and Nucl.Int.Length=45.7729

when you cal the Geant converter and call /material/g4/printMaterial all this dispalys all the materials and you get

Material:  Silicon    density:  2.329 g/cm3   RadL:   9.370 cm   Nucl.Int.Length:  45.680 cm 
                       Imean: 173.000 eV 

   --->  Element: SILICON (SI)   Z = 14.0   N =    28   A = 28.085 g/mole
         --->  Isotope:  SI28   Z = 14   N =  28   A =  27.98 g/mole   abundance: 92.230 %
         --->  Isotope:  SI29   Z = 14   N =  29   A =  28.98 g/mole   abundance:  4.683 %
         --->  Isotope:  SI30   Z = 14   N =  30   A =  29.97 g/mole   abundance:  3.087 %
          ElmMassFraction: 100.00 %  ElmAbundance 100.00 % 

As you can see RadL and Nucl.Int.Length are not the same.

Even when they use the same formula, they will not produce the same value, since Geant treats Silicon as composition of several isotopes and ROOT as a single isotope.

[petricm]

It's even worse, if you look at Beryllium which does not have isotopes and repeat the exercise from the previous comment you get

ROOT  -> RadL:  34.4683 cm   Nucl.Int.Length:  39.4488 cm 
GEANT -> RadL:  35.276  cm   Nucl.Int.Length:  39.413 cm 

So the treatment is different.

[petricm]

Continuation, ROOT computes RadL with this code

if (radlen>=0) {
    //taken grom Geant3 routine GSMATE
     const Double_t alr2av=1.39621E-03, al183=5.20948;
     fRadLen = fA/(alr2av*fDensity*fZ*(fZ +TGeoMaterial::ScreenFactor(fZ))*
            (al183-TMath::Log(fZ)/3-TGeoMaterial::Coulomb(fZ)));
 }

and Geant does this in G4Material.cc

void G4Material::ComputeRadiationLength()
{
  G4double radinv = 0.0 ;
  for (G4int i=0;i<fNumberOfElements;++i) {
     radinv += VecNbOfAtomsPerVolume[i]*((*theElementVector)[i]->GetfRadTsai());
  }
  fRadlen = (radinv <= 0.0 ? DBL_MAX : 1./radinv);
}

which calls the Tsai formula, by calling this GetfRadTsai in G4Element.h

void G4Element::ComputeLradTsaiFactor()
{
  //
  //  Compute Tsai's Expression for the Radiation Length
  //  (Phys Rev. D50 3-1 (1994) page 1254)

  static const G4double Lrad_light[]  = {5.31  , 4.79  , 4.74 ,  4.71} ;
  static const G4double Lprad_light[] = {6.144 , 5.621 , 5.805 , 5.924} ;

  const G4double logZ3 = std::log(fZeff)/3.;

  G4double Lrad, Lprad;
  G4int iz = (G4int)(fZeff+0.5) - 1 ;
  if (iz <= 3) { Lrad = Lrad_light[iz] ;  Lprad = Lprad_light[iz] ; }
    else { Lrad = std::log(184.15) - logZ3 ; Lprad = std::log(1194.) - 2*logZ3;}

  fRadTsai = 4*alpha_rcl2*fZeff*(fZeff*(Lrad-fCoulomb) + Lprad); 
}

This is a completely different formula!

[agheata]

Hi, If having G4-like radiation/interaction length is a requirement, one can supply pre-computed G4 values to TGeoMaterial using SetRadLen(-g4radlen, -g4intlen).

[MarkusFrankATcernch]

Hi Andrei, Thanks for still listening :-). The deeper question is: Why the difference in the formulae? Are both equally correct or do both represent best knowledge of a topic, which in general is not too well known? If one is generally accepted to be correct, it should be feasible to standardize on one - or do I miss something?

[petricm]

Just for completeness ROOT formula for Nucl.Int.Length

   // Compute interaction length using the same formula as in GEANT4
   if (intlen>=0) {
      const Double_t cm = 1.;
      const Double_t g = 6.2415e21; // [gram = 1E-3*joule*s*s/(m*m)]
      const Double_t amu = 1.03642688246781065e-02; // [MeV/c^2]
      const Double_t lambda0 = 35.*g/(cm*cm);  // [g/cm^2]
      Double_t nilinv = 0.0;
      TGeoElement *elem = GetElement();
      if (!elem) {
         Fatal("SetRadLen", "Element not found for material %s", GetName());
         return;
      }
      Double_t nbAtomsPerVolume = TMath::Na()*fDensity/elem->A();
      nilinv += nbAtomsPerVolume*TMath::Power(elem->Neff(), 0.6666667);
      nilinv *= amu/lambda0;
      fIntLen = (nilinv<=0) ? TGeoShape::Big() : (1./nilinv);
   }

vs. Geant formula

void G4Material::ComputeNuclearInterLength()
{
  static const G4double lambda0  = 35*CLHEP::g/CLHEP::cm2;
  static const G4double twothird = 2.0/3.0;
  G4double NILinv = 0.0;
  for (G4int i=0; i<fNumberOfElements; ++i) {
    G4int Z = G4lrint( (*theElementVector)[i]->GetZ());
    G4double A = (*theElementVector)[i]->GetN();
    if(1 == Z) {
      NILinv += VecNbOfAtomsPerVolume[i]*A;
    } else {
      NILinv += VecNbOfAtomsPerVolume[i]*G4Exp(twothird*G4Log(A));
    } 
  }
  NILinv *= amu/lambda0; 
  fNuclInterLen = (NILinv <= 0.0 ? DBL_MAX : 1./NILinv);
}

Here the formula is the same between ROOT and Geant, but there is still some small numeric difference for Beryllium 39.4488 vs 39.413.

@agheata wouldn't it make sense, since you already use the Geant4 formula for Nucl.Int.Length to also update the RadL formula which is currently from Geant3? Or are there any reasons that the Geant3 formula is used?

[gaede]

The computation of the radiation length in Geant4 follows the on from Tsai, quoted in the PDG: http://pdg.lbl.gov/2013/reviews/rpp2013-rev-passage-particles-matter.pdf (p18). So I would imagine that this is the 'state of the art'.

[agheata]

Hi, The choice of the radiation length formula has only historical reasons, ALICE was the first user and we wanted to have the numbers matching Geant3 simulations. For the correctness, I cannot say which one is better, probably @gaede is right and it is the one in Geant4. It would indeed make sense to have the values matching Geant4, I'm only worried about backward compatibility and don't want to suddenly screw other people simulation/reconstruction results...

Since we are talking about a layer in DD4HEP, one possibility would be to apply the G4 formulas in a geometry afterburner and patch the materials/mixtures as I suggested before. It is even possible to actually implement this in TGeoMaterial/Mixture like say: TGeoMaterial::AveragePropertiesG4 which would be user callable. What do you think? I would accept this as contribution (to make sure of the numerical matching ;-) )

[gaede]

Simply adding TGeoMaterial::AveragePropertiesG4 is probably not enough, as the method is called (as a 'side-effect') by other methods such as AddElement() and SetDensity(). If one wants both behaviors, we'd either need a flag at construction ( useGeant4MaterialProperies ) or add a complete class TGeoMaterialG4...

[gaede]

... or the users (we) will have to know when to call TGeoMaterial::AveragePropertiesG4().

[agheata]

A static material flag would do it, however, when reading from file the user will wonder which parameters were used. So the flag should be persistent (new data member + getter/setter in TGeoManager).

[petricm]

Few more discrepancies regarding the ROOT implementation Nucl.Int.Length

const Double_t g = 6.2415e21;

this is CLHEP::g and is

const Double_t g = 6.24150964712041664e+21;

This could be improved, as you are already taking the amu definition as in CLHEP. This only makes a tiny difference.

There is another difference here:

The function TMaht::Na() return 6.022 141 99e+23 whereas in Geant defines 6.022 141 79e+23/mole; PDG 2016 is 6.022 140 857(74) vs NIST 6.022 140 857

Taking into account this I still don't get where does the difference for Nucl.Int.Length come from.

[MarkusFrankATcernch]

Both TGeoMaterial::GetInLen() and TGepMaterial::GetRadLen() are virtual. We could sub-class TGeoMaterial and override these two methods, so that the G4 formulae are used. Such a change would be transparent to any other client code. Since all materials are created explicitly in the geometry instantiation from the elements themselves sub-classing should be simple. At least as a temporary measure this should be thinkable if this would be a feasible way.

[MarkusFrankATcernch]

Coming back to the previous contribution.... I made a small practical test. What would be required is to declare the function TGeoMaterial::SetRadLen(...) as virtual. Then overloading:

class Mixture : public TGeoMixture  {
public:
  /// Initializing constructor
  Mixture(const char* name, Int_t nel, Double_t rho = -1);
  /// Default destructor
  virtual ~Mixture();
  /// Set radiation length and interaction length properties
  virtual void     SetRadLen(Double_t radlen, Double_t intlen=0.);      
};

would allow to apply the proper formulae as wanted.... At the same time such a change should not be invasive and also should not affect existing behaviour. Andrei, do you think this would be a (temporary) possibility?

[agheata]

Markus, this would be indeed non-intrusive, but quite heavy as one would need to anyway trigger (with a flag) the creation of these new materials, manage them transiently (e.g replacing TGeoMaterial with these) than deal with their persistency (or make sure the conversion is done whenever the TGeoManager is loaded from a file. Why don't you (someone) provide me the method dealing with the G4 formulae for rad/int length calculation, then I'll plug them nicely into the TGeoMaterial/Mixture class making sure of the consistency for dynamic changes (e.g density) , persistency, .... All the user would have to do is to set a flag that would trigger the automatic conversion.

[petricm]

@agheata would you have time to come on Thursday to our developers meeting and we can discuss this in person?

[agheata]

Yes, normally I'm free, I've put it in my agenda

[agheata]

Sorry, somehow I got the notification too late from my agenda, seems to be too late. Can I meet one of you in person to discuss this?

[petricm]

@agheata we have decided that we will go with your proposed solution and that I will provide you the functions.

[agheata]

Thanks, let me know.

[petricm]

Resolved with merge of root-project/root#1302

To get the new formula you have to call ComputeDerivedQuantities() on TGeoMixture then GetRadLen() and GetIntLen() return the G4 formula, otherwise you get the old behavior.