CTrackFitter.C
//-----------------------------------------------------------------------------
// $Header: /asis/offline/ceres/cool/project/RCS/CTrackFitter.C,v 3.0 1996/10/02 09:39:57 voigt Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class CTrackFitter
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include <limits.h>
#include <math.h>
#include "rw/tpordvec.h"
#include "CTrackFitter.h"
#include "CMinuitWrap.h"
#include "CCoordinate.h"
#include "CElectronTrack.h"
#include "CTrackFitterSetup.h"
#include "CListIterator.h"
static RWTPtrOrderedVector<CLabXYZCoord> *measuredPoint; // info of all detectors
static RWTPtrOrderedVector<CLabXYZCoord> *measuredPointErrors;
static double zOfRich1MirrorHalf = 50; // cm
static double rich2RadiatorHalfLength = 100;
static double zOfDeflectionPlane = 110;
//
// this is the evaluation function for MINUIT to get the min. chi2
//
void trackChi2Fit(int*, double*, double*, double*, int*, void*);
CTrackFitter::CTrackFitter(const char* setupFile)
{
//
// Create and read setup
//
setup = new CTrackFitterSetup;
if (setupFile == 0) {
CString setupName = CString(C_DEFAULT_SETUP_PATH)+
CString(C_SETUPFILE_TRACKFITTER);
setup->read(setupName.data());
}
else
setup->read(setupFile);
//
// Store setup parmeter to local variables
//
nParams = MaxFitParameters;
fmin = 0;
flag = 0;
ierr = 0;
istat = 0;
for (int i=0; i<MaxFitParameters; i++) {
grad [i] = 0;
result[i] = 0;
val [i] = 0;
step [i] = 0;
}
nFitPoints = 0;
// 0123456789
strcpy(&names[0][0], "x0 ");
strcpy(&names[0][1], "y0 ");
strcpy(&names[0][2], "z0 ");
strcpy(&names[0][3], "thetaPre ");
strcpy(&names[0][4], "phiPre ");
strcpy(&names[0][5], "thetaPost ");
strcpy(&names[0][6], "phiPost ");
strcpy(&names[0][7], "deltaZ0 ");
mninit(5, 6, 42); // initialize MINUIT
}
CTrackFitter::~CTrackFitter()
{
delete setup;
}
CTrackFitter::CTrackFitter(const CTrackFitter &)
{
}
CTrackFitter & CTrackFitter::operator = (const CTrackFitter &)
{
return *this;
}
void CTrackFitter::listSetup(ostream& ost) const
{
ost << "\nSetup listing of track fitter" << endl;
CString line('=', 46);
ost << line << endl;
setup->list(ost);
}
void CTrackFitter::compareTracksAndFitResults(CList<CTrack>& tlist, ostream& ost)
{
CListIterator<CTrack> titer(tlist);
CTrack* track;
int i = 0;
ost << "\n Result of track fit" << endl;
CString line('=', 46);
ost << line << endl;
ost << "\n measured\t\t" << "fitted" << endl;
while (track = titer()) {
ost << "track no.\t" << i++ << endl;
ost << "starting point x:\t" << track->getVertex()->getX() << "\t" << val[0] << endl;
ost << "starting point y:\t" << track->getVertex()->getY() << "\t" << val[1] << endl;
ost << "starting point z:\t" << track->getVertex()->getZ() << "\t" << val[2] << endl;
ost << "theta pre field: \t" << track->getThetaSidc() << "\t" << val[3] << endl;
ost << "phi pre field: \t" << track->getPhiSidc() << "\t" << val[4] << endl;
ost << "theta post field:\t" << track->getClosestRich2Ring()->getPolarCenter().getTheta() << "\t" << val[5] << endl;
ost << "phi post field:\t" << track->getClosestRich2Ring()->getPolarCenter().getPhi() << "\t" << val[6] << endl;
ost << "delta z :\t" << 0 << "\t" << val[7] << endl;
}
}
CBoolean CTrackFitter::fitTracks(CList<CTrack>& tlist)
{
CListIterator<CTrack> titer(tlist);
CTrack *track;
while (track = titer()) fitOneTrack(*track);
return True;
}
CBoolean CTrackFitter::fitOneTrack(CTrack& track)
{
nFitPoints = 0;
measuredPoint->clearAndDestroy();
//
// feed track info into the static toplevel container
//
CLabXYZCoord point, deflectionPoint;
CPadCPadCoord padcPoint;
point.setX( track.getVertex()->getX() );
point.setY( track.getVertex()->getY() );
point.setZ( track.getVertex()->getZ() );
measuredPoint->append(&point);
point = track.getClosestSidc1Hit()->getCenterInXYZ();
measuredPoint->append(&point);
point = track.getClosestSidc1Hit()->getCenterInXYZ();
measuredPoint->append(&point);
//
// translate the measured angles in the rich detectors into pseudo cartesian coordinates
//
point.setX( cos(track.getClosestRich1Ring()->getPolarCenter().getPhi())
*tan(track.getClosestRich1Ring()->getPolarCenter().getTheta())
*(zOfRich1MirrorHalf + track.getVertex()->getZ()) );
point.setY( sin(track.getClosestRich1Ring()->getPolarCenter().getPhi())
*tan(track.getClosestRich1Ring()->getPolarCenter().getTheta())
*(zOfRich1MirrorHalf + track.getVertex()->getZ()) );
point.setZ( zOfRich1MirrorHalf + track.getVertex()->getZ() );
measuredPoint->append(&point);
deflectionPoint.setX( cos(track.getClosestRich1Ring()->getPolarCenter().getPhi())
*tan(track.getClosestRich1Ring()->getPolarCenter().getTheta())
*(zOfDeflectionPlane + track.getVertex()->getZ()) );
deflectionPoint.setY( sin(track.getClosestRich1Ring()->getPolarCenter().getPhi())
*tan(track.getClosestRich1Ring()->getPolarCenter().getTheta())
*(zOfDeflectionPlane + track.getVertex()->getZ()) );
deflectionPoint.setZ( zOfDeflectionPlane );
point.setX( cos(track.getClosestRich2Ring()->getPolarCenter().getPhi())
*tan(track.getClosestRich2Ring()->getPolarCenter().getTheta())
*rich2RadiatorHalfLength );
point.setY( sin(track.getClosestRich2Ring()->getPolarCenter().getPhi())
*tan(track.getClosestRich2Ring()->getPolarCenter().getTheta())
*rich2RadiatorHalfLength );
point.setZ( rich2RadiatorHalfLength );
point+=deflectionPoint;
measuredPoint->append(&point);
padcPoint.setX( track.getClosestPadCHit()->getX() );
padcPoint.setY( track.getClosestPadCHit()->getY() );
//transform(point,padcPoint);
measuredPoint->append(&point);
//
// fill in all errors >>>>>>>not finished
//
measuredPointErrors->clearAndDestroy();
point.setX(0.1);
point.setY(0.1);
point.setZ(0.1);
measuredPointErrors->append(&point);
measuredPointErrors->append(&point);
measuredPointErrors->append(&point);
measuredPointErrors->append(&point);
measuredPointErrors->append(&point);
measuredPointErrors->append(&point);
//
// Prepare starting values for MINUIT fit
//
val[0] = track.getVertex()->getX();
val[1] = track.getVertex()->getY();
val[2] = track.getVertex()->getZ();
val[3] = track.getClosestRich1Ring()->getPolarCenter().getTheta();
val[4] = track.getClosestRich1Ring()->getPolarCenter().getPhi();
val[5] = track.getClosestRich2Ring()->getPolarCenter().getTheta();
val[6] = track.getClosestRich2Ring()->getPolarCenter().getPhi();
val[7] = 0;
//
// Prepare starting step sizes for MINUIT fit
//
step[0] = 0.01;
step[1] = 0.01;
step[2] = 0.01;
step[3] = 0.01;
step[4] = 0.01;
step[5] = 0.01;
step[6] = 0.01;
step[7] = 0.01;
//
// Prepare parameter constraints for MINUIT fit
//
int i;
for (i=0; i<MaxFitParameters; i++) {
lowerBounds[i] = 0;
upperBounds[i] = 0;
}
lowerBounds[7] = -2; // constraints in the uncertainty of deflection plane
upperBounds[7] = 2;
//
// Prepare global steering for MINUIT fit
//
double argList;
argList = -1;
mnexcm(&trackChi2Fit, "SET PRInt" , &argList, 1, &ierr, (void *) 0); // no printout
argList = 0;
mnexcm(&trackChi2Fit, "SET NOWarn" , &argList, 1, &ierr, (void *) 0);
mnexcm(&trackChi2Fit, "CLEar" , &argList, 1, &ierr, (void *) 0);
mnexcm(&trackChi2Fit, "SET STRategy", &argList, 1, &ierr, (void *) 0);
//
// Feed starting values into MINUIT
//
for ( i=0; i<MaxFitParameters; i++) {
mnparm(i+1, names[i], val[i], step[i], lowerBounds[i], upperBounds[i], &ierr);
if (ierr != 0) {
cerr << "CTrackFitter::fitTrack():" << endl;
cerr << "\tERROR" << endl;
cerr << "\tmnparm() returned " << ierr << " for i = " << i << endl;
return False;
}
}
//
// Perform track fit
//
mnexcm (&trackChi2Fit, "MINImize", 0, 0, &ierr, (void *) 0);
if (ierr != 0) {
cerr << "CTrackFitter::fitTrack():" << endl;
cerr << "\tERROR" << endl;
cerr << "\tmnexcm() returned error code " << ierr << endl;
return False;
}
//
// Check fit
//
int idummy = 0;
double errdef = 0, fedm = 0;
mnstat (&fmin, &fedm, &errdef, &idummy, &idummy, &istat);
//
// Store fitted values
//
double dummy = 0;
mnpout(1, " ", &val[0], errdef, dummy, dummy, idummy);
mnpout(2, " ", &val[1], errdef, dummy, dummy, idummy);
mnpout(3, " ", &val[2], errdef, dummy, dummy, idummy);
mnpout(4, " ", &val[3], errdef, dummy, dummy, idummy);
mnpout(5, " ", &val[4], errdef, dummy, dummy, idummy);
mnpout(6, " ", &val[5], errdef, dummy, dummy, idummy);
mnpout(7, " ", &val[6], errdef, dummy, dummy, idummy);
mnpout(8, " ", &val[7], errdef, dummy, dummy, idummy);
return True;
}
void trackChi2Fit(int*, double [8], double *f, double qs[8], int*, void*)
{
*f=0;
CLabXYZCoord predictor;
for (int i=0;i<5;i++) {
if ( (*measuredPoint)(i)->getZ() < zOfDeflectionPlane ) {
predictor.setX( qs[0] + cos(qs[4])*tan(qs[3])*((*measuredPoint)(i)->getZ()+fabs(qs[2])) );
predictor.setY( qs[1] + sin(qs[4])*tan(qs[3])*((*measuredPoint)(i)->getZ()+fabs(qs[2])) );
} else {
predictor.setX( qs[0] + (cos(qs[4])*tan(qs[3])-cos(qs[6])*tan(qs[5]))
*(zOfDeflectionPlane+fabs(qs[2])+qs[7])
+ (cos(qs[6])*tan(qs[5])*((*measuredPoint)(i)->getZ()+fabs(qs[2])-qs[7])) );
predictor.setY( qs[1] + (sin(qs[4])*tan(qs[3])-sin(qs[6])*tan(qs[5]))
*(zOfDeflectionPlane+fabs(qs[2])+qs[7])
+ (sin(qs[6])*tan(qs[5])*((*measuredPoint)(i)->getZ()+fabs(qs[2])-qs[7])) );
}
predictor.setZ( (*measuredPoint)(i)->getZ() );
*f += (*(*measuredPoint)(i)-predictor) * (*(*measuredPoint)(i)-predictor)
/(*(*measuredPointErrors)(i) * *(*measuredPointErrors)(i));
}
}