#ifdef R__HAVE_CONFIG
#include "RConfigure.h"
#endif
#include "Riostream.h"
#include "TObjArray.h"
#include "TVirtualGeoPainter.h"
#include "TGeoManager.h"
#include "TGeoElement.h"
static const Int_t gMaxElem = 110;
static const Int_t gMaxLevel = 8;
static const Int_t gMaxDecay = 15;
static const char gElName[gMaxElem][3] = {
"H ","He","Li","Be","B ","C ","N ","O ","F ","Ne","Na","Mg",
"Al","Si","P ","S ","Cl","Ar","K ","Ca","Sc","Ti","V ","Cr",
"Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr",
"Rb","Sr","Y ","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd",
"In","Sn","Sb","Te","I ","Xe","Cs","Ba","La","Ce","Pr","Nd",
"Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf",
"Ta","W ","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po",
"At","Rn","Fr","Ra","Ac","Th","Pa","U ","Np","Pu","Am","Cm",
"Bk","Cf","Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs",
"Mt","Ds" };
static const char *gDecayName[gMaxDecay+1] = {
"2BetaMinus", "BetaMinus", "NeutronEm", "ProtonEm", "Alpha", "ECF",
"ElecCapt", "IsoTrans", "I", "SpontFiss", "2ProtonEm", "2NeutronEm",
"2Alpha", "Carbon12", "Carbon14", "Stable" };
static const Int_t gDecayDeltaA[gMaxDecay] = {
0, 0, -1, -1, -4,
-99, 0, 0, -99, -99,
-2, -2, -8, -12, -14 };
static const Int_t gDecayDeltaZ[gMaxDecay] = {
2, 1, 0, -1, -2,
-99, -1, 0, -99, -99,
-2, 0, -4, -6, -6 };
static const char gLevName[gMaxLevel]=" mnpqrs";
ClassImp(TGeoElement)
TGeoElement::TGeoElement()
{
SetDefined(kFALSE);
SetUsed(kFALSE);
fZ = 0;
fA = 0.0;
}
TGeoElement::TGeoElement(const char *name, const char *title, Int_t z, Double_t a)
:TNamed(name, title)
{
SetDefined(kFALSE);
SetUsed(kFALSE);
fZ = z;
fA = a;
}
ClassImp(TGeoElementRN)
TGeoElementRN::TGeoElementRN()
{
TObject::SetBit(kElementChecked,kFALSE);
fENDFcode = 0;
fIso = 0;
fLevel = 0;
fDeltaM = 0;
fHalfLife = 0;
fNatAbun = 0;
fTH_F = 0;
fTG_F = 0;
fTH_S = 0;
fTG_S = 0;
fStatus = 0;
fRatio = 0;
fDecays = 0;
}
TGeoElementRN::TGeoElementRN(Int_t A, Int_t Z, Int_t iso, Double_t level,
Double_t deltaM, Double_t halfLife, const char* JP,
Double_t natAbun, Double_t th_f, Double_t tg_f, Double_t th_s,
Double_t tg_s, Int_t status)
:TGeoElement("", JP, Z, A)
{
TObject::SetBit(kElementChecked,kFALSE);
fENDFcode = ENDF(A,Z,iso);
fIso = iso;
fLevel = level;
fDeltaM = deltaM;
fHalfLife = halfLife;
fTitle = JP;
if (!fTitle.Length()) fTitle = "?";
fNatAbun = natAbun;
fTH_F = th_f;
fTG_F = tg_f;
fTH_S = th_s;
fTG_S = tg_s;
fStatus = status;
fDecays = 0;
fRatio = 0;
MakeName(A,Z,iso);
if ((TMath::Abs(fHalfLife)<1.e-30) || fHalfLife<-1) Warning("ctor","Element %s has T1/2=%g [s]", fName.Data(), fHalfLife);
}
TGeoElementRN::TGeoElementRN(const TGeoElementRN& elem) : TGeoElement(elem),
fENDFcode(elem.fENDFcode),
fIso(elem.fIso),
fLevel(elem.fLevel),
fDeltaM(elem.fDeltaM),
fHalfLife(elem.fHalfLife),
fNatAbun(elem.fNatAbun),
fTH_F(elem.fTH_F),
fTG_F(elem.fTG_F),
fTH_S(elem.fTH_S),
fTG_S(elem.fTG_S),
fStatus(elem.fStatus),
fRatio(NULL),
fDecays(NULL)
{
Error("cpy ctor", "Not to use !");
}
TGeoElementRN::~TGeoElementRN()
{
if (fDecays) {
fDecays->Delete();
delete fDecays;
}
if (fRatio) delete fRatio;
}
TGeoElementRN &TGeoElementRN::operator=(const TGeoElementRN&)
{
Error("operator=", "Not to use !");
return *this;
}
void TGeoElementRN::AddDecay(Int_t decay, Int_t diso, Double_t branchingRatio, Double_t qValue)
{
if (branchingRatio<1E-20) {
TString decayName;
TGeoDecayChannel::DecayName(decay, decayName);
Warning("AddDecay", "Decay %s of %s has BR=0. Not added.", decayName.Data(),fName.Data());
return;
}
TGeoDecayChannel *dc = new TGeoDecayChannel(decay,diso,branchingRatio, qValue);
dc->SetParent(this);
if (!fDecays) fDecays = new TObjArray(5);
fDecays->Add(dc);
}
void TGeoElementRN::AddDecay(TGeoDecayChannel *dc)
{
dc->SetParent(this);
if (!fDecays) fDecays = new TObjArray(5);
fDecays->Add(dc);
}
Int_t TGeoElementRN::GetNdecays() const
{
if (!fDecays) return 0;
return fDecays->GetEntriesFast();
}
Bool_t TGeoElementRN::CheckDecays() const
{
if (TObject::TestBit(kElementChecked)) return kTRUE;
TObject *oelem = (TObject*)this;
TGeoDecayChannel *dc;
TGeoElementRN *elem;
TGeoElementTable *table = GetElementTable();
TString decayName;
if (!table) {
Error("CheckDecays", "Element table not present");
return kFALSE;
}
Bool_t resultOK = kTRUE;
if (!fDecays) {
oelem->SetBit(kElementChecked,kTRUE);
return resultOK;
}
Double_t br = 0.;
Int_t decayResult = 0;
TIter next(fDecays);
while ((dc=(TGeoDecayChannel*)next())) {
br += dc->BranchingRatio();
decayResult = DecayResult(dc);
if (decayResult) {
elem = table->GetElementRN(decayResult);
if (!elem) {
TGeoDecayChannel::DecayName(dc->Decay(),decayName);
Error("CheckDecays", "Element after decay %s of %s not found in DB", decayName.Data(),fName.Data());
return kFALSE;
}
dc->SetDaughter(elem);
resultOK = elem->CheckDecays();
}
}
if (TMath::Abs(br-100) > 1.E-3) {
Warning("CheckDecays", "BR for decays of element %s sum-up = %f", fName.Data(), br);
resultOK = kFALSE;
}
oelem->SetBit(kElementChecked,kTRUE);
return resultOK;
}
Int_t TGeoElementRN::DecayResult(TGeoDecayChannel *dc) const
{
Int_t da, dz, diso;
dc->DecayShift(da, dz, diso);
if (da == -99 || dz == -99) return 0;
return ENDF(Int_t(fA)+da,fZ+dz,fIso+diso);
}
void TGeoElementRN::FillPopulation(TObjArray *population, Double_t precision, Double_t factor)
{
TGeoElementRN *elem;
TGeoElemIter next(this, precision);
TGeoBatemanSol s(this);
s.Normalize(factor);
AddRatio(s);
if (!population->FindObject(this)) population->Add(this);
while ((elem=next())) {
TGeoBatemanSol ratio(next.GetBranch());
ratio.Normalize(factor);
elem->AddRatio(ratio);
if (!population->FindObject(elem)) population->Add(elem);
}
}
void TGeoElementRN::MakeName(Int_t a, Int_t z, Int_t iso)
{
fName = "";
if (z==0 && a==1) {
fName = "neutron";
return;
}
if (z>=1 && z<= gMaxElem) fName += Form("%3d-%s-",z,gElName[z-1]);
else fName = "?? -?? -";
if (a>=1 && a<=999) fName += Form("%3.3d",a);
else fName += "??";
if (iso>0 && iso<gMaxLevel) fName += Form("%c", gLevName[iso]);
fName.ReplaceAll(" ","");
}
void TGeoElementRN::Print(Option_t *option) const
{
printf("\n%-12s ",fName.Data());
printf("ENDF=%d; ",fENDFcode);
printf("A=%d; ",(Int_t)fA);
printf("Z=%d; ",fZ);
printf("Iso=%d; ",fIso);
printf("Level=%g[MeV]; ",fLevel);
printf("Dmass=%g[MeV]; ",fDeltaM);
if (fHalfLife>0) printf("Hlife=%g[s]\n",fHalfLife);
else printf("Hlife=INF\n");
printf("%13s"," ");
printf("J/P=%s; ",fTitle.Data());
printf("Abund=%g; ",fNatAbun);
printf("Htox=%g; ",fTH_F);
printf("Itox=%g; ",fTG_F);
printf("Stat=%d\n",fStatus);
if(!fDecays) return;
printf("Decay modes:\n");
TIter next(fDecays);
TGeoDecayChannel *dc;
while ((dc=(TGeoDecayChannel*)next())) dc->Print(option);
}
TGeoElementRN *TGeoElementRN::ReadElementRN(const char *line, Int_t &ndecays)
{
Int_t a,z,iso,status;
Double_t level, deltaM, halfLife, natAbun, th_f, tg_f, th_s, tg_s;
char name[20],jp[20];
sscanf(&line[0], "%s%d%d%d%lg%lg%lg%s%lg%lg%lg%lg%lg%d%d", name,&a,&z,&iso,&level,&deltaM,
&halfLife,jp,&natAbun,&th_f,&tg_f,&th_s,&tg_s,&status,&ndecays);
TGeoElementRN *elem = new TGeoElementRN(a,z,iso,level,deltaM,halfLife,
jp,natAbun,th_f,tg_f,th_s,tg_s,status);
return elem;
}
void TGeoElementRN::SavePrimitive(ostream &out, Option_t *option)
{
if (!strcmp(option,"h")) {
out << "#====================================================================================================================================" << endl;
out << "# Name A Z ISO LEV[MeV] DM[MeV] T1/2[s] J/P ABUND[%] HTOX ITOX HTOX ITOX STAT NDCY" << endl;
out << "#====================================================================================================================================" << endl;
}
out << setw(11) << fName.Data();
out << setw(5) << (Int_t)fA;
out << setw(5) << fZ;
out << setw(5) << fIso;
out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fLevel;
out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fDeltaM;
out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << fHalfLife;
out << setw(13) << fTitle.Data();
out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fNatAbun;
out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fTH_F;
out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fTG_F;
out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fTH_S;
out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fTG_S;
out << setw(5) << fStatus;
Int_t ndecays = 0;
if (fDecays) ndecays = fDecays->GetEntries();
out << setw(5) << ndecays;
out << endl;
if (fDecays) {
TIter next(fDecays);
TGeoDecayChannel *dc;
while ((dc=(TGeoDecayChannel*)next())) dc->SavePrimitive(out);
}
}
void TGeoElementRN::AddRatio(TGeoBatemanSol &ratio)
{
if (!fRatio) fRatio = new TGeoBatemanSol(ratio);
else *fRatio += ratio;
}
void TGeoElementRN::ResetRatio()
{
if (fRatio) {
delete fRatio;
fRatio = 0;
}
}
ClassImp(TGeoDecayChannel)
TGeoDecayChannel& TGeoDecayChannel::operator=(const TGeoDecayChannel& dc)
{
if(this!=&dc) {
TObject::operator=(dc);
fDecay = dc.fDecay;
fDiso = dc.fDiso;
fBranchingRatio = dc.fBranchingRatio;
fQvalue = dc.fQvalue;
fParent = dc.fParent;
fDaughter = dc.fDaughter;
}
return *this;
}
const char *TGeoDecayChannel::GetName() const
{
static TString name = "";
name = "";
if (!fDecay) return gDecayName[gMaxDecay];
for (Int_t i=0; i<gMaxDecay; i++) {
if (1<<i & fDecay) {
if (name.Length()) name += "+";
name += gDecayName[i];
}
}
return name.Data();
}
void TGeoDecayChannel::DecayName(UInt_t decay, TString &name)
{
if (!decay) {
name = gDecayName[gMaxDecay];
return;
}
name = "";
for (Int_t i=0; i<gMaxDecay; i++) {
if (1<<i & decay) {
if (name.Length()) name += "+";
name += gDecayName[i];
}
}
}
Int_t TGeoDecayChannel::GetIndex() const
{
return fParent->Decays()->IndexOf(this);
}
void TGeoDecayChannel::Print(Option_t *) const
{
TString name;
DecayName(fDecay, name);
printf("%-20s Diso: %3d BR: %9.3f%% Qval: %g\n", name.Data(),fDiso,fBranchingRatio,fQvalue);
}
TGeoDecayChannel *TGeoDecayChannel::ReadDecay(const char *line)
{
char name[80];
Int_t decay,diso;
Double_t branchingRatio, qValue;
sscanf(&line[0], "%s%d%d%lg%lg", name,&decay,&diso,&branchingRatio,&qValue);
TGeoDecayChannel *dc = new TGeoDecayChannel(decay,diso,branchingRatio,qValue);
return dc;
}
void TGeoDecayChannel::SavePrimitive(ostream &out, Option_t *)
{
TString decayName;
DecayName(fDecay, decayName);
out << setw(50) << decayName.Data();
out << setw(10) << fDecay;
out << setw(10) << fDiso;
out << setw(12) << setiosflags(ios::fixed) << setprecision(6) << fBranchingRatio;
out << setw(12) << setiosflags(ios::fixed) << setprecision(6) << fQvalue;
out << endl;
}
void TGeoDecayChannel::DecayShift(Int_t &dA, Int_t &dZ, Int_t &dI) const
{
dA=dZ=0;
dI=fDiso;
for(Int_t i=0; i<gMaxDecay; ++i) {
if(1<<i & fDecay) {
if(gDecayDeltaA[i] == -99 || gDecayDeltaZ[i] == -99 ) {
dA=dZ=-99;
return;
}
dA += gDecayDeltaA[i];
dZ += gDecayDeltaZ[i];
}
}
}
ClassImp(TGeoElemIter)
TGeoElemIter::TGeoElemIter(TGeoElementRN *top, Double_t limit)
: fTop(top), fElem(top), fBranch(0), fLevel(0), fLimitRatio(limit), fRatio(1.)
{
fBranch = new TObjArray(10);
}
TGeoElemIter::TGeoElemIter(const TGeoElemIter &iter)
:fTop(iter.fTop),
fElem(iter.fElem),
fBranch(0),
fLevel(iter.fLevel),
fLimitRatio(iter.fLimitRatio),
fRatio(iter.fRatio)
{
if (iter.fBranch) {
fBranch = new TObjArray(10);
for (Int_t i=0; i<fLevel; i++) fBranch->Add(iter.fBranch->At(i));
}
}
TGeoElemIter::~TGeoElemIter()
{
if (fBranch) delete fBranch;
}
TGeoElemIter &TGeoElemIter::operator=(const TGeoElemIter &iter)
{
if (&iter == this) return *this;
fTop = iter.fTop;
fElem = iter.fElem;
fLevel = iter.fLevel;
if (iter.fBranch) {
fBranch = new TObjArray(10);
for (Int_t i=0; i<fLevel; i++) fBranch->Add(iter.fBranch->At(i));
}
fLimitRatio = iter.fLimitRatio;
fRatio = iter.fRatio;
return *this;
}
TGeoElementRN *TGeoElemIter::operator()()
{
return Next();
}
TGeoElementRN *TGeoElemIter::Up()
{
TGeoDecayChannel *dc;
Int_t ind, nd;
while (fLevel) {
dc = (TGeoDecayChannel*)fBranch->At(fLevel-1);
ind = dc->GetIndex();
nd = dc->Parent()->GetNdecays();
fRatio /= 0.01*dc->BranchingRatio();
fElem = dc->Parent();
fBranch->RemoveAt(--fLevel);
ind++;
while (ind<nd) {
if (Down(ind++)) return (TGeoElementRN*)fElem;
}
}
fElem = NULL;
return NULL;
}
TGeoElementRN *TGeoElemIter::Down(Int_t ibranch)
{
TGeoDecayChannel *dc = (TGeoDecayChannel*)fElem->Decays()->At(ibranch);
if (!dc->Daughter()) return NULL;
Double_t br = 0.01*fRatio*dc->BranchingRatio();
if (br < fLimitRatio) return NULL;
fLevel++;
fRatio = br;
fBranch->Add(dc);
fElem = dc->Daughter();
return (TGeoElementRN*)fElem;
}
TGeoElementRN *TGeoElemIter::Next()
{
if (!fElem) return NULL;
Int_t nd = fElem->GetNdecays();
for (Int_t i=0; i<nd; i++) if (Down(i)) return (TGeoElementRN*)fElem;
return Up();
}
void TGeoElemIter::Print(Option_t * ) const
{
TGeoElementRN *elem;
TGeoDecayChannel *dc;
TString indent = "";
printf("=== Chain with %g %%\n", 100*fRatio);
for (Int_t i=0; i<fLevel; i++) {
dc = (TGeoDecayChannel*)fBranch->At(i);
elem = dc->Parent();
printf("%s%s (%g%% %s) T1/2=%g\n", indent.Data(), elem->GetName(),dc->BranchingRatio(),dc->GetName(),elem->HalfLife());
indent += " ";
if (i==fLevel-1) {
elem = dc->Daughter();
printf("%s%s\n", indent.Data(), elem->GetName());
}
}
}
ClassImp(TGeoElementTable)
TGeoElementTable *TGeoElement::GetElementTable() const
{
return gGeoManager->GetElementTable();
}
TGeoElementTable::TGeoElementTable()
{
fNelements = 0;
fList = 0;
fListRN = 0;
}
TGeoElementTable::TGeoElementTable(Int_t )
{
fNelements = 0;
fList = new TObjArray(128);
fListRN = 0;
BuildDefaultElements();
}
TGeoElementTable::TGeoElementTable(const TGeoElementTable& get) :
TObject(get),
fNelements(get.fNelements),
fList(get.fList),
fListRN(get.fListRN)
{
}
TGeoElementTable& TGeoElementTable::operator=(const TGeoElementTable& get)
{
if(this!=&get) {
TObject::operator=(get);
fNelements=get.fNelements;
fList=get.fList;
fListRN=get.fListRN;
}
return *this;
}
TGeoElementTable::~TGeoElementTable()
{
if (fList) {
fList->Delete();
delete fList;
}
if (fListRN) {
fListRN->Delete();
delete fListRN;
}
}
void TGeoElementTable::AddElement(const char *name, const char *title, Int_t z, Double_t a)
{
if (!fList) fList = new TObjArray(128);
fList->AddAt(new TGeoElement(name,title,z,a), fNelements++);
}
void TGeoElementTable::AddElementRN(TGeoElementRN *elem)
{
if (!fListRN) fListRN = new TObjArray(3600);
if (HasRNElements() && GetElementRN(elem->ENDFCode())) return;
fListRN->Add(elem);
fNelementsRN++;
fElementsRN.insert(ElementRNMap_t::value_type(elem->ENDFCode(), elem));
}
void TGeoElementTable::BuildDefaultElements()
{
if (HasDefaultElements()) return;
AddElement("VACUUM","VACUUM" ,0, 0.0);
AddElement("H" ,"HYDROGEN" ,1, 1.00794);
AddElement("HE" ,"HELIUM" ,2, 4.002602);
AddElement("LI" ,"LITHIUM" ,3, 6.941);
AddElement("BE" ,"BERYLLIUM" ,4, 9.01218);
AddElement("B" ,"BORON" ,5, 10.811);
AddElement("C" ,"CARBON" ,6 ,12.0107);
AddElement("N" ,"NITROGEN" ,7 ,14.00674);
AddElement("O" ,"OXYGEN" ,8 ,15.9994);
AddElement("F" ,"FLUORINE" ,9 ,18.9984032);
AddElement("NE" ,"NEON" ,10 ,20.1797);
AddElement("NA" ,"SODIUM" ,11 ,22.989770);
AddElement("MG" ,"MAGNESIUM" ,12 ,24.3050);
AddElement("AL" ,"ALUMINIUM" ,13 ,26.981538);
AddElement("SI" ,"SILICON" ,14 ,28.0855);
AddElement("P" ,"PHOSPHORUS" ,15 ,30.973761);
AddElement("S" ,"SULFUR" ,16 ,32.066);
AddElement("CL" ,"CHLORINE" ,17 ,35.4527);
AddElement("AR" ,"ARGON" ,18 ,39.948);
AddElement("K" ,"POTASSIUM" ,19 ,39.0983);
AddElement("CA" ,"CALCIUM" ,20 ,40.078);
AddElement("SC" ,"SCANDIUM" ,21 ,44.955910);
AddElement("TI" ,"TITANIUM" ,22 ,47.867);
AddElement("V" ,"VANADIUM" ,23 ,50.9415);
AddElement("CR" ,"CHROMIUM" ,24 ,51.9961);
AddElement("MN" ,"MANGANESE" ,25 ,54.938049);
AddElement("FE" ,"IRON" ,26 ,55.845);
AddElement("CO" ,"COBALT" ,27 ,58.933200);
AddElement("NI" ,"NICKEL" ,28 ,58.6934);
AddElement("CU" ,"COPPER" ,29 ,63.546);
AddElement("ZN" ,"ZINC" ,30 ,65.39);
AddElement("GA" ,"GALLIUM" ,31 ,69.723);
AddElement("GE" ,"GERMANIUM" ,32 ,72.61);
AddElement("AS" ,"ARSENIC" ,33 ,74.92160);
AddElement("SE" ,"SELENIUM" ,34 ,78.96);
AddElement("BR" ,"BROMINE" ,35 ,79.904);
AddElement("KR" ,"KRYPTON" ,36 ,83.80);
AddElement("RB" ,"RUBIDIUM" ,37 ,85.4678);
AddElement("SR" ,"STRONTIUM" ,38 ,87.62);
AddElement("Y" ,"YTTRIUM" ,39 ,88.90585);
AddElement("ZR" ,"ZIRCONIUM" ,40 ,91.224);
AddElement("NB" ,"NIOBIUM" ,41 ,92.90638);
AddElement("MO" ,"MOLYBDENUM" ,42 ,95.94);
AddElement("TC" ,"TECHNETIUM" ,43 ,98.0);
AddElement("RU" ,"RUTHENIUM" ,44 ,101.07);
AddElement("RH" ,"RHODIUM" ,45 ,102.90550);
AddElement("PD" ,"PALLADIUM" ,46 ,106.42);
AddElement("AG" ,"SILVER" ,47 ,107.8682);
AddElement("CD" ,"CADMIUM" ,48 ,112.411);
AddElement("IN" ,"INDIUM" ,49 ,114.818);
AddElement("SN" ,"TIN" ,50 ,118.710);
AddElement("SB" ,"ANTIMONY" ,51 ,121.760);
AddElement("TE" ,"TELLURIUM" ,52 ,127.60);
AddElement("I" ,"IODINE" ,53 ,126.90447);
AddElement("XE" ,"XENON" ,54 ,131.29);
AddElement("CS" ,"CESIUM" ,55 ,132.90545);
AddElement("BA" ,"BARIUM" ,56 ,137.327);
AddElement("LA" ,"LANTHANUM" ,57 ,138.9055);
AddElement("CE" ,"CERIUM" ,58 ,140.116);
AddElement("PR" ,"PRASEODYMIUM" ,59 ,140.90765);
AddElement("ND" ,"NEODYMIUM" ,60 ,144.24);
AddElement("PM" ,"PROMETHIUM" ,61 ,145.0);
AddElement("SM" ,"SAMARIUM" ,62 ,150.36);
AddElement("EU" ,"EUROPIUM" ,63 ,151.964);
AddElement("GD" ,"GADOLINIUM" ,64 ,157.25);
AddElement("TB" ,"TERBIUM" ,65 ,158.92534);
AddElement("DY" ,"DYSPROSIUM" ,66 ,162.50);
AddElement("HO" ,"HOLMIUM" ,67 ,164.93032);
AddElement("ER" ,"ERBIUM" ,68 ,167.26);
AddElement("TM" ,"THULIUM" ,69 ,168.93421);
AddElement("YB" ,"YTTERBIUM" ,70 ,173.04);
AddElement("LU" ,"LUTETIUM" ,71 ,174.967);
AddElement("HF" ,"HAFNIUM" ,72 ,178.49);
AddElement("TA" ,"TANTALUM" ,73 ,180.9479);
AddElement("W" ,"TUNGSTEN" ,74 ,183.84);
AddElement("RE" ,"RHENIUM" ,75 ,186.207);
AddElement("OS" ,"OSMIUM" ,76 ,190.23);
AddElement("IR" ,"IRIDIUM" ,77 ,192.217);
AddElement("PT" ,"PLATINUM" ,78 ,195.078);
AddElement("AU" ,"GOLD" ,79 ,196.96655);
AddElement("HG" ,"MERCURY" ,80 ,200.59);
AddElement("TL" ,"THALLIUM" ,81 ,204.3833);
AddElement("PB" ,"LEAD" ,82 ,207.2);
AddElement("BI" ,"BISMUTH" ,83 ,208.98038);
AddElement("PO" ,"POLONIUM" ,84 ,209.0);
AddElement("AT" ,"ASTATINE" ,85 ,210.0);
AddElement("RN" ,"RADON" ,86 ,222.0);
AddElement("FR" ,"FRANCIUM" ,87 ,223.0);
AddElement("RA" ,"RADIUM" ,88 ,226.0);
AddElement("AC" ,"ACTINIUM" ,89 ,227.0);
AddElement("TH" ,"THORIUM" ,90 ,232.0381);
AddElement("PA" ,"PROTACTINIUM" ,91 ,231.03588);
AddElement("U" ,"URANIUM" ,92 ,238.0289);
AddElement("NP" ,"NEPTUNIUM" ,93 ,237.0);
AddElement("PU" ,"PLUTONIUM" ,94 ,244.0);
AddElement("AM" ,"AMERICIUM" ,95 ,243.0);
AddElement("CM" ,"CURIUM" ,96 ,247.0);
AddElement("BK" ,"BERKELIUM" ,97 ,247.0);
AddElement("CF" ,"CALIFORNIUM",98 ,251.0);
AddElement("ES" ,"EINSTEINIUM",99 ,252.0);
AddElement("FM" ,"FERMIUM" ,100 ,257.0);
AddElement("MD" ,"MENDELEVIUM",101 ,258.0);
AddElement("NO" ,"NOBELIUM" ,102 ,259.0);
AddElement("LR" ,"LAWRENCIUM" ,103 ,262.0);
AddElement("RF" ,"RUTHERFORDIUM" ,104,261.0);
AddElement("DB" ,"DUBNIUM" ,105 ,262.0);
AddElement("SG" ,"SEABORGIUM" ,106 ,263.0);
AddElement("BH" ,"BOHRIUM" ,107 ,262.0);
AddElement("HS" ,"HASSIUM" ,108 ,265.0);
AddElement("MT" ,"MEITNERIUM" ,109 ,266.0);
AddElement("UUN" ,"UNUNNILIUM" ,110 ,269.0);
AddElement("UUU" ,"UNUNUNIUM" ,111 ,272.0);
AddElement("UUB" ,"UNUNBIUM" ,112 ,277.0);
TObject::SetBit(kETDefaultElements,kTRUE);
}
void TGeoElementTable::ImportElementsRN()
{
if (HasRNElements()) return;
TGeoElementRN *elem;
TString rnf;
#ifdef ROOTETCDIR
rnf.Form("%s/RadioNuclides.txt", ROOTETCDIR);
#else
rnf.Form("%s/etc/RadioNuclides.txt", gSystem->Getenv("ROOTSYS"));
#endif
FILE *fp = fopen(rnf, "r");
if (!fp) {
Error("ImportElementsRN","File RadioNuclides.txt not found");
return;
}
char line[150];
Int_t ndecays = 0;
Int_t i;
while (fgets(&line[0],140,fp)) {
if (line[0]=='#') continue;
elem = TGeoElementRN::ReadElementRN(line, ndecays);
for (i=0; i<ndecays; i++) {
fgets(&line[0],140,fp);
TGeoDecayChannel *dc = TGeoDecayChannel::ReadDecay(line);
elem->AddDecay(dc);
}
AddElementRN(elem);
}
TObject::SetBit(kETRNElements,kTRUE);
CheckTable();
}
Bool_t TGeoElementTable::CheckTable() const
{
if (!HasRNElements()) return HasDefaultElements();
TGeoElementRN *elem;
Bool_t result = kTRUE;
TIter next(fListRN);
while ((elem=(TGeoElementRN*)next())) {
if (!elem->CheckDecays()) result = kFALSE;
}
return result;
}
void TGeoElementTable::ExportElementsRN(const char *filename)
{
if (!HasRNElements()) return;
TString sname = filename;
if (!sname.Length()) sname = "RadioNuclides.txt";
ofstream out;
out.open(sname.Data(), ios::out);
if (!out.good()) {
Error("ExportElementsRN", "Cannot open file %s", sname.Data());
return;
}
TGeoElementRN *elem;
TIter next(fListRN);
Int_t i=0;
while ((elem=(TGeoElementRN*)next())) {
if ((i%48)==0) elem->SavePrimitive(out,"h");
else elem->SavePrimitive(out);
i++;
}
out.close();
}
TGeoElement *TGeoElementTable::FindElement(const char *name)
{
TString s(name);
s.ToUpper();
TGeoElement *elem;
elem = (TGeoElement*)fList->FindObject(s.Data());
if (elem) return elem;
TIter next(fList);
while ((elem=(TGeoElement*)next())) {
if (s == elem->GetTitle()) return elem;
}
return 0;
}
TGeoElementRN *TGeoElementTable::GetElementRN(Int_t ENDFcode) const
{
if (!HasRNElements()) {
TGeoElementTable *table = (TGeoElementTable*)this;
table->ImportElementsRN();
if (!fListRN) return 0;
}
ElementRNMap_t::const_iterator it = fElementsRN.find(ENDFcode);
if (it != fElementsRN.end()) return it->second;
return 0;
}
TGeoElementRN *TGeoElementTable::GetElementRN(Int_t a, Int_t z, Int_t iso) const
{
return GetElementRN(TGeoElementRN::ENDF(a,z,iso));
}
ClassImp(TGeoBatemanSol)
TGeoBatemanSol::TGeoBatemanSol(TGeoElementRN *elem)
:TObject(), TAttLine(), TAttFill(), TAttMarker(),
fElem(elem),
fElemTop(elem),
fCsize(10),
fNcoeff(0),
fFactor(1.),
fTmin(0.),
fTmax(0.),
fCoeff(NULL)
{
fCoeff = new BtCoef_t[fCsize];
fNcoeff = 1;
fCoeff[0].cn = 1.;
Double_t t12 = elem->HalfLife();
if (t12 == 0.) t12 = 1.e-30;
if (elem->Stable()) fCoeff[0].lambda = 0.;
else fCoeff[0].lambda = 1./t12;
}
TGeoBatemanSol::TGeoBatemanSol(const TObjArray *chain)
:TObject(), TAttLine(), TAttFill(), TAttMarker(),
fElem(NULL),
fElemTop(NULL),
fCsize(0),
fNcoeff(0),
fFactor(1.),
fTmin(0.),
fTmax(0.),
fCoeff(NULL)
{
TGeoDecayChannel *dc = (TGeoDecayChannel*)chain->At(0);
if (dc) fElemTop = dc->Parent();
dc = (TGeoDecayChannel*)chain->At(chain->GetEntriesFast()-1);
if (dc) {
fElem = dc->Daughter();
fCsize = chain->GetEntriesFast()+1;
fCoeff = new BtCoef_t[fCsize];
FindSolution(chain);
}
}
TGeoBatemanSol::TGeoBatemanSol(const TGeoBatemanSol& other)
:TObject(other), TAttLine(other), TAttFill(other), TAttMarker(other),
fElem(other.fElem),
fElemTop(other.fElemTop),
fCsize(other.fCsize),
fNcoeff(other.fNcoeff),
fFactor(other.fFactor),
fTmin(other.fTmin),
fTmax(other.fTmax),
fCoeff(NULL)
{
if (fCsize) {
fCoeff = new BtCoef_t[fCsize];
for (Int_t i=0; i<fNcoeff; i++) {
fCoeff[i].cn = other.fCoeff[i].cn;
fCoeff[i].lambda = other.fCoeff[i].lambda;
}
}
}
TGeoBatemanSol::~TGeoBatemanSol()
{
if (fCoeff) delete [] fCoeff;
}
TGeoBatemanSol& TGeoBatemanSol::operator=(const TGeoBatemanSol& other)
{
if (this == &other) return *this;
TObject::operator=(other);
TAttLine::operator=(other);
TAttFill::operator=(other);
TAttMarker::operator=(other);
fElem = other.fElem;
fElemTop = other.fElemTop;
if (fCoeff) delete [] fCoeff;
fCoeff = 0;
fCsize = other.fCsize;
fNcoeff = other.fNcoeff;
fFactor = other.fFactor;
fTmin = other.fTmin;
fTmax = other.fTmax;
if (fCsize) {
fCoeff = new BtCoef_t[fCsize];
for (Int_t i=0; i<fNcoeff; i++) {
fCoeff[i].cn = other.fCoeff[i].cn;
fCoeff[i].lambda = other.fCoeff[i].lambda;
}
}
return *this;
}
TGeoBatemanSol& TGeoBatemanSol::operator+=(const TGeoBatemanSol& other)
{
if (other.GetElement() != fElem) {
Error("operator+=", "Cannot add 2 solutions for different elements");
return *this;
}
Int_t i,j;
BtCoef_t *coeff = fCoeff;
Int_t ncoeff = fNcoeff + other.fNcoeff;
if (ncoeff > fCsize) {
fCsize = ncoeff;
coeff = new BtCoef_t[ncoeff];
for (i=0; i<fNcoeff; i++) {
coeff[i].cn = fCoeff[i].cn;
coeff[i].lambda = fCoeff[i].lambda;
}
delete [] fCoeff;
fCoeff = coeff;
}
ncoeff = fNcoeff;
for (j=0; j<other.fNcoeff; j++) {
for (i=0; i<fNcoeff; i++) {
if (coeff[i].lambda == other.fCoeff[j].lambda) {
coeff[i].cn += fFactor * other.fCoeff[j].cn;
break;
}
}
if (i == fNcoeff) {
coeff[ncoeff].cn = fFactor * other.fCoeff[j].cn;
coeff[ncoeff].lambda = other.fCoeff[j].lambda;
ncoeff++;
}
}
fNcoeff = ncoeff;
return *this;
}
Double_t TGeoBatemanSol::Concentration(Double_t time) const
{
Double_t conc = 0.;
for (Int_t i=0; i<fNcoeff; i++)
conc += fCoeff[i].cn * TMath::Exp(-fCoeff[i].lambda * time);
return conc;
}
void TGeoBatemanSol::Draw(Option_t *option)
{
if (!gGeoManager) return;
gGeoManager->GetGeomPainter()->DrawBatemanSol(this, option);
}
void TGeoBatemanSol::FindSolution(const TObjArray *array)
{
fNcoeff = 0;
if (!array || !array->GetEntriesFast()) return;
Int_t n = array->GetEntriesFast();
TGeoDecayChannel *dc = (TGeoDecayChannel*)array->At(n-1);
TGeoElementRN *elem = dc->Daughter();
if (elem != fElem) {
Error("FindSolution", "Last element in the list must be %s\n", fElem->GetName());
return;
}
Int_t i,j;
Int_t order = n+1;
if (!fCoeff) {
fCsize = order;
fCoeff = new BtCoef_t[fCsize];
}
if (fCsize < order) {
delete [] fCoeff;
fCsize = order;
fCoeff = new BtCoef_t[fCsize];
}
Double_t *lambda = new Double_t[order];
Double_t *br = new Double_t[n];
Double_t halflife;
for (i=0; i<n; i++) {
dc = (TGeoDecayChannel*)array->At(i);
elem = dc->Parent();
br[i] = 0.01 * dc->BranchingRatio();
halflife = elem->HalfLife();
if (halflife==0.) halflife = 1.e-30;
if (elem->Stable()) lambda[i] = 0.;
else lambda[i] = 1./halflife;
if (i==n-1) {
elem = dc->Daughter();
halflife = elem->HalfLife();
if (halflife==0.) halflife = 1.e-30;
if (elem->Stable()) lambda[n] = 0.;
else lambda[n] = 1./halflife;
}
}
for (i=0; i<order-1; i++) {
for (j=i+1; j<order; j++) {
if (lambda[j] == lambda[i]) lambda[j] += 0.001*lambda[j];
}
}
Double_t ain;
Double_t pdlambda, plambdabr=1.;
for (j=0; j<n; j++) plambdabr *= lambda[j]*br[j];
for (i=0; i<order; i++) {
pdlambda = 1.;
for (j=0; j<n+1; j++) {
if (j == i) continue;
pdlambda *= lambda[j] - lambda[i];
}
if (pdlambda == 0.) {
Error("FindSolution", "pdlambda=0 !!!");
delete [] lambda;
delete [] br;
return;
}
ain = plambdabr/pdlambda;
fCoeff[i].cn = ain;
fCoeff[i].lambda = lambda[i];
}
fNcoeff = order;
Normalize(fFactor);
}
void TGeoBatemanSol::Normalize(Double_t factor)
{
for (Int_t i=0; i<fNcoeff; i++) fCoeff[i].cn *= factor;
}
void TGeoBatemanSol::Print(Option_t * ) const
{
TString formula = Form("N[%s]/N[%s] = ", fElem->GetName(), fElemTop->GetName());
for (Int_t i=0; i<fNcoeff; i++) {
if (i == fNcoeff-1) formula += Form("%g*exp(-%g*t)", fCoeff[i].cn, fCoeff[i].lambda);
else formula += Form("%g*exp(-%g*t) + ", fCoeff[i].cn, fCoeff[i].lambda);
}
printf("%s\n", formula.Data());
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.