TMultiDimFit - parameter output makes no sense, nan and inf

ROOT
1

Hey Rooters,

I recently recognized that there exists a wonderful package in ROOT that allows one to make a fit for a more dimensional data set.

I tried it instantaneously, cause this was exactly what I needed for my analysis. Unfortunatly my code didn’t work that well and I don’t know why.

So here is my code together with the command lines that ROOT gave back to me…
Maybe one of you can help me…

code:
--------------------------------------------------------start
#include “Riostream.h”
#include “TROOT.h”
#include “TApplication.h”
#include “TCanvas.h”
#include “TH1.h”
#include “TSystem.h”
#include “TBrowser.h”
#include “TFile.h”
#include “TRandom.h”
#include “TMultiDimFit.h”
#include “TVectorD.h”
#include
#include
#include
#include
#include

void multifit() {

////open new file////
TFile *output = new TFile(“output.root”,“recreate”);

////prepare datasample////
TChain chainTP1(“T”);
chainTP1.Add(“wholedataTP1.root”);

const Double_t pi = 3.14159265; //define the number PI
const Double_t tp1 = 0.914;

Double_t en, sz, xc, yc, ximp, yimp, coszen, el, id;

Double_t dist; //distance from center of field of view
Double_t logsz; //logarithm of size to the base 10
Double_t logen; //logarith of energy*1000 to the base 10
Double_t imp; //impact distance
Double_t elrad; //elevation in radians
Double_t sinel; //sinus of elevation
Double_t logsinel;
Double_t x[3];

chainTP1.SetBranchAddress(“id”, &id);
chainTP1.SetBranchAddress(“en”, &en);
chainTP1.SetBranchAddress(“sz”, &sz);
chainTP1.SetBranchAddress(“xc”, &xc);
chainTP1.SetBranchAddress(“yc”, &yc);
chainTP1.SetBranchAddress(“ximp”, &ximp);
chainTP1.SetBranchAddress(“yimp”, &yimp);
chainTP1.SetBranchAddress(“el”, &el);

////create multidim fitter////
TMultiDimFit *multifit = new TMultiDimFit(3,TMultiDimFit::kMonomials,"");

Int_t mPowers[] = {5,2,1};
multifit->SetMaxPowers(mPowers);
multifit->SetMaxFunctions(500);
multifit->SetMaxStudy(500);
multifit->SetMaxTerms(15);
multifit->SetPowerLimit(3);
//multifit->SetMinAngle(45);
//multifit->SetMaxAngle(5);
multifit->SetMinRelativeError(0.1);

Double_t d,e;

////fitting////
Int_t nEntriesTP1 = (Int_t)chainTP1.GetEntries();
for (Int_t i=0; iAddRow(x,d,E=0);

}

multifit->Print(“s”);
multifit->MakeHistograms();
multifit->FindParameterization("");
multifit->Print(“rc”);

//TBrowser browser = new TBrowser(“browser”,“browsing multifit”);
//multifit->Browse(browser);

multifit->MakeCode();

output->Write();
output->Close();

//delete output;
//delete multifit;

}
----------------------------------------------------------------------end

and root gave me the following answer:
----------------------------------------------------------------------start
root [0] .x multifit.C
Warning: Automatic variable E is allocated multifit.C:90:
Sample statistics:

             D          1          2          3

Max: -311.2601 0 -314.2 -316
Min: -311.294 0 -inf -316
Mean: -311.2855 0 -inf -316
Function Sum Squares: 7.532e+10

Results of Parameterisation:

Total reduction of square residuals nan
Relative precision obtained: nan
Error obtained: nan
Multiple correlation coefficient: nan
Reduced Chi square over sample: nan
Maximum residual value: -1e+11
Minimum residual value: 1e+11
Estimated root mean square: nan
Maximum powers used: 5 2 1
Function codes of candidate functions.
1: considered, 2: too little contribution, 3: accepted.
3333333333 3333300000 0000000000 000000
Loop over candidates stopped because max number of terms reached

Coefficients:

