TMinuit

From: Nevzat Guler <nguler_at_jlab.org>
Date: Sun, 29 Jan 2006 01:43:13 -0500 (EST)

	Hello,
	I am having some problems with TMinuit. I took the example Ifit.C
and modified it for myself and incorporated into my program.
	In my root program I am looping over events in a tree. Therefore I

cannot send you a real working version of the program but I am sending the program itself. The ntuple is really big so I cannot send it. The program is below (also in attachment). In the attachment there is also a header file
for it and the "output" of the program. Currently it does not work. MIGRAD fails to make improvement and I keep getting the following error message:

Error: G__CallFunc::SetArgArray() must be initialized with 'G__CallFunc::SetFunc(G__ClassInfo* cls,char* fname,char* args,long* poffset)' first

I don't know if this error is causing the program not to work or if the function that I am fitting is terribly bad.

        Thank you for your inputs.

znevzat.C

---

#define znevzat_cxx
#include "TROOT.h"
#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TStopwatch.h"
#include "TTree.h"
#include "TFile.h"
#include "TNtuple.h"
#include "TChain.h"
#include "TStyle.h"
#include "TMath.h"
#include "TFormula.h"
#include "TH1.h"
#include "TF1.h"

gSystem->Load("libMinuit");
#include "TMinuit.h"

#include "znevzat.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

const Int_t ndata=20000;

Double_t deltaPHI[ndata];
Double_t deltaZ[ndata];
Double_t errorDZ[ndata];
Double_t deltaTHETA[ndata];
Double_t sigmaDZ[300];
Double_t thetabin[50];
Double_t mombin[50];
Double_t rad2deg=180.0/3.14159;

//______________________________________________________________________________

Double_t func(Double_t Df,Double_t *par) {
 Double_t value = (par[0]*Df+par[1]);
 return value;
}

//============================================================================================

void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag, Int_t ndata)
{

        Int_t i;

////calculate chisquare

	Double_t chisq=0;
	Double_t chi;
	for(i=0;i<ndata;i++) {
		chi=(deltaZ[i]-func(deltaPHI[i],par));
		chisq+=chi*chi;
	}

   f = chisq;
}

//=============================================================================================

void znevzat::Loop()
{

   if (fChain == 0) return;
   Int_t nentries = Int_t(fChain->GetEntriesFast());

//====================== DEFINE
THINGS===============================================

TH2D *hh = new TH2D("hh", "temporary histogram", 100, -1., 1.,100,-1.,1.);

for(Int_t j=0;j<ndata;j++){
deltaPHI[j]=deltaZ[j]=errorDZ[j]=deltaTHETA[j]=0.0;}
//==============================================================================

   Int_t nbytes = 0, nb = 0;
   //for (Int_t jentry=0; jentry<nentries;jentry++) {

   for (Int_t jentry=0; jentry<ndata;jentry++) {
//====================== FILL ARRAYS and HISTOGRAMS
=============================


if(npro==1&&abs(prz-gprz)<1.0&&abs(prf-gprf)<0.5&&abs(gprm-0.95)<0.1&&abs(57.*prt-35.)<10.){
	hh->Fill((prz-gprz)/prm,prf-gprf);
//	hh->Fill(gprt*rad2deg,gprm);
	}

	deltaPHI[jentry]=prf-gprf;
	deltaZ[jentry]=prz-gprz;
	deltaTHETA[jentry]=prt-gprt;
	errorDZ[jentry]=prz-rprz;

//===============================================================================
      Int_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;

   }

//====================== ANALYZE RESULTS
========================================

TCanvas *c1 = new TCanvas("c1","",0,0,600,400); gStyle->SetOptFit(0111);
gPad->SetGrid();

hh->Draw();
TProfile *prof = hh->ProfileX();

prof->SetLineColor(kRed);
prof->Draw("same");
prof->Draw("");
prof->Fit("pol1","R","",-0.6.,0.6);

TMinuit *gMinuit = new TMinuit(2);
gMinuit->SetFCN(fcn);

Double_t arglist[10];
Int_t ierflg = 0;
arglist[0] = 1;
gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);

// Set starting values and step sizes for parameters

   static Double_t vstart[2] = {0.0 , 0.0};

   static Double_t step[2] = {0.0001 , 0.0001};
   gMinuit->mnparm(0, "a1", vstart[0], step[0], 0,0,ierflg);
   gMinuit->mnparm(1, "a2", vstart[1], step[1], 0,0,ierflg);

// Now ready for minimization step

   arglist[0] = 100000000;
   arglist[1] = 0.001.;
   gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);

// Print results

   Double_t amin,edm,errdef;
   Int_t nvpar,nparx,icstat;
   gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);    gMinuit->mnprin(3,amin);
}

---

Received on Sun Jan 29 2006 - 07:43:28 MET