Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoMaterial.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 25/10/01
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 /** \class TGeoMaterial
13 \ingroup Geometry_classes
14 
15 Base class describing materials.
16 
17 \image html geom_material.jpg
18 */
19 
20 #include "Riostream.h"
21 #include "TMath.h"
22 #include "TObjArray.h"
23 #include "TStyle.h"
24 #include "TList.h"
25 #include "TGeoManager.h"
26 #include "TGeoExtension.h"
27 #include "TGeoMaterial.h"
28 #include "TGeoPhysicalConstants.h"
29 #include "TGeant4PhysicalConstants.h"
30 #include "TGDMLMatrix.h"
31 
32 // statics and globals
33 
34 ClassImp(TGeoMaterial);
35 
36 ////////////////////////////////////////////////////////////////////////////////
37 /// Default constructor
38 
39 TGeoMaterial::TGeoMaterial()
40  :TNamed(), TAttFill(),
41  fIndex(0),
42  fA(0.),
43  fZ(0.),
44  fDensity(0.),
45  fRadLen(0.),
46  fIntLen(0.),
47  fTemperature(0.),
48  fPressure(0.),
49  fState(kMatStateUndefined),
50  fShader(NULL),
51  fCerenkov(NULL),
52  fElement(NULL),
53  fUserExtension(0),
54  fFWExtension(0)
55 {
56  TGeoUnit::setUnitType(TGeoUnit::unitType()); // Ensure nobody changes the units afterwards
57  SetUsed(kFALSE);
58  fIndex = -1;
59  fTemperature = STP_temperature;
60  fPressure = STP_pressure;
61  fState = kMatStateUndefined;
62 }
63 
64 ////////////////////////////////////////////////////////////////////////////////
65 /// constructor
66 
67 TGeoMaterial::TGeoMaterial(const char *name)
68  :TNamed(name, ""), TAttFill(),
69  fIndex(0),
70  fA(0.),
71  fZ(0.),
72  fDensity(0.),
73  fRadLen(0.),
74  fIntLen(0.),
75  fTemperature(0.),
76  fPressure(0.),
77  fState(kMatStateUndefined),
78  fShader(NULL),
79  fCerenkov(NULL),
80  fElement(NULL),
81  fUserExtension(0),
82  fFWExtension(0)
83 {
84  TGeoUnit::setUnitType(TGeoUnit::unitType()); // Ensure nobody changes the units afterwards
85  fName = fName.Strip();
86  SetUsed(kFALSE);
87  fIndex = -1;
88  fTemperature = STP_temperature;
89  fPressure = STP_pressure;
90  fState = kMatStateUndefined;
91 
92  if (!gGeoManager) {
93  gGeoManager = new TGeoManager("Geometry", "default geometry");
94  }
95  gGeoManager->AddMaterial(this);
96 }
97 
98 ////////////////////////////////////////////////////////////////////////////////
99 /// constructor
100 
101 TGeoMaterial::TGeoMaterial(const char *name, Double_t a, Double_t z,
102  Double_t rho, Double_t radlen, Double_t intlen)
103  :TNamed(name, ""), TAttFill(),
104  fIndex(0),
105  fA(a),
106  fZ(z),
107  fDensity(rho),
108  fRadLen(0.),
109  fIntLen(0.),
110  fTemperature(0.),
111  fPressure(0.),
112  fState(kMatStateUndefined),
113  fShader(NULL),
114  fCerenkov(NULL),
115  fElement(NULL),
116  fUserExtension(0),
117  fFWExtension(0)
118 {
119  TGeoUnit::setUnitType(TGeoUnit::unitType()); // Ensure nobody changes the units afterwards
120  fName = fName.Strip();
121  SetUsed(kFALSE);
122  fIndex = -1;
123  fA = a;
124  fZ = z;
125  fDensity = rho;
126  fTemperature = STP_temperature;
127  fPressure = STP_pressure;
128  fState = kMatStateUndefined;
129  SetRadLen(radlen, intlen);
130  if (!gGeoManager) {
131  gGeoManager = new TGeoManager("Geometry", "default geometry");
132  }
133  if (fZ - Int_t(fZ) > 1E-3)
134  Warning("ctor", "Material %s defined with fractional Z=%f", GetName(), fZ);
135  if (GetElement()) GetElement()->SetUsed();
136  gGeoManager->AddMaterial(this);
137 }
138 
139 ////////////////////////////////////////////////////////////////////////////////
140 /// Constructor with state, temperature and pressure.
141 
142 TGeoMaterial::TGeoMaterial(const char *name, Double_t a, Double_t z, Double_t rho,
143  EGeoMaterialState state, Double_t temperature, Double_t pressure)
144  :TNamed(name, ""), TAttFill(),
145  fIndex(0),
146  fA(a),
147  fZ(z),
148  fDensity(rho),
149  fRadLen(0.),
150  fIntLen(0.),
151  fTemperature(temperature),
152  fPressure(pressure),
153  fState(state),
154  fShader(NULL),
155  fCerenkov(NULL),
156  fElement(NULL),
157  fUserExtension(0),
158  fFWExtension(0)
159 {
160  TGeoUnit::setUnitType(TGeoUnit::unitType()); // Ensure nobody changes the units afterwards
161  fName = fName.Strip();
162  SetUsed(kFALSE);
163  fIndex = -1;
164  SetRadLen(0,0);
165  if (!gGeoManager) {
166  gGeoManager = new TGeoManager("Geometry", "default geometry");
167  }
168  if (fZ - Int_t(fZ) > 1E-3)
169  Warning("ctor", "Material %s defined with fractional Z=%f", GetName(), fZ);
170  if (GetElement()) GetElement()->SetUsed();
171  gGeoManager->AddMaterial(this);
172 }
173 
174 ////////////////////////////////////////////////////////////////////////////////
175 /// constructor
176 
177 TGeoMaterial::TGeoMaterial(const char *name, TGeoElement *elem, Double_t rho)
178  :TNamed(name, ""), TAttFill(),
179  fIndex(0),
180  fA(0.),
181  fZ(0.),
182  fDensity(rho),
183  fRadLen(0.),
184  fIntLen(0.),
185  fTemperature(0.),
186  fPressure(0.),
187  fState(kMatStateUndefined),
188  fShader(NULL),
189  fCerenkov(NULL),
190  fElement(elem),
191  fUserExtension(0),
192  fFWExtension(0)
193 {
194  TGeoUnit::setUnitType(TGeoUnit::unitType()); // Ensure nobody changes the units afterwards
195  fName = fName.Strip();
196  SetUsed(kFALSE);
197  fIndex = -1;
198  fA = elem->A();
199  fZ = elem->Z();
200  SetRadLen(0,0);
201  fTemperature = STP_temperature;
202  fPressure = STP_pressure;
203  fState = kMatStateUndefined;
204  if (!gGeoManager) {
205  gGeoManager = new TGeoManager("Geometry", "default geometry");
206  }
207  if (fZ - Int_t(fZ) > 1E-3)
208  Warning("ctor", "Material %s defined with fractional Z=%f", GetName(), fZ);
209  if (GetElement()) GetElement()->SetUsed();
210  gGeoManager->AddMaterial(this);
211 }
212 
213 ////////////////////////////////////////////////////////////////////////////////
214 
215 TGeoMaterial::TGeoMaterial(const TGeoMaterial& gm) :
216  TNamed(gm),
217  TAttFill(gm),
218  fIndex(gm.fIndex),
219  fA(gm.fA),
220  fZ(gm.fZ),
221  fDensity(gm.fDensity),
222  fRadLen(gm.fRadLen),
223  fIntLen(gm.fIntLen),
224  fTemperature(gm.fTemperature),
225  fPressure(gm.fPressure),
226  fState(gm.fState),
227  fShader(gm.fShader),
228  fCerenkov(gm.fCerenkov),
229  fElement(gm.fElement),
230  fUserExtension(gm.fUserExtension->Grab()),
231  fFWExtension(gm.fFWExtension->Grab())
232 
233 {
234  //copy constructor
235  TGeoUnit::setUnitType(TGeoUnit::unitType()); // Ensure nobody changes the units afterwards
236  fProperties.SetOwner();
237  TIter next(&fProperties);
238  TNamed *property;
239  while ((property = (TNamed*)next())) fProperties.Add(new TNamed(*property));
240 }
241 
242 ////////////////////////////////////////////////////////////////////////////////
243 ///assignment operator
244 
245 TGeoMaterial& TGeoMaterial::operator=(const TGeoMaterial& gm)
246 {
247  if(this!=&gm) {
248  TNamed::operator=(gm);
249  TAttFill::operator=(gm);
250  fIndex=gm.fIndex;
251  fA=gm.fA;
252  fZ=gm.fZ;
253  fDensity=gm.fDensity;
254  fRadLen=gm.fRadLen;
255  fIntLen=gm.fIntLen;
256  fTemperature=gm.fTemperature;
257  fPressure=gm.fPressure;
258  fState=gm.fState;
259  fShader=gm.fShader;
260  fCerenkov=gm.fCerenkov;
261  fElement=gm.fElement;
262  fUserExtension = gm.fUserExtension->Grab();
263  fFWExtension = gm.fFWExtension->Grab();
264  fProperties.SetOwner();
265  TIter next(&fProperties);
266  TNamed *property;
267  while ((property = (TNamed*)next())) fProperties.Add(new TNamed(*property));
268  }
269  return *this;
270 }
271 
272 ////////////////////////////////////////////////////////////////////////////////
273 /// Destructor
274 
275 TGeoMaterial::~TGeoMaterial()
276 {
277  if (fUserExtension) {fUserExtension->Release(); fUserExtension=0;}
278  if (fFWExtension) {fFWExtension->Release(); fFWExtension=0;}
279 }
280 
281 ////////////////////////////////////////////////////////////////////////////////
282 /// Connect user-defined extension to the material. The material "grabs" a copy, so
283 /// the original object can be released by the producer. Release the previously
284 /// connected extension if any.
285 ///
286 /// NOTE: This interface is intended for user extensions and is guaranteed not
287 /// to be used by TGeo
288 
289 void TGeoMaterial::SetUserExtension(TGeoExtension *ext)
290 {
291  if (fUserExtension) fUserExtension->Release();
292  fUserExtension = 0;
293  if (ext) fUserExtension = ext->Grab();
294 }
295 
296 //_____________________________________________________________________________
297 const char *TGeoMaterial::GetPropertyRef(const char *property) const
298 {
299  // Find reference for a given property
300  TNamed *prop = (TNamed*)fProperties.FindObject(property);
301  return (prop) ? prop->GetTitle() : nullptr;
302 }
303 
304 //_____________________________________________________________________________
305 TGDMLMatrix *TGeoMaterial::GetProperty(const char *property) const
306 {
307  // Find reference for a given property
308  TNamed *prop = (TNamed*)fProperties.FindObject(property);
309  if ( !prop ) return nullptr;
310  return gGeoManager->GetGDMLMatrix(prop->GetTitle());
311 }
312 
313 //_____________________________________________________________________________
314 TGDMLMatrix *TGeoMaterial::GetProperty(Int_t i) const
315 {
316  // Find reference for a given property
317  TNamed *prop = (TNamed*)fProperties.At(i);
318  if ( !prop ) return nullptr;
319  return gGeoManager->GetGDMLMatrix(prop->GetTitle());
320 }
321 
322 //_____________________________________________________________________________
323 const char *TGeoMaterial::GetConstPropertyRef(const char *property) const
324 {
325  // Find reference for a given constant property
326  TNamed *prop = (TNamed*)fConstProperties.FindObject(property);
327  return (prop) ? prop->GetTitle() : nullptr;
328 }
329 
330 //_____________________________________________________________________________
331 Double_t TGeoMaterial::GetConstProperty(const char *property, Bool_t *err) const
332 {
333  // Find reference for a given constant property
334  TNamed *prop = (TNamed*)fConstProperties.FindObject(property);
335  if (!prop) {
336  if (err) *err = kTRUE;
337  return 0.;
338  }
339  return gGeoManager->GetProperty(prop->GetTitle(), err);
340 }
341 
342 //_____________________________________________________________________________
343 Double_t TGeoMaterial::GetConstProperty(Int_t i, Bool_t *err) const
344 {
345  // Find reference for a given constant property
346  TNamed *prop = (TNamed*)fConstProperties.At(i);
347  if (!prop) {
348  if (err) *err = kTRUE;
349  return 0.;
350  }
351  return gGeoManager->GetProperty(prop->GetTitle(), err);
352 }
353 
354 //_____________________________________________________________________________
355 bool TGeoMaterial::AddProperty(const char *property, const char *ref)
356 {
357  fProperties.SetOwner();
358  if (GetPropertyRef(property)) {
359  Error("AddProperty", "Property %s already added to material %s",
360  property, GetName());
361  return false;
362  }
363  fProperties.Add(new TNamed(property, ref));
364  return true;
365 }
366 
367 //_____________________________________________________________________________
368 bool TGeoMaterial::AddConstProperty(const char *property, const char *ref)
369 {
370  fConstProperties.SetOwner();
371  if (GetConstPropertyRef(property)) {
372  Error("AddConstProperty", "Constant property %s already added to material %s",
373  property, GetName());
374  return false;
375  }
376  fConstProperties.Add(new TNamed(property, ref));
377  return true;
378 }
379 
380 ////////////////////////////////////////////////////////////////////////////////
381 /// Connect framework defined extension to the material. The material "grabs" a copy,
382 /// so the original object can be released by the producer. Release the previously
383 /// connected extension if any.
384 ///
385 /// NOTE: This interface is intended for the use by TGeo and the users should
386 /// NOT connect extensions using this method
387 
388 void TGeoMaterial::SetFWExtension(TGeoExtension *ext)
389 {
390  if (fFWExtension) fFWExtension->Release();
391  fFWExtension = 0;
392  if (ext) fFWExtension = ext->Grab();
393 }
394 
395 ////////////////////////////////////////////////////////////////////////////////
396 /// Get a copy of the user extension pointer. The user must call Release() on
397 /// the copy pointer once this pointer is not needed anymore (equivalent to
398 /// delete() after calling new())
399 
400 TGeoExtension *TGeoMaterial::GrabUserExtension() const
401 {
402  if (fUserExtension) return fUserExtension->Grab();
403  return 0;
404 }
405 
406 ////////////////////////////////////////////////////////////////////////////////
407 /// Get a copy of the framework extension pointer. The user must call Release() on
408 /// the copy pointer once this pointer is not needed anymore (equivalent to
409 /// delete() after calling new())
410 
411 TGeoExtension *TGeoMaterial::GrabFWExtension() const
412 {
413  if (fFWExtension) return fFWExtension->Grab();
414  return 0;
415 }
416 
417 ////////////////////////////////////////////////////////////////////////////////
418 /// Provide a pointer name containing uid.
419 
420 char *TGeoMaterial::GetPointerName() const
421 {
422  static TString name;
423  name = TString::Format("pMat%d", GetUniqueID());
424  return (char*)name.Data();
425 }
426 
427 ////////////////////////////////////////////////////////////////////////////////
428 /// Set radiation/absorption lengths. If the values are negative, their absolute value
429 /// is taken, otherwise radlen is recomputed using G3 formula.
430 
431 void TGeoMaterial::SetRadLen(Double_t radlen, Double_t intlen)
432 {
433  fRadLen = TMath::Abs(radlen);
434  fIntLen = TMath::Abs(intlen);
435  // Check for vacuum
436  if (fA<0.9 || fZ<0.9) {
437  if (radlen<-1e5 || intlen<-1e-5) {
438  Error("SetRadLen","Material %s: user values taken for vacuum: radlen=%g or intlen=%g - too small", GetName(),fRadLen, fIntLen);
439  return;
440  }
441  // Ignore positive values and take big numbers
442  if (radlen>=0) fRadLen = 1.E30;
443  if (intlen>=0) fIntLen = 1.E30;
444  return;
445  }
446  TGeoUnit::UnitType typ = TGeoUnit::unitType();
447  // compute radlen systematically with G3 formula for a valid material
448  if ( typ == TGeoUnit::kTGeoUnits && radlen>=0 ) {
449  //taken grom Geant3 routine GSMATE
450  constexpr Double_t alr2av = 1.39621E-03*TGeoUnit::cm2;
451  constexpr Double_t al183 = 5.20948;
452  fRadLen = fA/(alr2av*fDensity*fZ*(fZ +TGeoMaterial::ScreenFactor(fZ))*
453  (al183-TMath::Log(fZ)/3-TGeoMaterial::Coulomb(fZ)));
454  fRadLen *= TGeoUnit::cm;
455  }
456  else if ( typ == TGeoUnit::kTGeant4Units && radlen>=0 ) {
457  //taken grom Geant3 routine GSMATE
458  constexpr Double_t alr2av = 1.39621E-03*TGeant4Unit::cm2;
459  constexpr Double_t al183 = 5.20948;
460  fRadLen = fA/(alr2av*fDensity*fZ*(fZ +TGeoMaterial::ScreenFactor(fZ))*
461  (al183-TMath::Log(fZ)/3-TGeoMaterial::Coulomb(fZ)));
462  fRadLen *= TGeant4Unit::cm;
463  }
464  // Compute interaction length using the same formula as in GEANT4
465  if ( typ == TGeoUnit::kTGeoUnits && intlen>=0 ) {
466  constexpr Double_t lambda0 = 35.*TGeoUnit::g/TGeoUnit::cm2; // [g/cm^2]
467  Double_t nilinv = 0.0;
468  TGeoElement *elem = GetElement();
469  if (!elem) {
470  Fatal("SetRadLen", "Element not found for material %s", GetName());
471  return;
472  }
473  Double_t nbAtomsPerVolume = TGeoUnit::Avogadro*fDensity/elem->A();
474  nilinv += nbAtomsPerVolume*TMath::Power(elem->Neff(), 0.6666667);
475  nilinv *= TGeoUnit::amu/lambda0;
476  fIntLen = (nilinv<=0) ? TGeoShape::Big() : (TGeoUnit::cm/nilinv);
477  }
478  else if ( typ == TGeoUnit::kTGeant4Units && intlen>=0 ) {
479  constexpr Double_t lambda0 = 35.*TGeant4Unit::g/TGeant4Unit::cm2; // [g/cm^2]
480  Double_t nilinv = 0.0;
481  TGeoElement *elem = GetElement();
482  if (!elem) {
483  Fatal("SetRadLen", "Element not found for material %s", GetName());
484  return;
485  }
486  Double_t nbAtomsPerVolume = TGeant4Unit::Avogadro*fDensity/elem->A();
487  nilinv += nbAtomsPerVolume*TMath::Power(elem->Neff(), 0.6666667);
488  nilinv *= TGeant4Unit::amu/lambda0;
489  fIntLen = (nilinv<=0) ? TGeoShape::Big() : (TGeant4Unit::cm/nilinv);
490  }
491 }
492 
493 ////////////////////////////////////////////////////////////////////////////////
494 /// static function
495 /// Compute Coulomb correction for pair production and Brem
496 /// REFERENCE : EGS MANUAL SLAC 210 - UC32 - JUNE 78
497 /// FORMULA 2.7.17
498 
499 Double_t TGeoMaterial::Coulomb(Double_t z)
500 {
501  Double_t az = TGeoUnit::unitType() == TGeoUnit::kTGeoUnits
502  ? TGeoUnit::fine_structure_const*z : TGeant4Unit::fine_structure_const*z;
503  Double_t az2 = az*az;
504  Double_t az4 = az2 * az2;
505  Double_t fp = ( 0.0083*az4 + 0.20206 + 1./(1.+az2) ) * az2;
506  Double_t fm = ( 0.0020*az4 + 0.0369 ) * az4;
507  return fp - fm;
508 }
509 
510 ////////////////////////////////////////////////////////////////////////////////
511 /// return true if the other material has the same physical properties
512 
513 Bool_t TGeoMaterial::IsEq(const TGeoMaterial *other) const
514 {
515  if (other==this) return kTRUE;
516  if (other->IsMixture()) return kFALSE;
517  if (TMath::Abs(fA-other->GetA())>1E-3) return kFALSE;
518  if (TMath::Abs(fZ-other->GetZ())>1E-3) return kFALSE;
519  if (TMath::Abs(fDensity-other->GetDensity())>1E-6) return kFALSE;
520  if (GetCerenkovProperties() != other->GetCerenkovProperties()) return kFALSE;
521 // if (fRadLen != other->GetRadLen()) return kFALSE;
522 // if (fIntLen != other->GetIntLen()) return kFALSE;
523  return kTRUE;
524 }
525 
526 ////////////////////////////////////////////////////////////////////////////////
527 /// print characteristics of this material
528 
529 void TGeoMaterial::Print(const Option_t * /*option*/) const
530 {
531  printf("Material %s %s A=%g Z=%g rho=%g radlen=%g intlen=%g index=%i\n", GetName(), GetTitle(),
532  fA,fZ,fDensity, fRadLen, fIntLen, fIndex);
533 }
534 
535 ////////////////////////////////////////////////////////////////////////////////
536 /// Save a primitive as a C++ statement(s) on output stream "out".
537 
538 void TGeoMaterial::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
539 {
540  if (TestBit(TGeoMaterial::kMatSavePrimitive)) return;
541  char *name = GetPointerName();
542  out << "// Material: " << GetName() << std::endl;
543  out << " a = " << fA << ";" << std::endl;
544  out << " z = " << fZ << ";" << std::endl;
545  out << " density = " << fDensity << ";" << std::endl;
546  out << " radl = " << fRadLen << ";" << std::endl;
547  out << " absl = " << fIntLen << ";" << std::endl;
548 
549  out << " " << name << " = new TGeoMaterial(\"" << GetName() << "\", a,z,density,radl,absl);" << std::endl;
550  out << " " << name << "->SetIndex(" << GetIndex() << ");" << std::endl;
551  SetBit(TGeoMaterial::kMatSavePrimitive);
552 }
553 
554 ////////////////////////////////////////////////////////////////////////////////
555 /// Get some default color related to this material.
556 
557 Int_t TGeoMaterial::GetDefaultColor() const
558 {
559  Int_t id = 1+ gGeoManager->GetListOfMaterials()->IndexOf(this);
560  return (2+id%6);
561 }
562 
563 ////////////////////////////////////////////////////////////////////////////////
564 /// Get a pointer to the element this material is made of.
565 
566 TGeoElement *TGeoMaterial::GetElement(Int_t) const
567 {
568  if (fElement) return fElement;
569  TGeoElementTable *table = gGeoManager->GetElementTable();
570  return table->GetElement(Int_t(fZ));
571 }
572 ////////////////////////////////////////////////////////////////////////////////
573 /// Single interface to get element properties.
574 
575 void TGeoMaterial::GetElementProp(Double_t &a, Double_t &z, Double_t &w, Int_t)
576 {
577  a = fA;
578  z = fZ;
579  w = 1.;
580 }
581 
582 ////////////////////////////////////////////////////////////////////////////////
583 /// Retrieve material index in the list of materials
584 
585 Int_t TGeoMaterial::GetIndex()
586 {
587  if (fIndex>=0) return fIndex;
588  TList *matlist = gGeoManager->GetListOfMaterials();
589  fIndex = matlist->IndexOf(this);
590  return fIndex;
591 }
592 
593 ////////////////////////////////////////////////////////////////////////////////
594 /// Create the material representing the decay product of this material at a
595 /// given time. The precision represent the minimum cumulative branching ratio for
596 /// which decay products are still taken into account.
597 
598 TGeoMaterial *TGeoMaterial::DecayMaterial(Double_t time, Double_t precision)
599 {
600  TObjArray *pop = new TObjArray();
601  if (!fElement || !fElement->IsRadioNuclide()) return this;
602  FillMaterialEvolution(pop, precision);
603  Int_t ncomp = pop->GetEntriesFast();
604  if (!ncomp) return this;
605  TGeoElementRN *el;
606  Double_t *weight = new Double_t[ncomp];
607  Double_t amed = 0.;
608  Int_t i;
609  for (i=0; i<ncomp; i++) {
610  el = (TGeoElementRN *)pop->At(i);
611  weight[i] = el->Ratio()->Concentration(time) * el->A();
612  amed += weight[i];
613  }
614  Double_t rho = fDensity*amed/fA;
615  TGeoMixture *mix = 0;
616  Int_t ncomp1 = ncomp;
617  for (i=0; i<ncomp; i++) {
618  if ((weight[i]/amed)<precision) {
619  amed -= weight[i];
620  ncomp1--;
621  }
622  }
623  if (ncomp1<2) {
624  el = (TGeoElementRN *)pop->At(0);
625  delete [] weight;
626  delete pop;
627  if (ncomp1==1) return new TGeoMaterial(TString::Format("%s-evol",GetName()), el, rho);
628  return NULL;
629  }
630  mix = new TGeoMixture(TString::Format("%s-evol",GetName()), ncomp, rho);
631  for (i=0; i<ncomp; i++) {
632  weight[i] /= amed;
633  if (weight[i]<precision) continue;
634  el = (TGeoElementRN *)pop->At(i);
635  mix->AddElement(el, weight[i]);
636  }
637  delete [] weight;
638  delete pop;
639  return mix;
640 }
641 
642 ////////////////////////////////////////////////////////////////////////////////
643 /// Fills a user array with all the elements deriving from the possible
644 /// decay of the top element composing the mixture. Each element contained
645 /// by <population> may be a radionuclide having a Bateman solution attached.
646 /// The precision represent the minimum cumulative branching ratio for
647 /// which decay products are still taken into account.
648 /// To visualize the time evolution of each decay product one can use:
649 /// ~~~ {.cpp}
650 /// TGeoElement *elem = population->At(index);
651 /// TGeoElementRN *elemrn = 0;
652 /// if (elem->IsRadioNuclide()) elemrn = (TGeoElementRN*)elem;
653 /// ~~~
654 /// One can get Ni/N1(t=0) at any moment of time. Ni is the number of atoms
655 /// of one of the decay products, N1(0) is the number of atoms of the top
656 /// element at t=0.
657 /// ~~~ {.cpp}
658 /// Double_t fraction_weight = elemrn->Ratio()->Concentration(time);
659 /// ~~~
660 /// One can also display the time evolution of the fractional weight:
661 /// ~~~ {.cpp}
662 /// elemrn->Ratio()->Draw(option);
663 /// ~~~
664 
665 void TGeoMaterial::FillMaterialEvolution(TObjArray *population, Double_t precision)
666 {
667  if (population->GetEntriesFast()) {
668  Error("FillMaterialEvolution", "Provide an empty array !");
669  return;
670  }
671  TGeoElementTable *table = gGeoManager->GetElementTable();
672  TGeoElement *elem;
673  TGeoElementRN *elemrn;
674  TIter next(table->GetElementsRN());
675  while ((elemrn=(TGeoElementRN*)next())) elemrn->ResetRatio();
676  elem = GetElement();
677  if (!elem) {
678  Fatal("FillMaterialEvolution", "Element not found for material %s", GetName());
679  return;
680  }
681  if (!elem->IsRadioNuclide()) {
682  population->Add(elem);
683  return;
684  }
685  elemrn = (TGeoElementRN*)elem;
686  elemrn->FillPopulation(population, precision);
687 }
688 
689 /** \class TGeoMixture
690 \ingroup Geometry_classes
691 
692 Mixtures of elements.
693 
694 */
695 
696 ClassImp(TGeoMixture);
697 
698 ////////////////////////////////////////////////////////////////////////////////
699 /// Default constructor
700 
701 TGeoMixture::TGeoMixture()
702 {
703  fNelements = 0;
704  fZmixture = 0;
705  fAmixture = 0;
706  fWeights = 0;
707  fNatoms = 0;
708  fVecNbOfAtomsPerVolume = 0;
709  fElements = 0;
710 }
711 
712 ////////////////////////////////////////////////////////////////////////////////
713 /// constructor
714 
715 TGeoMixture::TGeoMixture(const char *name, Int_t /*nel*/, Double_t rho)
716  :TGeoMaterial(name)
717 {
718  fZmixture = 0;
719  fAmixture = 0;
720  fWeights = 0;
721  fNelements = 0;
722  fNatoms = 0;
723  fVecNbOfAtomsPerVolume = 0;
724  fDensity = rho;
725  fElements = 0;
726  if (fDensity < 0) fDensity = 0.001;
727 }
728 
729 ////////////////////////////////////////////////////////////////////////////////
730 ///copy constructor
731 
732 TGeoMixture::TGeoMixture(const TGeoMixture& gm) :
733  TGeoMaterial(gm),
734  fNelements(gm.fNelements),
735  fZmixture(gm.fZmixture),
736  fAmixture(gm.fAmixture),
737  fWeights(gm.fWeights),
738  fNatoms(gm.fNatoms),
739  fVecNbOfAtomsPerVolume(gm.fVecNbOfAtomsPerVolume),
740  fElements(gm.fElements)
741 {
742 }
743 
744 ////////////////////////////////////////////////////////////////////////////////
745 ///assignment operator
746 
747 TGeoMixture& TGeoMixture::operator=(const TGeoMixture& gm)
748 {
749  if(this!=&gm) {
750  TGeoMaterial::operator=(gm);
751  fNelements=gm.fNelements;
752  fZmixture=gm.fZmixture;
753  fAmixture=gm.fAmixture;
754  fWeights=gm.fWeights;
755  fNatoms = gm.fNatoms;
756  fVecNbOfAtomsPerVolume = gm.fVecNbOfAtomsPerVolume;
757  fElements = gm.fElements;
758  }
759  return *this;
760 }
761 
762 ////////////////////////////////////////////////////////////////////////////////
763 /// Destructor
764 
765 TGeoMixture::~TGeoMixture()
766 {
767  if (fZmixture) delete[] fZmixture;
768  if (fAmixture) delete[] fAmixture;
769  if (fWeights) delete[] fWeights;
770  if (fNatoms) delete[] fNatoms;
771  if (fVecNbOfAtomsPerVolume) delete[] fVecNbOfAtomsPerVolume;
772  if (fElements) delete fElements;
773 }
774 
775 ////////////////////////////////////////////////////////////////////////////////
776 /// Compute effective A/Z and radiation length
777 
778 void TGeoMixture::AverageProperties()
779 {
780  TGeoUnit::UnitType typ = TGeoUnit::unitType();
781  const Double_t cm = (typ==TGeoUnit::kTGeoUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
782  const Double_t cm2 = (typ==TGeoUnit::kTGeoUnits) ? TGeoUnit::cm2 : TGeant4Unit::cm2;
783  const Double_t amu = (typ==TGeoUnit::kTGeoUnits) ? TGeoUnit::amu : TGeant4Unit::amu; // [MeV/c^2]
784  const Double_t gram = (typ==TGeoUnit::kTGeoUnits) ? TGeoUnit::gram : TGeant4Unit::gram;
785  const Double_t na = (typ==TGeoUnit::kTGeoUnits) ? TGeoUnit::Avogadro : TGeant4Unit::Avogadro;
786  const Double_t alr2av = 1.39621E-03 * cm2;
787  const Double_t al183 = 5.20948;
788  const Double_t lambda0 = 35.*gram/cm2; // [g/cm^2]
789  Double_t radinv = 0.0;
790  Double_t nilinv = 0.0;
791  Double_t nbAtomsPerVolume;
792  fA = 0;
793  fZ = 0;
794  for (Int_t j=0;j<fNelements;j++) {
795  if (fWeights[j] <= 0) continue;
796  fA += fWeights[j]*fAmixture[j];
797  fZ += fWeights[j]*fZmixture[j];
798  nbAtomsPerVolume = na*fDensity*fWeights[j]/GetElement(j)->A();
799  nilinv += nbAtomsPerVolume*TMath::Power(GetElement(j)->Neff(), 0.6666667);
800  Double_t zc = fZmixture[j];
801  Double_t alz = TMath::Log(zc)/3.;
802  Double_t xinv = zc*(zc+TGeoMaterial::ScreenFactor(zc))*
803  (al183-alz-TGeoMaterial::Coulomb(zc))/fAmixture[j];
804  radinv += xinv*fWeights[j];
805  }
806  radinv *= alr2av*fDensity;
807  if (radinv > 0) fRadLen = cm/radinv;
808  // Compute interaction length
809  nilinv *= amu/lambda0;
810  fIntLen = (nilinv<=0) ? TGeoShape::Big() : (cm/nilinv);
811 }
812 
813 ////////////////////////////////////////////////////////////////////////////////
814 /// add an element to the mixture using fraction by weight
815 /// Check if the element is already defined
816 
817 void TGeoMixture::AddElement(Double_t a, Double_t z, Double_t weight)
818 {
819  TGeoElementTable *table = gGeoManager->GetElementTable();
820  if (z<1 || z>table->GetNelements()-1)
821  Fatal("AddElement", "Cannot add element having Z=%d to mixture %s", (Int_t)z, GetName());
822  Int_t i;
823  for (i=0; i<fNelements; i++) {
824  if (TMath::Abs(z-fZmixture[i])<1.e-6 && TMath::Abs(a-fAmixture[i])<1.e-6) {
825  fWeights[i] += weight;
826  AverageProperties();
827  return;
828  }
829  }
830  if (!fNelements) {
831  fZmixture = new Double_t[1];
832  fAmixture = new Double_t[1];
833  fWeights = new Double_t[1];
834  } else {
835  Int_t nelements = fNelements+1;
836  Double_t *zmixture = new Double_t[nelements];
837  Double_t *amixture = new Double_t[nelements];
838  Double_t *weights = new Double_t[nelements];
839  for (Int_t j=0; j<fNelements; j++) {
840  zmixture[j] = fZmixture[j];
841  amixture[j] = fAmixture[j];
842  weights[j] = fWeights[j];
843  }
844  delete [] fZmixture;
845  delete [] fAmixture;
846  delete [] fWeights;
847  fZmixture = zmixture;
848  fAmixture = amixture;
849  fWeights = weights;
850  }
851 
852  fNelements++;
853  i = fNelements - 1;
854  fZmixture[i] = z;
855  fAmixture[i] = a;
856  fWeights[i] = weight;
857  if (z - Int_t(z) > 1E-3)
858  Warning("DefineElement", "Mixture %s has element defined with fractional Z=%f", GetName(), z);
859  GetElement(i)->SetDefined();
860  table->GetElement((Int_t)z)->SetDefined();
861 
862  //compute equivalent radiation length (taken from Geant3/GSMIXT)
863  AverageProperties();
864 }
865 
866 ////////////////////////////////////////////////////////////////////////////////
867 /// Define one component of the mixture as an existing material/mixture.
868 
869 void TGeoMixture::AddElement(TGeoMaterial *mat, Double_t weight)
870 {
871  TGeoElement *elnew, *elem;
872  Double_t a,z;
873  if (!mat->IsMixture()) {
874  elem = mat->GetBaseElement();
875  if (elem) {
876  AddElement(elem, weight);
877  } else {
878  a = mat->GetA();
879  z = mat->GetZ();
880  AddElement(a, z, weight);
881  }
882  return;
883  }
884  // The material is a mixture.
885  TGeoMixture *mix = (TGeoMixture*)mat;
886  Double_t wnew;
887  Int_t nelem = mix->GetNelements();
888  Bool_t elfound;
889  Int_t i,j;
890  // loop the elements of the daughter mixture
891  for (i=0; i<nelem; i++) {
892  elfound = kFALSE;
893  elnew = mix->GetElement(i);
894  if (!elnew) continue;
895  // check if we have the element already defined in the parent mixture
896  for (j=0; j<fNelements; j++) {
897  if (fWeights[j]<=0) continue;
898  elem = GetElement(j);
899  if (elem == elnew) {
900  // element found, compute new weight
901  fWeights[j] += weight * (mix->GetWmixt())[i];
902  elfound = kTRUE;
903  break;
904  }
905  }
906  if (elfound) continue;
907  // element not found, define it
908  wnew = weight * (mix->GetWmixt())[i];
909  AddElement(elnew, wnew);
910  }
911 }
912 
913 ////////////////////////////////////////////////////////////////////////////////
914 /// add an element to the mixture using fraction by weight
915 
916 void TGeoMixture::AddElement(TGeoElement *elem, Double_t weight)
917 {
918  TGeoElement *elemold;
919  TGeoElementTable *table = gGeoManager->GetElementTable();
920  if (!fElements) fElements = new TObjArray(128);
921  Bool_t exist = kFALSE;
922  // If previous elements were defined by A/Z, add corresponding TGeoElements
923  for (Int_t i=0; i<fNelements; i++) {
924  elemold = (TGeoElement*)fElements->At(i);
925  if (!elemold) fElements->AddAt(elemold = table->GetElement((Int_t)fZmixture[i]), i);
926  if (elemold == elem) exist = kTRUE;
927  }
928  if (!exist) fElements->AddAtAndExpand(elem, fNelements);
929  AddElement(elem->A(), elem->Z(), weight);
930 }
931 
932 ////////////////////////////////////////////////////////////////////////////////
933 /// Add a mixture element by number of atoms in the chemical formula.
934 
935 void TGeoMixture::AddElement(TGeoElement *elem, Int_t natoms)
936 {
937  Int_t i,j;
938  Double_t amol;
939  TGeoElement *elemold;
940  TGeoElementTable *table = gGeoManager->GetElementTable();
941  if (!fElements) fElements = new TObjArray(128);
942  // Check if the element is already defined
943  for (i=0; i<fNelements; i++) {
944  elemold = (TGeoElement*)fElements->At(i);
945  if (!elemold) fElements->AddAt(table->GetElement((Int_t)fZmixture[i]), i);
946  else if (elemold != elem) continue;
947  if ((elem==elemold) ||
948  (TMath::Abs(elem->Z()-fZmixture[i])<1.e-6 && TMath::Abs(elem->A()-fAmixture[i])<1.e-6)) {
949  fNatoms[i] += natoms;
950  amol = 0.;
951  for (j=0; j<fNelements; j++) amol += fAmixture[j]*fNatoms[j];
952  for (j=0; j<fNelements; j++) fWeights[j] = fNatoms[j]*fAmixture[j]/amol;
953  AverageProperties();
954  return;
955  }
956  }
957  // New element
958  if (!fNelements) {
959  fZmixture = new Double_t[1];
960  fAmixture = new Double_t[1];
961  fWeights = new Double_t[1];
962  fNatoms = new Int_t[1];
963  } else {
964  if (!fNatoms) {
965  Fatal("AddElement", "Cannot add element by natoms in mixture %s after defining elements by weight",
966  GetName());
967  return;
968  }
969  Int_t nelements = fNelements+1;
970  Double_t *zmixture = new Double_t[nelements];
971  Double_t *amixture = new Double_t[nelements];
972  Double_t *weights = new Double_t[nelements];
973  Int_t *nnatoms = new Int_t[nelements];
974  for (j=0; j<fNelements; j++) {
975  zmixture[j] = fZmixture[j];
976  amixture[j] = fAmixture[j];
977  weights[j] = fWeights[j];
978  nnatoms[j] = fNatoms[j];
979  }
980  delete [] fZmixture;
981  delete [] fAmixture;
982  delete [] fWeights;
983  delete [] fNatoms;
984  fZmixture = zmixture;
985  fAmixture = amixture;
986  fWeights = weights;
987  fNatoms = nnatoms;
988  }
989  fNelements++;
990  Int_t iel = fNelements-1;
991  fZmixture[iel] = elem->Z();
992  fAmixture[iel] = elem->A();
993  fNatoms[iel] = natoms;
994  fElements->AddAtAndExpand(elem, iel);
995  amol = 0.;
996  for (i=0; i<fNelements; i++) {
997  if (fNatoms[i]<=0) return;
998  amol += fAmixture[i]*fNatoms[i];
999  }
1000  for (i=0; i<fNelements; i++) fWeights[i] = fNatoms[i]*fAmixture[i]/amol;
1001  table->GetElement(elem->Z())->SetDefined();
1002  AverageProperties();
1003 }
1004 
1005 ////////////////////////////////////////////////////////////////////////////////
1006 /// Define the mixture element at index iel by number of atoms in the chemical formula.
1007 
1008 void TGeoMixture::DefineElement(Int_t /*iel*/, Int_t z, Int_t natoms)
1009 {
1010  TGeoElementTable *table = gGeoManager->GetElementTable();
1011  TGeoElement *elem = table->GetElement(z);
1012  if (!elem) {
1013  Fatal("DefineElement", "In mixture %s, element with Z=%i not found",GetName(),z);
1014  return;
1015  }
1016  AddElement(elem, natoms);
1017 }
1018 
1019 ////////////////////////////////////////////////////////////////////////////////
1020 /// Retrieve the pointer to the element corresponding to component I.
1021 
1022 TGeoElement *TGeoMixture::GetElement(Int_t i) const
1023 {
1024  if (i<0 || i>=fNelements) {
1025  Error("GetElement", "Mixture %s has only %d elements", GetName(), fNelements);
1026  return 0;
1027  }
1028  TGeoElement *elem = 0;
1029  if (fElements) elem = (TGeoElement*)fElements->At(i);
1030  if (elem) return elem;
1031  TGeoElementTable *table = gGeoManager->GetElementTable();
1032  return table->GetElement(Int_t(fZmixture[i]));
1033 }
1034 
1035 ////////////////////////////////////////////////////////////////////////////////
1036 /// Get specific activity (in Bq/gram) for the whole mixture (no argument) or
1037 /// for a given component.
1038 
1039 Double_t TGeoMixture::GetSpecificActivity(Int_t i) const
1040 {
1041  if (i>=0 && i<fNelements) return fWeights[i]*GetElement(i)->GetSpecificActivity();
1042  Double_t sa = 0;
1043  for (Int_t iel=0; iel<fNelements; iel++) {
1044  sa += fWeights[iel]*GetElement(iel)->GetSpecificActivity();
1045  }
1046  return sa;
1047 }
1048 
1049 ////////////////////////////////////////////////////////////////////////////////
1050 /// Return true if the other material has the same physical properties
1051 
1052 Bool_t TGeoMixture::IsEq(const TGeoMaterial *other) const
1053 {
1054  if (other->IsEqual(this)) return kTRUE;
1055  if (!other->IsMixture()) return kFALSE;
1056  TGeoMixture *mix = (TGeoMixture*)other;
1057  if (!mix) return kFALSE;
1058  if (fNelements != mix->GetNelements()) return kFALSE;
1059  if (TMath::Abs(fA-other->GetA())>1E-3) return kFALSE;
1060  if (TMath::Abs(fZ-other->GetZ())>1E-3) return kFALSE;
1061  if (TMath::Abs(fDensity-other->GetDensity())>1E-6) return kFALSE;
1062  if (GetCerenkovProperties() != other->GetCerenkovProperties()) return kFALSE;
1063 // if (fRadLen != other->GetRadLen()) return kFALSE;
1064 // if (fIntLen != other->GetIntLen()) return kFALSE;
1065  for (Int_t i=0; i<fNelements; i++) {
1066  if (TMath::Abs(fZmixture[i]-(mix->GetZmixt())[i])>1E-3) return kFALSE;
1067  if (TMath::Abs(fAmixture[i]-(mix->GetAmixt())[i])>1E-3) return kFALSE;
1068  if (TMath::Abs(fWeights[i]-(mix->GetWmixt())[i])>1E-3) return kFALSE;
1069  }
1070  return kTRUE;
1071 }
1072 
1073 ////////////////////////////////////////////////////////////////////////////////
1074 /// print characteristics of this material
1075 
1076 void TGeoMixture::Print(const Option_t * /*option*/) const
1077 {
1078  printf("Mixture %s %s Aeff=%g Zeff=%g rho=%g radlen=%g intlen=%g index=%i\n", GetName(), GetTitle(),
1079  fA,fZ,fDensity, fRadLen, fIntLen, fIndex);
1080  for (Int_t i=0; i<fNelements; i++) {
1081  if (fNatoms) printf(" Element #%i : %s Z=%6.2f A=%6.2f w=%6.3f natoms=%d\n", i, GetElement(i)->GetName(),fZmixture[i],
1082  fAmixture[i], fWeights[i], fNatoms[i]);
1083  else printf(" Element #%i : %s Z=%6.2f A=%6.2f w=%6.3f\n", i, GetElement(i)->GetName(),fZmixture[i],
1084  fAmixture[i], fWeights[i]);
1085  }
1086 }
1087 
1088 ////////////////////////////////////////////////////////////////////////////////
1089 /// Save a primitive as a C++ statement(s) on output stream "out".
1090 
1091 void TGeoMixture::SavePrimitive(std::ostream &out, Option_t * /*option*/ /*= ""*/)
1092 {
1093  if (TestBit(TGeoMaterial::kMatSavePrimitive)) return;
1094  char *name = GetPointerName();
1095  out << "// Mixture: " << GetName() << std::endl;
1096  out << " nel = " << fNelements << ";" << std::endl;
1097  out << " density = " << fDensity << ";" << std::endl;
1098  out << " " << name << " = new TGeoMixture(\"" << GetName() << "\", nel,density);" << std::endl;
1099  for (Int_t i=0; i<fNelements; i++) {
1100  TGeoElement *el = GetElement(i);
1101  out << " a = " << fAmixture[i] << "; z = "<< fZmixture[i] << "; w = " << fWeights[i] << "; // " << el->GetName() << std::endl;
1102  out << " " << name << "->DefineElement(" << i << ",a,z,w);" << std::endl;
1103  }
1104  out << " " << name << "->SetIndex(" << GetIndex() << ");" << std::endl;
1105  SetBit(TGeoMaterial::kMatSavePrimitive);
1106 }
1107 
1108 ////////////////////////////////////////////////////////////////////////////////
1109 /// Create the mixture representing the decay product of this material at a
1110 /// given time. The precision represent the minimum cumulative branching ratio for
1111 /// which decay products are still taken into account.
1112 
1113 TGeoMaterial *TGeoMixture::DecayMaterial(Double_t time, Double_t precision)
1114 {
1115  TObjArray *pop = new TObjArray();
1116  FillMaterialEvolution(pop, precision);
1117  Int_t ncomp = pop->GetEntriesFast();
1118  if (!ncomp) return this;
1119  TGeoElement *elem;
1120  TGeoElementRN *el;
1121  Double_t *weight = new Double_t[ncomp];
1122  Double_t amed = 0.;
1123  Int_t i, j;
1124  for (i=0; i<ncomp; i++) {
1125  elem = (TGeoElement *)pop->At(i);
1126  if (!elem->IsRadioNuclide()) {
1127  j = fElements->IndexOf(elem);
1128  weight[i] = fWeights[j]*fAmixture[0]/fWeights[0];
1129  } else {
1130  el = (TGeoElementRN*)elem;
1131  weight[i] = el->Ratio()->Concentration(time) * el->A();
1132  }
1133  amed += weight[i];
1134  }
1135  Double_t rho = fDensity * fWeights[0] * amed/fAmixture[0];
1136  TGeoMixture *mix = 0;
1137  Int_t ncomp1 = ncomp;
1138  for (i=0; i<ncomp; i++) {
1139  if ((weight[i]/amed)<precision) {
1140  amed -= weight[i];
1141  ncomp1--;
1142  }
1143  }
1144  if (ncomp1<2) {
1145  el = (TGeoElementRN *)pop->At(0);
1146  delete [] weight;
1147  delete pop;
1148  if (ncomp1==1) return new TGeoMaterial(TString::Format("%s-evol",GetName()), el, rho);
1149  return NULL;
1150  }
1151  mix = new TGeoMixture(TString::Format("%s-evol",GetName()), ncomp, rho);
1152  for (i=0; i<ncomp; i++) {
1153  weight[i] /= amed;
1154  if (weight[i]<precision) continue;
1155  el = (TGeoElementRN *)pop->At(i);
1156  mix->AddElement(el, weight[i]);
1157  }
1158  delete [] weight;
1159  delete pop;
1160  return mix;
1161 }
1162 
1163 ////////////////////////////////////////////////////////////////////////////////
1164 /// Fills a user array with all the elements deriving from the possible
1165 /// decay of the top elements composing the mixture. Each element contained
1166 /// by <population> may be a radionuclide having a Bateman solution attached.
1167 /// The precision represent the minimum cumulative branching ratio for
1168 /// which decay products are still taken into account.
1169 /// To visualize the time evolution of each decay product one can use:
1170 /// ~~~ {.cpp}
1171 /// TGeoElement *elem = population->At(index);
1172 /// TGeoElementRN *elemrn = 0;
1173 /// if (elem->IsRadioNuclide()) elemrn = (TGeoElementRN*)elem;
1174 /// ~~~
1175 /// One can get Ni/N1(t=0) at any moment of time. Ni is the number of atoms
1176 /// of one of the decay products, N1(0) is the number of atoms of the first top
1177 /// element at t=0.
1178 /// ~~~ {.cpp}
1179 /// Double_t fraction_weight = elemrn->Ratio()->Concentration(time);
1180 /// ~~~
1181 /// One can also display the time evolution of the fractional weight:
1182 /// ~~~ {.cpp}
1183 /// elemrn->Ratio()->Draw(option);
1184 /// ~~~
1185 
1186 void TGeoMixture::FillMaterialEvolution(TObjArray *population, Double_t precision)
1187 {
1188  if (population->GetEntriesFast()) {
1189  Error("FillMaterialEvolution", "Provide an empty array !");
1190  return;
1191  }
1192  TGeoElementTable *table = gGeoManager->GetElementTable();
1193  TGeoElement *elem;
1194  TGeoElementRN *elemrn;
1195  TIter next(table->GetElementsRN());
1196  while ((elemrn=(TGeoElementRN*)next())) elemrn->ResetRatio();
1197  Double_t factor;
1198  for (Int_t i=0; i<fNelements; i++) {
1199  elem = GetElement(i);
1200  if (!elem->IsRadioNuclide()) {
1201  population->Add(elem);
1202  continue;
1203  }
1204  elemrn = (TGeoElementRN*)elem;
1205  factor = fWeights[i]*fAmixture[0]/(fWeights[0]*fAmixture[i]);
1206  elemrn->FillPopulation(population, precision, factor);
1207  }
1208 }
1209 
1210 ////////////////////////////////////////////////////////////////////////////////
1211 /// static function
1212 /// Compute screening factor for pair production and Bremsstrahlung
1213 /// REFERENCE : EGS MANUAL SLAC 210 - UC32 - JUNE 78
1214 /// FORMULA 2.7.22
1215 
1216 Double_t TGeoMaterial::ScreenFactor(Double_t z)
1217 {
1218  const Double_t al183= 5.20948 , al1440 = 7.27239;
1219  Double_t alz = TMath::Log(z)/3.;
1220  Double_t factor = (al1440 - 2*alz) / (al183 - alz - TGeoMaterial::Coulomb(z));
1221  return factor;
1222 }
1223 
1224 ////////////////////////////////////////////////////////////////////////////////
1225 /// Compute Derived Quantities as in Geant4
1226 
1227 void TGeoMixture::ComputeDerivedQuantities()
1228 {
1229  const Double_t Na = (TGeoUnit::unitType()==TGeoUnit::kTGeoUnits)
1230  ? TGeoUnit::Avogadro : TGeant4Unit::Avogadro;
1231 
1232  if ( fVecNbOfAtomsPerVolume ) delete [] fVecNbOfAtomsPerVolume;
1233 
1234  fVecNbOfAtomsPerVolume = new Double_t[fNelements];
1235 
1236  // Formula taken from G4Material.cxx L312
1237  for (Int_t i=0; i<fNelements; ++i) {
1238  fVecNbOfAtomsPerVolume[i] = Na*fDensity*fWeights[i]/((TGeoElement*)fElements->At(i))->A();
1239  }
1240  ComputeRadiationLength();
1241  ComputeNuclearInterLength();
1242 }
1243 
1244 
1245 ////////////////////////////////////////////////////////////////////////////////
1246 /// Compute Radiation Length based on Geant4 formula
1247 
1248 void TGeoMixture::ComputeRadiationLength()
1249 {
1250  // Formula taken from G4Material.cxx L556
1251  const Double_t cm = (TGeoUnit::unitType()==TGeoUnit::kTGeoUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
1252  Double_t radinv = 0.0 ;
1253  for (Int_t i=0;i<fNelements;++i) {
1254  radinv += fVecNbOfAtomsPerVolume[i]*((TGeoElement*)fElements->At(i))->GetfRadTsai();
1255  }
1256  fRadLen = (radinv <= 0.0 ? DBL_MAX : cm/radinv);
1257 }
1258 
1259 ////////////////////////////////////////////////////////////////////////////////
1260 /// Compute Nuclear Interaction Length based on Geant4 formula
1261 void TGeoMixture::ComputeNuclearInterLength()
1262 {
1263  // Formula taken from G4Material.cxx L567
1264  TGeoUnit::UnitType typ = TGeoUnit::unitType();
1265  const Double_t g = (typ==TGeoUnit::kTGeoUnits) ? TGeoUnit::g : TGeant4Unit::g;
1266  const Double_t cm = (typ==TGeoUnit::kTGeoUnits) ? TGeoUnit::cm : TGeant4Unit::cm;
1267  const Double_t amu = (typ==TGeoUnit::kTGeoUnits) ? TGeoUnit::amu : TGeant4Unit::amu;
1268  const Double_t lambda0 = 35*g/(cm*cm);
1269  const Double_t twothird = 2.0/3.0;
1270  Double_t NILinv = 0.0;
1271  for (Int_t i=0; i<fNelements; ++i) {
1272  Int_t Z = static_cast<Int_t>(((TGeoElement*)fElements->At(i))->Z()+0.5);
1273  Double_t A = ((TGeoElement*)fElements->At(i))->Neff();
1274  if(1 == Z) {
1275  NILinv += fVecNbOfAtomsPerVolume[i]*A;
1276  } else {
1277  NILinv += fVecNbOfAtomsPerVolume[i]*TMath::Exp(twothird*TMath::Log(A));
1278  }
1279  }
1280  NILinv *= amu/lambda0;
1281  fIntLen = (NILinv <= 0.0 ? DBL_MAX : cm/NILinv);
1282 }