Dear All,
We are trying to fit a function that consists of other functions, however this function is not continuous. Therefore we put some constraints such that there are the appropriate gaps. Then we fit and it almost works, however the two parameters that we are interested in stay zero with infinite uncertainty. Moreover, if we try to set those parameters before the fit, the fit fails miserably.
The code looks like this (see below), and you need the txt file that we added in order to get the data. Noise_scan_100events_1.txt (673 Bytes)
Now the question is, why does it fail? When it fits the four functions separately the fits look great (see plot that is produced by the code) but when we fit them all together it ends up looking like a sum of four straight lines. The reason why we want to do this is because in reality we want to constrain the four fits of the individual lines in such a way that par[2] of f0 == par[2] of f2 and par[2] of f1 == par[2] of f3.
Thanks a lot in advance,
Nikkie
#include <iostream>
#include <fstream>
#include "TCanvas.h"
#include "TMultiGraph.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TLegend.h"
#include "TF1.h"
#include "TAxis.h"
double _x, _y0, _y1, _y2, _y3, X0[25], X1[25], X2[25], X3[25], Y0[25], Y1[25], Y2[25], Y3[25];
double xMin = 0, xMax = 250;
int nMeasurements = 25;
double errX[25], errY[25];
double noiseFunction (double *x, double *par) {
return par[0]+par[1]*((x[0]<par[2])*sqrt(par[2]/x[0]) + (x[0]>=par[2]));
}//end noiseFunction
void getData (const char *fileName = "Noise_scan_100events_1.txt") {
std::ifstream in_file1;
in_file1.open(fileName);
for (int iMeasurement=0; iMeasurement<nMeasurements; iMeasurement++) {
in_file1 >> _x >> _y0 >> _y1 >> _y2 >> _y3;
X0[iMeasurement] = _x;
X1[iMeasurement] = _x + 250;
X2[iMeasurement] = _x + 500;
X3[iMeasurement] = _x + 750;
Y0[iMeasurement] = _y0;
Y1[iMeasurement] = _y1;
Y2[iMeasurement] = _y2;
Y3[iMeasurement] = _y3;
errX[iMeasurement] = 0.1;
errY[iMeasurement] = 0.1;
}//end for i
in_file1.close();
}//end void getData ()
void biasVoltageFit () {
getData ();
TCanvas *c0 = new TCanvas();
TMultiGraph *mg0 = new TMultiGraph("mg1","Noise scan");
mg0->SetTitle("; #bf{Voltage (-V)}; Noise");
TGraphErrors *gr0 = new TGraphErrors(nMeasurements, X0, Y0, errX, errY);
gr0->SetLineColor(2);
gr0->SetLineWidth(0);
gr0->SetMarkerColor(2);
gr0->SetMarkerStyle(20);
gr0->SetMarkerSize(1.5);
TGraphErrors *gr1 = new TGraphErrors(nMeasurements, X1, Y1, errX, errY);
gr1->SetMarkerStyle(23);
gr1->SetMarkerSize(1.5);
gr1->SetMarkerColor(6);
gr1->SetLineColor(4);
gr1->SetLineWidth(0);
TGraphErrors *gr2 = new TGraphErrors(nMeasurements, X2, Y2, errX, errY);
gr2->SetMarkerStyle(21);
gr2->SetMarkerSize(1.5);
gr2->SetMarkerColor(4);
gr2->SetLineColor(4);
gr2->SetLineWidth(0);
TGraphErrors *gr3 = new TGraphErrors(nMeasurements, X3, Y3, errX, errY);
gr3->SetLineColor(3);
gr3->SetLineWidth(0);
gr3->SetMarkerColor(3);
gr3->SetMarkerStyle(22);
gr3->SetMarkerSize(1.5);
mg0->Add(gr0);
mg0->Add(gr1);
mg0->Add(gr2);
mg0->Add(gr3);
mg0->Draw("AP");
TLegend *l0= new TLegend(0.85,0.85,0.6,0.65);
l0->AddEntry(gr0,"Fe0CBC0_even","p");
l0->AddEntry(gr1,"Fe0CBC0_odd","p");
l0->AddEntry(gr2,"Fe1CBC1_even","p");
l0->AddEntry(gr3,"Fe1CBC1_odd","p");
l0->Draw();
TF1 *f0 = new TF1("f0", "[0]+[1]*((x<[2])*sqrt([2]/x) + (x>=[2]))", 10, 250);
TF1 *f1 = new TF1("f1", "[0]+[1]*(((x-250)<[2])*sqrt(([2])/(x-250)) + ((x-250)>=[2]))", 260, 500);
TF1 *f2 = new TF1("f2", "[0]+[1]*(((x-500)<[2])*sqrt([2]/(x-500)) + ((x-500)>=[2]))", 510, 750);
TF1 *f3 = new TF1("f3", "[0]+[1]*(((x-750)<[2])*sqrt([2]/(x-750)) + ((x-750)>=[2]))", 760, 1000);
TF1 *total = new TF1("total", "(([0]+[1]*((x<[2])*sqrt([2]/x) + (x>=[2])))*(x>10)*(x<=250) + ([3]+[4]*(((x-250)<[9])*sqrt([9]/(x-250)) + ((x-250)>=[9])))*(x>260)*(x<=500) + ([5]+[6]*(((x-500)<[2])*sqrt([2]/(x-500)) + ((x-500)>=[2])))*(x>510)*(x<=750) + ([7]+[8]*(((x-750)<[9])*sqrt([9]/(x-750)) + ((x-750)>=[9])))*(x>760)*(x<=1000))", 10, 250);
total->SetParLimits(0, 0.01, 10);
total->SetParLimits(1, 0.01, 10);
//total->SetParLimits(2, 40, 100);
total->SetParLimits(3, 0.01, 10);
total->SetParLimits(4, 0.01, 10);
total->SetParLimits(5, 0.01, 10);
total->SetParLimits(6, 0.01, 10);
total->SetParLimits(7, 0.01, 10);
total->SetParLimits(8, 0.01, 10);
//total->SetParLimits(9, 40, 100);
total->Print();
f0->SetParLimits(0, 0, 10);
f0->SetParLimits(1, 0, 10);
f0->SetParLimits(2, 40, 100);
f0->SetParameter(0, 1);
f0->SetParameter(1, 1.5);
f0->SetParameter(2, 75);
f1->SetParLimits(0, 0, 10);
f1->SetParLimits(1, 0, 10);
f1->SetParLimits(2, 40, 100);
f1->SetParameter(0, 1);
f1->SetParameter(1, 1.5);
f1->SetParameter(2, 75);
f2->SetParLimits(0, 0, 10);
f2->SetParLimits(1, 0, 10);
f2->SetParLimits(2, 40, 100);
f2->SetParameter(0, 1);
f2->SetParameter(1, 1.5);
f2->SetParameter(2, 75);
f3->SetParLimits(0, 0, 10);
f3->SetParLimits(1, 0, 10);
f3->SetParLimits(2, 40, 100);
f3->SetParameter(0, 1);
f3->SetParameter(1, 1.5);
f3->SetParameter(2, 75);
mg0->Fit(f0, "R");
mg0->Fit(f1, "R+");
mg0->Fit(f2, "R+");
mg0->Fit(f3, "R+");
total->SetParameter(0, f0->GetParameter(0));
total->SetParameter(1, f0->GetParameter(1));
//total->SetParameter(2, 64);
total->SetParameter(3, f1->GetParameter(0));
total->SetParameter(4, f1->GetParameter(1));
total->SetParameter(5, f2->GetParameter(0));
total->SetParameter(6, f2->GetParameter(1));
total->SetParameter(7, f3->GetParameter(0));
total->SetParameter(8, f3->GetParameter(1));
//total->SetParameter(9, 64);
mg0->Fit(total, "R+");
}//end void biasVoltageFit ()