Hi,
Is there a simple public example of using complex numbers in ROOT? I am trying to fit a cross-section for two Breit-Wigner’s plus a constant term. I did the single Breit-Wigner using real numbers, and that works fine, but the code would be clearer using complex numbers. But, when I try to run it, ROOT goes into an infinite loop. This is ROOT 5.28 - old, but probably not ancient. I could update if that would help. The code follows; I have commented out the parts that pertain to the second Breit-Wigner, so it should be equivalent to the old version, which worked.
A pointer to any mistakes, and/or to a simple ROOT example using complex numbers would be appreciated. I searched the ROOT documentation, but the examples I found were pretty heavy-duty.
Thanks.
Spencer Klein, Berkeley
double fitf(Double_t *v, Double_t *par){
#include <TComplex.h>
/* This reads histograms from Bill's TTree and fits for a rho, direct pipi and eventually an omega signal.
April 9, 2015*/
// This is the fit function for the rho
// v[0] is the dipion mass
// par is the array containing the parameters being fitted
// par[0] is the rho mass
// par[1] is the rho width
// par[2] is 'A;
//par[3] is 'B;
// par[4] is the constant background
// par[5] is the linear (in mass)background
// par[6] is reserved for the omega
double mpipi=v[0];
double mass2=mpipi*mpipi;
double polybg=par[4]+par[5]*mpipi;
// To add the omega, switch to complex numbers
double rhomass=par[0];
double rhomass2=rhomass*rhomass;
double rhowidth=par[1];
double pionmass=0.13957;
double pionmass2=pionmass*pionmass;
// rho mass-dependent width
double val=pow((mass2-4*pionmass2)/(rhomass2-4*pionmass2),1.5);
double mdwidth=rhowidth*(rhomass/mpipi)*val; // mass dependent width
// now omega, with mass-dependent width
// This is probably too narrow, since it neglects detector broadening
double omegamass=0.7872;
double omegamass2=omegamass*omegamass;
double omegawidth=0.085;
double omegaval= val=pow((mass2-4*pionmass2)/(omegamass2-4*pionmass2),1.5);
double omegamdwidth=omegawidth*(omegamass/mpipi)*val;
// try it with complex numbers
std::complex<double> ix=(0,1);
std::complex<double> denom=(mass2-rhomass2)+ix*rhomass*mdwidth;
std::complex<double> aterm=par[2]*sqrt(mpipi*rhomass*mdwidth)/denom;
std::complex<double> bterm=par[3];
// std::complex<double> omegadenom=(mass2-omegamass2)+ix*omegamass*omegamdwidth;
std::complex<double>cterm=0.;
// std::complex<double> cterm=par[6]*sqrt(mpipi*omegamass*omegamdwidth)/omegadenom;
std::complex<double> total=aterm+bterm+cterm;
double amp=abs(total);
double fitval=amp+polybg;
return fitval;
}
Hi Spencer,
could you share the code in its entirety and some of the data to run on?
Cheers,
Danilo
The code is used to fit a histogram. The pertinent parts are:
TH1F *Tnet = (TH1F*)rhosig->Clone("Tnet");
// Now do the fit
TF1 *func = new TF1("func",fitf,0.5,1.3,7);
// This fit will do much better with some initial guesses
func->SetParameters(0.768,0.151,25.,-100,0.,0.,0.);
func->SetParLimits(0,0.74,0.78);
func->SetParLimits(1,0.1,0.2);
Tnet->Fit("func");
This all worked OK when func has 6 free parameters (no omega==par[6]), and was done using all real numbers.
Spencer
Hi Spencer,
The usage of complex numbers for ROOT in this case can hardly be a problem as it is confined in your fitting function.
What do you mean with infinite loop? Can you share a reproducer?
Cheers,
Danilo
When I run this code with an all-real version of the fitting code, I get a nice fit that reproduces the data. When I try the complex version, ROOT runs indefinitely (I’ve waited ~ 1 minute, while the all-real version fits in under a second), occupying 25% of the CPU on my 4-processor machine (i.e one CPU is fully engaged). I do not want to post the histogram file I’m using, since it involves private STAR data, but I’d be happy to send it to you privately. That, plus the code I’ve already posted should be enough to reproduce the problem.
I have verified that I can call the fit function in a simple macro:
TF1 *func = new TF1("func",fitf,0.5,1.3,7);
func->SetParameters(0.768,0.151,25.,-100,0.,0.,0.);
cout<< func->Eval(1.0)<<endl;
and it returns a reasonable answer, so I am also a bit confused. Does ROOT have a mechanism to detect that a fit it not converging, and stop execution? Is it possible to give it a fit function that flails indefinitely, without producing any output, or are there limits to ensure convergence, and/or giving up?
Spencer
Hi,
this looks like a non converging fit.
It is hard to say without a reproducer - it is of course understood that there might be some privacy issues involved, but the reproducer of the bug must not include the accurate numbers coming from star.
Did you check that the implementation which involves complex numbers gives similar results to the one with real numbers in the 0.5-1.3 range? Are there any discontinuities?
Did you try to print the values of the parameters during the fit in the two cases to see when they started diverging?
Hi again:
I don’t think that it’s a non-converging fit. I played with the code a bit more, so it just does a linear fit (N=a+b*mass) to the histogram, with the other parameters fixed. When I commented out all of the complex number code in fitf.C (i.e. from ix=(0,1) to double amp =abs(total)), then the code ran in under a second. When I uncommented the code, so that the calculations were done, but their results were not used (i.e. setting amp=0), then the code ran and returned theidentical result, but took 3 minutes (same dedicated single CPU…). So, maybe the problem is that the complex evaluation might work, but takes an unacceptably long time. I could understand the complex version being a few times slower, but 100X slower isn’t acceptable. If the full, 6-7 parameter fit takes 15 seconds in the double version, this would be 30 minutes with complex numbers - too long to be useful. So, unless you have a better idea, I’ll go back to using doubles.
Thanks.
Spencer
Hi Spencer,
I am sorry but without a reproducer and/or a comparison of the values of the two functions for the values in the interval of interest, I would not know how to help you further.
Danilo
#include <cmath>
#include <complex>
double fitf(double *v, double *par) {
/* This reads histograms from Bill's TTree and fits for a rho,
direct pipi and eventually an omega signal.
April 9, 2015*/
// This is the fit function for the rho
// v[0] is the dipion mass
// par is the array containing the parameters being fitted
// par[0] is the rho mass
// par[1] is the rho width
// par[2] is 'A;
// par[3] is 'B;
// par[4] is the constant background
// par[5] is the linear (in mass)background
// par[6] is reserved for the omega
double mpipi=v[0];
double mass2=mpipi*mpipi;
double polybg=par[4]+par[5]*mpipi;
// To add the omega, switch to complex numbers
double rhomass=par[0];
double rhomass2=rhomass*rhomass;
double rhowidth=par[1];
double pionmass=0.13957;
double pionmass2=pionmass*pionmass;
// rho mass-dependent width
double val=std::pow((mass2-4*pionmass2)/(rhomass2-4*pionmass2),1.5);
double mdwidth=rhowidth*(rhomass/mpipi)*val; // mass dependent width
// now omega, with mass-dependent width
// This is probably too narrow, since it neglects detector broadening
double omegamass=0.7872;
double omegamass2=omegamass*omegamass;
double omegawidth=0.085;
double omegaval=std::pow((mass2-4*pionmass2)/(omegamass2-4*pionmass2),1.5);
double omegamdwidth=omegawidth*(omegamass/mpipi)*omegaval;
// try it with complex numbers
std::complex<double> ix=std::complex<double>(0,1);
std::complex<double> denom=(mass2-rhomass2)+ix*rhomass*mdwidth;
std::complex<double> aterm=par[2]*std::sqrt(mpipi*rhomass*mdwidth)/denom;
std::complex<double> bterm=par[3];
// std::complex<double> omegadenom=(mass2-omegamass2)+ix*omegamass*omegamdwidth;
std::complex<double> cterm=0.;
// std::complex<double> cterm=par[6]*std::sqrt(mpipi*omegamass*omegamdwidth)/omegadenom;
std::complex<double> total=aterm+bterm+cterm;
double amp=std::abs(total);
double fitval=amp+polybg;
return fitval;
}
Hi,
Is there a clear ROOT example showing how to efficiently use complex numbers in a simple way?
I tried changing my fitf.C to include the headers you listed:
double fitf(Double_t *v, Double_t *par){
#include <cmath>
#include <complex>
etc.
However, when I run it, for each call, I get the messages:
Note: File “cmath” already loaded
Note: File “complex” already loaded
This repeated loading does not seem to lead to efficient execution. For reference, I get this when I type
.L fitf.C
.x omegafit.C
using code as in my earlier mesage.
I also get these notes when I load fitf.C and then execute a simple macro:
{
gROOT->Reset();
#include <TMath.h>
#include <iostream>
using namespace std;
cout <<"Starting program."<<endl;
TF1 *func = new TF1("func",fitf,0.5,1.3,7);
func->SetParameters(0.768,0.151,25.,-100,0.,0.,0.);
cout<< func->Eval(1.0)<<endl;
}
Spencer
In my previous post here, you can see where you’re supposed to put “#include” directives.