Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
80void DrawPopulation(TObjArray *vect, TCanvas *can, Double_t tmin=0.,
81 Double_t tmax=0., Bool_t logx=kFALSE);
82
83void RadioNuclides()
84{
85 TGeoManager *geom = new TGeoManager("","");
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
259void 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}
int Int_t
Definition RtypesCore.h:45
const Bool_t kFALSE
Definition RtypesCore.h:92
bool Bool_t
Definition RtypesCore.h:63
double Double_t
Definition RtypesCore.h:59
const Bool_t kTRUE
Definition RtypesCore.h:91
R__EXTERN TGeoManager * gGeoManager
char * Form(const char *fmt,...)
Draw all kinds of Arrows.
Definition TArrow.h:29
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition TAttFill.h:37
virtual void SetLineWidth(Width_t lwidth)
Set the line width.
Definition TAttLine.h:43
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition TAttLine.h:40
virtual void SetTextAlign(Short_t align=11)
Set the text alignment.
Definition TAttText.h:41
virtual void SetTextColor(Color_t tcolor=1)
Set the text color.
Definition TAttText.h:43
virtual void SetTextSize(Float_t tsize=1)
Set the text size.
Definition TAttText.h:46
The Canvas class.
Definition TCanvas.h:23
1-Dim function class
Definition TF1.h:213
virtual void SetTitle(const char *title="")
Set function title if title has the form "fffffff;xxxx;yyyy", it is assumed that the function title i...
Definition TF1.cxx:3562
virtual void Draw(Option_t *option="")
Draw the solution of Bateman equation versus time.
void SetRange(Double_t tmin=0., Double_t tmax=0.)
TGeoElementRN * GetElement() const
Class representing a radionuclidevoid TGeoManager::SetDefaultRootUnits() { if ( fgDefaultUnits == kRo...
TGeoBatemanSol * Ratio() const
Table of elements.
TGeoElementRN * GetElementRN(Int_t ENDFcode) const
Retrieve a radionuclide by ENDF code.
Base class for chemical elements.
Definition TGeoElement.h:37
virtual Bool_t IsRadioNuclide() const
Definition TGeoElement.h:87
The manager class for any TGeo geometry.
Definition TGeoManager.h:45
TGeoElementTable * GetElementTable()
Returns material table. Creates it if not existing.
Base class describing materials.
virtual void Print(const Option_t *option="") const
print characteristics of this material
virtual TGeoMaterial * DecayMaterial(Double_t time, Double_t precision=0.001)
Create the material representing the decay product of this material at a given time.
virtual void FillMaterialEvolution(TObjArray *population, Double_t precision=0.001)
Fills a user array with all the elements deriving from the possible decay of the top element composin...
Mixtures of elements.
void AddElement(Double_t a, Double_t z, Double_t weight)
add an element to the mixture using fraction by weight Check if the element is already defined
virtual void Print(const Option_t *option="") const
print characteristics of this material
virtual void FillMaterialEvolution(TObjArray *population, Double_t precision=0.001)
Fills a user array with all the elements deriving from the possible decay of the top elements composi...
virtual TGeoMaterial * DecayMaterial(Double_t time, Double_t precision=0.001)
Create the mixture representing the decay product of this material at a given time.
To draw Mathematical Formula.
Definition TLatex.h:18
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47
An array of TObjects.
Definition TObjArray.h:37
Int_t GetEntriesFast() const
Definition TObjArray.h:64
virtual void Clear(Option_t *option="")
Remove all objects from the array.
TObject * At(Int_t idx) const
Definition TObjArray.h:166
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition TObject.cxx:197
void SetLogy(Int_t value=1) override
Set Lin/Log scale for Y.
Definition TPad.cxx:5947
TObject * FindObject(const char *name) const override
Search if object named name is inside this pad or in pads inside this pad.
Definition TPad.cxx:2625
void SetLogx(Int_t value=1) override
Set Lin/Log scale for X.
Definition TPad.cxx:5933
const char * GetName() const override
Returns name of object.
Definition TPad.h:254
A Pave (see TPave) with text, lines or/and boxes inside.
Definition TPaveText.h:21
virtual TText * AddText(Double_t x1, Double_t y1, const char *label)
Add a new Text line to this pavetext at given coordinates.
virtual void Draw(Option_t *option="")
Draw this pavetext with its current attributes.
TPaveText * pt
return c1
Definition legend1.C:41
const Int_t n
Definition legend1.C:16
return c2
Definition legend2.C:14