I’m trying to use Lorenzo Montana’s code to fit 3D data into a line but I get too many errors. The first part of this code plots my data. The second part does the fitting. Any suggestions would be highly appreciated.
#define Example_cxx
#include "Example.h"
#include <TH2.h>
#include <TH3F.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TPolyLine3D.h>
#include <TH2F.h>
#include <TCanvas.h>
#include <TLegend.h>
void Example::Loop()
{
TCanvas *c1=new TCanvas("c1","c1",800,600);
// 3D fit
TH3F* hd3 = new TH3F("h3", "h2", 100, 0, 70, 100, 0, 70, 100, 0, 70);
// ### These functions check if the ROOT file is loaded properly ###
if (fChain == 0) return;
// ### This gets the total number of entries (e.g. events) stored in the ROOT file
// ### and event is one simulated particle passing through the liquid argon
Long64_t nentries = fChain->GetEntriesFast();
// ### These are variables we won't use, but ROOT defines for us
Long64_t nbytes = 0, nb = 0;
// ### This for loop loops over all the events stored in the ROOT file
// ### jentry is the current event number under analyis
//for (Long64_t jentry=0; jentry<nentries;jentry++)
std::cout<<" Total_entries "<<nentries<<std::endl;
for (Long64_t jentry=11; jentry<21;jentry++)
{
// ### Load the information from the tree for this event (e.g. jentry)
Long64_t ientry = LoadTree(jentry);
// ### We will break the code if the returned value of the object in memory is < 0
if (ientry < 0) break;
// #### store the current file information
nb = fChain->GetEntry(jentry); nbytes += nb;
// ### Code starts here ###
if(maxTrackIDE < 1){continue;}
//std::cout<<maxTrackIDE<<endl;
// ### Loop over the IDE's ###
//TPolyLine3D *l1 = new TPolyLine3D(maxTrackIDE);
//std::cout<<"max_track_ide"<<maxTrackIDE<<std::endl;
for(int ide = 1; ide < maxTrackIDE; ide++)
{
//std::cout<<ide<<endl;
// ### This one didn't work ###
//l1->SetPoint(ide,IDEPos[ide][0],IDEPos[ide][1],IDEPos[ide][2]);
// ### This one works fine ###
hd3->Fill(IDEPos[ide][0],IDEPos[ide][1],IDEPos[ide][2]);
std::cout<<"X: "<<IDEPos[ide][0]<<std::endl;
std::cout<<"Y: "<<IDEPos[ide][1]<<std::endl;
std::cout<<"Z: "<<IDEPos[ide][2]<<std::endl;
//hXposition->Fill(IDEPos[ide][0]);
}//<---End ide loop
//l1->SetLineColor(kRed);
//l1->Draw();
//delete l1;
// Step 1: ### Loop over all the track IDE's for one event and store the X, Y, Z position
// ### into a container
}
hd3->Draw();
c1->Update();
}
// define the parametric line equation
void line(double t, const double *p, double &x, double &y, double &z) {
// a parametric line is define from 6 parameters but 4 are independent
// x0,y0,z0,z1,y1,z1 which are the coordinates of two points on the line
// can choose z0 = 0 if line not parallel to x-y plane and z1 = 1;
x = p[0] + p[1]*t;
y = p[2] + p[3]*t;
z = t;
}
bool first = true;
// function Object to be minimized
struct SumDistance2 {
// the TGraph is a data member of the object
TGraph2D *fGraph;
SumDistance2(TGraph2D *g) : fGraph(g) {}
// calculate distance line-point
double distance2(double X,double Y,double Z, const double *p) {
// distance line point is D= | (xp-x0) cross ux |
// where ux is direction of line and x0 is a point in the line (like t = 0)
ROOT::Math::XYZVector xp(x,y,z);
ROOT::Math::XYZVector x0(p[0], p[2], 0. );
ROOT::Math::XYZVector x1(p[0] + p[1], p[2] + p[3], 1. );
ROOT::Math::XYZVector u = (x1-x0).Unit();
double d2 = ((xp-x0).Cross(u)).Mag2();
return d2;
}
// implementation of the function to be minimized
double operator() (const double *par) {
assert(fGraph != 0);
double * x = fGraph->GetX();
double * y = fGraph->GetY();
double * z = fGraph->GetZ();
int npoints = fGraph->GetN();
double sum = 0;
for (int i = 0; i < npoints; ++i) {
double d = distance2(x[i],y[i],z[i],par);
sum += d;
}
if (first) {
std::cout << "Total Initial distance square = " << sum << std::endl;
}
first = false;
return sum;
}
};
Int_t line3Dfit()
{
gStyle->SetOptStat(0);
gStyle->SetOptFit();
//double e = 0.1;
Int_t nd = 10000;
// double xmin = 0; double ymin = 0;
// double xmax = 10; double ymax = 10;
TGraph2D * gr = new TGraph2D();
// Fill the 2D graph
double p0[4] = {10,20,1,2};
// generate graph with the 3d points
for (Int_t N=0; N<nd; N++) {
double x,y,z = 0;
// Generate a random number
double t = gRandom->Uniform(0,10);
line(t,p0,x,y,z);
double err = 1;
// do a gaussian smearing around the points in all coordinates
x += gRandom->Gaus(0,err);
y += gRandom->Gaus(0,err);
z += gRandom->Gaus(0,err);
gr->SetPoint(N,x,y,z);
//dt->SetPointError(N,0,0,err);
}
// fit the graph now
ROOT::Fit::Fitter fitter;
// make the functor objet
SumDistance2 sdist(gr);
ROOT::Math::Functor fcn(sdist,4);
// set the function and the initial parameter values
double pStart[4] = {1,1,1,1};
fitter.SetFCN(fcn,pStart);
// set step sizes different than default ones (0.3 times parameter values)
for (int i = 0; i < 4; ++i) fitter.Config().ParSettings(i).SetStepSize(0.01);
bool ok = fitter.FitFCN();
if (!ok) {
Error("line3Dfit","Line3D Fit failed");
return 1;
}
const ROOT::Fit::FitResult & result = fitter.Result();
std::cout << "Total final distance square " << result.MinFcnValue() << std::endl;
result.Print(std::cout);
gr->Draw("p0");
// get fit parameters
const double * parFit = result.GetParams();
// draw the fitted line
int n = 1000;
double t0 = 0;
double dt = 10;
TPolyLine3D *l = new TPolyLine3D(n);
for (int i = 0; i <n;++i) {
double t = t0+ dt*i/n;
double x,y,z;
line(t,parFit,x,y,z);
l->SetPoint(i,x,y,z);
}
l->SetLineColor(kRed);
l->Draw("same");
// draw original line
TPolyLine3D *l0 = new TPolyLine3D(n);
for (int i = 0; i <n;++i) {
double t = t0+ dt*i/n;
double x,y,z;
line(t,p0,x,y,z);
l0->SetPoint(i,x,y,z);
}
l0->SetLineColor(kBlue);
l0->Draw("same");
return 0;
}
int main() {
return line3Dfit();
}
//hXposition->Draw();
//<---This is the end of the Loop() function/COde
Please read tips for efficient and successful posting and posting code
_ROOT Version:_ Not Provided
_Platform:_ Not Provided
_Compiler:_ Not Provided
___