Value Error Powers

0 nan nan 0 0 0
1 nan nan 1 0 0
2 nan nan 2 0 0
3 nan nan 0 1 0
4 nan nan 3 0 0
5 nan nan 1 1 0
6 nan nan 4 0 0
7 nan nan 2 1 0
8 nan nan 0 0 1
9 nan nan 5 0 0
10 nan nan 0 2 0
11 nan nan 3 1 0
12 nan nan 1 0 1
13 nan nan 1 2 0
14 nan nan 4 1 0

root [1]
-------------------------------------------------------------------end

You see, there are a lot of "not a number"s and "infinity"s where there shouldn’t be…
I’m glad about every piece of help here…
Thanks a lot,

Chris

2

I cannot test your code that requires an input file.
However, I see a compilation error at
for (Int_t i=0; iAddRow(x,d,E=0); //this line is wrong
You should
#include "TChain.h"
and remove
#include “TMultiDimFit” (already TMultiDimFit.h)
I suggest to do
root > .L multifit.C+

Rene

3

Rene,

thanks for the quick reply. Indeed, the line with the for loop you mentioned would be a big coding error. But this is just a error that occured because of bad “copy’n’paste”-ing into the forum. It should be

////fitting////
Int_t nEntriesTP1 = (Int_t)chainTP1.GetEntries();
for (Int_t i=0; iAddRow(x,d,E=0);

}

the rest should be the same as mentioned the message befor. I’m going to include the header you mentioned and to remove the . I hope this works.

Thanks a lot,
Chris

4

again an error of “copy’n’paste”…
but this time it will work…

Double_t d,e;

////fitting////
Int_t nEntriesTP1 = (Int_t)chainTP1.GetEntries();
for (Int_t i=0; i
chainTP1.GetEntry(i);

dist = sqrt(xcxc + ycyc);
logsz = log10(sz/tp1);
elrad = (el/180.)pi;
sinel = sin(elrad);
logsinel = log10(sinel);
logen = log10(1000en);

x[0] = dist;
x[1] = logsz;
x[2] = logsinel;

d = logen;

multifit->AddRow(x,d,E=0);

}

multifit->Print(“s”);
multifit->MakeHistograms();
multifit->FindParameterization(“”);
multifit->Print(“rc”);

ok, that’s it…
thanks,
Chris

5

Hi again,

I’m going to attach a data file, so you can try by your self. The file is not a good representative but it’s the smallest one I have. I also attach a modified c-file including my code.

thanks,
Chris
multifit3.C (2.4 KB)
wholedataTP1E4.root (915 KB)

6

Your Tree variables are floats , not doubles!!
Change the line
Double_t en, sz, xc, yc, ximp, yimp, el, id;
to
Float_t en, sz, xc, yc, ximp, yimp, el, id;

Then, you have cases where your variable sz = 0 and you do
logsz = log10(sz/tp1);
This is very clearly shown by the output of TMuktiDimFit.

Rene

7

Rene,

of course - size=0 makes no sense. I already had trouble with that obvious case some time ago. Thanks for the reminder. I tried to follow the example in the reference guide and I thought the doubles were necessary. But I didn’t keep in mind that I already declared my variables as floats.

This should do it,
thanks,
Chris

8

Ok, this did the job. Impressive package!! But there is one thing. When I look into that MDF.C the powers within this file are not the same than the powers ROOT gives me onto the screen. They are always one integer bigger. I show you:

This I get from ROOT onto my screen:

Coefficients:

Value Error Powers

0 0.07638 0.004339 0 0 0
1 0.6954 0.01052 0 1 0
2 -0.05692 0.01763 1 0 0
3 -0.03708 0.04539 2 0 0
4 0.8 0.08127 3 0 0
5 -0.08459 0.00617 0 0 1
6 -0.1615 0.0393 1 1 0
7 1.211 0.08197 4 0 0
8 -0.4356 0.09412 2 1 0
9 0.1755 0.01426 0 2 0
10 0.9 0.1115 3 1 0
11 1.506 0.1815 4 1 0
12 0.354 0.1018 5 0 0
13 -0.02523 0.01405 0 1 1
14 0.002813 0.01556 1 0 1
15 0.06147 0.0541 1 1 1
16 0.1078 0.1237 4 2 0
17 0.06847 0.06925 5 0 1
18 0.05803 0.06148 1 2 1
19 -0.03643 0.2082 2 2 1
20 0.01925 0.3022 4 2 1
21 -0.1857 0.222 5 2 1
22 0.1174 0.1063 4 0 1
23 0.03153 0.06974 2 1 1
24 0.00348 0.02508 0 2 1
25 -0.04514 0.1454 3 1 1
26 -0.05866 0.06294 2 0 1
27 0.06201 0.1514 3 2 0
28 -0.02339 0.06081 1 2 0

And this I read from the MDF.C:

// Assignment to coefficients vector.
static double gCoefficient[] = {
0.0763792,
0.695381,
-0.056921,
-0.0370823,
0.800029,
-0.0845872,
-0.16145,
1.21051,
-0.435628,
0.1755,
0.900025,
1.50641,
0.354029,
-0.0252297,
0.00281289,
0.0614668,
0.107766,
0.0684715,
0.0580282,
-0.0364296,
0.019249,
-0.185716,
0.117379,
0.0315298,
0.00348003,
-0.0451418,
-0.0586608,
0.0620083,
-0.0233934
};

// Assignment to powers vector.
// The powers are stored row-wise, that is
// p_ij = gPower[i * NVariables + j];
static int gPower[] = {
1, 1, 1,
1, 2, 1,
2, 1, 1,
3, 1, 1,
4, 1, 1,
1, 1, 2,
2, 2, 1,
5, 1, 1,
3, 2, 1,
1, 3, 1,
4, 2, 1,
5, 2, 1,
6, 1, 1,
1, 2, 2,
2, 1, 2,
2, 2, 2,
5, 3, 1,
6, 1, 2,
2, 3, 2,
3, 3, 2,
5, 3, 2,
6, 3, 2,
5, 1, 2,
3, 2, 2,
1, 3, 2,
4, 2, 2,
3, 1, 2,
4, 3, 1,
2, 3, 1
};

You can see, the powers in the MDF.C are always exactly one integer bigger than what ROOT gives out onto screen.

Is that supposed to be that way?

Thanks,
Chris

9

to avoid confusion, as you can see I’m using three variables. It’s hard to distinguish the powers in my pasted expression for the ROOT output at the first look. Sorry for that. But it should be obvious what I mean.

cheers,
Chris

10

Hi Chris,

You say

That’s how it’s supposed to be. If you look at how the terms are actually calculated in the produced script you’ll see why. The reason it was done like this was to simplify the code :slight_smile: You see, the same code is reused for Chebyshev and Legendre polynomials, and there it makes a difference.

As you’ll see in the code, to calculate the term (before multiplying by the coefficient), we do

switch(power) { case 1: r = 1; break; case 2: r = v; break; default: p2 = v; for (k = 3; k <= power; k++) { p3 = p2 * v; p1 = p2; p2 = p3; } r = p3; }

That is, if power==1 (that is x_i^0), then we simply set the term to be 1. If power==2 (that is x_i^1), then the term is set to x_i. For higher order terms, and algorithm that ensures too large numbers are not multiplied is used. Suppose power==5 (x_i^4), p2=x_i initially, and the loop goes from 3 to 5 - that is 3 iterations where p2 is updated with a factor x_i, meaning we calculate exactly x_i^4.

So there’s meaning to the madness.

I hope that clarifies things for you.

Thank you for the kind words.

Yours,

Christian

11

Hi,

thanks for the quick reply, that answered my question. These words are not just kind. I try to fit a function onto a multidim dataset since weeks/months now and just recently discovered that there exists this package. If you want I give you feedback how well my fit and the fit made by your package agrees.

cheers,
Chris

12

to avoid more confusion. I mean I’m really impressed. No kidding…