ROOT_528-00b_version: tutorials/geom/RadioNuclides.C Source File

00001 void RadioNuclides()
00002 {
00003 // Macro that demonstrates usage of radioactive elements/materials/mixtures
00004 // with TGeo package.
00005 //
00006 // A radionuclide (TGeoElementRN) derives from the class TGeoElement and
00007 // provides additional information related to its radioactive properties and
00008 // decay modes.
00009 //
00010 // The radionuclides table is loaded on demand by any call:
00011 //    TGeoElementRN *TGeoElementTable::GetElementRN(Int_t atomic_number, 
00012 //                                                  Int_t atomic_charge, 
00013 //                                                  Int_t isomeric_number)
00014 // The isomeric number is optional and the default value is 0.
00015 //
00016 // To create a radioactive material based on a radionuclide, one should use the
00017 // constructor:
00018 //    TGeoMaterial(const char *name, TGeoElement *elem, Double_t density)
00019 // To create a radioactive mixture, one can use radionuclides as well as stable
00020 // elements:
00021 //    TGeoMixture(const char *name, Int_t nelements, Double_t density);
00022 //    TGeoMixture::AddElement(TGeoElement *elem, Double_t weight_fraction);
00023 // Once defined, one can retrieve the time evolution for the radioactive 
00024 // materials/mixtures by using one of the 2 methods:
00025 //
00026 //    void TGeoMaterial::FillMaterialEvolution(TObjArray *population,
00027 //                                             Double_t   precision=0.001)
00028 // To use this method, one has to provide an empty TObjArray object that will
00029 // be filled with all elements coming from the decay chain of the initial 
00030 // radionuclides contained by the material/mixture. The precision represent the
00031 // cumulative branching ratio for which decay products are still considered.
00032 // The POPULATION list may contain stable elements as well as radionuclides,
00033 // depending on the initial elements. To test if an element is a radionuclide:
00034 //    Bool_t TGeoElement::IsRadioNuclide() const
00035 // All radionuclides in the output population list have attached objects that
00036 // represent the time evolution of their fraction of nuclei with respect to the
00037 // top radionuclide in the decay chain. These objects (Bateman solutions) can be
00038 // retrieved and drawn:
00039 //    TGeoBatemanSol *TGeoElementRN::Ratio();
00040 //    void TGeoBatemanSol::Draw();
00041 //
00042 // Another method allows to create the evolution of a given radioactive 
00043 // material/mixture at a given moment in time:
00044 //    TGeoMaterial::DecayMaterial(Double_t time, Double_t precision=0.001)
00045 // The method will create the mixture that result from the decay of a initial
00046 // material/mixture at TIME, while all resulting elements having a fractional
00047 // weight less than PRECISION are excluded.
00048 //Author: Mihaela Gheata
00049 
00050    TGeoManager *geom = new TGeoManager("","");
00051    TGeoElementTable *table = gGeoManager->GetElementTable();
00052    TGeoElementRN *c14 = table->GetElementRN(14,6);
00053    TGeoElementRN *el1 = table->GetElementRN(53,20);
00054    TGeoElementRN *el2 = table->GetElementRN(78,38);
00055    // Radioactive material
00056    TGeoMaterial *mat = new TGeoMaterial("C14", c14, 1.3);
00057    printf("___________________________________________________________\n");
00058    printf("Radioactive material:\n");
00059    mat->Print();
00060    Double_t time = 1.5e11; // seconds
00061    TGeoMaterial *decaymat = mat->DecayMaterial(time);
00062    printf("Radioactive material evolution after %g years:\n", time/3.1536e7);
00063    decaymat->Print();
00064    //Radioactive mixture
00065    TGeoMixture *mix = new TGeoMixture("mix", 2, 7.3);
00066    mix->AddElement(el1, 0.35);
00067    mix->AddElement(el2, 0.65);
00068    printf("___________________________________________________________\n");
00069    printf("Radioactive mixture:\n");
00070    mix->Print();
00071    time = 1000.;
00072    decaymat = mix->DecayMaterial(time);
00073    printf("Radioactive mixture evolution after %g seconds:\n", time);
00074    decaymat->Print();
00075    TObjArray *vect = new TObjArray();
00076    TCanvas *c1 = new TCanvas("c1","C14 decay", 800,600);
00077    c1->SetGrid();
00078    mat->FillMaterialEvolution(vect);
00079    DrawPopulation(vect, c1, 0, 1.4e12);
00080    TLatex *tex = new TLatex(8.35e11,0.564871,"C_{N^{14}_{7}}");
00081    tex->SetTextSize(0.0388601);
00082    tex->SetLineWidth(2);
00083    tex->Draw();
00084    tex = new TLatex(3.33e11,0.0620678,"C_{C^{14}_{6}}");
00085    tex->SetTextSize(0.0388601);
00086    tex->SetLineWidth(2);
00087    tex->Draw();
00088    tex = new TLatex(9.4e11,0.098,"C_{X}=#frac{N_{X}(t)}{N_{0}(t=0)}=\
00089    #sum_{j}#alpha_{j}e^{-#lambda_{j}t}");
00090    tex->SetTextSize(0.0388601);
00091    tex->SetLineWidth(2);
00092    tex->Draw();
00093    TPaveText *pt = new TPaveText(2.6903e+11,0.0042727,1.11791e+12,0.0325138,"br");
00094    pt->SetFillColor(5);
00095    pt->SetTextAlign(12);
00096    pt->SetTextColor(4);
00097    text = pt->AddText("Time evolution of a population of radionuclides.");
00098    text = pt->AddText("The concentration of a nuclide X represent the  ");
00099    text = pt->AddText("ratio between the number of X nuclei and the    ");
00100    text = pt->AddText("number of nuclei of the top element of the decay");
00101    text = pt->AddText("from which X derives from at T=0.               ");
00102    pt->Draw();
00103    c1->Modified();
00104    vect->Clear();
00105    TCanvas *c2 = new TCanvas("c2","Mixture decay", 1000,800);
00106    c2->SetGrid();
00107    mix->FillMaterialEvolution(vect);
00108    DrawPopulation(vect, c2, 0.01, 1000., kTRUE);
00109    tex = new TLatex(0.019,0.861,"C_{Ca^{53}_{20}}");
00110    tex->SetTextSize(0.0388601);
00111    tex->SetTextColor(1);
00112    tex->Draw();
00113    tex = new TLatex(0.0311,0.078064,"C_{Sc^{52}_{21}}");
00114    tex->SetTextSize(0.0388601);
00115    tex->SetTextColor(2);
00116    tex->Draw();
00117    tex = new TLatex(0.1337,0.010208,"C_{Ti^{52}_{22}}");
00118    tex->SetTextSize(0.0388601);
00119    tex->SetTextColor(3);
00120    tex->Draw();
00121    tex = new TLatex(1.54158,0.00229644,"C_{V^{52}_{23}}");
00122    tex->SetTextSize(0.0388601);
00123    tex->SetTextColor(4);
00124    tex->Draw();
00125    tex = new TLatex(25.0522,0.00135315,"C_{Cr^{52}_{24}}");
00126    tex->SetTextSize(0.0388601);
00127    tex->SetTextColor(5);
00128    tex->Draw();
00129    tex = new TLatex(0.1056,0.5429,"C_{Sc^{53}_{21}}");
00130    tex->SetTextSize(0.0388601);
00131    tex->SetTextColor(6);
00132    tex->Draw();
00133    tex = new TLatex(0.411,0.1044,"C_{Ti^{53}_{22}}");
00134    tex->SetTextSize(0.0388601);
00135    tex->SetTextColor(7);
00136    tex->Draw();
00137    tex = new TLatex(2.93358,0.0139452,"C_{V^{53}_{23}}");
00138    tex->SetTextSize(0.0388601);
00139    tex->SetTextColor(8);
00140    tex->Draw();
00141    tex = new TLatex(10.6235,0.00440327,"C_{Cr^{53}_{24}}");
00142    tex->SetTextSize(0.0388601);
00143    tex->SetTextColor(9);
00144    tex->Draw();
00145    tex = new TLatex(15.6288,0.782976,"C_{Sr^{78}_{38}}");
00146    tex->SetTextSize(0.0388601);
00147    tex->SetTextColor(1);
00148    tex->Draw();
00149    tex = new TLatex(20.2162,0.141779,"C_{Rb^{78}_{37}}");
00150    tex->SetTextSize(0.0388601);
00151    tex->SetTextColor(2);
00152    tex->Draw();
00153    tex = new TLatex(32.4055,0.0302101,"C_{Kr^{78}_{36}}");
00154    tex->SetTextSize(0.0388601);
00155    tex->SetTextColor(3);
00156    tex->Draw();
00157    tex = new TLatex(117.,1.52,"C_{X}=#frac{N_{X}(t)}{N_{0}(t=0)}=#sum_{j}\
00158    #alpha_{j}e^{-#lambda_{j}t}");
00159    tex->SetTextSize(0.03);
00160    tex->SetLineWidth(2);
00161    tex->Draw();
00162    TArrow *arrow = new TArrow(0.0235313,0.74106,0.0385371,0.115648,0.02,">");
00163    arrow->SetFillColor(1);
00164    arrow->SetFillStyle(1001);
00165    arrow->SetLineWidth(2);
00166    arrow->SetAngle(30);
00167    arrow->Draw();
00168    arrow = new TArrow(0.0543138,0.0586338,0.136594,0.0146596,0.02,">");
00169    arrow->SetFillColor(1);
00170    arrow->SetFillStyle(1001);
00171    arrow->SetLineWidth(2);
00172    arrow->SetAngle(30);
00173    arrow->Draw();
00174    arrow = new TArrow(0.31528,0.00722919,1.29852,0.00306079,0.02,">");
00175    arrow->SetFillColor(1);
00176    arrow->SetFillStyle(1001);
00177    arrow->SetLineWidth(2);
00178    arrow->SetAngle(30);
00179    arrow->Draw();
00180    arrow = new TArrow(4.13457,0.00201942,22.5047,0.00155182,0.02,">");
00181    arrow->SetFillColor(1);
00182    arrow->SetFillStyle(1001);
00183    arrow->SetLineWidth(2);
00184    arrow->SetAngle(30);
00185    arrow->Draw();
00186    arrow = new TArrow(0.0543138,0.761893,0.0928479,0.67253,0.02,">");
00187    arrow->SetFillColor(1);
00188    arrow->SetFillStyle(1001);
00189    arrow->SetLineWidth(2);
00190    arrow->SetAngle(30);
00191    arrow->Draw();
00192    arrow = new TArrow(0.238566,0.375717,0.416662,0.154727,0.02,">");
00193    arrow->SetFillColor(1);
00194    arrow->SetFillStyle(1001);
00195    arrow->SetLineWidth(2);
00196    arrow->SetAngle(30);
00197    arrow->Draw();
00198    arrow = new TArrow(0.653714,0.074215,2.41863,0.0213142,0.02,">");
00199    arrow->SetFillColor(1);
00200    arrow->SetFillStyle(1001);
00201    arrow->SetLineWidth(2);
00202    arrow->SetAngle(30);
00203    arrow->Draw();
00204    arrow = new TArrow(5.58256,0.00953882,10.6235,0.00629343,0.02,">");
00205    arrow->SetFillColor(1);
00206    arrow->SetFillStyle(1001);
00207    arrow->SetLineWidth(2);
00208    arrow->SetAngle(30);
00209    arrow->Draw();
00210    arrow = new TArrow(22.0271,0.601935,22.9926,0.218812,0.02,">");
00211    arrow->SetFillColor(1);
00212    arrow->SetFillStyle(1001);
00213    arrow->SetLineWidth(2);
00214    arrow->SetAngle(30);
00215    arrow->Draw();
00216    arrow = new TArrow(27.2962,0.102084,36.8557,0.045686,0.02,">");
00217    arrow->SetFillColor(1);
00218    arrow->SetFillStyle(1001);
00219    arrow->SetLineWidth(2);
00220    arrow->SetAngle(30);
00221    arrow->Draw();
00222 }
00223 
00224 void DrawPopulation(TObjArray *vect, TCanvas *can, Double_t tmin=0., 
00225    Double_t tmax=0., Bool_t logx=kFALSE)
00226 {
00227    Int_t n = vect->GetEntriesFast();
00228    TGeoElementRN *elem;
00229    TGeoBatemanSol *sol;
00230    can->SetLogy();
00231    
00232    if (logx) can->SetLogx();
00233    
00234 
00235    for (Int_t i=0; i<n; i++) {
00236       TGeoElement *el = (TGeoElement*)vect->At(i);
00237       if (!el->IsRadioNuclide()) continue;
00238       TGeoElementRN *elem = (TGeoElementRN *)el;
00239       TGeoBatemanSol *sol = elem->Ratio();
00240       if (sol) {
00241          sol->SetLineColor(1+(i%9));
00242          sol->SetLineWidth(2);
00243          if (tmax>0.) sol->SetRange(tmin,tmax);
00244          if (i==0) {
00245             sol->Draw();
00246             TF1 *func = (TF1*)can->FindObject(
00247                Form("conc%s",sol->GetElement()->GetName()));
00248             if (func) {
00249                if (!strcmp(can->GetName(),"c1")) func->SetTitle(
00250                   "Concentration of C14 derived elements;time[s];Ni/N0(C14)");
00251                else func->SetTitle(
00252                   "Concentration of elements derived from mixture Ca53+Sr78;\
00253                   time[s];Ni/N0(Ca53)");
00254             }   
00255          }   
00256          else      sol->Draw("SAME");
00257       }
00258    }      
00259 }