TGraphSmooth and TSpline + Eval failure in loop

ROOT
1

Dear all,

I am experiencing a strange problem.

I made a routine and in this I make a graph of time vs position, I smooth it and then Fit with a Spline. Finally I use Eval to determine at which time a certain position is reached.
The routine works fine.

Then I inserted (copied and pasted) the routine in a more general program. The routine is in a loop with
for(int i=1,i<10000,i++). No error is given, everything works fine, but if I look at the results of “Eval”, 5 times out of 10000 are not correct.

The point is that, if I ask the program to print in a file the “t” and “position” variable for the particular “i” for which I realized the failure happened, and read the file, plot it and fit with the first routine, I obtain the result I should obtain (different from the one given by the second program).

Has anyone already experienced this problem, pls?

Thanks,
regards,
Francesca

2

Hi,

This is strange, maybe is some memory corruption happening in the full program. To understand it, one would need the program to reproduce this result.

Lorenzo

3

Hello.

 In attachment are the files.

I “tarred” them so that the structure is maintained. The file Plot.C is instead a stand alone routine.
Makefile compiles what is in Roottalk.tar.gz.
To launch the program after compilation:

./step13 -1 -300000 90 0.0005 0.010 0.00075

The program calculates the space charge effects felt by an electron (-1 in the line above) generated between two electrodes. The distance between the 2 electrodes is 10 cm (it goes from y=-5 cm to y = +5 cm) and the electric field is -300000 V/m (second argument in the line to launch the routine). A proton beam of energy 90 MeV (third argument in the line to launch the routine) is travelling along the z direction passing through the middle of an ionization chamber. The proton beam has the following widths:
sigma_x =0.0005 m, sigma_y = 0.01 m, sigma_z= 0.00075 m (given in the line to launch the program).
The electron was generated between the two electrodes, at a position:

  • x sampled with a Gaussian distribution around 0 with sigma = sigma_x
  • y sampled with a Gaussian distribution around 0 with sigma = sigma_y
  • z sampled with Uniformly between -2.5 and 2.5 mm

Once the electron has been tracked all along the y direction (up to y = +5 cm), a TTree is Filled (but not saved up to the end of the program) and another electron is generated. The number of electrons generated is fixed by the variable “test” in Step13.C.
For this particular case (Ey = -300000 V/m, Proton_Energy = 90 MeV, sigma_x =0.0005 m, sigma_y = 0.01 m, sigma_z= 0.00075 m) the first time the program fails on my laptop is at test = 651.

The program uses an adaptive Runge-Kutta routine (partly copied from MATLAB) to determine for some times the position x,y,z of the electron (and since I am not a great programmer is rather slow…).
An array of times (“time”) is given as input, but the program does not really use it.
It somehow overwrites it and copy the new times chosen by itself in the array “vremia”.

Since the steps in y are defined by the “dt” chosen by program (and therefore the electron in the last step can cross the electrode passing for example from y=4.9 cm to y=5.1 cm) the program stops tracking the electron when the condition y_electron>5 cm is fulfilled. At this point a Spline is used to fit the data and find at which moment y=5 cm was reached.

Finally, in order not to save too many data, only the positions at the times = “timen” are saved.
They are obtain by fitting the spline at the times “timen”.
When the program fails, it writes on a Data.dat file:

“BEFORE” , vremia, y <-- for all the track, but the last step (since for the last step of “vremia” is overwritten with the value of the fit of the spline at y=5cm)
"BEFORE " ,variable_for_checking_t, variable_for_checking_t <— for the last step (variable_for_checking_t contains the vremia[last_step] value before it was overwritten)
"PROCESSED " vremia[length], distance <----- last step. vremia[length] is obtained by the fit of the spline
"AFTER ", timen, y <— from the spline. And they are all correct, in the sense that for that particular t, you’d obtain that y. So I really do not understand why the step “LAST” fails instead.
Have I made a mistake somewhere, pls?

If you open the Data.dat file, copy the data in the lines “BEFORE” to another file (ex. test.dat) and use the Plot.C (after cutting away the word “BEFORE” from every line) to fit them with a spline and ask it to give you the time when the electron reached the position y = +5cm, a value different from the one obtained in Data.dat, "entry = “PROCESSED” is obtained (and this value makes sense).

I am sorry the program is such a mess.

Thanks,
cheers,
Francesca
Spline.C (1.96 KB)
Roottalk.tar.gz (12.8 KB)

4

Dear all,

it as a mistake in the code.

If I replace

grout = gs->Approx(g,“linear”,length);

with

grout = gs->Approx(g,“linear”);

the problem is solved.

Cheers,
Francesca

5

Hello.

I thought it was solved, but now the problem occurs at test = 270.
Could anyone give a look, pls?

