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