Nested fit parameters

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 :slight_smile:
Cheers,
Nikkie


_ROOT Version:6
_Platform: CentOS CERN 8
_Compiler: gcc


Hello @NikkieD,

I am inviting @moneta to this post. I think he may be able to answer your question.

Cheers,
J.

Hi,

I have not understood the problem. When you fix the parameter 0 , 'Temperature, it is fixed for the TF1 object you are considering, in your case for the function, critCurrent. Now since critCurrent is implemented by two sub-functions,reduceFieldandreduceTemperature` , it will be fixed for those functions two.
If you want to create separate TF1 objects for those sub-functions, you would need to fix the parameter for those objects as well

Lorenzo

Hi @moneta

Thanks for you reply. So I’m indeed fixing parameter 0, temperature, and I thought it would be fixed for the other two function that are used to build critCurrent, but when I print out the values of each function, it looks like the temperature, when inside the other two functions is actually being optimized.

Here I have a sample of the output, that I get by using the std::cout statements in each funciton:

4.8425 reducedField
4.8425 reducedTemperature
4.2 reducedTemperature
4.2 criticalCurrent
4.8775 reducedField
4.8775 reducedTemperature
4.2 reducedTemperature
4.2 criticalCurrent
4.9125 reducedField
4.9125 reducedTemperature
4.2 reducedTemperature
4.2 criticalCurrent
4.9475 reducedField
4.9475 reducedTemperature
4.2 reducedTemperature
4.2 criticalCurrent
4.9825 reducedField
4.9825 reducedTemperature
4.2 reducedTemperature

If they would be fixed, I would expect the temperature always to be printed as 4.2, instead you see it converging in the reduced field and reduced temperature functions.

Is this what is really happening, or is this just an artifact of me putting the std::cout print statement in the wrong place?

BTW if you comment those print statements out, you’ll get the fit result at the bottom of your screen. Otherwise It might be hidden in the print statements, at least for me…

Cheers,
Nikkie

double b = reducedField(bfield, par);

1 Like

Thanks a lot, I’m sorry that was stupid of me :slight_smile:

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.