//
//-----------------------------------------------------------------
-#include <TMinuit.h>
+#include <TVirtualFitter.h>
#include <TGeoMatrix.h>
#include "AliAlignObj.h"
#include "AliTrackPointArray.h"
#include "AliTrackResidualsChi2.h"
+
+
+void TrackResidualsChi2Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag);
+
+
ClassImp(AliTrackResidualsChi2)
+
//______________________________________________________________________________
-void Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
+void TrackResidualsChi2Fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
// This function is called by minuit
// The corresponding member method is called
// using SetObjectFit/GetObjectFit methods of TMinuit
- AliTrackResidualsChi2* dummy = (AliTrackResidualsChi2 *)gMinuit->GetObjectFit();
+ AliTrackResidualsChi2* dummy = (AliTrackResidualsChi2 *)TVirtualFitter::GetFitter()->GetObjectFit();
dummy->Chi2(npar, gin, f, par, iflag);
}
{
// Implementation of Chi2 based minimization
// of track residuala sum
- TMinuit *gMinuit = new TMinuit(6); //initialize TMinuit
- gMinuit->SetObjectFit(this);
- gMinuit->SetFCN(Fcn);
-
+ //
+ // CHOLM: Modified to use the TVirtualFitter interface only. This
+ // makes it possible other fitters than TMinuit (say TMinuit2,
+ // TFumili) by simply changing ones configuration file, or setting
+ // it up in a script. Much more robust and flexible.
Double_t arglist[10];
Int_t ierflg = 0;
+ TVirtualFitter *fitter = TVirtualFitter::Fitter(this,6); //initialize TMinuit
+ arglist[0] = -1;
+ ierflg = fitter->ExecuteCommand("SET PRINT", arglist, 1);
+
+ fitter->SetFCN(TrackResidualsChi2Fcn);
arglist[0] = 1;
- gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
+ ierflg = fitter->ExecuteCommand("SET ERR", arglist ,1);
// Set starting values and step sizes for parameters
- Double_t tr[3],rot[3];
- fAlignObj->GetPars(tr,rot);
- static Double_t step[6] = {0,0,0,0,0,0};
- gMinuit->mnparm(0, "x", tr[0], step[0], 0,0,ierflg);
- gMinuit->mnparm(1, "y", tr[1], step[1], 0,0,ierflg);
- gMinuit->mnparm(2, "z", tr[2], step[2], 0,0,ierflg);
- gMinuit->mnparm(3, "psi", rot[0], step[3], 0,0,ierflg);
- gMinuit->mnparm(4, "theta", rot[1], step[4], 0,0,ierflg);
- gMinuit->mnparm(5, "phi", rot[2], step[5], 0,0,ierflg);
+ Double_t pars[6] = {0,0,0,0,0,0};
+ for(Int_t i=0;i<6;i++) if(fBFixed[i]) pars[i]=fFixed[i];
+ Double_t step[6] = {0.0001,0.0001,0.0001,0.0001,0.0001,0.0001};
+ ierflg = fitter->SetParameter(0, "dx", pars[0], step[0], 0,0);
+ ierflg = fitter->SetParameter(1, "dy", pars[1], step[1], 0,0);
+ ierflg = fitter->SetParameter(2, "dz", pars[2], step[2], 0,0);
+ ierflg = fitter->SetParameter(3, "psi", pars[3], step[3], 0,0);
+ ierflg = fitter->SetParameter(4, "theta", pars[4], step[4], 0,0);
+ ierflg = fitter->SetParameter(5, "phi", pars[5], step[5], 0,0);
+
+ // Fix parameters
+ if(fBFixed[0]) {printf("Fixing dx=%f\n",pars[0]); fitter->FixParameter(0);}
+ if(fBFixed[1]) {printf("Fixing dy=%f\n",pars[1]); fitter->FixParameter(1);}
+ if(fBFixed[2]) {printf("Fixing dz=%f\n",pars[2]); fitter->FixParameter(2);}
+ if(fBFixed[3]) {printf("Fixing psi=%f\n",pars[3]); fitter->FixParameter(3);}
+ if(fBFixed[4]) {printf("Fixing theta=%f\n",pars[4]); fitter->FixParameter(4);}
+ if(fBFixed[5]) {printf("Fixing phi=%f\n",pars[5]); fitter->FixParameter(5);}
// Now ready for minimization step
arglist[0] = 500;
arglist[1] = 1.;
- gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
+ fitter->ExecuteCommand("MIGRAD", arglist ,2);
// Print results
Double_t amin,edm,errdef;
Int_t nvpar,nparx,icstat;
- gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
- gMinuit->mnprin(3,amin);
-
+ icstat = fitter->GetStats(amin,edm,errdef,nvpar,nparx);
+
+ //Construct the covariance matrix for AlignObj
+ Double_t *cov=fitter->GetCovarianceMatrix();
+ Int_t unfixedparam=6;
+ Int_t fixedparamat[6]={0,0,0,0,0,0};
+ if(fBFixed[0]==kTRUE){
+ unfixedparam--;
+ fixedparamat[0]=1;
+ }
+ for(Int_t j=1;j<6;j++){
+ if(fBFixed[j]==kTRUE){
+ unfixedparam--;
+ fixedparamat[j]=fixedparamat[j-1]+1;
+ }
+ else fixedparamat[j]=fixedparamat[j-1];
+ }
+
+ Double_t cov2[36];
+ for(Int_t i=0;i<6;i++){
+ for(Int_t j=0;j<6;j++){
+ if(fBFixed[i]==kTRUE||fBFixed[j]==kTRUE){
+ cov2[i+6*j]=0.;
+ }
+ else cov2[i+6*j]=cov[i-fixedparamat[i]+6*(j-fixedparamat[j])];
+
+ }
+ }
+
+ Double_t covmatrarray[21];
+ for(Int_t j=0;j<6;j++){
+ for(Int_t i=j;i<6;i++){
+ covmatrarray[i*(i+1)/2+j]=cov2[i+6*j];
+ }
+ }
+
+ // printf("covar 2: %.10f \n %.10f ; %.10f ; \n; %.10f; %.10f; %.10f; \n %.10f; %.10f; %.10f; %.10f; \n %.10f ; %.10f ; %.10f ; %.10f ; %.10f \n %.10f ; %.10f ; %.10f ; %.10f ; %.10f ; %.10f ;\n",covmatrarray[0],covmatrarray[1],covmatrarray[2],covmatrarray[3],covmatrarray[4],covmatrarray[5],covmatrarray[6],covmatrarray[7],covmatrarray[8],covmatrarray[9],covmatrarray[10],covmatrarray[11],covmatrarray[12],covmatrarray[13],covmatrarray[14],covmatrarray[15],covmatrarray[16],covmatrarray[17],covmatrarray[18],covmatrarray[19],covmatrarray[20]);
+
+
+ fAlignObj->SetCorrMatrix(covmatrarray);
+ fChi2 = amin; fNdf -= nvpar;
+
return kTRUE;
}
Double_t chi2 = 0;
fAlignObj->SetPars(par[0],par[1],par[2],par[3],par[4],par[5]);
- TGeoHMatrix m;
- fAlignObj->GetMatrix(m);
- Double_t *rot = m.GetRotationMatrix();
- Double_t *tr = m.GetTranslation();
AliTrackPoint p1,p2;
- Double_t dR[3];
- Float_t xyzvol[3],xyztr[3];
- Float_t covvol[6],covtr[6];
+
+ Bool_t count = kFALSE;
+ if (fNdf == 0) count = kTRUE;
for (Int_t itrack = 0; itrack < fLast; itrack++) {
if (!fVolArray[itrack] || !fTrackArray[itrack]) continue;
for (Int_t ipoint = 0; ipoint < fVolArray[itrack]->GetNPoints(); ipoint++) {
fVolArray[itrack]->GetPoint(p1,ipoint);
+ fAlignObj->Transform(p1);
fTrackArray[itrack]->GetPoint(p2,ipoint);
- p1.GetXYZ(xyzvol,covvol);
- p2.GetXYZ(xyztr,covtr);
-
- for (Int_t i = 0; i < 3; i++)
- dR[i] = tr[i]
- + xyzvol[0]*rot[3*i]
- + xyzvol[1]*rot[3*i+1]
- + xyzvol[2]*rot[3*i+2]
- - xyztr[i];
- chi2 += dR[0]*dR[0]+dR[1]*dR[1]+dR[2]*dR[2];
+ Float_t residual = p2.GetResidual(p1,kTRUE);
+ chi2 += residual;
+ if (count) fNdf += 3;
}
}
f = chi2;
}
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff