CVertexFitter.C
//-----------------------------------------------------------------------------
// $Header: CVertexFitter.C,v 3.3 97/04/21 13:04:45 messer Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class CVertexFitter.
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include <limits.h>
#include <math.h>
#include "CVertexFitter.h"
#include "CMinuitWrap.h"
#include "CMCDigiHit.h"
static const CSortedList<CSidcHit> *hits1; // hits SiDC1
static const CSortedList<CSidcHit> *hits2; // hits SiDC2
static RWTValOrderedVector<int> indexListLower; // lower Phi range SiDC1
static RWTValOrderedVector<int> indexListUpper; // upper Phi range SiDC2
static double zOfSidc1 = 0;
static double zOfTarget = 0;
static double vertexTwoSigma2 = 0;
void robustChi2(int*, double*, double*, double*, int*, void*); // this is for MINUIT to get the min. chi2
void CVertexFitter::fillVertex(CVertex& vtx)
{
vtx.setFittedZ(val[0]);
vtx.setZ(getZofNearestTargetDisk());
vtx.setX(val[1]);
vtx.setY(val[2]);
vtx.setChi2(fmin);
vtx.setNTracks(nFitPoints);
vtx.setRC(istat);
}
double CVertexFitter::getZofNearestTargetDisk()
{
//
// return z positions as expected from pair values
//
double z = val[0];
if (- 8.95 <= z && z < - 8.65) return (- 8.79);
else if (- 9.25 <= z && z < - 8.95) return (- 9.11);
else if (- 9.55 <= z && z < - 9.25) return (- 9.40);
else if (- 9.85 <= z && z < - 9.55) return (- 9.71);
else if (-10.15 <= z && z < - 9.85) return (-10.01);
else if (-10.45 <= z && z < -10.15) return (-10.32);
else if (-10.75 <= z && z < -10.45) return (-10.62);
else if (-11.05 <= z && z < -10.75) return (-10.92);
else return z; // bad z-position: return fitted value...
}
void CVertexFitter::listSetup(ostream& ost) const
{
ost << "\nSetup listing of vertex fitter" << endl;
CString line('=', 46);
ost << line << endl;
setup->list(ost);
}
CVertexFitter::CVertexFitter(const char* setupFile)
{
//
// Create and read setup
//
setup = new CVertexFitterSetup;
if (setupFile == 0) {
CString setupName = CString(C_DEFAULT_SETUP_PATH)+
CString(C_SETUPFILE_VERTEXFITTER);
setup->read(setupName.data());
}
else
setup->read(setupFile);
//
// Store setup parmeter to local variables
//
zOfSidc1 = setup->getZOfSidc1();
zOfTarget = setup->getZOfTarget();
vertexTwoSigma2 = setup->getVertexTwoSigma2();
vertexDeltaPhiRange = setup->getVertexDeltaPhiRange() * ToRadian;
vertexScanSteps = setup->getVertexScanSteps();
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;
strcpy(&names[0][0], "x ");
strcpy(&names[1][0], "y ");
strcpy(&names[2][0], "z ");
mninit(5, 6, 42); // initialize MINUIT
}
CVertexFitter::~CVertexFitter() { delete setup; }
CBoolean CVertexFitter::makeVertex(CSidc& sd1, CSidc& sd2)
{
hits1 = &(sd1.getHits()); // getHits() returns reference
hits2 = &(sd2.getHits());
if ( hits1->length() == 0 ) return False;
if ( hits2->length() == 0 ) return False;
indexListLower.clear(); // reset the index lists
indexListUpper.clear();
indexListLower.resize(hits2->length()); // and resize them
indexListUpper.resize (hits2->length());
nFitPoints = 0;
//
// Create list of indices for upper and lower phi range of SiDC2 hits
// A slice of +/- 'vertexDeltaPhiRange' around -Pi/Pi is left out in order
// to avoid pileup in the indices lists.
//
int i2, i1;
CAngle phi1, phi2;
int minIndex, maxIndex;
for (i2=0; i2<hits2->length(); i2++) {
minIndex = INT_MAX;
maxIndex = 0;
phi2 = (*hits2)(i2)->getCenter().getPhi();
if (fabs(double(phi2)) < Pi - vertexDeltaPhiRange) {
for (i1=0; i1<hits1->length(); i1++) {
phi1 = (*hits1)(i1)->getCenter().getPhi();
if (fabs(double(phi2 - phi1)) < vertexDeltaPhiRange ) {
if (i1 < minIndex) minIndex = i1;
if (i1 > maxIndex) maxIndex = i1;
}
}
}
indexListLower(i2) = minIndex;
indexListUpper(i2) = maxIndex;
if (maxIndex != 0) nFitPoints++;
}
//
// Prepare starting values for MINUIT fit
//
val[0] = findDisk();
val[1] = 0;
val[2] = 0;
step[0] = 0.05;
step[1] = 0.01;
step[2] = 0.01;
double argList;
argList = -1;
mnexcm(&robustChi2, "SET PRInt" , &argList, 1, &ierr, (void *) 0); // no printout
argList = 0;
mnexcm(&robustChi2, "SET NOWarn" , &argList, 1, &ierr, (void *) 0);
mnexcm(&robustChi2, "CLEar" , &argList, 1, &ierr, (void *) 0);
mnexcm(&robustChi2, "SET STRategy", &argList, 1, &ierr, (void *) 0);
//
// Feed starting values into MINUIT
//
for (int i=0; i<MaxFitParameters; i++) {
mnparm(i+1, names[i], val[i], step[i], 0, 0, &ierr);
if (ierr != 0) {
cerr << "CVertexFitter::makeVertex():" << endl;
cerr << "\tERROR" << endl;
cerr << "\tmnparm() returned " << ierr << " for i = " << i << endl;
return False;
}
}
//
// Perform vertex fit
//
mnexcm (&robustChi2, "MINImize", 0, 0, &ierr, (void *) 0);
if (ierr != 0) {
cerr << "CVertexFitter::makeVertex():" << 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 x, y, z
//
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);
return True;
}
double CVertexFitter::findDisk()
{
int i;
double dz, vtxZ;
double zstart = -3.0;
for (i=0; i<MaxFitParameters; i++) result[i] = 0;
fmin = FLT_MAX;
result[0] = zOfTarget + zstart;
dz = 2.*fabs(zstart)/vertexScanSteps;
double f;
for (i=0; i<vertexScanSteps; i++) {
robustChi2(&nParams, grad, &f, result, &flag, (void*) 0);
if (f < fmin) {
fmin = f;
vtxZ = result[0];
}
result[0] += dz;
}
return vtxZ;
}
void robustChi2(int*, double [3], double *f, double qs[3], int*, void*)
{
int ih2, ih1;
double scale, x, y, arg;
*f = 0;
for (ih2=0; ih2<hits2->length(); ih2++) { // loop sorted list of sidc-2 hits
scale = (zOfSidc1 - (*hits2)(ih2)->centerInXYZ.getZ()) /
((*hits2)(ih2)->centerInXYZ.getZ() - qs[0]);
x = (*hits2)(ih2)->centerInXYZ.getX() +
scale * ((*hits2)(ih2)->centerInXYZ.getX() - qs[1]);
y = (*hits2)(ih2)->centerInXYZ.getY() +
scale * ((*hits2)(ih2)->centerInXYZ.getY() - qs[2]);
for (ih1=indexListLower(ih2); ih1<=indexListUpper(ih2); ih1++) {
arg = ( ((*hits1)(ih1)->centerInXYZ.getX() - x)*
((*hits1)(ih1)->centerInXYZ.getX() - x) +
((*hits1)(ih1)->centerInXYZ.getY() - y)*
((*hits1)(ih1)->centerInXYZ.getY() - y) );
arg /= vertexTwoSigma2;
if (arg < 15) *f -= exp(-arg); // prevent underflows
}
}
}