00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "RConfigure.h"
00023
00024 #include "Riostream.h"
00025
00026 #include "TSystem.h"
00027 #include "TObjArray.h"
00028 #include "TVirtualGeoPainter.h"
00029 #include "TGeoManager.h"
00030 #include "TGeoElement.h"
00031 #include "TMath.h"
00032
00033
00034 static const Int_t gMaxElem = 110;
00035 static const Int_t gMaxLevel = 8;
00036 static const Int_t gMaxDecay = 15;
00037
00038 static const char gElName[gMaxElem][3] = {
00039 "H ","He","Li","Be","B ","C ","N ","O ","F ","Ne","Na","Mg",
00040 "Al","Si","P ","S ","Cl","Ar","K ","Ca","Sc","Ti","V ","Cr",
00041 "Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr",
00042 "Rb","Sr","Y ","Zr","Nb","Mo","Tc","Ru","Rh","Pd","Ag","Cd",
00043 "In","Sn","Sb","Te","I ","Xe","Cs","Ba","La","Ce","Pr","Nd",
00044 "Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf",
00045 "Ta","W ","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po",
00046 "At","Rn","Fr","Ra","Ac","Th","Pa","U ","Np","Pu","Am","Cm",
00047 "Bk","Cf","Es","Fm","Md","No","Lr","Rf","Db","Sg","Bh","Hs",
00048 "Mt","Ds" };
00049
00050 static const char *gDecayName[gMaxDecay+1] = {
00051 "2BetaMinus", "BetaMinus", "NeutronEm", "ProtonEm", "Alpha", "ECF",
00052 "ElecCapt", "IsoTrans", "I", "SpontFiss", "2ProtonEm", "2NeutronEm",
00053 "2Alpha", "Carbon12", "Carbon14", "Stable" };
00054
00055 static const Int_t gDecayDeltaA[gMaxDecay] = {
00056 0, 0, -1, -1, -4,
00057 -99, 0, 0, -99, -99,
00058 -2, -2, -8, -12, -14 };
00059
00060 static const Int_t gDecayDeltaZ[gMaxDecay] = {
00061 2, 1, 0, -1, -2,
00062 -99, -1, 0, -99, -99,
00063 -2, 0, -4, -6, -6 };
00064 static const char gLevName[gMaxLevel]=" mnpqrs";
00065
00066 ClassImp(TGeoElement)
00067
00068
00069 TGeoElement::TGeoElement()
00070 {
00071
00072 SetDefined(kFALSE);
00073 SetUsed(kFALSE);
00074 fZ = 0;
00075 fN = 0;
00076 fNisotopes = 0;
00077 fA = 0.0;
00078 fIsotopes = NULL;
00079 fAbundances = NULL;
00080 }
00081
00082
00083 TGeoElement::TGeoElement(const char *name, const char *title, Int_t z, Double_t a)
00084 :TNamed(name, title)
00085 {
00086
00087 SetDefined(kFALSE);
00088 SetUsed(kFALSE);
00089 fZ = z;
00090 fN = Int_t(a);
00091 fNisotopes = 0;
00092 fA = a;
00093 fIsotopes = NULL;
00094 fAbundances = NULL;
00095 }
00096
00097
00098 TGeoElement::TGeoElement(const char *name, const char *title, Int_t nisotopes)
00099 :TNamed(name, title)
00100 {
00101
00102 SetDefined(kFALSE);
00103 SetUsed(kFALSE);
00104 fZ = 0;
00105 fN = 0;
00106 fNisotopes = nisotopes;
00107 fA = 0.0;
00108 fIsotopes = new TObjArray(nisotopes);
00109 fAbundances = new Double_t[nisotopes];
00110 }
00111
00112
00113 TGeoElement::TGeoElement(const char *name, const char *title, Int_t z, Int_t n, Double_t a)
00114 :TNamed(name, title)
00115 {
00116
00117 SetDefined(kFALSE);
00118 SetUsed(kFALSE);
00119 fZ = z;
00120 fN = n;
00121 fNisotopes = 0;
00122 fA = a;
00123 fIsotopes = NULL;
00124 fAbundances = NULL;
00125 }
00126
00127
00128 void TGeoElement::Print(Option_t *option) const
00129 {
00130
00131 printf("Element: %s Z=%d N=%f A=%f [g/mole]\n", GetName(), fZ,Neff(),fA);
00132 if (HasIsotopes()) {
00133 for (Int_t i=0; i<fNisotopes; i++) {
00134 TGeoIsotope *iso = GetIsotope(i);
00135 printf("=>Isotope %s, abundance=%f :\n", iso->GetName(), fAbundances[i]);
00136 iso->Print(option);
00137 }
00138 }
00139 }
00140
00141
00142 TGeoElementTable *TGeoElement::GetElementTable()
00143 {
00144
00145 if (!gGeoManager) {
00146 ::Error("TGeoElementTable::GetElementTable", "Create a geometry manager first");
00147 return NULL;
00148 }
00149 return gGeoManager->GetElementTable();
00150 }
00151
00152
00153 void TGeoElement::AddIsotope(TGeoIsotope *isotope, Double_t relativeAbundance)
00154 {
00155
00156 if (!fIsotopes) {
00157 Fatal("AddIsotope", "Cannot add isotopes to normal elements. Use constructor with number of isotopes.");
00158 return;
00159 }
00160 Int_t ncurrent = 0;
00161 TGeoIsotope *isocrt;
00162 for (ncurrent=0; ncurrent<fNisotopes; ncurrent++)
00163 if (!fIsotopes->At(ncurrent)) break;
00164 if (ncurrent==fNisotopes) {
00165 Error("AddIsotope", "All %d isotopes of element %s already defined", fNisotopes, GetName());
00166 return;
00167 }
00168
00169 if ((fZ!=0) && (isotope->GetZ()!=fZ)) {
00170 Fatal("AddIsotope", "Trying to add isotope %s with different Z to the same element %s",
00171 isotope->GetName(), GetName());
00172 return;
00173 } else {
00174 fZ = isotope->GetZ();
00175 }
00176 fIsotopes->Add(isotope);
00177 fAbundances[ncurrent] = relativeAbundance;
00178 if (ncurrent==fNisotopes-1) {
00179 Double_t weight = 0.0;
00180 Double_t aeff = 0.0;
00181 Double_t neff = 0.0;
00182 for (Int_t i=0; i<fNisotopes; i++) {
00183 isocrt = (TGeoIsotope*)fIsotopes->At(i);
00184 aeff += fAbundances[i]*isocrt->GetA();
00185 neff += fAbundances[i]*isocrt->GetN();
00186 weight += fAbundances[i];
00187 }
00188 aeff /= weight;
00189 neff /= weight;
00190 fN = (Int_t)neff;
00191 fA = aeff;
00192 }
00193 }
00194
00195
00196 Double_t TGeoElement::Neff() const
00197 {
00198
00199 if (!fNisotopes) return fN;
00200 TGeoIsotope *isocrt;
00201 Double_t weight = 0.0;
00202 Double_t neff = 0.0;
00203 for (Int_t i=0; i<fNisotopes; i++) {
00204 isocrt = (TGeoIsotope*)fIsotopes->At(i);
00205 neff += fAbundances[i]*isocrt->GetN();
00206 weight += fAbundances[i];
00207 }
00208 neff /= weight;
00209 return neff;
00210 }
00211
00212
00213 TGeoIsotope *TGeoElement::GetIsotope(Int_t i) const
00214 {
00215
00216 if (i>=0 && i<fNisotopes) {
00217 return (TGeoIsotope*)fIsotopes->At(i);
00218 }
00219 return NULL;
00220 }
00221
00222
00223 Double_t TGeoElement::GetRelativeAbundance(Int_t i) const
00224 {
00225
00226 if (i>=0 && i<fNisotopes) return fAbundances[i];
00227 return 0.0;
00228 }
00229
00230 ClassImp(TGeoIsotope)
00231
00232
00233 TGeoIsotope::TGeoIsotope()
00234 :TNamed(),
00235 fZ(0),
00236 fN(0),
00237 fA(0)
00238 {
00239
00240 }
00241
00242
00243 TGeoIsotope::TGeoIsotope(const char *name, Int_t z, Int_t n, Double_t a)
00244 :TNamed(name,""),
00245 fZ(z),
00246 fN(n),
00247 fA(a)
00248 {
00249
00250 if (z<1) Fatal("ctor", "Not allowed Z=%d (<1) for isotope: %s", z,name);
00251 if (n<z) Fatal("ctor", "Not allowed Z=%d < N=%d for isotope: %s", z,n,name);
00252 TGeoElement::GetElementTable()->AddIsotope(this);
00253 }
00254
00255
00256 TGeoIsotope *TGeoIsotope::FindIsotope(const char *name)
00257 {
00258
00259 TGeoElementTable *elTable = TGeoElement::GetElementTable();
00260 return elTable->FindIsotope(name);
00261 }
00262
00263
00264 void TGeoIsotope::Print(Option_t *) const
00265 {
00266
00267 printf("Isotope: %s Z=%d N=%d A=%f [g/mole]\n", GetName(), fZ,fN,fA);
00268 }
00269
00270 ClassImp(TGeoElementRN)
00271
00272
00273 TGeoElementRN::TGeoElementRN()
00274 {
00275
00276 TObject::SetBit(kElementChecked,kFALSE);
00277 fENDFcode = 0;
00278 fIso = 0;
00279 fLevel = 0;
00280 fDeltaM = 0;
00281 fHalfLife = 0;
00282 fNatAbun = 0;
00283 fTH_F = 0;
00284 fTG_F = 0;
00285 fTH_S = 0;
00286 fTG_S = 0;
00287 fStatus = 0;
00288 fRatio = 0;
00289 fDecays = 0;
00290 }
00291
00292
00293 TGeoElementRN::TGeoElementRN(Int_t aa, Int_t zz, Int_t iso, Double_t level,
00294 Double_t deltaM, Double_t halfLife, const char* JP,
00295 Double_t natAbun, Double_t th_f, Double_t tg_f, Double_t th_s,
00296 Double_t tg_s, Int_t status)
00297 :TGeoElement("", JP, zz, aa)
00298 {
00299
00300 TObject::SetBit(kElementChecked,kFALSE);
00301 fENDFcode = ENDF(aa,zz,iso);
00302 fIso = iso;
00303 fLevel = level;
00304 fDeltaM = deltaM;
00305 fHalfLife = halfLife;
00306 fTitle = JP;
00307 if (!fTitle.Length()) fTitle = "?";
00308 fNatAbun = natAbun;
00309 fTH_F = th_f;
00310 fTG_F = tg_f;
00311 fTH_S = th_s;
00312 fTG_S = tg_s;
00313 fStatus = status;
00314 fDecays = 0;
00315 fRatio = 0;
00316 MakeName(aa,zz,iso);
00317 if ((TMath::Abs(fHalfLife)<1.e-30) || fHalfLife<-1) Warning("ctor","Element %s has T1/2=%g [s]", fName.Data(), fHalfLife);
00318 }
00319
00320
00321 TGeoElementRN::~TGeoElementRN()
00322 {
00323
00324 if (fDecays) {
00325 fDecays->Delete();
00326 delete fDecays;
00327 }
00328 if (fRatio) delete fRatio;
00329 }
00330
00331
00332 void TGeoElementRN::AddDecay(Int_t decay, Int_t diso, Double_t branchingRatio, Double_t qValue)
00333 {
00334
00335 if (branchingRatio<1E-20) {
00336 TString decayName;
00337 TGeoDecayChannel::DecayName(decay, decayName);
00338 Warning("AddDecay", "Decay %s of %s has BR=0. Not added.", decayName.Data(),fName.Data());
00339 return;
00340 }
00341 TGeoDecayChannel *dc = new TGeoDecayChannel(decay,diso,branchingRatio, qValue);
00342 dc->SetParent(this);
00343 if (!fDecays) fDecays = new TObjArray(5);
00344 fDecays->Add(dc);
00345 }
00346
00347
00348 void TGeoElementRN::AddDecay(TGeoDecayChannel *dc)
00349 {
00350
00351 dc->SetParent(this);
00352 if (!fDecays) fDecays = new TObjArray(5);
00353 fDecays->Add(dc);
00354 }
00355
00356
00357 Int_t TGeoElementRN::GetNdecays() const
00358 {
00359
00360 if (!fDecays) return 0;
00361 return fDecays->GetEntriesFast();
00362 }
00363
00364
00365 Bool_t TGeoElementRN::CheckDecays() const
00366 {
00367
00368 if (TObject::TestBit(kElementChecked)) return kTRUE;
00369 TObject *oelem = (TObject*)this;
00370 TGeoDecayChannel *dc;
00371 TGeoElementRN *elem;
00372 TGeoElementTable *table = GetElementTable();
00373 TString decayName;
00374 if (!table) {
00375 Error("CheckDecays", "Element table not present");
00376 return kFALSE;
00377 }
00378 Bool_t resultOK = kTRUE;
00379 if (!fDecays) {
00380 oelem->SetBit(kElementChecked,kTRUE);
00381 return resultOK;
00382 }
00383 Double_t br = 0.;
00384 Int_t decayResult = 0;
00385 TIter next(fDecays);
00386 while ((dc=(TGeoDecayChannel*)next())) {
00387 br += dc->BranchingRatio();
00388 decayResult = DecayResult(dc);
00389 if (decayResult) {
00390 elem = table->GetElementRN(decayResult);
00391 if (!elem) {
00392 TGeoDecayChannel::DecayName(dc->Decay(),decayName);
00393 Error("CheckDecays", "Element after decay %s of %s not found in DB", decayName.Data(),fName.Data());
00394 return kFALSE;
00395 }
00396 dc->SetDaughter(elem);
00397 resultOK = elem->CheckDecays();
00398 }
00399 }
00400 if (TMath::Abs(br-100) > 1.E-3) {
00401 Warning("CheckDecays", "BR for decays of element %s sum-up = %f", fName.Data(), br);
00402 resultOK = kFALSE;
00403 }
00404 oelem->SetBit(kElementChecked,kTRUE);
00405 return resultOK;
00406 }
00407
00408
00409 Int_t TGeoElementRN::DecayResult(TGeoDecayChannel *dc) const
00410 {
00411
00412 Int_t da, dz, diso;
00413 dc->DecayShift(da, dz, diso);
00414 if (da == -99 || dz == -99) return 0;
00415 return ENDF(Int_t(fA)+da,fZ+dz,fIso+diso);
00416 }
00417
00418
00419 void TGeoElementRN::FillPopulation(TObjArray *population, Double_t precision, Double_t factor)
00420 {
00421
00422
00423
00424
00425
00426
00427 TGeoElementRN *elem;
00428 TGeoElemIter next(this, precision);
00429 TGeoBatemanSol s(this);
00430 s.Normalize(factor);
00431 AddRatio(s);
00432 if (!population->FindObject(this)) population->Add(this);
00433 while ((elem=next())) {
00434 TGeoBatemanSol ratio(next.GetBranch());
00435 ratio.Normalize(factor);
00436 elem->AddRatio(ratio);
00437 if (!population->FindObject(elem)) population->Add(elem);
00438 }
00439 }
00440
00441
00442 void TGeoElementRN::MakeName(Int_t a, Int_t z, Int_t iso)
00443 {
00444
00445 fName = "";
00446 if (z==0 && a==1) {
00447 fName = "neutron";
00448 return;
00449 }
00450 if (z>=1 && z<= gMaxElem) fName += TString::Format("%3d-%s-",z,gElName[z-1]);
00451 else fName = "?? -?? -";
00452 if (a>=1 && a<=999) fName += TString::Format("%3.3d",a);
00453 else fName += "??";
00454 if (iso>0 && iso<gMaxLevel) fName += TString::Format("%c", gLevName[iso]);
00455 fName.ReplaceAll(" ","");
00456 }
00457
00458
00459 void TGeoElementRN::Print(Option_t *option) const
00460 {
00461
00462 printf("\n%-12s ",fName.Data());
00463 printf("ENDF=%d; ",fENDFcode);
00464 printf("A=%d; ",(Int_t)fA);
00465 printf("Z=%d; ",fZ);
00466 printf("Iso=%d; ",fIso);
00467 printf("Level=%g[MeV]; ",fLevel);
00468 printf("Dmass=%g[MeV]; ",fDeltaM);
00469 if (fHalfLife>0) printf("Hlife=%g[s]\n",fHalfLife);
00470 else printf("Hlife=INF\n");
00471 printf("%13s"," ");
00472 printf("J/P=%s; ",fTitle.Data());
00473 printf("Abund=%g; ",fNatAbun);
00474 printf("Htox=%g; ",fTH_F);
00475 printf("Itox=%g; ",fTG_F);
00476 printf("Stat=%d\n",fStatus);
00477 if(!fDecays) return;
00478 printf("Decay modes:\n");
00479 TIter next(fDecays);
00480 TGeoDecayChannel *dc;
00481 while ((dc=(TGeoDecayChannel*)next())) dc->Print(option);
00482 }
00483
00484
00485 TGeoElementRN *TGeoElementRN::ReadElementRN(const char *line, Int_t &ndecays)
00486 {
00487
00488 Int_t a,z,iso,status;
00489 Double_t level, deltaM, halfLife, natAbun, th_f, tg_f, th_s, tg_s;
00490 char name[20],jp[20];
00491 sscanf(&line[0], "%s%d%d%d%lg%lg%lg%s%lg%lg%lg%lg%lg%d%d", name,&a,&z,&iso,&level,&deltaM,
00492 &halfLife,jp,&natAbun,&th_f,&tg_f,&th_s,&tg_s,&status,&ndecays);
00493 TGeoElementRN *elem = new TGeoElementRN(a,z,iso,level,deltaM,halfLife,
00494 jp,natAbun,th_f,tg_f,th_s,tg_s,status);
00495 return elem;
00496 }
00497
00498
00499 void TGeoElementRN::SavePrimitive(ostream &out, Option_t *option)
00500 {
00501
00502 if (!strcmp(option,"h")) {
00503
00504 out << "#====================================================================================================================================" << endl;
00505 out << "# Name A Z ISO LEV[MeV] DM[MeV] T1/2[s] J/P ABUND[%] HTOX ITOX HTOX ITOX STAT NDCY" << endl;
00506 out << "#====================================================================================================================================" << endl;
00507 }
00508 out << setw(11) << fName.Data();
00509 out << setw(5) << (Int_t)fA;
00510 out << setw(5) << fZ;
00511 out << setw(5) << fIso;
00512 out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fLevel;
00513 out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fDeltaM;
00514 out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << fHalfLife;
00515 out << setw(13) << fTitle.Data();
00516 out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fNatAbun;
00517 out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fTH_F;
00518 out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fTG_F;
00519 out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fTH_S;
00520 out << setw(10) << setiosflags(ios::fixed) << setprecision(5) << fTG_S;
00521 out << setw(5) << fStatus;
00522 Int_t ndecays = 0;
00523 if (fDecays) ndecays = fDecays->GetEntries();
00524 out << setw(5) << ndecays;
00525 out << endl;
00526 if (fDecays) {
00527 TIter next(fDecays);
00528 TGeoDecayChannel *dc;
00529 while ((dc=(TGeoDecayChannel*)next())) dc->SavePrimitive(out);
00530 }
00531 }
00532
00533
00534 void TGeoElementRN::AddRatio(TGeoBatemanSol &ratio)
00535 {
00536
00537 if (!fRatio) fRatio = new TGeoBatemanSol(ratio);
00538 else *fRatio += ratio;
00539 }
00540
00541
00542 void TGeoElementRN::ResetRatio()
00543 {
00544
00545 if (fRatio) {
00546 delete fRatio;
00547 fRatio = 0;
00548 }
00549 }
00550
00551 ClassImp(TGeoDecayChannel)
00552
00553
00554 TGeoDecayChannel& TGeoDecayChannel::operator=(const TGeoDecayChannel& dc)
00555 {
00556
00557
00558 if(this!=&dc) {
00559 TObject::operator=(dc);
00560 fDecay = dc.fDecay;
00561 fDiso = dc.fDiso;
00562 fBranchingRatio = dc.fBranchingRatio;
00563 fQvalue = dc.fQvalue;
00564 fParent = dc.fParent;
00565 fDaughter = dc.fDaughter;
00566 }
00567 return *this;
00568 }
00569
00570
00571 const char *TGeoDecayChannel::GetName() const
00572 {
00573
00574 static TString name = "";
00575 name = "";
00576 if (!fDecay) return gDecayName[gMaxDecay];
00577 for (Int_t i=0; i<gMaxDecay; i++) {
00578 if (1<<i & fDecay) {
00579 if (name.Length()) name += "+";
00580 name += gDecayName[i];
00581 }
00582 }
00583 return name.Data();
00584 }
00585
00586
00587 void TGeoDecayChannel::DecayName(UInt_t decay, TString &name)
00588 {
00589
00590 if (!decay) {
00591 name = gDecayName[gMaxDecay];
00592 return;
00593 }
00594 name = "";
00595 for (Int_t i=0; i<gMaxDecay; i++) {
00596 if (1<<i & decay) {
00597 if (name.Length()) name += "+";
00598 name += gDecayName[i];
00599 }
00600 }
00601 }
00602
00603
00604 Int_t TGeoDecayChannel::GetIndex() const
00605 {
00606
00607 return fParent->Decays()->IndexOf(this);
00608 }
00609
00610
00611 void TGeoDecayChannel::Print(Option_t *) const
00612 {
00613
00614 TString name;
00615 DecayName(fDecay, name);
00616 printf("%-20s Diso: %3d BR: %9.3f%% Qval: %g\n", name.Data(),fDiso,fBranchingRatio,fQvalue);
00617 }
00618
00619
00620 TGeoDecayChannel *TGeoDecayChannel::ReadDecay(const char *line)
00621 {
00622
00623 char name[80];
00624 Int_t decay,diso;
00625 Double_t branchingRatio, qValue;
00626 sscanf(&line[0], "%s%d%d%lg%lg", name,&decay,&diso,&branchingRatio,&qValue);
00627 TGeoDecayChannel *dc = new TGeoDecayChannel(decay,diso,branchingRatio,qValue);
00628 return dc;
00629 }
00630
00631
00632 void TGeoDecayChannel::SavePrimitive(ostream &out, Option_t *)
00633 {
00634
00635 TString decayName;
00636 DecayName(fDecay, decayName);
00637 out << setw(50) << decayName.Data();
00638 out << setw(10) << fDecay;
00639 out << setw(10) << fDiso;
00640 out << setw(12) << setiosflags(ios::fixed) << setprecision(6) << fBranchingRatio;
00641 out << setw(12) << setiosflags(ios::fixed) << setprecision(6) << fQvalue;
00642 out << endl;
00643 }
00644
00645
00646 void TGeoDecayChannel::DecayShift(Int_t &dA, Int_t &dZ, Int_t &dI) const
00647 {
00648
00649 dA=dZ=0;
00650 dI=fDiso;
00651 for(Int_t i=0; i<gMaxDecay; ++i) {
00652 if(1<<i & fDecay) {
00653 if(gDecayDeltaA[i] == -99 || gDecayDeltaZ[i] == -99 ) {
00654 dA=dZ=-99;
00655 return;
00656 }
00657 dA += gDecayDeltaA[i];
00658 dZ += gDecayDeltaZ[i];
00659 }
00660 }
00661 }
00662
00663 ClassImp(TGeoElemIter)
00664
00665
00666 TGeoElemIter::TGeoElemIter(TGeoElementRN *top, Double_t limit)
00667 : fTop(top), fElem(top), fBranch(0), fLevel(0), fLimitRatio(limit), fRatio(1.)
00668 {
00669
00670 fBranch = new TObjArray(10);
00671 }
00672
00673
00674 TGeoElemIter::TGeoElemIter(const TGeoElemIter &iter)
00675 :fTop(iter.fTop),
00676 fElem(iter.fElem),
00677 fBranch(0),
00678 fLevel(iter.fLevel),
00679 fLimitRatio(iter.fLimitRatio),
00680 fRatio(iter.fRatio)
00681 {
00682
00683 if (iter.fBranch) {
00684 fBranch = new TObjArray(10);
00685 for (Int_t i=0; i<fLevel; i++) fBranch->Add(iter.fBranch->At(i));
00686 }
00687 }
00688
00689
00690 TGeoElemIter::~TGeoElemIter()
00691 {
00692
00693 if (fBranch) delete fBranch;
00694 }
00695
00696
00697 TGeoElemIter &TGeoElemIter::operator=(const TGeoElemIter &iter)
00698 {
00699
00700 if (&iter == this) return *this;
00701 fTop = iter.fTop;
00702 fElem = iter.fElem;
00703 fLevel = iter.fLevel;
00704 if (iter.fBranch) {
00705 fBranch = new TObjArray(10);
00706 for (Int_t i=0; i<fLevel; i++) fBranch->Add(iter.fBranch->At(i));
00707 }
00708 fLimitRatio = iter.fLimitRatio;
00709 fRatio = iter.fRatio;
00710 return *this;
00711 }
00712
00713
00714 TGeoElementRN *TGeoElemIter::operator()()
00715 {
00716
00717 return Next();
00718 }
00719
00720
00721 TGeoElementRN *TGeoElemIter::Up()
00722 {
00723
00724 TGeoDecayChannel *dc;
00725 Int_t ind, nd;
00726 while (fLevel) {
00727
00728 dc = (TGeoDecayChannel*)fBranch->At(fLevel-1);
00729 ind = dc->GetIndex();
00730 nd = dc->Parent()->GetNdecays();
00731 fRatio /= 0.01*dc->BranchingRatio();
00732 fElem = dc->Parent();
00733 fBranch->RemoveAt(--fLevel);
00734 ind++;
00735 while (ind<nd) {
00736 if (Down(ind++)) return (TGeoElementRN*)fElem;
00737 }
00738 }
00739 fElem = NULL;
00740 return NULL;
00741 }
00742
00743
00744 TGeoElementRN *TGeoElemIter::Down(Int_t ibranch)
00745 {
00746
00747
00748 TGeoDecayChannel *dc = (TGeoDecayChannel*)fElem->Decays()->At(ibranch);
00749 if (!dc->Daughter()) return NULL;
00750 Double_t br = 0.01*fRatio*dc->BranchingRatio();
00751 if (br < fLimitRatio) return NULL;
00752 fLevel++;
00753 fRatio = br;
00754 fBranch->Add(dc);
00755 fElem = dc->Daughter();
00756 return (TGeoElementRN*)fElem;
00757 }
00758
00759
00760 TGeoElementRN *TGeoElemIter::Next()
00761 {
00762
00763 if (!fElem) return NULL;
00764
00765 Int_t nd = fElem->GetNdecays();
00766 for (Int_t i=0; i<nd; i++) if (Down(i)) return (TGeoElementRN*)fElem;
00767 return Up();
00768 }
00769
00770
00771 void TGeoElemIter::Print(Option_t * ) const
00772 {
00773
00774 TGeoElementRN *elem;
00775 TGeoDecayChannel *dc;
00776 TString indent = "";
00777 printf("=== Chain with %g %%\n", 100*fRatio);
00778 for (Int_t i=0; i<fLevel; i++) {
00779 dc = (TGeoDecayChannel*)fBranch->At(i);
00780 elem = dc->Parent();
00781 printf("%s%s (%g%% %s) T1/2=%g\n", indent.Data(), elem->GetName(),dc->BranchingRatio(),dc->GetName(),elem->HalfLife());
00782 indent += " ";
00783 if (i==fLevel-1) {
00784 elem = dc->Daughter();
00785 printf("%s%s\n", indent.Data(), elem->GetName());
00786 }
00787 }
00788 }
00789
00790 ClassImp(TGeoElementTable)
00791
00792
00793 TGeoElementTable::TGeoElementTable()
00794 {
00795
00796 fNelements = 0;
00797 fNelementsRN = 0;
00798 fNisotopes = 0;
00799 fList = 0;
00800 fListRN = 0;
00801 fIsotopes = 0;
00802 }
00803
00804
00805 TGeoElementTable::TGeoElementTable(Int_t )
00806 {
00807
00808 fNelements = 0;
00809 fNelementsRN = 0;
00810 fNisotopes = 0;
00811 fList = new TObjArray(128);
00812 fListRN = 0;
00813 fIsotopes = 0;
00814 BuildDefaultElements();
00815
00816 }
00817
00818
00819 TGeoElementTable::TGeoElementTable(const TGeoElementTable& get) :
00820 TObject(get),
00821 fNelements(get.fNelements),
00822 fNelementsRN(get.fNelementsRN),
00823 fNisotopes(get.fNisotopes),
00824 fList(get.fList),
00825 fListRN(get.fListRN),
00826 fIsotopes(0)
00827 {
00828
00829 }
00830
00831
00832 TGeoElementTable& TGeoElementTable::operator=(const TGeoElementTable& get)
00833 {
00834
00835 if(this!=&get) {
00836 TObject::operator=(get);
00837 fNelements=get.fNelements;
00838 fNelementsRN=get.fNelementsRN;
00839 fNisotopes=get.fNisotopes;
00840 fList=get.fList;
00841 fListRN=get.fListRN;
00842 fIsotopes = 0;
00843 }
00844 return *this;
00845 }
00846
00847
00848 TGeoElementTable::~TGeoElementTable()
00849 {
00850
00851 if (fList) {
00852 fList->Delete();
00853 delete fList;
00854 }
00855 if (fListRN) {
00856 fListRN->Delete();
00857 delete fListRN;
00858 }
00859 if (fIsotopes) {
00860 fIsotopes->Delete();
00861 delete fIsotopes;
00862 }
00863 }
00864
00865
00866 void TGeoElementTable::AddElement(const char *name, const char *title, Int_t z, Double_t a)
00867 {
00868
00869 if (!fList) fList = new TObjArray(128);
00870 fList->AddAtAndExpand(new TGeoElement(name,title,z,a), fNelements++);
00871 }
00872
00873
00874 void TGeoElementTable::AddElement(const char *name, const char *title, Int_t z, Int_t n, Double_t a)
00875 {
00876
00877 if (!fList) fList = new TObjArray(128);
00878 fList->AddAtAndExpand(new TGeoElement(name,title,z,n,a), fNelements++);
00879 }
00880
00881
00882 void TGeoElementTable::AddElementRN(TGeoElementRN *elem)
00883 {
00884
00885 if (!fListRN) fListRN = new TObjArray(3600);
00886 if (HasRNElements() && GetElementRN(elem->ENDFCode())) return;
00887
00888 fListRN->Add(elem);
00889 fNelementsRN++;
00890 fElementsRN.insert(ElementRNMap_t::value_type(elem->ENDFCode(), elem));
00891 }
00892
00893
00894 void TGeoElementTable::AddIsotope(TGeoIsotope *isotope)
00895 {
00896
00897 if (FindIsotope(isotope->GetName())) {
00898 Error("AddIsotope", "Isotope with the same name: %s already in table. Not adding.",isotope->GetName());
00899 return;
00900 }
00901 if (!fIsotopes) fIsotopes = new TObjArray();
00902 fIsotopes->Add(isotope);
00903 }
00904
00905
00906 void TGeoElementTable::BuildDefaultElements()
00907 {
00908
00909 if (HasDefaultElements()) return;
00910 AddElement("VACUUM","VACUUM" ,0, 0, 0.0);
00911 AddElement("H" ,"HYDROGEN" ,1, 1, 1.00794);
00912 AddElement("HE" ,"HELIUM" ,2, 4, 4.002602);
00913 AddElement("LI" ,"LITHIUM" ,3, 7, 6.941);
00914 AddElement("BE" ,"BERYLLIUM" ,4, 9, 9.01218);
00915 AddElement("B" ,"BORON" ,5, 11, 10.811);
00916 AddElement("C" ,"CARBON" ,6, 12, 12.0107);
00917 AddElement("N" ,"NITROGEN" ,7, 14, 14.00674);
00918 AddElement("O" ,"OXYGEN" ,8, 16, 15.9994);
00919 AddElement("F" ,"FLUORINE" ,9, 19, 18.9984032);
00920 AddElement("NE" ,"NEON" ,10, 20, 20.1797);
00921 AddElement("NA" ,"SODIUM" ,11, 23, 22.989770);
00922 AddElement("MG" ,"MAGNESIUM" ,12, 24, 24.3050);
00923 AddElement("AL" ,"ALUMINIUM" ,13, 27, 26.981538);
00924 AddElement("SI" ,"SILICON" ,14, 28, 28.0855);
00925 AddElement("P" ,"PHOSPHORUS" ,15, 31, 30.973761);
00926 AddElement("S" ,"SULFUR" ,16, 32, 32.066);
00927 AddElement("CL" ,"CHLORINE" ,17, 35, 35.4527);
00928 AddElement("AR" ,"ARGON" ,18, 40, 39.948);
00929 AddElement("K" ,"POTASSIUM" ,19, 39, 39.0983);
00930 AddElement("CA" ,"CALCIUM" ,20, 40, 40.078);
00931 AddElement("SC" ,"SCANDIUM" ,21, 45, 44.955910);
00932 AddElement("TI" ,"TITANIUM" ,22, 48, 47.867);
00933 AddElement("V" ,"VANADIUM" ,23, 51, 50.9415);
00934 AddElement("CR" ,"CHROMIUM" ,24, 52, 51.9961);
00935 AddElement("MN" ,"MANGANESE" ,25, 55, 54.938049);
00936 AddElement("FE" ,"IRON" ,26, 56, 55.845);
00937 AddElement("CO" ,"COBALT" ,27, 59, 58.933200);
00938 AddElement("NI" ,"NICKEL" ,28, 59, 58.6934);
00939 AddElement("CU" ,"COPPER" ,29, 64, 63.546);
00940 AddElement("ZN" ,"ZINC" ,30, 65, 65.39);
00941 AddElement("GA" ,"GALLIUM" ,31, 70, 69.723);
00942 AddElement("GE" ,"GERMANIUM" ,32, 73, 72.61);
00943 AddElement("AS" ,"ARSENIC" ,33, 75, 74.92160);
00944 AddElement("SE" ,"SELENIUM" ,34, 79, 78.96);
00945 AddElement("BR" ,"BROMINE" ,35, 80, 79.904);
00946 AddElement("KR" ,"KRYPTON" ,36, 84, 83.80);
00947 AddElement("RB" ,"RUBIDIUM" ,37, 85, 85.4678);
00948 AddElement("SR" ,"STRONTIUM" ,38, 88, 87.62);
00949 AddElement("Y" ,"YTTRIUM" ,39, 89, 88.90585);
00950 AddElement("ZR" ,"ZIRCONIUM" ,40, 91, 91.224);
00951 AddElement("NB" ,"NIOBIUM" ,41, 93, 92.90638);
00952 AddElement("MO" ,"MOLYBDENUM" ,42, 96, 95.94);
00953 AddElement("TC" ,"TECHNETIUM" ,43, 98, 98.0);
00954 AddElement("RU" ,"RUTHENIUM" ,44, 101, 101.07);
00955 AddElement("RH" ,"RHODIUM" ,45, 103, 102.90550);
00956 AddElement("PD" ,"PALLADIUM" ,46, 106, 106.42);
00957 AddElement("AG" ,"SILVER" ,47, 108, 107.8682);
00958 AddElement("CD" ,"CADMIUM" ,48, 112, 112.411);
00959 AddElement("IN" ,"INDIUM" ,49, 115, 114.818);
00960 AddElement("SN" ,"TIN" ,50, 119, 118.710);
00961 AddElement("SB" ,"ANTIMONY" ,51, 122, 121.760);
00962 AddElement("TE" ,"TELLURIUM" ,52, 128, 127.60);
00963 AddElement("I" ,"IODINE" ,53, 127, 126.90447);
00964 AddElement("XE" ,"XENON" ,54, 131, 131.29);
00965 AddElement("CS" ,"CESIUM" ,55, 133, 132.90545);
00966 AddElement("BA" ,"BARIUM" ,56, 137, 137.327);
00967 AddElement("LA" ,"LANTHANUM" ,57, 139, 138.9055);
00968 AddElement("CE" ,"CERIUM" ,58, 140, 140.116);
00969 AddElement("PR" ,"PRASEODYMIUM" ,59, 141, 140.90765);
00970 AddElement("ND" ,"NEODYMIUM" ,60, 144, 144.24);
00971 AddElement("PM" ,"PROMETHIUM" ,61, 145, 145.0);
00972 AddElement("SM" ,"SAMARIUM" ,62, 150, 150.36);
00973 AddElement("EU" ,"EUROPIUM" ,63, 152, 151.964);
00974 AddElement("GD" ,"GADOLINIUM" ,64, 157, 157.25);
00975 AddElement("TB" ,"TERBIUM" ,65, 159, 158.92534);
00976 AddElement("DY" ,"DYSPROSIUM" ,66, 162, 162.50);
00977 AddElement("HO" ,"HOLMIUM" ,67, 165, 164.93032);
00978 AddElement("ER" ,"ERBIUM" ,68, 167, 167.26);
00979 AddElement("TM" ,"THULIUM" ,69, 169, 168.93421);
00980 AddElement("YB" ,"YTTERBIUM" ,70, 173, 173.04);
00981 AddElement("LU" ,"LUTETIUM" ,71, 175, 174.967);
00982 AddElement("HF" ,"HAFNIUM" ,72, 178, 178.49);
00983 AddElement("TA" ,"TANTALUM" ,73, 181, 180.9479);
00984 AddElement("W" ,"TUNGSTEN" ,74, 184, 183.84);
00985 AddElement("RE" ,"RHENIUM" ,75, 186, 186.207);
00986 AddElement("OS" ,"OSMIUM" ,76, 190, 190.23);
00987 AddElement("IR" ,"IRIDIUM" ,77, 192, 192.217);
00988 AddElement("PT" ,"PLATINUM" ,78, 195, 195.078);
00989 AddElement("AU" ,"GOLD" ,79, 197, 196.96655);
00990 AddElement("HG" ,"MERCURY" ,80, 200, 200.59);
00991 AddElement("TL" ,"THALLIUM" ,81, 204, 204.3833);
00992 AddElement("PB" ,"LEAD" ,82, 207, 207.2);
00993 AddElement("BI" ,"BISMUTH" ,83, 209, 208.98038);
00994 AddElement("PO" ,"POLONIUM" ,84, 209, 209.0);
00995 AddElement("AT" ,"ASTATINE" ,85, 210, 210.0);
00996 AddElement("RN" ,"RADON" ,86, 222, 222.0);
00997 AddElement("FR" ,"FRANCIUM" ,87, 223, 223.0);
00998 AddElement("RA" ,"RADIUM" ,88, 226, 226.0);
00999 AddElement("AC" ,"ACTINIUM" ,89, 227, 227.0);
01000 AddElement("TH" ,"THORIUM" ,90, 232, 232.0381);
01001 AddElement("PA" ,"PROTACTINIUM" ,91, 231, 231.03588);
01002 AddElement("U" ,"URANIUM" ,92, 238, 238.0289);
01003 AddElement("NP" ,"NEPTUNIUM" ,93, 237, 237.0);
01004 AddElement("PU" ,"PLUTONIUM" ,94, 244, 244.0);
01005 AddElement("AM" ,"AMERICIUM" ,95, 243, 243.0);
01006 AddElement("CM" ,"CURIUM" ,96, 247, 247.0);
01007 AddElement("BK" ,"BERKELIUM" ,97, 247, 247.0);
01008 AddElement("CF" ,"CALIFORNIUM",98, 251, 251.0);
01009 AddElement("ES" ,"EINSTEINIUM",99, 252, 252.0);
01010 AddElement("FM" ,"FERMIUM" ,100, 257, 257.0);
01011 AddElement("MD" ,"MENDELEVIUM",101, 258, 258.0);
01012 AddElement("NO" ,"NOBELIUM" ,102, 259, 259.0);
01013 AddElement("LR" ,"LAWRENCIUM" ,103, 262, 262.0);
01014 AddElement("RF" ,"RUTHERFORDIUM",104, 261, 261.0);
01015 AddElement("DB" ,"DUBNIUM" ,105, 262, 262.0);
01016 AddElement("SG" ,"SEABORGIUM" ,106, 263, 263.0);
01017 AddElement("BH" ,"BOHRIUM" ,107, 262, 262.0);
01018 AddElement("HS" ,"HASSIUM" ,108, 265, 265.0);
01019 AddElement("MT" ,"MEITNERIUM" ,109, 266, 266.0);
01020 AddElement("UUN" ,"UNUNNILIUM" ,110, 269, 269.0);
01021 AddElement("UUU" ,"UNUNUNIUM" ,111, 272, 272.0);
01022 AddElement("UUB" ,"UNUNBIUM" ,112, 277, 277.0);
01023
01024 TObject::SetBit(kETDefaultElements,kTRUE);
01025 }
01026
01027
01028 void TGeoElementTable::ImportElementsRN()
01029 {
01030
01031 if (HasRNElements()) return;
01032 TGeoElementRN *elem;
01033 TString rnf;
01034 #ifdef ROOTETCDIR
01035 rnf.Form("%s/RadioNuclides.txt", ROOTETCDIR);
01036 #else
01037 rnf.Form("%s/etc/RadioNuclides.txt", gSystem->Getenv("ROOTSYS"));
01038 #endif
01039 FILE *fp = fopen(rnf, "r");
01040 if (!fp) {
01041 Error("ImportElementsRN","File RadioNuclides.txt not found");
01042 return;
01043 }
01044 char line[150];
01045 Int_t ndecays = 0;
01046 Int_t i;
01047 while (fgets(&line[0],140,fp)) {
01048 if (line[0]=='#') continue;
01049 elem = TGeoElementRN::ReadElementRN(line, ndecays);
01050 for (i=0; i<ndecays; i++) {
01051 if (!fgets(&line[0],140,fp)) {
01052 Error("ImportElementsRN", "Error parsing RadioNuclides.txt file");
01053 fclose(fp);
01054 return;
01055 }
01056 TGeoDecayChannel *dc = TGeoDecayChannel::ReadDecay(line);
01057 elem->AddDecay(dc);
01058 }
01059 AddElementRN(elem);
01060
01061 }
01062 TObject::SetBit(kETRNElements,kTRUE);
01063 CheckTable();
01064 fclose(fp);
01065 }
01066
01067
01068 Bool_t TGeoElementTable::CheckTable() const
01069 {
01070
01071 if (!HasRNElements()) return HasDefaultElements();
01072 TGeoElementRN *elem;
01073 Bool_t result = kTRUE;
01074 TIter next(fListRN);
01075 while ((elem=(TGeoElementRN*)next())) {
01076 if (!elem->CheckDecays()) result = kFALSE;
01077 }
01078 return result;
01079 }
01080
01081
01082 void TGeoElementTable::ExportElementsRN(const char *filename)
01083 {
01084
01085 if (!HasRNElements()) return;
01086 TString sname = filename;
01087 if (!sname.Length()) sname = "RadioNuclides.txt";
01088 ofstream out;
01089 out.open(sname.Data(), ios::out);
01090 if (!out.good()) {
01091 Error("ExportElementsRN", "Cannot open file %s", sname.Data());
01092 return;
01093 }
01094
01095 TGeoElementRN *elem;
01096 TIter next(fListRN);
01097 Int_t i=0;
01098 while ((elem=(TGeoElementRN*)next())) {
01099 if ((i%48)==0) elem->SavePrimitive(out,"h");
01100 else elem->SavePrimitive(out);
01101 i++;
01102 }
01103 out.close();
01104 }
01105
01106
01107 TGeoElement *TGeoElementTable::FindElement(const char *name) const
01108 {
01109
01110 TString s(name);
01111 s.ToUpper();
01112 TGeoElement *elem;
01113 elem = (TGeoElement*)fList->FindObject(s.Data());
01114 if (elem) return elem;
01115
01116 TIter next(fList);
01117 while ((elem=(TGeoElement*)next())) {
01118 if (s == elem->GetTitle()) return elem;
01119 }
01120 return 0;
01121 }
01122
01123
01124 TGeoIsotope *TGeoElementTable::FindIsotope(const char *name) const
01125 {
01126
01127 if (!fIsotopes) return NULL;
01128 return (TGeoIsotope*)fIsotopes->FindObject(name);
01129 }
01130
01131
01132 TGeoElementRN *TGeoElementTable::GetElementRN(Int_t ENDFcode) const
01133 {
01134
01135 if (!HasRNElements()) {
01136 TGeoElementTable *table = (TGeoElementTable*)this;
01137 table->ImportElementsRN();
01138 if (!fListRN) return 0;
01139 }
01140 ElementRNMap_t::const_iterator it = fElementsRN.find(ENDFcode);
01141 if (it != fElementsRN.end()) return it->second;
01142 return 0;
01143 }
01144
01145
01146 TGeoElementRN *TGeoElementTable::GetElementRN(Int_t a, Int_t z, Int_t iso) const
01147 {
01148
01149 return GetElementRN(TGeoElementRN::ENDF(a,z,iso));
01150 }
01151
01152 ClassImp(TGeoBatemanSol)
01153
01154
01155 TGeoBatemanSol::TGeoBatemanSol(TGeoElementRN *elem)
01156 :TObject(), TAttLine(), TAttFill(), TAttMarker(),
01157 fElem(elem),
01158 fElemTop(elem),
01159 fCsize(10),
01160 fNcoeff(0),
01161 fFactor(1.),
01162 fTmin(0.),
01163 fTmax(0.),
01164 fCoeff(NULL)
01165 {
01166
01167 fCoeff = new BtCoef_t[fCsize];
01168 fNcoeff = 1;
01169 fCoeff[0].cn = 1.;
01170 Double_t t12 = elem->HalfLife();
01171 if (t12 == 0.) t12 = 1.e-30;
01172 if (elem->Stable()) fCoeff[0].lambda = 0.;
01173 else fCoeff[0].lambda = TMath::Log(2.)/t12;
01174 }
01175
01176
01177 TGeoBatemanSol::TGeoBatemanSol(const TObjArray *chain)
01178 :TObject(), TAttLine(), TAttFill(), TAttMarker(),
01179 fElem(NULL),
01180 fElemTop(NULL),
01181 fCsize(0),
01182 fNcoeff(0),
01183 fFactor(1.),
01184 fTmin(0.),
01185 fTmax(0.),
01186 fCoeff(NULL)
01187 {
01188
01189 TGeoDecayChannel *dc = (TGeoDecayChannel*)chain->At(0);
01190 if (dc) fElemTop = dc->Parent();
01191 dc = (TGeoDecayChannel*)chain->At(chain->GetEntriesFast()-1);
01192 if (dc) {
01193 fElem = dc->Daughter();
01194 fCsize = chain->GetEntriesFast()+1;
01195 fCoeff = new BtCoef_t[fCsize];
01196 FindSolution(chain);
01197 }
01198 }
01199
01200
01201 TGeoBatemanSol::TGeoBatemanSol(const TGeoBatemanSol& other)
01202 :TObject(other), TAttLine(other), TAttFill(other), TAttMarker(other),
01203 fElem(other.fElem),
01204 fElemTop(other.fElemTop),
01205 fCsize(other.fCsize),
01206 fNcoeff(other.fNcoeff),
01207 fFactor(other.fFactor),
01208 fTmin(other.fTmin),
01209 fTmax(other.fTmax),
01210 fCoeff(NULL)
01211 {
01212
01213 if (fCsize) {
01214 fCoeff = new BtCoef_t[fCsize];
01215 for (Int_t i=0; i<fNcoeff; i++) {
01216 fCoeff[i].cn = other.fCoeff[i].cn;
01217 fCoeff[i].lambda = other.fCoeff[i].lambda;
01218 }
01219 }
01220 }
01221
01222
01223 TGeoBatemanSol::~TGeoBatemanSol()
01224 {
01225
01226 if (fCoeff) delete [] fCoeff;
01227 }
01228
01229
01230 TGeoBatemanSol& TGeoBatemanSol::operator=(const TGeoBatemanSol& other)
01231 {
01232
01233 if (this == &other) return *this;
01234 TObject::operator=(other);
01235 TAttLine::operator=(other);
01236 TAttFill::operator=(other);
01237 TAttMarker::operator=(other);
01238 fElem = other.fElem;
01239 fElemTop = other.fElemTop;
01240 if (fCoeff) delete [] fCoeff;
01241 fCoeff = 0;
01242 fCsize = other.fCsize;
01243 fNcoeff = other.fNcoeff;
01244 fFactor = other.fFactor;
01245 fTmin = other.fTmin;
01246 fTmax = other.fTmax;
01247 if (fCsize) {
01248 fCoeff = new BtCoef_t[fCsize];
01249 for (Int_t i=0; i<fNcoeff; i++) {
01250 fCoeff[i].cn = other.fCoeff[i].cn;
01251 fCoeff[i].lambda = other.fCoeff[i].lambda;
01252 }
01253 }
01254 return *this;
01255 }
01256
01257
01258 TGeoBatemanSol& TGeoBatemanSol::operator+=(const TGeoBatemanSol& other)
01259 {
01260
01261 if (other.GetElement() != fElem) {
01262 Error("operator+=", "Cannot add 2 solutions for different elements");
01263 return *this;
01264 }
01265 Int_t i,j;
01266 BtCoef_t *coeff = fCoeff;
01267 Int_t ncoeff = fNcoeff + other.fNcoeff;
01268 if (ncoeff > fCsize) {
01269 fCsize = ncoeff;
01270 coeff = new BtCoef_t[ncoeff];
01271 for (i=0; i<fNcoeff; i++) {
01272 coeff[i].cn = fCoeff[i].cn;
01273 coeff[i].lambda = fCoeff[i].lambda;
01274 }
01275 delete [] fCoeff;
01276 fCoeff = coeff;
01277 }
01278 ncoeff = fNcoeff;
01279 for (j=0; j<other.fNcoeff; j++) {
01280 for (i=0; i<fNcoeff; i++) {
01281 if (coeff[i].lambda == other.fCoeff[j].lambda) {
01282 coeff[i].cn += fFactor * other.fCoeff[j].cn;
01283 break;
01284 }
01285 }
01286 if (i == fNcoeff) {
01287 coeff[ncoeff].cn = fFactor * other.fCoeff[j].cn;
01288 coeff[ncoeff].lambda = other.fCoeff[j].lambda;
01289 ncoeff++;
01290 }
01291 }
01292 fNcoeff = ncoeff;
01293 return *this;
01294 }
01295
01296 Double_t TGeoBatemanSol::Concentration(Double_t time) const
01297 {
01298
01299 Double_t conc = 0.;
01300 for (Int_t i=0; i<fNcoeff; i++)
01301 conc += fCoeff[i].cn * TMath::Exp(-fCoeff[i].lambda * time);
01302 return conc;
01303 }
01304
01305
01306 void TGeoBatemanSol::Draw(Option_t *option)
01307 {
01308
01309 if (!gGeoManager) return;
01310 gGeoManager->GetGeomPainter()->DrawBatemanSol(this, option);
01311 }
01312
01313
01314 void TGeoBatemanSol::FindSolution(const TObjArray *array)
01315 {
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325 fNcoeff = 0;
01326 if (!array || !array->GetEntriesFast()) return;
01327 Int_t n = array->GetEntriesFast();
01328 TGeoDecayChannel *dc = (TGeoDecayChannel*)array->At(n-1);
01329 TGeoElementRN *elem = dc->Daughter();
01330 if (elem != fElem) {
01331 Error("FindSolution", "Last element in the list must be %s\n", fElem->GetName());
01332 return;
01333 }
01334 Int_t i,j;
01335 Int_t order = n+1;
01336 if (!fCoeff) {
01337 fCsize = order;
01338 fCoeff = new BtCoef_t[fCsize];
01339 }
01340 if (fCsize < order) {
01341 delete [] fCoeff;
01342 fCsize = order;
01343 fCoeff = new BtCoef_t[fCsize];
01344 }
01345
01346 Double_t *lambda = new Double_t[order];
01347 Double_t *br = new Double_t[n];
01348 Double_t halflife;
01349 for (i=0; i<n; i++) {
01350 dc = (TGeoDecayChannel*)array->At(i);
01351 elem = dc->Parent();
01352 br[i] = 0.01 * dc->BranchingRatio();
01353 halflife = elem->HalfLife();
01354 if (halflife==0.) halflife = 1.e-30;
01355 if (elem->Stable()) lambda[i] = 0.;
01356 else lambda[i] = TMath::Log(2.)/halflife;
01357 if (i==n-1) {
01358 elem = dc->Daughter();
01359 halflife = elem->HalfLife();
01360 if (halflife==0.) halflife = 1.e-30;
01361 if (elem->Stable()) lambda[n] = 0.;
01362 else lambda[n] = TMath::Log(2.)/halflife;
01363 }
01364 }
01365
01366 for (i=0; i<order-1; i++) {
01367 for (j=i+1; j<order; j++) {
01368 if (lambda[j] == lambda[i]) lambda[j] += 0.001*lambda[j];
01369 }
01370 }
01371 Double_t ain;
01372 Double_t pdlambda, plambdabr=1.;
01373 for (j=0; j<n; j++) plambdabr *= lambda[j]*br[j];
01374 for (i=0; i<order; i++) {
01375 pdlambda = 1.;
01376 for (j=0; j<n+1; j++) {
01377 if (j == i) continue;
01378 pdlambda *= lambda[j] - lambda[i];
01379 }
01380 if (pdlambda == 0.) {
01381 Error("FindSolution", "pdlambda=0 !!!");
01382 delete [] lambda;
01383 delete [] br;
01384 return;
01385 }
01386 ain = plambdabr/pdlambda;
01387 fCoeff[i].cn = ain;
01388 fCoeff[i].lambda = lambda[i];
01389 }
01390 fNcoeff = order;
01391 Normalize(fFactor);
01392 delete [] lambda;
01393 delete [] br;
01394 }
01395
01396
01397 void TGeoBatemanSol::Normalize(Double_t factor)
01398 {
01399
01400 for (Int_t i=0; i<fNcoeff; i++) fCoeff[i].cn *= factor;
01401 }
01402
01403
01404 void TGeoBatemanSol::Print(Option_t * ) const
01405 {
01406
01407 TString formula;
01408 formula.Form("N[%s]/N[%s] = ", fElem->GetName(), fElemTop->GetName());
01409 for (Int_t i=0; i<fNcoeff; i++) {
01410 if (i == fNcoeff-1) formula += TString::Format("%g*exp(-%g*t)", fCoeff[i].cn, fCoeff[i].lambda);
01411 else formula += TString::Format("%g*exp(-%g*t) + ", fCoeff[i].cn, fCoeff[i].lambda);
01412 }
01413 printf("%s\n", formula.Data());
01414 }
01415
