ClusterFitter3DMinuit.h

Go to the documentation of this file.

#pragma once
 
#include <TH3D.h>
#include <TFile.h>
#include <TMinuit.h>
#include <TFitResult.h>
 
#include "Line.h"
#include "AMSTrdMCTrack.h"
 
#include <vector>
#include <string>
 
//
#include <Math/Functor.h>
#include <TPolyLine3D.h>
#include <Math/Vector3D.h>
#include <Fit/Fitter.h>
 
using namespace ROOT::Math;
 
using std::cerr;
using std::endl;
 
 
// Global variables to store data points for Minuit function
TH3D* gHistogram = nullptr;
 
 
namespace ClusterFitter3DMinuit
{
    
    // Objective function to be minimized
    void objectiveFunction(int& nDim, double* gout, double& result, double params[], int flag)
    {
        result = 0.0;
        for (int ix = 1; ix <= gHistogram->GetNbinsX(); ++ix) {
            for (int iy = 1; iy <= gHistogram->GetNbinsY(); ++iy) {
                for (int iz = 1; iz <= gHistogram->GetNbinsZ(); ++iz) {
                    if (gHistogram->GetBinContent(ix, iy, iz) > 0)
                    {
                        double coords[3];
                        coords[0] = gHistogram->GetXaxis()->GetBinCenter(ix);
                        coords[1] = gHistogram->GetYaxis()->GetBinCenter(iy);
                        coords[2] = gHistogram->GetZaxis()->GetBinCenter(iz);
                        result += Line::distanceFromLine3D(coords, params);
                    }
                }
            }
        }
    }
 
 
    Line::FitLine3DMinuit fitClustersToLine3D(const AMSTrdMCTrack::Track &track, const int& eventID, TFile * file)
    {
        Line::FitLine3DMinuit result;
        result.trkID = track.trkID;
        result.g3PID = track.g3PID;
        
        TH3D * histogram = AMSTrdMCTrack::setClustersInTH3D(track);
        gHistogram = histogram;
        TMinuit minuit(6);
        minuit.SetFCN(objectiveFunction);
 
        // Initialize parameters
        // (Int_t parNo, const char *name, Double_t initVal, Double_t initErr, Double_t lowerLimit, Double_t upperLimit)
        minuit.DefineParameter(0, "px",   0, 0.1, -100, 100);   // positions
        minuit.DefineParameter(1, "py",   0, 0.1, -100, 100);
        minuit.DefineParameter(2, "pz", 100, 0.1,   90, 140);
        minuit.DefineParameter(3, "vx",   1, 0.1, -1000, 1000); // directions
        minuit.DefineParameter(4, "vy",   0, 0.1, -1000, 1000);
        minuit.DefineParameter(5, "vz",   0, 0.1, -1000, 1000);
 
        // Minimize
        minuit.Migrad();
        
        // Extract parameters and errors
        for (int i = 0; i < 6; i++) {
            double paramValue, paramError;
            minuit.GetParameter(i, paramValue, paramError);
            result.fitParameters.push_back(paramValue);
            result.fitErrors.push_back(paramError);
        }
        
        // Optionally, get the chi-square and degrees of freedom
        // result.chi2 = minuit.GetChisquare();
        // result.ndf = minuit.GetNDF();
 
        // Check if the fit was successful
        result.fitSuccess = (minuit.GetStatus() == 0);  // Assuming status 0 indicates a successful fit
        
        printf("\n\nFor track.trkID = %d and g3PID: %d: minuit.GetStatus() = %d\n\n\n", track.trkID, track.g3PID, minuit.GetStatus());
        
        // if (file && result.fitResult)
        {
            file->cd();
            histogram->Write(Form("graph2D_%d_%d_%d_%d", eventID, track.trkID, track.g3PID, track.clusters.size()));
        }
    
        delete histogram;
        return result;
    }
}