Hi,
I have an oscillation (simple cosine wave) that I try to fit. When I use a chi^2 minimization, the fit gives me what I expect. When I switch to log-likelihood, the fit fails.
Any ideas of what I may be doing wrong?
Cheers,
–Christos
// histogram constants
const float LOW = 0.0;
const float HI = 3500.0;
const int Nbins = 100;
const float bin_size = (HI - LOW)/Nbins;
// lifetime, mass difference, speed of light
const float ctau = 461.7; // in micrometers
const float DeltaM = 0.5; // in ps^{-1}
const float speed_of_light = 299.792; // in microns/ps
// proper time distribution for unmixed events
inline Double_t pure_unm(Double_t * x, Double_t * par)
{
/* par[0]: # of events, par[1]: c*lifetime, par[2]: delta-m */
return (par[0]*bin_size/(2*par[1]))*TMath::Exp(-x[0]/par[1])*
(1 + TMath::Cos(par[2]*x[0]/speed_of_light));
}
// proper time distribution for mixed events
inline Double_t pure_mix(Double_t * x, Double_t * par)
{
/* par[0]: # of events, par[1]: c*lifetime, par[2]: delta-m */
return (par[0]*bin_size/(2*par[1]))*TMath::Exp(-x[0]/par[1])*
(1 - TMath::Cos(par[2]*x[0]/speed_of_light));
}
// proper time distribution for asymmetry
inline Double_t pure_asym(Double_t * x, Double_t * par)
{
/* par[0]: amplitude(=1), par[1]: delta-m */
return (par[0])*TMath::Cos(par[1]*x[0]/speed_of_light);
}
int fit_asym(void)
{
// setup histograms
TH1F * h21 = new TH1F("unm","unmixed", Nbins, LOW, HI);
TH1F * h22 = new TH1F("mix","mixed", Nbins, LOW, HI);
h21->GetXaxis()->SetTitle("c_{ }#tau (#mum)");
h22->GetXaxis()->SetTitle("c_{ }#tau (#mum)");
// setup distributions
f21 = new TF1("funm", pure_unm, LOW, HI, 3);
f22 = new TF1("fmix", pure_mix, LOW, HI, 3);
f23 = new TF1("fasym", pure_asym, LOW, HI, 2);
f21->SetParameters(1, ctau, DeltaM);
f22->SetParameters(1, ctau, DeltaM);
// fraction of mixed events
float chi_d = f22->Integral(LOW, HI)/
(f21->Integral(LOW, HI) + f22->Integral(LOW, HI));
// fill histograms
for(int i = 0; i != 100000; ++i)
{
if(gRandom->Rndm() < chi_d)
h22->Fill(f22->GetRandom());
else
h21->Fill(f21->GetRandom());
}
// get asymmetry
h23 = new TH1F(* (TH1F*) h21->GetAsymmetry(h22) );
h23->SetNameTitle("asym","asymmetry");
h23->GetXaxis()->SetTitle("c_{ }#tau (#mum)");
TH1F * h10 = (TH1F *) (h23->Clone());
TH1F * h20 = (TH1F *) (h23->Clone());
// fit - plot histograms
TCanvas * c1 = new TCanvas();
c1->Divide(1, 2);
c1->cd(1);
f23->SetParameters(1, DeltaM);
f23->SetParNames("Amplitude", "#Deltam");
// chi2 minimization
h10->Fit(f23, "RIV");
c1->cd(2);
f23->SetParameters(1, DeltaM);
f23->SetParNames("Amplitude", "#Deltam");
// log-likelihood minimization
h20->Fit(f23, "RIVLL");
return 0;
}