RadioNuclides.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_geom
/// Macro that demonstrates usage of radioactive elements/materials/mixtures
/// with TGeo package.
///
/// A radionuclide (TGeoElementRN) derives from the class TGeoElement and
/// provides additional information related to its radioactive properties and
/// decay modes.
///
/// The radionuclides table is loaded on demand by any call:
///
/// ~~~{.cpp}
///    TGeoElementRN *TGeoElementTable::GetElementRN(Int_t atomic_number,
///                                                  Int_t atomic_charge,
///                                                  Int_t isomeric_number)
/// ~~~
///
/// The isomeric number is optional and the default value is 0.
///
/// To create a radioactive material based on a radionuclide, one should use the
/// constructor:
///
/// ~~~{.cpp}
///    TGeoMaterial(const char *name, TGeoElement *elem, Double_t density)
/// ~~~
///
/// To create a radioactive mixture, one can use radionuclides as well as stable
/// elements:
///
/// ~~~{.cpp}
///    TGeoMixture(const char *name, Int_t nelements, Double_t density);
///    TGeoMixture::AddElement(TGeoElement *elem, Double_t weight_fraction);
/// ~~~
///
/// Once defined, one can retrieve the time evolution for the radioactive
/// materials/mixtures by using one of the 2 methods:
///
/// ~~~{.cpp}
///    void TGeoMaterial::FillMaterialEvolution(TObjArray *population,
///                                             Double_t   precision=0.001)
/// ~~~
///
/// To use this method, one has to provide an empty TObjArray object that will
/// be filled with all elements coming from the decay chain of the initial
/// radionuclides contained by the material/mixture. The precision represent the
/// cumulative branching ratio for which decay products are still considered.
/// The POPULATION list may contain stable elements as well as radionuclides,
/// depending on the initial elements. To test if an element is a radionuclide:
///
/// ~~~{.cpp}
///    Bool_t TGeoElement::IsRadioNuclide() const
/// ~~~
///
/// All radionuclides in the output population list have attached objects that
/// represent the time evolution of their fraction of nuclei with respect to the
/// top radionuclide in the decay chain. These objects (Bateman solutions) can be
/// retrieved and drawn:
///
/// ~~~{.cpp}
///    TGeoBatemanSol *TGeoElementRN::Ratio();
///    void TGeoBatemanSol::Draw();
/// ~~~
///
/// Another method allows to create the evolution of a given radioactive
/// material/mixture at a given moment in time:
///
/// ~~~{.cpp}
///    TGeoMaterial::DecayMaterial(Double_t time, Double_t precision=0.001)
/// ~~~
///
/// The method will create the mixture that result from the decay of a initial
/// material/mixture at TIME, while all resulting elements having a fractional
/// weight less than PRECISION are excluded.
///
/// \macro_image
/// \macro_code
///
/// \author Mihaela Gheata
 
void DrawPopulation(TObjArray *vect, TCanvas *can, Double_t tmin = 0., Double_t tmax = 0., Bool_t logx = kFALSE);
 
