Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
RadioNuclides.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_geom
3 /// Macro that demonstrates usage of radioactive elements/materials/mixtures
4 /// with TGeo package.
5 ///
6 /// A radionuclide (TGeoElementRN) derives from the class TGeoElement and
7 /// provides additional information related to its radioactive properties and
8 /// decay modes.
9 ///
10 /// The radionuclides table is loaded on demand by any call:
11 ///
12 /// ~~~{.cpp}
13 /// TGeoElementRN *TGeoElementTable::GetElementRN(Int_t atomic_number,
14 /// Int_t atomic_charge,
15 /// Int_t isomeric_number)
16 /// ~~~
17 ///
18 /// The isomeric number is optional and the default value is 0.
19 ///
20 /// To create a radioactive material based on a radionuclide, one should use the
21 /// constructor:
22 ///
23 /// ~~~{.cpp}
24 /// TGeoMaterial(const char *name, TGeoElement *elem, Double_t density)
25 /// ~~~
26 ///
27 /// To create a radioactive mixture, one can use radionuclides as well as stable
28 /// elements:
29 ///
30 /// ~~~{.cpp}
31 /// TGeoMixture(const char *name, Int_t nelements, Double_t density);
32 /// TGeoMixture::AddElement(TGeoElement *elem, Double_t weight_fraction);
33 /// ~~~
34 ///
35 /// Once defined, one can retrieve the time evolution for the radioactive
36 /// materials/mixtures by using one of the 2 methods:
37 ///
38 /// ~~~{.cpp}
39 /// void TGeoMaterial::FillMaterialEvolution(TObjArray *population,
40 /// Double_t precision=0.001)
41 /// ~~~
42 ///
43 /// To use this method, one has to provide an empty TObjArray object that will
44 /// be filled with all elements coming from the decay chain of the initial
45 /// radionuclides contained by the material/mixture. The precision represent the
46 /// cumulative branching ratio for which decay products are still considered.
47 /// The POPULATION list may contain stable elements as well as radionuclides,
48 /// depending on the initial elements. To test if an element is a radionuclide:
49 ///
50 /// ~~~{.cpp}
51 /// Bool_t TGeoElement::IsRadioNuclide() const
52 /// ~~~
53 ///
54 /// All radionuclides in the output population list have attached objects that
55 /// represent the time evolution of their fraction of nuclei with respect to the
56 /// top radionuclide in the decay chain. These objects (Bateman solutions) can be
57 /// retrieved and drawn:
58 ///
59 /// ~~~{.cpp}
60 /// TGeoBatemanSol *TGeoElementRN::Ratio();
61 /// void TGeoBatemanSol::Draw();
62 /// ~~~
63 ///
64 /// Another method allows to create the evolution of a given radioactive
65 /// material/mixture at a given moment in time:
66 ///
67 /// ~~~{.cpp}
68 /// TGeoMaterial::DecayMaterial(Double_t time, Double_t precision=0.001)
69 /// ~~~
70 ///
71 /// The method will create the mixture that result from the decay of a initial
72 /// material/mixture at TIME, while all resulting elements having a fractional
73 /// weight less than PRECISION are excluded.
74 ///
75 /// \macro_image
76 /// \macro_code
77 ///
78 /// \author Mihaela Gheata
79 
80 void DrawPopulation(TObjArray *vect, TCanvas *can, Double_t tmin=0.,
81  Double_t tmax=0., Bool_t logx=kFALSE);
82 
83 void RadioNuclides()
84 {
85  TGeoManager *geom = new TGeoManager("","");
86  TGeoElementTable *table = gGeoManager->GetElementTable();
87  TGeoElementRN *c14 = table->GetElementRN(14,6);
88  TGeoElementRN *el1 = table->GetElementRN(53,20);
89  TGeoElementRN *el2 = table->GetElementRN(78,38);
90  // Radioactive material
91  TGeoMaterial *mat = new TGeoMaterial("C14", c14, 1.3);
92  printf("___________________________________________________________\n");
93  printf("Radioactive material:\n");
94  mat->Print();
95  Double_t time = 1.5e11; // seconds
96  TGeoMaterial *decaymat = mat->DecayMaterial(time);
97  printf("Radioactive material evolution after %g years:\n", time/3.1536e7);
98  decaymat->Print();
99  //Radioactive mixture
100  TGeoMixture *mix = new TGeoMixture("mix", 2, 7.3);
101  mix->AddElement(el1, 0.35);
102  mix->AddElement(el2, 0.65);
103  printf("___________________________________________________________\n");
104  printf("Radioactive mixture:\n");
105  mix->Print();
106  time = 1000.;
107  decaymat = mix->DecayMaterial(time);
108  printf("Radioactive mixture evolution after %g seconds:\n", time);
109  decaymat->Print();
110  TObjArray *vect = new TObjArray();
111  TCanvas *c1 = new TCanvas("c1","C14 decay", 800,600);
112  c1->SetGrid();
113  mat->FillMaterialEvolution(vect);
114  DrawPopulation(vect, c1, 0, 1.4e12);
115  TLatex *tex = new TLatex(8.35e11,0.564871,"C_{N^{14}_{7}}");
116  tex->SetTextSize(0.0388601);
117  tex->SetLineWidth(2);
118  tex->Draw();
119  tex = new TLatex(3.33e11,0.0620678,"C_{C^{14}_{6}}");
120  tex->SetTextSize(0.0388601);
121  tex->SetLineWidth(2);
122  tex->Draw();
123  tex = new TLatex(9.4e11,0.098,"C_{X}=#frac{N_{X}(t)}{N_{0}(t=0)}=\
124  #sum_{j}#alpha_{j}e^{-#lambda_{j}t}");
125  tex->SetTextSize(0.0388601);
126  tex->SetLineWidth(2);
127  tex->Draw();
128  TPaveText *pt = new TPaveText(2.6903e+11,0.0042727,1.11791e+12,0.0325138,"br");
129  pt->SetFillColor(5);
130  pt->SetTextAlign(12);
131  pt->SetTextColor(4);
132  pt->AddText("Time evolution of a population of radionuclides.");
133  pt->AddText("The concentration of a nuclide X represent the ");
134  pt->AddText("ratio between the number of X nuclei and the ");
135  pt->AddText("number of nuclei of the top element of the decay");
136  pt->AddText("from which X derives from at T=0. ");
137  pt->Draw();
138  c1->Modified();
139  vect->Clear();
140  TCanvas *c2 = new TCanvas("c2","Mixture decay", 1000,800);
141  c2->SetGrid();
142  mix->FillMaterialEvolution(vect);
143  DrawPopulation(vect, c2, 0.01, 1000., kTRUE);
144  tex = new TLatex(0.019,0.861,"C_{Ca^{53}_{20}}");
145  tex->SetTextSize(0.0388601);
146  tex->SetTextColor(1);
147  tex->Draw();
148  tex = new TLatex(0.0311,0.078064,"C_{Sc^{52}_{21}}");
149  tex->SetTextSize(0.0388601);
150  tex->SetTextColor(2);
151  tex->Draw();
152  tex = new TLatex(0.1337,0.010208,"C_{Ti^{52}_{22}}");
153  tex->SetTextSize(0.0388601);
154  tex->SetTextColor(3);
155  tex->Draw();
156  tex = new TLatex(1.54158,0.00229644,"C_{V^{52}_{23}}");
157  tex->SetTextSize(0.0388601);
158  tex->SetTextColor(4);
159  tex->Draw();
160  tex = new TLatex(25.0522,0.00135315,"C_{Cr^{52}_{24}}");
161  tex->SetTextSize(0.0388601);
162  tex->SetTextColor(5);
163  tex->Draw();
164  tex = new TLatex(0.1056,0.5429,"C_{Sc^{53}_{21}}");
165  tex->SetTextSize(0.0388601);
166  tex->SetTextColor(6);
167  tex->Draw();
168  tex = new TLatex(0.411,0.1044,"C_{Ti^{53}_{22}}");
169  tex->SetTextSize(0.0388601);
170  tex->SetTextColor(7);
171  tex->Draw();
172  tex = new TLatex(2.93358,0.0139452,"C_{V^{53}_{23}}");
173  tex->SetTextSize(0.0388601);
174  tex->SetTextColor(8);
175  tex->Draw();
176  tex = new TLatex(10.6235,0.00440327,"C_{Cr^{53}_{24}}");
177  tex->SetTextSize(0.0388601);
178  tex->SetTextColor(9);
179  tex->Draw();
180  tex = new TLatex(15.6288,0.782976,"C_{Sr^{78}_{38}}");
181  tex->SetTextSize(0.0388601);
182  tex->SetTextColor(1);
183  tex->Draw();
184  tex = new TLatex(20.2162,0.141779,"C_{Rb^{78}_{37}}");
185  tex->SetTextSize(0.0388601);
186  tex->SetTextColor(2);
187  tex->Draw();
188  tex = new TLatex(32.4055,0.0302101,"C_{Kr^{78}_{36}}");
189  tex->SetTextSize(0.0388601);
190  tex->SetTextColor(3);
191  tex->Draw();
192  tex = new TLatex(117.,1.52,"C_{X}=#frac{N_{X}(t)}{N_{0}(t=0)}=#sum_{j}\
193  #alpha_{j}e^{-#lambda_{j}t}");
194  tex->SetTextSize(0.03);
195  tex->SetLineWidth(2);
196  tex->Draw();
197  TArrow *arrow = new TArrow(0.0235313,0.74106,0.0385371,0.115648,0.02,">");
198  arrow->SetFillColor(1);
199  arrow->SetFillStyle(1001);
200  arrow->SetLineWidth(2);
201  arrow->SetAngle(30);
202  arrow->Draw();
203  arrow = new TArrow(0.0543138,0.0586338,0.136594,0.0146596,0.02,">");
204  arrow->SetFillColor(1);
205  arrow->SetFillStyle(1001);
206  arrow->SetLineWidth(2);
207  arrow->SetAngle(30);
208  arrow->Draw();
209  arrow = new TArrow(0.31528,0.00722919,1.29852,0.00306079,0.02,">");
210  arrow->SetFillColor(1);
211  arrow->SetFillStyle(1001);
212  arrow->SetLineWidth(2);
213  arrow->SetAngle(30);
214  arrow->Draw();
215  arrow = new TArrow(4.13457,0.00201942,22.5047,0.00155182,0.02,">");
216  arrow->SetFillColor(1);
217  arrow->SetFillStyle(1001);
218  arrow->SetLineWidth(2);
219  arrow->SetAngle(30);
220  arrow->Draw();
221  arrow = new TArrow(0.0543138,0.761893,0.0928479,0.67253,0.02,">");
222  arrow->SetFillColor(1);
223  arrow->SetFillStyle(1001);
224  arrow->SetLineWidth(2);
225  arrow->SetAngle(30);
226  arrow->Draw();
227  arrow = new TArrow(0.238566,0.375717,0.416662,0.154727,0.02,">");
228  arrow->SetFillColor(1);
229  arrow->SetFillStyle(1001);
230  arrow->SetLineWidth(2);
231  arrow->SetAngle(30);
232  arrow->Draw();
233  arrow = new TArrow(0.653714,0.074215,2.41863,0.0213142,0.02,">");
234  arrow->SetFillColor(1);
235  arrow->SetFillStyle(1001);
236  arrow->SetLineWidth(2);
237  arrow->SetAngle(30);
238  arrow->Draw();
239  arrow = new TArrow(5.58256,0.00953882,10.6235,0.00629343,0.02,">");
240  arrow->SetFillColor(1);
241  arrow->SetFillStyle(1001);
242  arrow->SetLineWidth(2);
243  arrow->SetAngle(30);
244  arrow->Draw();
245  arrow = new TArrow(22.0271,0.601935,22.9926,0.218812,0.02,">");
246  arrow->SetFillColor(1);
247  arrow->SetFillStyle(1001);
248  arrow->SetLineWidth(2);
249  arrow->SetAngle(30);
250  arrow->Draw();
251  arrow = new TArrow(27.2962,0.102084,36.8557,0.045686,0.02,">");
252  arrow->SetFillColor(1);
253  arrow->SetFillStyle(1001);
254  arrow->SetLineWidth(2);
255  arrow->SetAngle(30);
256  arrow->Draw();
257 }
258 
259 void DrawPopulation(TObjArray *vect, TCanvas *can, Double_t tmin,
260  Double_t tmax, Bool_t logx)
261 {
262  Int_t n = vect->GetEntriesFast();
263  TGeoElementRN *elem;
264  TGeoBatemanSol *sol;
265  can->SetLogy();
266 
267  if (logx) can->SetLogx();
268 
269 
270  for (Int_t i=0; i<n; i++) {
271  TGeoElement *el = (TGeoElement*)vect->At(i);
272  if (!el->IsRadioNuclide()) continue;
273  TGeoElementRN *elem = (TGeoElementRN *)el;
274  TGeoBatemanSol *sol = elem->Ratio();
275  if (sol) {
276  sol->SetLineColor(1+(i%9));
277  sol->SetLineWidth(2);
278  if (tmax>0.) sol->SetRange(tmin,tmax);
279  if (i==0) {
280  sol->Draw();
281  TF1 *func = (TF1*)can->FindObject(
282  Form("conc%s",sol->GetElement()->GetName()));
283  if (func) {
284  if (!strcmp(can->GetName(),"c1")) func->SetTitle(
285  "Concentration of C14 derived elements;time[s];Ni/N0(C14)");
286  else func->SetTitle(
287  "Concentration of elements derived from mixture Ca53+Sr78;\
288  time[s];Ni/N0(Ca53)");
289  }
290  }
291  else sol->Draw("SAME");
292  }
293  }
294 }