How to fit a process with memory

Newbie
1

Hello,
I am trying to build a tool to fit measurement data to a mathematical model. The goal is to characterize defects on irradiated silicon sensors.

The measurement consists in recording the electrical current through the detector while changing its temperature. So, from my measurement I obtain three arrays of data: (time (t), Temperature (T), Current (I))

The current is mathematically described as:

I[k] = n_0[k] x e_0[k] + ... + n_i[k] x e_i[k] + ... + n_N[k] x e_N[k]

where N is the number of defect levels in my irradiated sensor, and for each level i, n_i[k] and e_i[k] are described as follows:

e_i[k] = p0_i * T[k]^2 * F(T[k]) * exp (-p1_i / (K * T[k]) )

n_i[k] = p2_i * exp ( -e_i[k] * C) * n_i[k-1]

where F(T) is a tabulated function, and K and C are known constants. p0_i , p1_i and p2_i are the parameters that I am trying to fit.

I do not know if such a complicated function can be implemented and fitted in ROOT somehow. I am especially concerned by the fact that the measured current depends not only on the instant temperature but on the history (n_i[k] depends on n_i[k-1]).

I am quite new in ROOT so I am still trying to understand which set of functions I should use (or, even, if this is feasible at all). For the moment I am quite lost, so any hint will be greatly appreciated.

Thanks!

Isidre

2

This doesn’t seem too daunting for root. In the end the whole function I[k] is just a 1-dimensional function of k, with 3 parameters. All the complicated calculations would happen internally to the function, and you only need to tell ROOT that it’s a 1D function with 3 parameters.

I would suggest reading up about fitting and minimization in ROOT, and going through some of the tutorial examples:

https://root.cern.ch/numerical-algorithms
https://root.cern.ch/fitting
https://root.cern.ch/root/html/tutorials/fit/index.html

That’ll be way more information than you need, so focus on the information about 1D functions.

1 Like
3

Hello Isidre,

another way to go forward is RooFit. It’s designed to implement complex statistical models.

The easy part is writing down the exponentials, coefficients and summations, and fitting these to data. I can point you to plenty of tutorials that show how to do this.
What will be difficult (but probably solvable) is the tabulated function. Is it a C++ function that can be called?
What even more difficult is the dependence on the previous entry in the dataset. Is it really just a double / float that needs to be remembered across function evaluations? One thing that comes to mind is to write a function with an evil hack, i.e. it remembers the previous entry.

Let me know what’s needed.

PS:
I added three “`” before and after the formulae in your post, so they get typeset in a font that’s more easy to read.

4

Hello,
Thank you so much to both for your answers.

First of all, I just realized an error in my original post, the parameter p2_i is actually the initial value of n_i[k], and should not appear as multiplication factor in the original formula. It should be instead:

n_i[k] = exp ( -e_i[k] * C) * n_i[k-1] 
with n_i[0] = p2_i

And n_i[k] can be rewritten to depend on the integral of e_n[k] since the origin of the measurement. That is:

n_i[k] = exp ( -e_i[k] * C) * n_i[k-1] =  exp ( -e_i[k] * C - e_i[k-1]*C) * n_i[k-2] =  p2_i * exp ( -(e_i[k] + ...e_i[1])*C)

So, I think the problem boils down to having access to the whole time and temperature vector from inside the fitting function, so that I can integrate e_i[k]. What I did as first try is to create my own class and declare the time and temperature vectors as private data members. The same goes for the tabulated function, which I can pre-compute using the temperature information. The fitting function is implemented in the class, so it can access both the temperature and the tabulated function without need for passing them as arguments to the function.

In summary, my code looks something like this:

class MyClass{
private:
   vector<double> Temp, Ftab, time;

public:

  MyClass(vector<double> time vector<double> Temp, vector<double> Ftab){
       this->time = time;
       this->Temp = Temp;
       this->Ftab = Ftab;
    }

   double operator()(Double_t *x, Double_t *par){
        [...]
        for(k=1; time.at(k)<=x[0] ;k++){
            // Here I calculate integral from t0 to current time (x[0]), using Ftab and Temp as well.
        }
        [...]
   }
}

int main(){

   [...]
   MyClass fitObject(time, Temp, Ftab)
   TF1 *fitFunc = new TF1("fitFunc", fitObject, t_0, t_end, N_parameters); // This is now my fitting function
  [...]
}

This seems to work now. I tried with some simulated data and only one defect level and the fit converged to the right values.

I will see what happens for data with more levels. If I get into trouble again I will give it a try with RooFit.

Thanks again!

Isidre