Dear All,
I’m trying to fit a complicated function by making the function in stages. Then I want to fix one parameter, because that one was measured to be a certain value. However, this parameter appears in more than one of the functions that are used to build the complete fit function, and then I see that for the “final fit” the parameter is indeed fixed, but for the other functions it is being fitted/optimized.
Is there a way to fix this parameter for all functions?
My script is the following:
// c++ libraries
#include <iostream>
#include <cmath>
// Command line parser library
#include <tclap/CmdLine.h>
// Root librabries
#include "TF1.h"
#include "TH1.h"
#include "TGraph.h"
#include "TFitResult.h"
#include "TTree.h"
#include "TLine.h"
#include "TLatex.h"
#include "TPad.h"
#include "TCanvas.h"
// Fit settings
const double load = 3.9/2.35; // Load in [T / kA]
const double criticalTemperature = 9.2; // Critical temperature of NbTi in [K]
const double upperCriticalField = 14.4; // Maximum upper critical field at T=0K [T]
const double n = 1.7; // Some exponent
const double C0 = 37.7; // Normalization constant
const double fAlpha = 0.89; // Fit parameter
const double fBeta = 1.1; // Fit parameter
const double fGamma = 2.09; // Fit parameter
const double preFactorINR = 3.6; // Pre-factor of INR to scale to Bottura fit
// Plot settings
const int loadLineColor = 861; // Color of the load line
const int critSurfColor = 801; // Color of the crit surf line
const int lineWidth = 2; // Line width for plots
TF1* critCurrent;
TF1* loadline;
// Load line
double calcLoadLine ( double* x, double* par ) {
// if you want current: par[0] = load
// else if you want Bfield: par[0] = 1/load
return x[0] * par[0];
} // end double loadLine ()
// Reduced temperature, input for critical surface
double reducedTemperature ( double* par ) {
// par[0] = temperature
std::cout << par[0] << " reducedTemperature" << std::endl;
return par[0] / criticalTemperature;
}
// Reduced field, input for critical surface
double reducedField ( double* bfield, double* par ) {
// par[0] = temperature [K]
std::cout << par[0] << " reducedField" << std::endl;
double bC2 = upperCriticalField * ( 1 - std::pow(reducedTemperature(par), n));
return bfield[0] / bC2;
}
// Critical surface, I[B]
double criticalCurrent ( double* bfield, double* par ) {
// par[0] = temperature [K]
// par[1] = Norm (C0 has been removed)
// par[2] = fAlpha
// par[3] = fBeta
// par[4] = fGamma
std::cout << par[0] << " criticalCurrent" << std::endl;
double b = reducedField(par, bfield);
double t = reducedTemperature(par);
return par[1] * (1/bfield[0]) * (std::pow(b, par[2])) * (std::pow((1 - b), par[3])) * (std::pow((1 - std::pow(t, n)), par[4]));
}
// Intersect function
double findIntersect ( double* bfield, double* par ) {
return TMath::Abs(critCurrent->EvalPar(bfield, par) - loadline->EvalPar(bfield, par));
}
TGraph* measuredCriticalCurrent( int sample ) {
// Sample/I:Bext/D:B_ec5/D:Ic_ec5/D:n_ec5/D:B_ec4/D:Ic_ec4/D:n_ec4/D:B_rho/D:IC_rho/D
TTree *tree = new TTree();
tree->ReadFile("melcData.txt", "", '\t');
tree->Draw("Ic_ec5:B_ec5", Form("Sample==%d", sample));
auto gtemp = (TGraph*)gPad->GetPrimitive("Graph");
// Make a nicer graph
gtemp->SetTitle(";B [T]; I[B]");
gtemp->SetMarkerStyle(23);
gtemp->GetXaxis()->SetLimits(0., 10.);
gtemp->SetMinimum(0.);
gtemp->SetMaximum(10000.);
return gtemp;
}
// Quick test script for the ROOT interpreter
void operatingCurrentField () {
double temperature = 4.2;
// Get the data from the Ic measurements:
TGraph* sample24Hist = measuredCriticalCurrent(18);
// Now I want to plot this
//TCanvas *c0 = new TCanvas();
//sample24Hist->Draw("ap");
// And fit this
critCurrent = new TF1("critCurr", criticalCurrent, 0.5, 15., 5);
critCurrent->SetLineColor(critSurfColor);
critCurrent->SetLineWidth(lineWidth);
critCurrent->SetParName(0, "Temperature");
critCurrent->SetParName(1, "Normalization");
critCurrent->SetParName(2, "Alpha");
critCurrent->SetParName(3, "Beta");
critCurrent->SetParName(4, "Gamma");
critCurrent->FixParameter(0, temperature);
//critCurrent->SetParameter("Temperature", 4.2);
//critCurrent->SetParLimits(0, 3.5, 4.7);
critCurrent->SetParameter(1, 4000.*37);
//critCurrent->SetParLimits(1, 3000., 5000.);
// critCurrent->SetParLimits(2, 0., 10.);
// critCurrent->SetParLimits(3, 0., 20.);
// critCurrent->SetParLimits(4, 0., 30.);
// critCurrent->FixParameter(2, fAlpha);
// critCurrent->FixParameter(3, fBeta);
// critCurrent->FixParameter(4, fGamma);
// critCurrent->SetParameter("Alpha", fAlpha);
// critCurrent->SetParameter("Beta", fBeta);
// critCurrent->SetParameter("Gamma", fGamma);
TFitResultPtr fitR = sample24Hist->Fit(critCurrent, "SR", "", 1., 5.);
fitR->Print("V");
// double preFactor = preFactorINR;
// // Make the critical surface function
} // end void operatingCurrentField
The parameter that I wanted to fix was temperature, parameter 0 in all the functions.
The data file is given here: melcData.txt (1.6 KB)
I hope this is possible to do
Cheers,
Nikkie
_ROOT Version:6
_Platform: CentOS CERN 8
_Compiler: gcc