void RadioNuclides()
{
   TGeoManager *geom = new TGeoManager("", "");
   TGeoElementTable *table = gGeoManager->GetElementTable();
   TGeoElementRN *c14 = table->GetElementRN(14, 6);
   TGeoElementRN *el1 = table->GetElementRN(53, 20);
   TGeoElementRN *el2 = table->GetElementRN(78, 38);
   // Radioactive material
   TGeoMaterial *mat = new TGeoMaterial("C14", c14, 1.3);
   printf("___________________________________________________________\n");
   printf("Radioactive material:\n");
   mat->Print();
   Double_t time = 1.5e11; // seconds
   TGeoMaterial *decaymat = mat->DecayMaterial(time);
   printf("Radioactive material evolution after %g years:\n", time / 3.1536e7);
   decaymat->Print();
   // Radioactive mixture
   TGeoMixture *mix = new TGeoMixture("mix", 2, 7.3);
   mix->AddElement(el1, 0.35);
   mix->AddElement(el2, 0.65);
   printf("___________________________________________________________\n");
   printf("Radioactive mixture:\n");
   mix->Print();
   time = 1000.;
   decaymat = mix->DecayMaterial(time);
   printf("Radioactive mixture evolution after %g seconds:\n", time);
   decaymat->Print();
   TObjArray *vect = new TObjArray();
   TCanvas *c1 = new TCanvas("c1", "C14 decay", 800, 600);
   c1->SetGrid();
   mat->FillMaterialEvolution(vect);
   DrawPopulation(vect, c1, 0, 1.4e12);
   TLatex *tex = new TLatex(8.35e11, 0.564871, "C_{N^{14}_{7}}");
   tex->SetTextSize(0.0388601);
   tex->SetLineWidth(2);
   tex->Draw();
   tex = new TLatex(3.33e11, 0.0620678, "C_{C^{14}_{6}}");
   tex->SetTextSize(0.0388601);
   tex->SetLineWidth(2);
   tex->Draw();
   tex = new TLatex(9.4e11, 0.098, "C_{X}=#frac{N_{X}(t)}{N_{0}(t=0)}=\
   #sum_{j}#alpha_{j}e^{-#lambda_{j}t}");
   tex->SetTextSize(0.0388601);
   tex->SetLineWidth(2);
   tex->Draw();
   TPaveText *pt = new TPaveText(2.6903e+11, 0.0042727, 1.11791e+12, 0.0325138, "br");
   pt->SetFillColor(5);
   pt->SetTextAlign(12);
   pt->SetTextColor(4);
   pt->AddText("Time evolution of a population of radionuclides.");
   pt->AddText("The concentration of a nuclide X represent the  ");
   pt->AddText("ratio between the number of X nuclei and the    ");
   pt->AddText("number of nuclei of the top element of the decay");
   pt->AddText("from which X derives from at T=0.               ");
   pt->Draw();
   c1->Modified();
   vect->Clear();
   TCanvas *c2 = new TCanvas("c2", "Mixture decay", 1000, 800);
   c2->SetGrid();
   mix->FillMaterialEvolution(vect);
   DrawPopulation(vect, c2, 0.01, 1000., kTRUE);
   tex = new TLatex(0.019, 0.861, "C_{Ca^{53}_{20}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(1);
   tex->Draw();
   tex = new TLatex(0.0311, 0.078064, "C_{Sc^{52}_{21}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(2);
   tex->Draw();
   tex = new TLatex(0.1337, 0.010208, "C_{Ti^{52}_{22}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(3);
   tex->Draw();
   tex = new TLatex(1.54158, 0.00229644, "C_{V^{52}_{23}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(4);
   tex->Draw();
   tex = new TLatex(25.0522, 0.00135315, "C_{Cr^{52}_{24}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(5);
   tex->Draw();
   tex = new TLatex(0.1056, 0.5429, "C_{Sc^{53}_{21}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(6);
   tex->Draw();
   tex = new TLatex(0.411, 0.1044, "C_{Ti^{53}_{22}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(7);
   tex->Draw();
   tex = new TLatex(2.93358, 0.0139452, "C_{V^{53}_{23}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(8);
   tex->Draw();
   tex = new TLatex(10.6235, 0.00440327, "C_{Cr^{53}_{24}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(9);
   tex->Draw();
   tex = new TLatex(15.6288, 0.782976, "C_{Sr^{78}_{38}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(1);
   tex->Draw();
   tex = new TLatex(20.2162, 0.141779, "C_{Rb^{78}_{37}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(2);
   tex->Draw();
   tex = new TLatex(32.4055, 0.0302101, "C_{Kr^{78}_{36}}");
   tex->SetTextSize(0.0388601);
   tex->SetTextColor(3);
   tex->Draw();
   tex = new TLatex(117., 1.52, "C_{X}=#frac{N_{X}(t)}{N_{0}(t=0)}=#sum_{j}\
   #alpha_{j}e^{-#lambda_{j}t}");
   tex->SetTextSize(0.03);
   tex->SetLineWidth(2);
   tex->Draw();
   TArrow *arrow = new TArrow(0.0235313, 0.74106, 0.0385371, 0.115648, 0.02, ">");
   arrow->SetFillColor(1);
   arrow->SetFillStyle(1001);
   arrow->SetLineWidth(2);
   arrow->SetAngle(30);
   arrow->Draw();
   arrow = new TArrow(0.0543138, 0.0586338, 0.136594, 0.0146596, 0.02, ">");
   arrow->SetFillColor(1);
   arrow->SetFillStyle(1001);
   arrow->SetLineWidth(2);
   arrow->SetAngle(30);
   arrow->Draw();
   arrow = new TArrow(0.31528, 0.00722919, 1.29852, 0.00306079, 0.02, ">");
   arrow->SetFillColor(1);
   arrow->SetFillStyle(1001);
   arrow->SetLineWidth(2);
   arrow->SetAngle(30);
   arrow->Draw();
   arrow = new TArrow(4.13457, 0.00201942, 22.5047, 0.00155182, 0.02, ">");
   arrow->SetFillColor(1);
   arrow->SetFillStyle(1001);
   arrow->SetLineWidth(2);
   arrow->SetAngle(30);
   arrow->Draw();
   arrow = new TArrow(0.0543138, 0.761893, 0.0928479, 0.67253, 0.02, ">");
   arrow->SetFillColor(1);
   arrow->SetFillStyle(1001);
   arrow->SetLineWidth(2);
   arrow->SetAngle(30);
   arrow->Draw();
   arrow = new TArrow(0.238566, 0.375717, 0.416662, 0.154727, 0.02, ">");
   arrow->SetFillColor(1);
   arrow->SetFillStyle(1001);
   arrow->SetLineWidth(2);
   arrow->SetAngle(30);
   arrow->Draw();
   arrow = new TArrow(0.653714, 0.074215, 2.41863, 0.0213142, 0.02, ">");
   arrow->SetFillColor(1);
   arrow->SetFillStyle(1001);
   arrow->SetLineWidth(2);
   arrow->SetAngle(30);
   arrow->Draw();
   arrow = new TArrow(5.58256, 0.00953882, 10.6235, 0.00629343, 0.02, ">");
   arrow->SetFillColor(1);
   arrow->SetFillStyle(1001);
   arrow->SetLineWidth(2);
   arrow->SetAngle(30);
   arrow->Draw();
   arrow = new TArrow(22.0271, 0.601935, 22.9926, 0.218812, 0.02, ">");
   arrow->SetFillColor(1);
   arrow->SetFillStyle(1001);
   arrow->SetLineWidth(2);
   arrow->SetAngle(30);
   arrow->Draw();
   arrow = new TArrow(27.2962, 0.102084, 36.8557, 0.045686, 0.02, ">");
   arrow->SetFillColor(1);
   arrow->SetFillStyle(1001);
   arrow->SetLineWidth(2);
   arrow->SetAngle(30);
   arrow->Draw();
}
 
void DrawPopulation(TObjArray *vect, TCanvas *can, Double_t tmin, Double_t tmax, Bool_t logx)
{
   Int_t n = vect->GetEntriesFast();
   TGeoElementRN *elem;
   TGeoBatemanSol *sol;
   can->SetLogy();
 
   if (logx)
      can->SetLogx();
 
   for (Int_t i = 0; i < n; i++) {
      TGeoElement *el = (TGeoElement *)vect->At(i);
      if (!el->IsRadioNuclide())
         continue;
      TGeoElementRN *elem = (TGeoElementRN *)el;
      TGeoBatemanSol *sol = elem->Ratio();
      if (sol) {
         sol->SetLineColor(1 + (i % 9));
         sol->SetLineWidth(2);
         if (tmax > 0.)
            sol->SetRange(tmin, tmax);
         if (i == 0) {
            sol->Draw();
            TF1 *func = (TF1 *)can->FindObject(Form("conc%s", sol->GetElement()->GetName()));
            if (func) {
               if (!strcmp(can->GetName(), "c1"))
                  func->SetTitle("Concentration of C14 derived elements;time[s];Ni/N0(C14)");
               else
                  func->SetTitle("Concentration of elements derived from mixture Ca53+Sr78;\
                  time[s];Ni/N0(Ca53)");
            }
         } else
            sol->Draw("SAME");
      }
   }
}