Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TGeoElement.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 17/06/04
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 TGeoElement
13 \ingroup Geometry_classes
14 Base class for chemical elements
15 */
16 
17 /** \class TGeoElementRN
18 \ingroup Geometry_classes
19 Class representing a radionuclide
20 */
21 
22 /** \class TGeoElemIter
23 \ingroup Geometry_classes
24 Iterator for decay branches
25 */
26 
27 /** \class TGeoDecayChannel
28 \ingroup Geometry_classes
29 A decay channel for a radionuclide
30 */
31 
32 /** \class TGeoElementTable
33 \ingroup Geometry_classes
34 Table of elements
35 */
36 
37 #include "RConfigure.h"
38 
39 #include "Riostream.h"
40 
41 #include "TSystem.h"
42 #include "TROOT.h"
43 #include "TObjArray.h"
44 #include "TVirtualGeoPainter.h"
45 #include "TGeoManager.h"
46 #include "TGeoElement.h"
47 #include "TMath.h"
48 #include "TGeoPhysicalConstants.h"
49 #include "TGeant4PhysicalConstants.h"
50 
51 // statics and globals
52 static const Int_t gMaxElem = 110;
53 static const Int_t gMaxLevel = 8;
54 static const Int_t gMaxDecay = 15;
55 
56 static const char gElName[gMaxElem][3] = {
57  "H ","He","Li","Be","B ","C ","N ","O ","F ","Ne","Na","Mg",
58  "Al","Si","P ","S ","Cl","Ar","K ","Ca","Sc","Ti","V ","Cr",
59  "Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr",
60  "Rb","Sr","Y ","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd",
61  "In","Sn","Sb","Te","I ","Xe","Cs","Ba","La","Ce","Pr","Nd",
62  "Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf",
63  "Ta","W ","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po",
64  "At","Rn","Fr","Ra","Ac","Th","Pa","U ","Np","Pu","Am","Cm",
65  "Bk","Cf","Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs",
66  "Mt","Ds" };
67 
68 static const char *gDecayName[gMaxDecay+1] = {
69  "2BetaMinus", "BetaMinus", "NeutronEm", "ProtonEm", "Alpha", "ECF",
70  "ElecCapt", "IsoTrans", "I", "SpontFiss", "2ProtonEm", "2NeutronEm",
71  "2Alpha", "Carbon12", "Carbon14", "Stable" };
72 
73 static const Int_t gDecayDeltaA[gMaxDecay] = {
74  0, 0, -1, -1, -4,
75  -99, 0, 0, -99, -99,
76  -2, -2, -8, -12, -14 };
77 
78 static const Int_t gDecayDeltaZ[gMaxDecay] = {
79  2, 1, 0, -1, -2,
80  -99, -1, 0, -99, -99,
81  -2, 0, -4, -6, -6 };
82 static const char gLevName[gMaxLevel]=" mnpqrs";
83 
84 ClassImp(TGeoElement);
85 
86 ////////////////////////////////////////////////////////////////////////////////
87 /// Default constructor
88 
89 TGeoElement::TGeoElement()
90 {
91  TGeoUnit::setUnitType(TGeoUnit::unitType()); // Ensure nobody changes the units afterwards
92  SetDefined(kFALSE);
93  SetUsed(kFALSE);
94  fZ = 0;
95  fN = 0;
96  fNisotopes = 0;
97  fA = 0.0;
98  fIsotopes = NULL;
99  fAbundances = NULL;
100 }
101 
102 ////////////////////////////////////////////////////////////////////////////////
103 /// Obsolete constructor
104 
105 TGeoElement::TGeoElement(const char *name, const char *title, Int_t z, Double_t a)
106  :TNamed(name, title)
107 {
108  TGeoUnit::setUnitType(TGeoUnit::unitType()); // Ensure nobody changes the units afterwards
109  SetDefined(kFALSE);
110  SetUsed(kFALSE);
111  fZ = z;
112  fN = Int_t(a);
113  fNisotopes = 0;
114  fA = a;
115  fIsotopes = NULL;
116  fAbundances = NULL;
117  ComputeDerivedQuantities();
118 }
119 
120 ////////////////////////////////////////////////////////////////////////////////
121 /// Element having isotopes.
122 
123 TGeoElement::TGeoElement(const char *name, const char *title, Int_t nisotopes)
124  :TNamed(name, title)
125 {
126  TGeoUnit::setUnitType(TGeoUnit::unitType()); // Ensure nobody changes the units afterwards
127  SetDefined(kFALSE);
128  SetUsed(kFALSE);
129  fZ = 0;
130  fN = 0;
131  fNisotopes = nisotopes;
132  fA = 0.0;
133  fIsotopes = new TObjArray(nisotopes);
134  fAbundances = new Double_t[nisotopes];
135 }
136 
137 ////////////////////////////////////////////////////////////////////////////////
138 /// Constructor
139 
140 TGeoElement::TGeoElement(const char *name, const char *title, Int_t z, Int_t n, Double_t a)
141  :TNamed(name, title)
142 {
143  TGeoUnit::setUnitType(TGeoUnit::unitType()); // Ensure nobody changes the units afterwards
144  SetDefined(kFALSE);
145  SetUsed(kFALSE);
146  fZ = z;
147  fN = n;
148  fNisotopes = 0;
149  fA = a;
150  fIsotopes = NULL;
151  fAbundances = NULL;
152  ComputeDerivedQuantities();
153 }
154 ////////////////////////////////////////////////////////////////////////////////
155 /// Calculate properties for an atomic number
156 
157 void TGeoElement::ComputeDerivedQuantities()
158 {
159  // Radiation Length
160  ComputeCoulombFactor();
161  ComputeLradTsaiFactor();
162 }
163 ////////////////////////////////////////////////////////////////////////////////
164 /// Compute Coulomb correction factor (Phys Rev. D50 3-1 (1994) page 1254)
165 
166 void TGeoElement::ComputeCoulombFactor()
167 {
168  static constexpr Double_t k1 = 0.0083 , k2 = 0.20206 ,k3 = 0.0020 , k4 = 0.0369;
169  Double_t fsc = TGeoUnit::unitType() == TGeoUnit::kTGeoUnits
170  ? TGeoUnit::fine_structure_const : TGeant4Unit::fine_structure_const;
171  Double_t az2 = (fsc*fZ)*(fsc*fZ);
172  Double_t az4 = az2 * az2;
173 
174  fCoulomb = (k1*az4 + k2 + 1./(1.+az2))*az2 - (k3*az4 + k4)*az4;
175 }
176 ////////////////////////////////////////////////////////////////////////////////
177 /// Compute Tsai's Expression for the Radiation Length (Phys Rev. D50 3-1 (1994) page 1254)
178 
179 void TGeoElement::ComputeLradTsaiFactor()
180 {
181  static constexpr Double_t Lrad_light[] = {5.31 , 4.79 , 4.74 , 4.71} ;
182  static constexpr Double_t Lprad_light[] = {6.144 , 5.621 , 5.805 , 5.924} ;
183 
184  fRadTsai = 0.0;
185  if (fZ == 0) return;
186  const Double_t logZ3 = TMath::Log(fZ)/3.;
187 
188  Double_t Lrad, Lprad;
189  Double_t alpha_rcl2 = TGeoUnit::unitType() == TGeoUnit::kTGeoUnits
190  ? TGeoUnit::alpha_rcl2 : TGeant4Unit::alpha_rcl2;
191  Int_t iz = static_cast<Int_t>(fZ+0.5) - 1 ; // The static cast comes from G4lrint
192  static const Double_t log184 = TMath::Log(184.15);
193  static const Double_t log1194 = TMath::Log(1194.);
194  if (iz <= 3) { Lrad = Lrad_light[iz] ; Lprad = Lprad_light[iz] ; }
195  else { Lrad = log184 - logZ3 ; Lprad = log1194 - 2*logZ3;}
196 
197  fRadTsai = 4*alpha_rcl2*fZ*(fZ*(Lrad-fCoulomb) + Lprad);
198 }
199 ////////////////////////////////////////////////////////////////////////////////
200 /// Print this isotope
201 
202 void TGeoElement::Print(Option_t *option) const
203 {
204  printf("Element: %s Z=%d N=%f A=%f [g/mole]\n", GetName(), fZ,Neff(),fA);
205  if (HasIsotopes()) {
206  for (Int_t i=0; i<fNisotopes; i++) {
207  TGeoIsotope *iso = GetIsotope(i);
208  printf("=>Isotope %s, abundance=%f :\n", iso->GetName(), fAbundances[i]);
209  iso->Print(option);
210  }
211  }
212 }
213 
214 ////////////////////////////////////////////////////////////////////////////////
215 /// Returns pointer to the table.
216 
217 TGeoElementTable *TGeoElement::GetElementTable()
218 {
219  if (!gGeoManager) {
220  ::Error("TGeoElementTable::GetElementTable", "Create a geometry manager first");
221  return NULL;
222  }
223  return gGeoManager->GetElementTable();
224 }
225 
226 ////////////////////////////////////////////////////////////////////////////////
227 /// Add an isotope for this element. All isotopes have to be isotopes of the same element.
228 
229 void TGeoElement::AddIsotope(TGeoIsotope *isotope, Double_t relativeAbundance)
230 {
231  if (!fIsotopes) {
232  Fatal("AddIsotope", "Cannot add isotopes to normal elements. Use constructor with number of isotopes.");
233  return;
234  }
235  Int_t ncurrent = 0;
236  TGeoIsotope *isocrt;
237  for (ncurrent=0; ncurrent<fNisotopes; ncurrent++)
238  if (!fIsotopes->At(ncurrent)) break;
239  if (ncurrent==fNisotopes) {
240  Error("AddIsotope", "All %d isotopes of element %s already defined", fNisotopes, GetName());
241  return;
242  }
243  // Check Z of the new isotope
244  if ((fZ!=0) && (isotope->GetZ()!=fZ)) {
245  Fatal("AddIsotope", "Trying to add isotope %s with different Z to the same element %s",
246  isotope->GetName(), GetName());
247  return;
248  } else {
249  fZ = isotope->GetZ();
250  }
251  fIsotopes->Add(isotope);
252  fAbundances[ncurrent] = relativeAbundance;
253  if (ncurrent==fNisotopes-1) {
254  Double_t weight = 0.0;
255  Double_t aeff = 0.0;
256  Double_t neff = 0.0;
257  for (Int_t i=0; i<fNisotopes; i++) {
258  isocrt = (TGeoIsotope*)fIsotopes->At(i);
259  aeff += fAbundances[i]*isocrt->GetA();
260  neff += fAbundances[i]*isocrt->GetN();
261  weight += fAbundances[i];
262  }
263  aeff /= weight;
264  neff /= weight;
265  fN = (Int_t)neff;
266  fA = aeff;
267  }
268  ComputeDerivedQuantities();
269 }
270 
271 ////////////////////////////////////////////////////////////////////////////////
272 /// Returns effective number of nucleons.
273 
274 Double_t TGeoElement::Neff() const
275 {
276  if (!fNisotopes) return fN;
277  TGeoIsotope *isocrt;
278  Double_t weight = 0.0;
279  Double_t neff = 0.0;
280  for (Int_t i=0; i<fNisotopes; i++) {
281  isocrt = (TGeoIsotope*)fIsotopes->At(i);
282  neff += fAbundances[i]*isocrt->GetN();
283  weight += fAbundances[i];
284  }
285  neff /= weight;
286  return neff;
287 }
288 
289 ////////////////////////////////////////////////////////////////////////////////
290 /// Return i-th isotope in the element.
291 
292 TGeoIsotope *TGeoElement::GetIsotope(Int_t i) const
293 {
294  if (i>=0 && i<fNisotopes) {
295  return (TGeoIsotope*)fIsotopes->At(i);
296  }
297  return NULL;
298 }
299 
300 ////////////////////////////////////////////////////////////////////////////////
301 /// Return relative abundance of i-th isotope in this element.
302 
303 Double_t TGeoElement::GetRelativeAbundance(Int_t i) const
304 {
305  if (i>=0 && i<fNisotopes) return fAbundances[i];
306  return 0.0;
307 }
308 
309 ClassImp(TGeoIsotope);
310 
311 ////////////////////////////////////////////////////////////////////////////////
312 /// Dummy I/O constructor
313 
314 TGeoIsotope::TGeoIsotope()
315  :TNamed(),
316  fZ(0),
317  fN(0),
318  fA(0)
319 {
320 }
321 
322 ////////////////////////////////////////////////////////////////////////////////
323 /// Constructor
324 
325 TGeoIsotope::TGeoIsotope(const char *name, Int_t z, Int_t n, Double_t a)
326  :TNamed(name,""),
327  fZ(z),
328  fN(n),
329  fA(a)
330 {
331  if (z<1) Fatal("ctor", "Not allowed Z=%d (<1) for isotope: %s", z,name);
332  if (n<z) Fatal("ctor", "Not allowed Z=%d < N=%d for isotope: %s", z,n,name);
333  TGeoElement::GetElementTable()->AddIsotope(this);
334 }
335 
336 ////////////////////////////////////////////////////////////////////////////////
337 /// Find existing isotope by name.
338 
339 TGeoIsotope *TGeoIsotope::FindIsotope(const char *name)
340 {
341  TGeoElementTable *elTable = TGeoElement::GetElementTable();
342  if (!elTable) return 0;
343  return elTable->FindIsotope(name);
344 }
345 
346 ////////////////////////////////////////////////////////////////////////////////
347 /// Print this isotope
348 
349 void TGeoIsotope::Print(Option_t *) const
350 {
351  printf("Isotope: %s Z=%d N=%d A=%f [g/mole]\n", GetName(), fZ,fN,fA);
352 }
353 
354 ClassImp(TGeoElementRN);
355 
356 ////////////////////////////////////////////////////////////////////////////////
357 /// Default constructor
358 
359 TGeoElementRN::TGeoElementRN()
360 {
361  TObject::SetBit(kElementChecked,kFALSE);
362  fENDFcode = 0;
363  fIso = 0;
364  fLevel = 0;
365  fDeltaM = 0;
366  fHalfLife = 0;
367  fNatAbun = 0;
368  fTH_F = 0;
369  fTG_F = 0;
370  fTH_S = 0;
371  fTG_S = 0;
372  fStatus = 0;
373  fRatio = 0;
374  fDecays = 0;
375 }
376 
377 ////////////////////////////////////////////////////////////////////////////////
378 /// Constructor.
379 
380 TGeoElementRN::TGeoElementRN(Int_t aa, Int_t zz, Int_t iso, Double_t level,
381  Double_t deltaM, Double_t halfLife, const char* JP,
382  Double_t natAbun, Double_t th_f, Double_t tg_f, Double_t th_s,
383  Double_t tg_s, Int_t status)
384  :TGeoElement("", JP, zz, aa)
385 {
386  TObject::SetBit(kElementChecked,kFALSE);
387  fENDFcode = ENDF(aa,zz,iso);
388  fIso = iso;
389  fLevel = level;
390  fDeltaM = deltaM;
391  fHalfLife = halfLife;
392  fTitle = JP;
393  if (!fTitle.Length()) fTitle = "?";
394  fNatAbun = natAbun;
395  fTH_F = th_f;
396  fTG_F = tg_f;
397  fTH_S = th_s;
398  fTG_S = tg_s;
399  fStatus = status;
400  fDecays = 0;
401  fRatio = 0;
402  MakeName(aa,zz,iso);
403  if ((TMath::Abs(fHalfLife)<1.e-30) || fHalfLife<-1) Warning("ctor","Element %s has T1/2=%g [s]", fName.Data(), fHalfLife);
404 }
405 
406 ////////////////////////////////////////////////////////////////////////////////
407 /// Destructor.
408 
409 TGeoElementRN::~TGeoElementRN()
410 {
411  if (fDecays) {
412  fDecays->Delete();
413  delete fDecays;
414  }
415  if (fRatio) delete fRatio;
416 }
417 
418 ////////////////////////////////////////////////////////////////////////////////
419 /// Adds a decay mode for this element.
420 
421 void TGeoElementRN::AddDecay(Int_t decay, Int_t diso, Double_t branchingRatio, Double_t qValue)
422 {
423  if (branchingRatio<1E-20) {
424  TString decayName;
425  TGeoDecayChannel::DecayName(decay, decayName);
426  Warning("AddDecay", "Decay %s of %s has BR=0. Not added.", decayName.Data(),fName.Data());
427  return;
428  }
429  TGeoDecayChannel *dc = new TGeoDecayChannel(decay,diso,branchingRatio, qValue);
430  dc->SetParent(this);
431  if (!fDecays) fDecays = new TObjArray(5);
432  fDecays->Add(dc);
433 }
434 
435 ////////////////////////////////////////////////////////////////////////////////
436 /// Adds a decay channel to the list of decays.
437 
438 void TGeoElementRN::AddDecay(TGeoDecayChannel *dc)
439 {
440  dc->SetParent(this);
441  if (!fDecays) fDecays = new TObjArray(5);
442  fDecays->Add(dc);
443 }
444 
445 ////////////////////////////////////////////////////////////////////////////////
446 /// Get number of decay channels of this element.
447 
448 Int_t TGeoElementRN::GetNdecays() const
449 {
450  if (!fDecays) return 0;
451  return fDecays->GetEntriesFast();
452 }
453 
454 ////////////////////////////////////////////////////////////////////////////////
455 /// Get the activity in Bq of a gram of material made from this element.
456 
457 Double_t TGeoElementRN::GetSpecificActivity() const
458 {
459  static const Double_t ln2 = TMath::Log(2.);
460  Double_t sa = (fHalfLife>0 && fA>0)?(ln2*TMath::Na()/fHalfLife/fA):0.;
461  return sa;
462 }
463 
464 ////////////////////////////////////////////////////////////////////////////////
465 /// Check if all decay chain of the element is OK.
466 
467 Bool_t TGeoElementRN::CheckDecays() const
468 {
469  if (TObject::TestBit(kElementChecked)) return kTRUE;
470  TObject *oelem = (TObject*)this;
471  TGeoDecayChannel *dc;
472  TGeoElementRN *elem;
473  TGeoElementTable *table = GetElementTable();
474  TString decayName;
475  if (!table) {
476  Error("CheckDecays", "Element table not present");
477  return kFALSE;
478  }
479  Bool_t resultOK = kTRUE;
480  if (!fDecays) {
481  oelem->SetBit(kElementChecked,kTRUE);
482  return resultOK;
483  }
484  Double_t br = 0.;
485  Int_t decayResult = 0;
486  TIter next(fDecays);
487  while ((dc=(TGeoDecayChannel*)next())) {
488  br += dc->BranchingRatio();
489  decayResult = DecayResult(dc);
490  if (decayResult) {
491  elem = table->GetElementRN(decayResult);
492  if (!elem) {
493  TGeoDecayChannel::DecayName(dc->Decay(),decayName);
494  Error("CheckDecays", "Element after decay %s of %s not found in DB", decayName.Data(),fName.Data());
495  return kFALSE;
496  }
497  dc->SetDaughter(elem);
498  resultOK = elem->CheckDecays();
499  }
500  }
501  if (TMath::Abs(br-100) > 1.E-3) {
502  Warning("CheckDecays", "BR for decays of element %s sum-up = %f", fName.Data(), br);
503  resultOK = kFALSE;
504  }
505  oelem->SetBit(kElementChecked,kTRUE);
506  return resultOK;
507 }
508 
509 ////////////////////////////////////////////////////////////////////////////////
510 /// Returns ENDF code of decay result.
511 
512 Int_t TGeoElementRN::DecayResult(TGeoDecayChannel *dc) const
513 {
514  Int_t da, dz, diso;
515  dc->DecayShift(da, dz, diso);
516  if (da == -99 || dz == -99) return 0;
517  return ENDF(Int_t(fA)+da,fZ+dz,fIso+diso);
518 }
519 
520 ////////////////////////////////////////////////////////////////////////////////
521 /// Fills the input array with the set of RN elements resulting from the decay of
522 /// this one. All element in the list will contain the time evolution of their
523 /// proportion by number with respect to this element. The proportion can be
524 /// retrieved via the method TGeoElementRN::Ratio().
525 /// The precision represent the minimum cumulative branching ratio for
526 /// which decay products are still taken into account.
527 
528 void TGeoElementRN::FillPopulation(TObjArray *population, Double_t precision, Double_t factor)
529 {
530  TGeoElementRN *elem;
531  TGeoElemIter next(this, precision);
532  TGeoBatemanSol s(this);
533  s.Normalize(factor);
534  AddRatio(s);
535  if (!population->FindObject(this)) population->Add(this);
536  while ((elem=next())) {
537  TGeoBatemanSol ratio(next.GetBranch());
538  ratio.Normalize(factor);
539  elem->AddRatio(ratio);
540  if (!population->FindObject(elem)) population->Add(elem);
541  }
542 }
543 
544 ////////////////////////////////////////////////////////////////////////////////
545 /// Generate a default name for the element.
546 
547 void TGeoElementRN::MakeName(Int_t a, Int_t z, Int_t iso)
548 {
549  fName = "";
550  if (z==0 && a==1) {
551  fName = "neutron";
552  return;
553  }
554  if (z>=1 && z<= gMaxElem) fName += TString::Format("%3d-%s-",z,gElName[z-1]);
555  else fName = "?? -?? -";
556  if (a>=1 && a<=999) fName += TString::Format("%3.3d",a);
557  else fName += "??";
558  if (iso>0 && iso<gMaxLevel) fName += TString::Format("%c", gLevName[iso]);
559  fName.ReplaceAll(" ","");
560 }
561 
562 ////////////////////////////////////////////////////////////////////////////////
563 /// Print info about the element;
564 
565 void TGeoElementRN::Print(Option_t *option) const
566 {
567  printf("\n%-12s ",fName.Data());
568  printf("ENDF=%d; ",fENDFcode);
569  printf("A=%d; ",(Int_t)fA);
570  printf("Z=%d; ",fZ);
571  printf("Iso=%d; ",fIso);
572  printf("Level=%g[MeV]; ",fLevel);
573  printf("Dmass=%g[MeV]; ",fDeltaM);
574  if (fHalfLife>0) printf("Hlife=%g[s]\n",fHalfLife);
575  else printf("Hlife=INF\n");
576  printf("%13s"," ");
577  printf("J/P=%s; ",fTitle.Data());
578  printf("Abund=%g; ",fNatAbun);
579  printf("Htox=%g; ",fTH_F);
580  printf("Itox=%g; ",fTG_F);
581  printf("Stat=%d\n",fStatus);
582  if(!fDecays) return;
583  printf("Decay modes:\n");
584  TIter next(fDecays);
585  TGeoDecayChannel *dc;
586  while ((dc=(TGeoDecayChannel*)next())) dc->Print(option);
587 }
588 
589 ////////////////////////////////////////////////////////////////////////////////
590 /// Create element from line record.
591 
592 TGeoElementRN *TGeoElementRN::ReadElementRN(const char *line, Int_t &ndecays)
593 {
594  Int_t a,z,iso,status;
595  Double_t level, deltaM, halfLife, natAbun, th_f, tg_f, th_s, tg_s;
596  char name[20],jp[20];
597  sscanf(&line[0], "%s%d%d%d%lg%lg%lg%s%lg%lg%lg%lg%lg%d%d", name,&a,&z,&iso,&level,&deltaM,
598  &halfLife,jp,&natAbun,&th_f,&tg_f,&th_s,&tg_s,&status,&ndecays);
599  TGeoElementRN *elem = new TGeoElementRN(a,z,iso,level,deltaM,halfLife,
600  jp,natAbun,th_f,tg_f,th_s,tg_s,status);
601  return elem;
602 }
603 
604 ////////////////////////////////////////////////////////////////////////////////
605 /// Save primitive for RN elements.
606 
607 void TGeoElementRN::SavePrimitive(std::ostream &out, Option_t *option)
608 {
609  if (!strcmp(option,"h")) {
610  // print a header if requested
611  out << "#====================================================================================================================================" << std::endl;
612  out << "# Name A Z ISO LEV[MeV] DM[MeV] T1/2[s] J/P ABUND[%] HTOX ITOX HTOX ITOX STAT NDCY" << std::endl;
613  out << "#====================================================================================================================================" << std::endl;
614  }
615  out << std::setw(11) << fName.Data();
616  out << std::setw(5) << (Int_t)fA;
617  out << std::setw(5) << fZ;
618  out << std::setw(5) << fIso;
619  out << std::setw(10) << std::setiosflags(std::ios::fixed) << std::setprecision(5) << fLevel;
620  out << std::setw(10) << std::setiosflags(std::ios::fixed) << std::setprecision(5) << fDeltaM;
621  out << std::setw(10) << std::setiosflags(std::ios::scientific) << std::setprecision(3) << fHalfLife;
622  out << std::setw(13) << fTitle.Data();
623  out << std::setw(10) << std::setiosflags(std::ios::fixed) << std::setprecision(5) << fNatAbun;
624  out << std::setw(10) << std::setiosflags(std::ios::fixed) << std::setprecision(5) << fTH_F;
625  out << std::setw(10) << std::setiosflags(std::ios::fixed) << std::setprecision(5) << fTG_F;
626  out << std::setw(10) << std::setiosflags(std::ios::fixed) << std::setprecision(5) << fTH_S;
627  out << std::setw(10) << std::setiosflags(std::ios::fixed) << std::setprecision(5) << fTG_S;
628  out << std::setw(5) << fStatus;
629  Int_t ndecays = 0;
630  if (fDecays) ndecays = fDecays->GetEntries();
631  out << std::setw(5) << ndecays;
632  out << std::endl;
633  if (fDecays) {
634  TIter next(fDecays);
635  TGeoDecayChannel *dc;
636  while ((dc=(TGeoDecayChannel*)next())) dc->SavePrimitive(out);
637  }
638 }
639 
640 ////////////////////////////////////////////////////////////////////////////////
641 /// Adds a proportion ratio to the existing one.
642 
643 void TGeoElementRN::AddRatio(TGeoBatemanSol &ratio)
644 {
645  if (!fRatio) fRatio = new TGeoBatemanSol(ratio);
646  else *fRatio += ratio;
647 }
648 
649 ////////////////////////////////////////////////////////////////////////////////
650 /// Clears the existing ratio.
651 
652 void TGeoElementRN::ResetRatio()
653 {
654  if (fRatio) {
655  delete fRatio;
656  fRatio = 0;
657  }
658 }
659 
660 ClassImp(TGeoDecayChannel);
661 
662 ////////////////////////////////////////////////////////////////////////////////
663 /// Assignment.
664 ///assignment operator
665 
666 TGeoDecayChannel& TGeoDecayChannel::operator=(const TGeoDecayChannel& dc)
667 {
668  if(this!=&dc) {
669  TObject::operator=(dc);
670  fDecay = dc.fDecay;
671  fDiso = dc.fDiso;
672  fBranchingRatio = dc.fBranchingRatio;
673  fQvalue = dc.fQvalue;
674  fParent = dc.fParent;
675  fDaughter = dc.fDaughter;
676  }
677  return *this;
678 }
679 
680 ////////////////////////////////////////////////////////////////////////////////
681 /// Returns name of decay.
682 
683 const char *TGeoDecayChannel::GetName() const
684 {
685  static TString name = "";
686  name = "";
687  if (!fDecay) return gDecayName[gMaxDecay];
688  for (Int_t i=0; i<gMaxDecay; i++) {
689  if (1<<i & fDecay) {
690  if (name.Length()) name += "+";
691  name += gDecayName[i];
692  }
693  }
694  return name.Data();
695 }
696 
697 ////////////////////////////////////////////////////////////////////////////////
698 /// Returns decay name.
699 
700 void TGeoDecayChannel::DecayName(UInt_t decay, TString &name)
701 {
702  if (!decay) {
703  name = gDecayName[gMaxDecay];
704  return;
705  }
706  name = "";
707  for (Int_t i=0; i<gMaxDecay; i++) {
708  if (1<<i & decay) {
709  if (name.Length()) name += "+";
710  name += gDecayName[i];
711  }
712  }
713 }
714 
715 ////////////////////////////////////////////////////////////////////////////////
716 /// Get index of this channel in the list of decays of the parent nuclide.
717 
718 Int_t TGeoDecayChannel::GetIndex() const
719 {
720  return fParent->Decays()->IndexOf(this);
721 }
722 
723 ////////////////////////////////////////////////////////////////////////////////
724 /// Prints decay info.
725 
726 void TGeoDecayChannel::Print(Option_t *) const
727 {
728  TString name;
729  DecayName(fDecay, name);
730  printf("%-20s Diso: %3d BR: %9.3f%% Qval: %g\n", name.Data(),fDiso,fBranchingRatio,fQvalue);
731 }
732 
733 ////////////////////////////////////////////////////////////////////////////////
734 /// Create element from line record.
735 
736 TGeoDecayChannel *TGeoDecayChannel::ReadDecay(const char *line)
737 {
738  char name[80];
739  Int_t decay,diso;
740  Double_t branchingRatio, qValue;
741  sscanf(&line[0], "%s%d%d%lg%lg", name,&decay,&diso,&branchingRatio,&qValue);
742  TGeoDecayChannel *dc = new TGeoDecayChannel(decay,diso,branchingRatio,qValue);
743  return dc;
744 }
745 
746 ////////////////////////////////////////////////////////////////////////////////
747 /// Save primitive for decays.
748 
749 void TGeoDecayChannel::SavePrimitive(std::ostream &out, Option_t *)
750 {
751  TString decayName;
752  DecayName(fDecay, decayName);
753  out << std::setw(50) << decayName.Data();
754  out << std::setw(10) << fDecay;
755  out << std::setw(10) << fDiso;
756  out << std::setw(12) << std::setiosflags(std::ios::fixed) << std::setprecision(6) << fBranchingRatio;
757  out << std::setw(12) << std::setiosflags(std::ios::fixed) << std::setprecision(6) << fQvalue;
758  out << std::endl;
759 }
760 
761 ////////////////////////////////////////////////////////////////////////////////
762 /// Returns variation in A, Z and Iso after decay.
763 
764 void TGeoDecayChannel::DecayShift(Int_t &dA, Int_t &dZ, Int_t &dI) const
765 {
766  dA=dZ=0;
767  dI=fDiso;
768  for(Int_t i=0; i<gMaxDecay; ++i) {
769  if(1<<i & fDecay) {
770  if(gDecayDeltaA[i] == -99 || gDecayDeltaZ[i] == -99 ) {
771  dA=dZ=-99;
772  return;
773  }
774  dA += gDecayDeltaA[i];
775  dZ += gDecayDeltaZ[i];
776  }
777  }
778 }
779 
780 ClassImp(TGeoElemIter);
781 
782 ////////////////////////////////////////////////////////////////////////////////
783 /// Default constructor.
784 
785 TGeoElemIter::TGeoElemIter(TGeoElementRN *top, Double_t limit)
786  : fTop(top), fElem(top), fBranch(0), fLevel(0), fLimitRatio(limit), fRatio(1.)
787 {
788  fBranch = new TObjArray(10);
789 }
790 
791 ////////////////////////////////////////////////////////////////////////////////
792 /// Copy ctor.
793 
794 TGeoElemIter::TGeoElemIter(const TGeoElemIter &iter)
795  :fTop(iter.fTop),
796  fElem(iter.fElem),
797  fBranch(0),
798  fLevel(iter.fLevel),
799  fLimitRatio(iter.fLimitRatio),
800  fRatio(iter.fRatio)
801 {
802  if (iter.fBranch) {
803  fBranch = new TObjArray(10);
804  for (Int_t i=0; i<fLevel; i++) fBranch->Add(iter.fBranch->At(i));
805  }
806 }
807 
808 ////////////////////////////////////////////////////////////////////////////////
809 /// Destructor.
810 
811 TGeoElemIter::~TGeoElemIter()
812 {
813  if (fBranch) delete fBranch;
814 }
815 
816 ////////////////////////////////////////////////////////////////////////////////
817 /// Assignment.
818 
819 TGeoElemIter &TGeoElemIter::operator=(const TGeoElemIter &iter)
820 {
821  if (&iter == this) return *this;
822  fTop = iter.fTop;
823  fElem = iter.fElem;
824  fLevel = iter.fLevel;
825  if (iter.fBranch) {
826  fBranch = new TObjArray(10);
827  for (Int_t i=0; i<fLevel; i++) fBranch->Add(iter.fBranch->At(i));
828  }
829  fLimitRatio = iter.fLimitRatio;
830  fRatio = iter.fRatio;
831  return *this;
832 }
833 
834 ////////////////////////////////////////////////////////////////////////////////
835 /// () operator.
836 
837 TGeoElementRN *TGeoElemIter::operator()()
838 {
839  return Next();
840 }
841 
842 ////////////////////////////////////////////////////////////////////////////////
843 /// Go upwards from the current location until the next branching, then down.
844 
845 TGeoElementRN *TGeoElemIter::Up()
846 {
847  TGeoDecayChannel *dc;
848  Int_t ind, nd;
849  while (fLevel) {
850  // Current decay channel
851  dc = (TGeoDecayChannel*)fBranch->At(fLevel-1);
852  ind = dc->GetIndex();
853  nd = dc->Parent()->GetNdecays();
854  fRatio /= 0.01*dc->BranchingRatio();
855  fElem = dc->Parent();
856  fBranch->RemoveAt(--fLevel);
857  ind++;
858  while (ind<nd) {
859  if (Down(ind++)) return (TGeoElementRN*)fElem;
860  }
861  }
862  fElem = NULL;
863  return NULL;
864 }
865 
866 ////////////////////////////////////////////////////////////////////////////////
867 /// Go downwards from current level via ibranch as low in the tree as possible.
868 /// Return value flags if the operation was successful.
869 
870 TGeoElementRN *TGeoElemIter::Down(Int_t ibranch)
871 {
872  TGeoDecayChannel *dc = (TGeoDecayChannel*)fElem->Decays()->At(ibranch);
873  if (!dc->Daughter()) return NULL;
874  Double_t br = 0.01*fRatio*dc->BranchingRatio();
875  if (br < fLimitRatio) return NULL;
876  fLevel++;
877  fRatio = br;
878  fBranch->Add(dc);
879  fElem = dc->Daughter();
880  return (TGeoElementRN*)fElem;
881 }
882 
883 ////////////////////////////////////////////////////////////////////////////////
884 /// Return next element.
885 
886 TGeoElementRN *TGeoElemIter::Next()
887 {
888  if (!fElem) return NULL;
889  // Check if this is the first iteration.
890  Int_t nd = fElem->GetNdecays();
891  for (Int_t i=0; i<nd; i++) if (Down(i)) return (TGeoElementRN*)fElem;
892  return Up();
893 }
894 
895 ////////////////////////////////////////////////////////////////////////////////
896 /// Print info about the current decay branch.
897 
898 void TGeoElemIter::Print(Option_t * /*option*/) const
899 {
900  TGeoElementRN *elem;
901  TGeoDecayChannel *dc;
902  TString indent = "";
903  printf("=== Chain with %g %%\n", 100*fRatio);
904  for (Int_t i=0; i<fLevel; i++) {
905  dc = (TGeoDecayChannel*)fBranch->At(i);
906  elem = dc->Parent();
907  printf("%s%s (%g%% %s) T1/2=%g\n", indent.Data(), elem->GetName(),dc->BranchingRatio(),dc->GetName(),elem->HalfLife());
908  indent += " ";
909  if (i==fLevel-1) {
910  elem = dc->Daughter();
911  printf("%s%s\n", indent.Data(), elem->GetName());
912  }
913  }
914 }
915 
916 ClassImp(TGeoElementTable);
917 
918 ////////////////////////////////////////////////////////////////////////////////
919 /// default constructor
920 
921 TGeoElementTable::TGeoElementTable()
922 {
923  fNelements = 0;
924  fNelementsRN = 0;
925  fNisotopes = 0;
926  fList = 0;
927  fListRN = 0;
928  fIsotopes = 0;
929 }
930 
931 ////////////////////////////////////////////////////////////////////////////////
932 /// constructor
933 
934 TGeoElementTable::TGeoElementTable(Int_t /*nelements*/)
935 {
936  fNelements = 0;
937  fNelementsRN = 0;
938  fNisotopes = 0;
939  fList = new TObjArray(128);
940  fListRN = 0;
941  fIsotopes = 0;
942  BuildDefaultElements();
943 // BuildElementsRN();
944 }
945 
946 ////////////////////////////////////////////////////////////////////////////////
947 ///copy constructor
948 
949 TGeoElementTable::TGeoElementTable(const TGeoElementTable& get) :
950  TObject(get),
951  fNelements(get.fNelements),
952  fNelementsRN(get.fNelementsRN),
953  fNisotopes(get.fNisotopes),
954  fList(get.fList),
955  fListRN(get.fListRN),
956  fIsotopes(0)
957 {
958 }
959 
960 ////////////////////////////////////////////////////////////////////////////////
961 ///assignment operator
962 
963 TGeoElementTable& TGeoElementTable::operator=(const TGeoElementTable& get)
964 {
965  if(this!=&get) {
966  TObject::operator=(get);
967  fNelements=get.fNelements;
968  fNelementsRN=get.fNelementsRN;
969  fNisotopes=get.fNisotopes;
970  fList=get.fList;
971  fListRN=get.fListRN;
972  fIsotopes = 0;
973  }
974  return *this;
975 }
976 
977 ////////////////////////////////////////////////////////////////////////////////
978 /// destructor
979 
980 TGeoElementTable::~TGeoElementTable()
981 {
982  if (fList) {
983  fList->Delete();
984  delete fList;
985  }
986  if (fListRN) {
987  fListRN->Delete();
988  delete fListRN;
989  }
990  if (fIsotopes) {
991  fIsotopes->Delete();
992  delete fIsotopes;
993  }
994 }
995 
996 ////////////////////////////////////////////////////////////////////////////////
997 /// Add an element to the table. Obsolete.
998 
999 void TGeoElementTable::AddElement(const char *name, const char *title, Int_t z, Double_t a)
1000 {
1001  if (!fList) fList = new TObjArray(128);
1002  fList->AddAtAndExpand(new TGeoElement(name,title,z,a), fNelements++);
1003 }
1004 
1005 ////////////////////////////////////////////////////////////////////////////////
1006 /// Add an element to the table.
1007 
1008 void TGeoElementTable::AddElement(const char *name, const char *title, Int_t z, Int_t n, Double_t a)
1009 {
1010  if (!fList) fList = new TObjArray(128);
1011  fList->AddAtAndExpand(new TGeoElement(name,title,z,n,a), fNelements++);
1012 }
1013 
1014 ////////////////////////////////////////////////////////////////////////////////
1015 /// Add a custom element to the table.
1016 
1017 void TGeoElementTable::AddElement(TGeoElement *elem)
1018 {
1019  if (!fList) fList = new TObjArray(128);
1020  TGeoElement *orig = FindElement(elem->GetName());
1021  if (orig) {
1022  Error("AddElement", "Found element with same name: %s (%s). Cannot add to table.",
1023  orig->GetName(), orig->GetTitle());
1024  return;
1025  }
1026  fList->AddAtAndExpand(elem, fNelements++);
1027 }
1028 
1029 ////////////////////////////////////////////////////////////////////////////////
1030 /// Add a radionuclide to the table and map it.
1031 
1032 void TGeoElementTable::AddElementRN(TGeoElementRN *elem)
1033 {
1034  if (!fListRN) fListRN = new TObjArray(3600);
1035  if (HasRNElements() && GetElementRN(elem->ENDFCode())) return;
1036 // elem->Print();
1037  fListRN->Add(elem);
1038  fNelementsRN++;
1039  fElementsRN.insert(ElementRNMap_t::value_type(elem->ENDFCode(), elem));
1040 }
1041 
1042 ////////////////////////////////////////////////////////////////////////////////
1043 /// Add isotope to the table.
1044 
1045 void TGeoElementTable::AddIsotope(TGeoIsotope *isotope)
1046 {
1047  if (FindIsotope(isotope->GetName())) {
1048  Error("AddIsotope", "Isotope with the same name: %s already in table. Not adding.",isotope->GetName());
1049  return;
1050  }
1051  if (!fIsotopes) fIsotopes = new TObjArray();
1052  fIsotopes->Add(isotope);
1053 }
1054 
1055 ////////////////////////////////////////////////////////////////////////////////
1056 /// Creates the default element table
1057 
1058 void TGeoElementTable::BuildDefaultElements()
1059 {
1060  if (HasDefaultElements()) return;
1061  AddElement("VACUUM","VACUUM" ,0, 0, 0.0);
1062  AddElement("H" ,"HYDROGEN" ,1, 1, 1.00794);
1063  AddElement("HE" ,"HELIUM" ,2, 4, 4.002602);
1064  AddElement("LI" ,"LITHIUM" ,3, 7, 6.941);
1065  AddElement("BE" ,"BERYLLIUM" ,4, 9, 9.01218);
1066  AddElement("B" ,"BORON" ,5, 11, 10.811);
1067  AddElement("C" ,"CARBON" ,6, 12, 12.0107);
1068  AddElement("N" ,"NITROGEN" ,7, 14, 14.00674);
1069  AddElement("O" ,"OXYGEN" ,8, 16, 15.9994);
1070  AddElement("F" ,"FLUORINE" ,9, 19, 18.9984032);
1071  AddElement("NE" ,"NEON" ,10, 20, 20.1797);
1072  AddElement("NA" ,"SODIUM" ,11, 23, 22.989770);
1073  AddElement("MG" ,"MAGNESIUM" ,12, 24, 24.3050);
1074  AddElement("AL" ,"ALUMINIUM" ,13, 27, 26.981538);
1075  AddElement("SI" ,"SILICON" ,14, 28, 28.0855);
1076  AddElement("P" ,"PHOSPHORUS" ,15, 31, 30.973761);
1077  AddElement("S" ,"SULFUR" ,16, 32, 32.066);
1078  AddElement("CL" ,"CHLORINE" ,17, 35, 35.4527);
1079  AddElement("AR" ,"ARGON" ,18, 40, 39.948);
1080  AddElement("K" ,"POTASSIUM" ,19, 39, 39.0983);
1081  AddElement("CA" ,"CALCIUM" ,20, 40, 40.078);
1082  AddElement("SC" ,"SCANDIUM" ,21, 45, 44.955910);
1083  AddElement("TI" ,"TITANIUM" ,22, 48, 47.867);
1084  AddElement("V" ,"VANADIUM" ,23, 51, 50.9415);
1085  AddElement("CR" ,"CHROMIUM" ,24, 52, 51.9961);
1086  AddElement("MN" ,"MANGANESE" ,25, 55, 54.938049);
1087  AddElement("FE" ,"IRON" ,26, 56, 55.845);
1088  AddElement("CO" ,"COBALT" ,27, 59, 58.933200);
1089  AddElement("NI" ,"NICKEL" ,28, 59, 58.6934);
1090  AddElement("CU" ,"COPPER" ,29, 64, 63.546);
1091  AddElement("ZN" ,"ZINC" ,30, 65, 65.39);
1092  AddElement("GA" ,"GALLIUM" ,31, 70, 69.723);
1093  AddElement("GE" ,"GERMANIUM" ,32, 73, 72.61);
1094  AddElement("AS" ,"ARSENIC" ,33, 75, 74.92160);
1095  AddElement("SE" ,"SELENIUM" ,34, 79, 78.96);
1096  AddElement("BR" ,"BROMINE" ,35, 80, 79.904);
1097  AddElement("KR" ,"KRYPTON" ,36, 84, 83.80);
1098  AddElement("RB" ,"RUBIDIUM" ,37, 85, 85.4678);
1099  AddElement("SR" ,"STRONTIUM" ,38, 88, 87.62);
1100  AddElement("Y" ,"YTTRIUM" ,39, 89, 88.90585);
1101  AddElement("ZR" ,"ZIRCONIUM" ,40, 91, 91.224);
1102  AddElement("NB" ,"NIOBIUM" ,41, 93, 92.90638);
1103  AddElement("MO" ,"MOLYBDENUM" ,42, 96, 95.94);
1104  AddElement("TC" ,"TECHNETIUM" ,43, 98, 98.0);
1105  AddElement("RU" ,"RUTHENIUM" ,44, 101, 101.07);
1106  AddElement("RH" ,"RHODIUM" ,45, 103, 102.90550);
1107  AddElement("PD" ,"PALLADIUM" ,46, 106, 106.42);
1108  AddElement("AG" ,"SILVER" ,47, 108, 107.8682);
1109  AddElement("CD" ,"CADMIUM" ,48, 112, 112.411);
1110  AddElement("IN" ,"INDIUM" ,49, 115, 114.818);
1111  AddElement("SN" ,"TIN" ,50, 119, 118.710);
1112  AddElement("SB" ,"ANTIMONY" ,51, 122, 121.760);
1113  AddElement("TE" ,"TELLURIUM" ,52, 128, 127.60);
1114  AddElement("I" ,"IODINE" ,53, 127, 126.90447);
1115  AddElement("XE" ,"XENON" ,54, 131, 131.29);
1116  AddElement("CS" ,"CESIUM" ,55, 133, 132.90545);
1117  AddElement("BA" ,"BARIUM" ,56, 137, 137.327);
1118  AddElement("LA" ,"LANTHANUM" ,57, 139, 138.9055);
1119  AddElement("CE" ,"CERIUM" ,58, 140, 140.116);
1120  AddElement("PR" ,"PRASEODYMIUM" ,59, 141, 140.90765);
1121  AddElement("ND" ,"NEODYMIUM" ,60, 144, 144.24);
1122  AddElement("PM" ,"PROMETHIUM" ,61, 145, 145.0);
1123  AddElement("SM" ,"SAMARIUM" ,62, 150, 150.36);
1124  AddElement("EU" ,"EUROPIUM" ,63, 152, 151.964);
1125  AddElement("GD" ,"GADOLINIUM" ,64, 157, 157.25);
1126  AddElement("TB" ,"TERBIUM" ,65, 159, 158.92534);
1127  AddElement("DY" ,"DYSPROSIUM" ,66, 162, 162.50);
1128  AddElement("HO" ,"HOLMIUM" ,67, 165, 164.93032);
1129  AddElement("ER" ,"ERBIUM" ,68, 167, 167.26);
1130  AddElement("TM" ,"THULIUM" ,69, 169, 168.93421);
1131  AddElement("YB" ,"YTTERBIUM" ,70, 173, 173.04);
1132  AddElement("LU" ,"LUTETIUM" ,71, 175, 174.967);
1133  AddElement("HF" ,"HAFNIUM" ,72, 178, 178.49);
1134  AddElement("TA" ,"TANTALUM" ,73, 181, 180.9479);
1135  AddElement("W" ,"TUNGSTEN" ,74, 184, 183.84);
1136  AddElement("RE" ,"RHENIUM" ,75, 186, 186.207);
1137  AddElement("OS" ,"OSMIUM" ,76, 190, 190.23);
1138  AddElement("IR" ,"IRIDIUM" ,77, 192, 192.217);
1139  AddElement("PT" ,"PLATINUM" ,78, 195, 195.078);
1140  AddElement("AU" ,"GOLD" ,79, 197, 196.96655);
1141  AddElement("HG" ,"MERCURY" ,80, 200, 200.59);
1142  AddElement("TL" ,"THALLIUM" ,81, 204, 204.3833);
1143  AddElement("PB" ,"LEAD" ,82, 207, 207.2);
1144  AddElement("BI" ,"BISMUTH" ,83, 209, 208.98038);
1145  AddElement("PO" ,"POLONIUM" ,84, 209, 209.0);
1146  AddElement("AT" ,"ASTATINE" ,85, 210, 210.0);
1147  AddElement("RN" ,"RADON" ,86, 222, 222.0);
1148  AddElement("FR" ,"FRANCIUM" ,87, 223, 223.0);
1149  AddElement("RA" ,"RADIUM" ,88, 226, 226.0);
1150  AddElement("AC" ,"ACTINIUM" ,89, 227, 227.0);
1151  AddElement("TH" ,"THORIUM" ,90, 232, 232.0381);
1152  AddElement("PA" ,"PROTACTINIUM" ,91, 231, 231.03588);
1153  AddElement("U" ,"URANIUM" ,92, 238, 238.0289);
1154  AddElement("NP" ,"NEPTUNIUM" ,93, 237, 237.0);
1155  AddElement("PU" ,"PLUTONIUM" ,94, 244, 244.0);
1156  AddElement("AM" ,"AMERICIUM" ,95, 243, 243.0);
1157  AddElement("CM" ,"CURIUM" ,96, 247, 247.0);
1158  AddElement("BK" ,"BERKELIUM" ,97, 247, 247.0);
1159  AddElement("CF" ,"CALIFORNIUM",98, 251, 251.0);
1160  AddElement("ES" ,"EINSTEINIUM",99, 252, 252.0);
1161  AddElement("FM" ,"FERMIUM" ,100, 257, 257.0);
1162  AddElement("MD" ,"MENDELEVIUM",101, 258, 258.0);
1163  AddElement("NO" ,"NOBELIUM" ,102, 259, 259.0);
1164  AddElement("LR" ,"LAWRENCIUM" ,103, 262, 262.0);
1165  AddElement("RF" ,"RUTHERFORDIUM",104, 261, 261.0);
1166  AddElement("DB" ,"DUBNIUM" ,105, 262, 262.0);
1167  AddElement("SG" ,"SEABORGIUM" ,106, 263, 263.0);
1168  AddElement("BH" ,"BOHRIUM" ,107, 262, 262.0);
1169  AddElement("HS" ,"HASSIUM" ,108, 265, 265.0);
1170  AddElement("MT" ,"MEITNERIUM" ,109, 266, 266.0);
1171  AddElement("UUN" ,"UNUNNILIUM" ,110, 269, 269.0);
1172  AddElement("UUU" ,"UNUNUNIUM" ,111, 272, 272.0);
1173  AddElement("UUB" ,"UNUNBIUM" ,112, 277, 277.0);
1174 
1175  TObject::SetBit(kETDefaultElements,kTRUE);
1176 }
1177 
1178 ////////////////////////////////////////////////////////////////////////////////
1179 /// Creates the list of radionuclides.
1180 
1181 void TGeoElementTable::ImportElementsRN()
1182 {
1183  if (HasRNElements()) return;
1184  TGeoElementRN *elem;
1185  TString rnf = "RadioNuclides.txt";
1186  gSystem->PrependPathName(TROOT::GetEtcDir(), rnf);
1187  FILE *fp = fopen(rnf, "r");
1188  if (!fp) {
1189  Error("ImportElementsRN","File RadioNuclides.txt not found");
1190  return;
1191  }
1192  char line[150];
1193  Int_t ndecays = 0;
1194  Int_t i;
1195  while (fgets(&line[0],140,fp)) {
1196  if (line[0]=='#') continue;
1197  elem = TGeoElementRN::ReadElementRN(line, ndecays);
1198  for (i=0; i<ndecays; i++) {
1199  if (!fgets(&line[0],140,fp)) {
1200  Error("ImportElementsRN", "Error parsing RadioNuclides.txt file");
1201  fclose(fp);
1202  return;
1203  }
1204  TGeoDecayChannel *dc = TGeoDecayChannel::ReadDecay(line);
1205  elem->AddDecay(dc);
1206  }
1207  AddElementRN(elem);
1208 // elem->Print();
1209  }
1210  TObject::SetBit(kETRNElements,kTRUE);
1211  CheckTable();
1212  fclose(fp);
1213 }
1214 
1215 ////////////////////////////////////////////////////////////////////////////////
1216 /// Checks status of element table.
1217 
1218 Bool_t TGeoElementTable::CheckTable() const
1219 {
1220  if (!HasRNElements()) return HasDefaultElements();
1221  TGeoElementRN *elem;
1222  Bool_t result = kTRUE;
1223  TIter next(fListRN);
1224  while ((elem=(TGeoElementRN*)next())) {
1225  if (!elem->CheckDecays()) result = kFALSE;
1226  }
1227  return result;
1228 }
1229 
1230 ////////////////////////////////////////////////////////////////////////////////
1231 /// Export radionuclides in a file.
1232 
1233 void TGeoElementTable::ExportElementsRN(const char *filename)
1234 {
1235  if (!HasRNElements()) return;
1236  TString sname = filename;
1237  if (!sname.Length()) sname = "RadioNuclides.txt";
1238  std::ofstream out;
1239  out.open(sname.Data(), std::ios::out);
1240  if (!out.good()) {
1241  Error("ExportElementsRN", "Cannot open file %s", sname.Data());
1242  return;
1243  }
1244 
1245  TGeoElementRN *elem;
1246  TIter next(fListRN);
1247  Int_t i=0;
1248  while ((elem=(TGeoElementRN*)next())) {
1249  if ((i%48)==0) elem->SavePrimitive(out,"h");
1250  else elem->SavePrimitive(out);
1251  i++;
1252  }
1253  out.close();
1254 }
1255 
1256 ////////////////////////////////////////////////////////////////////////////////
1257 /// Search an element by symbol or full name
1258 /// Exact matching
1259 
1260 TGeoElement *TGeoElementTable::FindElement(const char *name) const
1261 {
1262  TGeoElement *elem;
1263  elem = (TGeoElement*)fList->FindObject(name);
1264  if (elem) return elem;
1265  // Search case insensitive by element name
1266  TString s(name);
1267  s.ToUpper();
1268  elem = (TGeoElement*)fList->FindObject(s.Data());
1269  if (elem) return elem;
1270  // Search by full name
1271  TIter next(fList);
1272  while ((elem=(TGeoElement*)next())) {
1273  if (s == elem->GetTitle()) return elem;
1274  }
1275  return 0;
1276 }
1277 
1278 ////////////////////////////////////////////////////////////////////////////////
1279 /// Find existing isotope by name. Not optimized for a big number of isotopes.
1280 
1281 TGeoIsotope *TGeoElementTable::FindIsotope(const char *name) const
1282 {
1283  if (!fIsotopes) return NULL;
1284  return (TGeoIsotope*)fIsotopes->FindObject(name);
1285 }
1286 
1287 ////////////////////////////////////////////////////////////////////////////////
1288 /// Retrieve a radionuclide by ENDF code.
1289 
1290 TGeoElementRN *TGeoElementTable::GetElementRN(Int_t ENDFcode) const
1291 {
1292  if (!HasRNElements()) {
1293  TGeoElementTable *table = (TGeoElementTable*)this;
1294  table->ImportElementsRN();
1295  if (!fListRN) return 0;
1296  }
1297  ElementRNMap_t::const_iterator it = fElementsRN.find(ENDFcode);
1298  if (it != fElementsRN.end()) return it->second;
1299  return 0;
1300 }
1301 
1302 ////////////////////////////////////////////////////////////////////////////////
1303 /// Retrieve a radionuclide by a, z, and isomeric state.
1304 
1305 TGeoElementRN *TGeoElementTable::GetElementRN(Int_t a, Int_t z, Int_t iso) const
1306 {
1307  return GetElementRN(TGeoElementRN::ENDF(a,z,iso));
1308 }
1309 
1310 ////////////////////////////////////////////////////////////////////////////////
1311 /// Print table of elements. The accepted options are:
1312 /// "" - prints everything by default
1313 /// "D" - prints default elements only
1314 /// "I" - prints isotopes
1315 /// "R" - prints radio-nuclides only if imported
1316 /// "U" - prints user-defined elements only
1317 
1318 void TGeoElementTable::Print(Option_t *option) const
1319 {
1320  TString opt(option);
1321  opt.ToUpper();
1322  Int_t induser = HasDefaultElements() ? 113 : 0;
1323  // Default elements
1324  if (opt=="" || opt=="D") {
1325  if (induser) printf("================\nDefault elements\n================\n");
1326  for (Int_t iel=0; iel<induser; ++iel) fList->At(iel)->Print();
1327  }
1328  // Isotopes
1329  if (opt=="" || opt=="I") {
1330  if (fIsotopes) {
1331  printf("================\nIsotopes\n================\n");
1332  fIsotopes->Print();
1333  }
1334  }
1335  // Radio-nuclides
1336  if (opt=="" || opt=="R") {
1337  if (HasRNElements()) {
1338  printf("================\nRadio-nuclides\n================\n");
1339  fListRN->Print();
1340  }
1341  }
1342  // User-defined elements
1343  if (opt=="" || opt=="U") {
1344  if (fNelements>induser) printf("================\nUser elements\n================\n");
1345  for (Int_t iel=induser; iel<fNelements; ++iel) fList->At(iel)->Print();
1346  }
1347 }
1348 
1349 ClassImp(TGeoBatemanSol);
1350 
1351 ////////////////////////////////////////////////////////////////////////////////
1352 /// Default ctor.
1353 
1354 TGeoBatemanSol::TGeoBatemanSol(TGeoElementRN *elem)
1355  :TObject(), TAttLine(), TAttFill(), TAttMarker(),
1356  fElem(elem),
1357  fElemTop(elem),
1358  fCsize(10),
1359  fNcoeff(0),
1360  fFactor(1.),
1361  fTmin(0.),
1362  fTmax(0.),
1363  fCoeff(NULL)
1364 {
1365  fCoeff = new BtCoef_t[fCsize];
1366  fNcoeff = 1;
1367  fCoeff[0].cn = 1.;
1368  Double_t t12 = elem->HalfLife();
1369  if (t12 == 0.) t12 = 1.e-30;
1370  if (elem->Stable()) fCoeff[0].lambda = 0.;
1371  else fCoeff[0].lambda = TMath::Log(2.)/t12;
1372 }
1373 
1374 ////////////////////////////////////////////////////////////////////////////////
1375 /// Default ctor.
1376 
1377 TGeoBatemanSol::TGeoBatemanSol(const TObjArray *chain)
1378  :TObject(), TAttLine(), TAttFill(), TAttMarker(),
1379  fElem(NULL),
1380  fElemTop(NULL),
1381  fCsize(0),
1382  fNcoeff(0),
1383  fFactor(1.),
1384  fTmin(0.),
1385  fTmax(0.),
1386  fCoeff(NULL)
1387 {
1388  TGeoDecayChannel *dc = (TGeoDecayChannel*)chain->At(0);
1389  if (dc) fElemTop = dc->Parent();
1390  dc = (TGeoDecayChannel*)chain->At(chain->GetEntriesFast()-1);
1391  if (dc) {
1392  fElem = dc->Daughter();
1393  fCsize = chain->GetEntriesFast()+1;
1394  fCoeff = new BtCoef_t[fCsize];
1395  FindSolution(chain);
1396  }
1397 }
1398 
1399 ////////////////////////////////////////////////////////////////////////////////
1400 /// Copy constructor.
1401 
1402 TGeoBatemanSol::TGeoBatemanSol(const TGeoBatemanSol& other)
1403  :TObject(other), TAttLine(other), TAttFill(other), TAttMarker(other),
1404  fElem(other.fElem),
1405  fElemTop(other.fElemTop),
1406  fCsize(other.fCsize),
1407  fNcoeff(other.fNcoeff),
1408  fFactor(other.fFactor),
1409  fTmin(other.fTmin),
1410  fTmax(other.fTmax),
1411  fCoeff(NULL)
1412 {
1413  if (fCsize) {
1414  fCoeff = new BtCoef_t[fCsize];
1415  for (Int_t i=0; i<fNcoeff; i++) {
1416  fCoeff[i].cn = other.fCoeff[i].cn;
1417  fCoeff[i].lambda = other.fCoeff[i].lambda;
1418  }
1419  }
1420 }
1421 
1422 ////////////////////////////////////////////////////////////////////////////////
1423 /// Destructor.
1424 
1425 TGeoBatemanSol::~TGeoBatemanSol()
1426 {
1427  if (fCoeff) delete [] fCoeff;
1428 }
1429 
1430 ////////////////////////////////////////////////////////////////////////////////
1431 /// Assignment.
1432 
1433 TGeoBatemanSol& TGeoBatemanSol::operator=(const TGeoBatemanSol& other)
1434 {
1435  if (this == &other) return *this;
1436  TObject::operator=(other);
1437  TAttLine::operator=(other);
1438  TAttFill::operator=(other);
1439  TAttMarker::operator=(other);
1440  fElem = other.fElem;
1441  fElemTop = other.fElemTop;
1442  if (fCoeff) delete [] fCoeff;
1443  fCoeff = 0;
1444  fCsize = other.fCsize;
1445  fNcoeff = other.fNcoeff;
1446  fFactor = other.fFactor;
1447  fTmin = other.fTmin;
1448  fTmax = other.fTmax;
1449  if (fCsize) {
1450  fCoeff = new BtCoef_t[fCsize];
1451  for (Int_t i=0; i<fNcoeff; i++) {
1452  fCoeff[i].cn = other.fCoeff[i].cn;
1453  fCoeff[i].lambda = other.fCoeff[i].lambda;
1454  }
1455  }
1456  return *this;
1457 }
1458 
1459 ////////////////////////////////////////////////////////////////////////////////
1460 /// Addition of other solution.
1461 
1462 TGeoBatemanSol& TGeoBatemanSol::operator+=(const TGeoBatemanSol& other)
1463 {
1464  if (other.GetElement() != fElem) {
1465  Error("operator+=", "Cannot add 2 solutions for different elements");
1466  return *this;
1467  }
1468  Int_t i,j;
1469  BtCoef_t *coeff = fCoeff;
1470  Int_t ncoeff = fNcoeff + other.fNcoeff;
1471  if (ncoeff > fCsize) {
1472  fCsize = ncoeff;
1473  coeff = new BtCoef_t[ncoeff];
1474  for (i=0; i<fNcoeff; i++) {
1475  coeff[i].cn = fCoeff[i].cn;
1476  coeff[i].lambda = fCoeff[i].lambda;
1477  }
1478  delete [] fCoeff;
1479  fCoeff = coeff;
1480  }
1481  ncoeff = fNcoeff;
1482  for (j=0; j<other.fNcoeff; j++) {
1483  for (i=0; i<fNcoeff; i++) {
1484  if (coeff[i].lambda == other.fCoeff[j].lambda) {
1485  coeff[i].cn += fFactor * other.fCoeff[j].cn;
1486  break;
1487  }
1488  }
1489  if (i == fNcoeff) {
1490  coeff[ncoeff].cn = fFactor * other.fCoeff[j].cn;
1491  coeff[ncoeff].lambda = other.fCoeff[j].lambda;
1492  ncoeff++;
1493  }
1494  }
1495  fNcoeff = ncoeff;
1496  return *this;
1497 }
1498 ////////////////////////////////////////////////////////////////////////////////
1499 /// Find concentration of the element at a given time.
1500 
1501 Double_t TGeoBatemanSol::Concentration(Double_t time) const
1502 {
1503  Double_t conc = 0.;
1504  for (Int_t i=0; i<fNcoeff; i++)
1505  conc += fCoeff[i].cn * TMath::Exp(-fCoeff[i].lambda * time);
1506  return conc;
1507 }
1508 
1509 ////////////////////////////////////////////////////////////////////////////////
1510 /// Draw the solution of Bateman equation versus time.
1511 
1512 void TGeoBatemanSol::Draw(Option_t *option)
1513 {
1514  if (!gGeoManager) return;
1515  gGeoManager->GetGeomPainter()->DrawBatemanSol(this, option);
1516 }
1517 
1518 ////////////////////////////////////////////////////////////////////////////////
1519 /// Find the solution for the Bateman equations corresponding to the decay
1520 /// chain described by an array ending with element X.
1521 /// A->B->...->X
1522 /// Cn = SUM [Ain * exp(-LMBDi*t)];
1523 /// Cn - concentration Nx/Na
1524 /// n - order of X in chain (A->B->X => n=3)
1525 /// LMBDi - decay constant for element of order i in the chain
1526 /// Ain = LMBD1*...*LMBD(n-1) * br1*...*br(n-1)/(LMBD1-LMBDi)...(LMBDn-LMBDi)
1527 /// bri - branching ratio for decay Ei->Ei+1
1528 
1529 void TGeoBatemanSol::FindSolution(const TObjArray *array)
1530 {
1531  fNcoeff = 0;
1532  if (!array || !array->GetEntriesFast()) return;
1533  Int_t n = array->GetEntriesFast();
1534  TGeoDecayChannel *dc = (TGeoDecayChannel*)array->At(n-1);
1535  TGeoElementRN *elem = dc->Daughter();
1536  if (elem != fElem) {
1537  Error("FindSolution", "Last element in the list must be %s\n", fElem->GetName());
1538  return;
1539  }
1540  Int_t i,j;
1541  Int_t order = n+1;
1542  if (!fCoeff) {
1543  fCsize = order;
1544  fCoeff = new BtCoef_t[fCsize];
1545  }
1546  if (fCsize < order) {
1547  delete [] fCoeff;
1548  fCsize = order;
1549  fCoeff = new BtCoef_t[fCsize];
1550  }
1551 
1552  Double_t *lambda = new Double_t[order];
1553  Double_t *br = new Double_t[n];
1554  Double_t halflife;
1555  for (i=0; i<n; i++) {
1556  dc = (TGeoDecayChannel*)array->At(i);
1557  elem = dc->Parent();
1558  br[i] = 0.01 * dc->BranchingRatio();
1559  halflife = elem->HalfLife();
1560  if (halflife==0.) halflife = 1.e-30;
1561  if (elem->Stable()) lambda[i] = 0.;
1562  else lambda[i] = TMath::Log(2.)/halflife;
1563  if (i==n-1) {
1564  elem = dc->Daughter();
1565  halflife = elem->HalfLife();
1566  if (halflife==0.) halflife = 1.e-30;
1567  if (elem->Stable()) lambda[n] = 0.;
1568  else lambda[n] = TMath::Log(2.)/halflife;
1569  }
1570  }
1571  // Check if we have equal lambdas
1572  for (i=0; i<order-1; i++) {
1573  for (j=i+1; j<order; j++) {
1574  if (lambda[j] == lambda[i]) lambda[j] += 0.001*lambda[j];
1575  }
1576  }
1577  Double_t ain;
1578  Double_t pdlambda, plambdabr=1.;
1579  for (j=0; j<n; j++) plambdabr *= lambda[j]*br[j];
1580  for (i=0; i<order; i++) {
1581  pdlambda = 1.;
1582  for (j=0; j<n+1; j++) {
1583  if (j == i) continue;
1584  pdlambda *= lambda[j] - lambda[i];
1585  }
1586  if (pdlambda == 0.) {
1587  Error("FindSolution", "pdlambda=0 !!!");
1588  delete [] lambda;
1589  delete [] br;
1590  return;
1591  }
1592  ain = plambdabr/pdlambda;
1593  fCoeff[i].cn = ain;
1594  fCoeff[i].lambda = lambda[i];
1595  }
1596  fNcoeff = order;
1597  Normalize(fFactor);
1598  delete [] lambda;
1599  delete [] br;
1600 }
1601 
1602 ////////////////////////////////////////////////////////////////////////////////
1603 /// Normalize all coefficients with a given factor.
1604 
1605 void TGeoBatemanSol::Normalize(Double_t factor)
1606 {
1607  for (Int_t i=0; i<fNcoeff; i++) fCoeff[i].cn *= factor;
1608 }
1609 
1610 ////////////////////////////////////////////////////////////////////////////////
1611 /// Print concentration evolution.
1612 
1613 void TGeoBatemanSol::Print(Option_t * /*option*/) const
1614 {
1615  TString formula;
1616  formula.Form("N[%s]/N[%s] = ", fElem->GetName(), fElemTop->GetName());
1617  for (Int_t i=0; i<fNcoeff; i++) {
1618  if (i == fNcoeff-1) formula += TString::Format("%g*exp(-%g*t)", fCoeff[i].cn, fCoeff[i].lambda);
1619  else formula += TString::Format("%g*exp(-%g*t) + ", fCoeff[i].cn, fCoeff[i].lambda);
1620  }
1621  printf("%s\n", formula.Data());
1622 }
1623