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