Moreover I do not understand why, if I try to delete gs, gs1, gs2,… to free memory (they are pointers defined with “new”) the program crushes.

Thanks,
regards,
cheers,
Francesca

6

Try the attached Makefile.

You have many places (in several source code files) where you use the operator new[] but you do not use the operator delete[] afterwards. You should fix all these places.

In the “Step13TTree.c”, I fixed several “most obvious” places only (again, you need to go through all another reported compiler warnings and fix them yourself): // ... TFile *f_res = 0; // ... const int n = 87;//105; // ... // ... for (unsigned int i=0;i<history_tm.size();i++) // ... for (unsigned int i=0;i<history_tm.size();i++) fprintf(fff, " AFTER %d %e %e \n", test, history_tm.at(i), history_ps_y.at(i));} // ... // ... if (f_res) f_res->Write(); // ...
Also, you create a “memory leak” when you have something like this: TGraph *grout = new TGraph(length); TGraphSmooth *gs = new TGraphSmooth("g"); grout = gs->Approx(g,"linear",length); You should have: TGraphSmooth *gs = new TGraphSmooth("g"); TGraph *grout = gs->Approx(g, "linear", length); Finally, you should not “delete grout;”. You should only “delete gs;” (and it will automatically delete “grout”, which it “owns”).

When you fix everything, we’ll see what happens (I was able to execute “./step13 -1 -300000 90 0.0005 0.010 0.00075” up to “particle nb 660” with one “EROOR FOR particle 651”).
Makefile.txt (864 Bytes)

7

Hello.

Thanks a lot. I run all the week end the code up to 10000 particles.
I get 32 failure of the Spline, the first of which at test = 269.
Could you help me to fix my code so that there are no failures, please?

Thanks,
regards,
Francesca

8

Replace the line: else (history_tm.at(history_tm.size()-1) = s->Eval(distance)); with: else if(s->Eval(distance)>history_tm.at(history_tm.size()-2)) history_tm.at(history_tm.size()-1) = s->Eval(distance); else std::cout << " ... possible problems @ " << test << " ..." << std::endl;

9

Hello.

Thanks a lot.
I did it, but still sometimes still fails thought the times are in increasing order.
I went around just saying that, if the time value estimated by using “Eval” is such that when “reversing” the graph (i.e. you use the time you got to evaluate the distance) you get a “y” value fulfilling the condition:
absolute_value(y-distance) >0.001

then test =test-1 (i.e. instead of running 10000 events, the program will run 10001 events and will discard the bad one).
But still I do not understand why from time to time Eval fails in my code…

Any help is appreciated.

Cheers,
Francesca

10

Try to dump your graphs and splines for these “problem cases” and analyse them “manually”.
My guess is that in these cases you do not get an “injective function” (especially around the maximal “distance” value) and then the “reversing” misbehaves.
For example, I tip that for all these “problem cases” you have:
s->Eval(distance) < history_tm.at(history_tm.size()-2) < history_tm.at(history_tm.size()-1)
or at least (I tried to “hide” this case by the patch from my previous post here):
history_tm.at(history_tm.size()-2) < s->Eval(distance) < history_tm.at(history_tm.size()-1)
while in general one would like to have:
history_tm.at(history_tm.size()-2) < history_tm.at(history_tm.size()-1) <= s->Eval(distance)

11

Hello.

I thank for your answer.

I’ve been dumping the problematic cases these days, but often it looks as “Eval” in my program fails with no reason.

For example, in the attachment is the case test=270.
Here the lines labelled “BEFORE” report the times (3rd column) and y positions (4th column) of the electron.
The line labelled “PROCESSED” is the result of Eval, therefore the time is evaluated at y = 0.05 cm.
The result (9.008819e-11 s) is wrong.

I attach another file: test.dat. This is nothing else that all lines labelled with “BEFORE” in Data.dat.
If I feed the routine “Spline.C” (loaded in my second post) with test.dat, this will give a result
t of the order of 10^-9 s (as it should be) and not the 9.008819e-11 s as reported at the line labelled with “PROCESS” in the Data.dat file.

I do not understand what’s going on, because the stand alone routine fits the same data as the program.
Even the code to do it is the same (in the program the code used in Spline.C is put in a loop), therefore if one fails, the other should fail too.

Are there other sources of memory leak in my program apart those ones already pointed out, pls?

Thanks,
cheers,
Francesca
test.dat (10.1 KB)
Data.dat (16.8 KB)

12

So, I have finally understood why your “splines” sometimes misbehave around the maximal “distance” value.
Modify ALL your “Approx” calls adding four parameters in the end: TGraph *grout = gs->Approx(g, "linear", length, 0, 0, 0, 2); // rule 2

13

Thanks!!!

Fra