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