Activation from TGeoBatemanSol?

I’m fiddling with the Bateman solver included in the Geometry package (TGeoBatemanSol) for modelling the sequential and/or parallel decay chains given an initial isotope concentration, i.e. something very similar to the ROOT example geom/RadioNuclides.C.

Whereas TGeoBatemanSol::Draw() draws the temporal development of the concentration coefficients cn(t) (actually cn(t)/ctop(0)), I would like to draw the temporal development of the activation as caused by each radioisotope, A_n(t) with decay rate lambda_n. Since the activation component and concentration are intimately related as A_n(t) = lambda_n * cn(t), my request is actually just an individual weighting of of the solutions already found from the solver. It seems to me there is no way of scaling the cn’s individually (only by an overall factor). Is this true? If so, is there any clever workaround to get the activations?

Below is a (non-complete) snippet of my code for reference.

  Double_t lambda, cn;
  for (Int_t i=0; i<n; i++) { // Element loop
    TGeoElement *el = (TGeoElement*)vect->At(i);
    if (!el->IsRadioNuclide()) continue;
    TGeoElementRN* elem = (TGeoElementRN *)el;
    TGeoBatemanSol* sol = elem->Ratio();
    if (sol) {
      if (tmax>0.) sol->SetRange(tmin,tmax);
      const Int_t npx = sol->GetNcoeff();
      for (Int_t j = 0; j < npx; j++){ // Branch loop
      	if (!j) sol->Normalize(Norm/(volume*dens));
        // Would like to normalize by lambda, i.e. show activity:
        // cn = lambda_n*cn.

Maybe an activation plot could be implemented as a draw option? Unfortunately, I don’t think my programming skills are at a sufficient level… :slight_smile:



This looks feasible. Activity was just the last step missing in the implementation. I will try to have a look once back from vacation.


Woaw, sounds terrific! I think the implementation should be somewhat straightforward and and it would be a very powerful tool compared to the implementation effort.

It just came to me that it would be wonderful, if a method / draw option would allow one to plot the combined activation (A = A_1 + A_2 +…+ A_N) as a function of time along with each radioisotope component’s contribution. Again, such an implmentation should be straightforward, now that you’re getting your hands dirty anyway.:wink:

Best regards,