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