CRobustVertexFitter.C
//-----------------------------------------------------------------------------
// $Header: /asis/offline/ceres/cool/project/RCS/CRobustVertexFitter.C,v 1.9 1997/08/19 10:16:57 messer Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class CRobustVertexFitter.
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include <limits.h>
#include <math.h>
#include "CRobustVertexFitter.h"
#include "CCyclicIterator.h"
#include "CElectronTrack.h"
#include "CSortedListIterator.h"
void CRobustVertexFitter::fillVertex(CVertex& vtx)
{
vtx.setFittedZ(center.getZ());
if (setup->getFixTargetDiscs())
vtx.setZ(getZOfNearestTargetDisk());
else
vtx.setZ(center.getZ());
vtx.setX(center.getX());
vtx.setY(center.getY());
vtx.setChi2(variance);
vtx.setNTracks(nFitPoints);
vtx.setRC(0);
}
double CRobustVertexFitter::getZOfNearestTargetDisk()
{
double z = center.getZ();
//
// Get the z position as expected from the target positions given by the setup values.
// The fixed target position for one disc is returned, if the fitted z value is within
// half the distance to the previous or next disc.
// The first and last disc in the list have to be treated individually. There we allow
// a distance of 1.6 mm from the predicted position.
// If no match is found, the fitted z position is returned.
//
int numberOfDiscs = setup->getNumberOfTargetDiscs();
double zfix, znext, zmin, zmax, zprevious;
for (int i=0; i<numberOfDiscs; i++) {
zfix = setup->getZOfTargetDiscs(i);
znext = i<numberOfDiscs-1 ? setup->getZOfTargetDiscs(i+1) : setup->getZOfTargetDiscs(i)-0.32;
zprevious = i>0 ? setup->getZOfTargetDiscs(i-1) : setup->getZOfTargetDiscs(i)+0.32;
zmin = zfix - (zfix - znext)/2.;
zmax = zfix + (zprevious - zfix)/2.;
//
// As we do not know in which order the discs were given,
// we simply check if z lies inbetween the borders.
//
if (zmin <= z && z < zmax) return (zfix);
}
return z; // bad z-position: return fitted value...
}
void CRobustVertexFitter::listSetup(ostream& ost) const
{
ost << "\nSetup listing of vertex fitter" << endl;
CString line('=', 46);
ost << line << endl;
setup->list(ost);
}
CRobustVertexFitter::CRobustVertexFitter(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 parameter to local variables
//
vertexDeltaPhiRange = CAngle(setup->getVertexDeltaPhiRange() * ToRadian);
nFitPoints = 0;
//
// initialize vertex
//
center.setX(0);
center.setY(0);
center.setZ(0);
variance = 0;
}
CRobustVertexFitter::~CRobustVertexFitter() { delete setup; }
CBoolean CRobustVertexFitter::makeVertex(CSidc& sd1, CSidc& sd2)
{
//
// Fill list of hits in both sidcs which are within the required radial range.
//
sidc1Hits.clear();
sidc2Hits.clear();
int i;
for (i=0; i<sd1.getHits().length(); i++) {
if (sd1.getHits()(i)->getCenter().getRadius() > setup->getRMinSidc1() &&
sd1.getHits()(i)->getCenter().getRadius() < setup->getRMaxSidc1()) {
sidc1Hits.append(sd1.getHits()(i));
}
}
for (i=0; i<sd2.getHits().length(); i++) {
if (sd2.getHits()(i)->getCenter().getRadius() > setup->getRMinSidc2() &&
sd2.getHits()(i)->getCenter().getRadius() < setup->getRMaxSidc2()) {
sidc2Hits.append(sd2.getHits()(i));
}
}
//
// still some hits left?
//
if ( sidc1Hits.length() == 0 ) return True;
if ( sidc2Hits.length() == 0 ) return True;
indexListLower.clear(); // reset the index lists
indexListUpper.clear();
indexListLower.resize(sidc2Hits.length()); // and resize them
indexListUpper.resize(sidc2Hits.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.
//
CAngle phi2;
CSidcHit testHit;
CLabCylinderCoord testCoord;
for (int i2=0; i2<sidc2Hits.length(); i2++) {
indexListLower(i2) = 0;
indexListUpper(i2) = 0;
phi2 = sidc2Hits(i2)->getCenter().getPhi();
if (fabs(double(phi2)) < Pi - double(vertexDeltaPhiRange)) {
testCoord.setPhi(phi2-vertexDeltaPhiRange);
testHit.setCenter(testCoord);
indexListLower(i2) = sidc1Hits.getIndexOfBestMatch(testHit);
testCoord.setPhi(phi2+vertexDeltaPhiRange);
testHit.setCenter(testCoord);
indexListUpper(i2) = sidc1Hits.getIndexOfBestMatch(testHit);
}
if (indexListUpper(i2) != 0) nFitPoints++;
}
//
// perform vertex fit
//
double A, B, C,D, E, F, G, H, K, L;
double ARG, w, varsqold;
double xold, yold, zold, dx, dy, dz, arg1;
const double c1 = 3;
int ihit2;
CLabXYZCoord hit1, hit2;
const double finalPrecision = 0.005; // in cm
CCyclicIterator<CSidcHit, CSortedList<CSidcHit> > sidc1HitIterator(sidc1Hits);
//
// set initial values
//
double x = 0.0;
double y = 0.0;
double z = setup->getZOfTargetCenter();
double varsq = 0.03;
for (int n=0; n<20; n++) {
A = B = C = D = E = F = G = L = 0.;
for (ihit2=0; ihit2<sidc2Hits.length(); ihit2++) {
if (indexListUpper(ihit2) == 0) continue;
hit2 = sidc2Hits(ihit2)->getCenterInXYZ();
sidc1HitIterator.setPosition(indexListLower(ihit2));
while(sidc1HitIterator.getPosition()-1 != indexListUpper(ihit2)) {
hit1 = sidc1HitIterator.next()->getCenterInXYZ();
H = (hit1.getX()-hit2.getX()) / (hit2.getZ()-hit1.getZ());
K = (hit1.getY()-hit2.getY()) / (hit2.getZ()-hit1.getZ());
ARG = (x-hit1.getX()+H*(z-hit1.getZ())) * (x-hit1.getX()+H*(z-hit1.getZ()))
+ (y-hit1.getY()+K*(z-hit1.getZ())) * (y-hit1.getY()+K*(z-hit1.getZ()));
arg1 = sqrt(ARG);
if(ARG <= c1*c1*varsq)
w = (1 - ARG/(c1*c1*varsq)) * (1 - ARG/(c1*c1*varsq));
else
w=0.0;
A += w;
B += w * H;
C += w * hit1.getX();
D += w * K;
E += w * hit1.getY();
G += w * (hit1.getX()*(hit2.getX()-hit1.getX()) +
hit1.getY()*(hit2.getY()-hit1.getY()))/
(hit2.getZ()-hit1.getZ());
F += w * ((hit2.getX()-hit1.getX()) * (hit2.getX()-hit1.getX()) +
(hit2.getY()-hit1.getY()) * (hit2.getY()-hit1.getY())) /
((hit2.getZ()-hit1.getZ()) * (hit2.getZ()-hit1.getZ()));
L += w * ARG;
}
}
xold = x;
yold = y;
zold = z;
varsqold = varsq;
if((F*A-D*D-B*B) == 0.0 || A == 0.0) {
clog << "CRobustVertexFitter::makeVertex()\n";
clog << "\tWARNING\n";
clog << "\tpossible floating exeption\n";
return False;
}
double z0sd1 = sd1.getSetup()->getZPosition();
z = (z0sd1*F*A-G*A-D*(E+z0sd1*D)-B*(C+z0sd1*B))/(F*A-D*D-B*B);
x = (C+(z0sd1-z)*B)/A;
y = (E+(z0sd1-z)*D)/A;
varsq = L/A;
dx = x - xold;
dy = y - yold;
dz = z - zold;
//
// leave iteration loop if we are good enough
//
double maxdxy = fabs(dx) > fabs(dy) ? dx : dy;
if(fabs(maxdxy) < finalPrecision && abs(dz) < finalPrecision) break;
}
//
// memorize final values
//
center.setX(x);
center.setY(y);
center.setZ(z);
variance = varsq;
return True;
}
void CRobustVertexFitter::refineVertex(CVertex& vertex, CSortedList<CElectronTrack>& trackList)
{
CSortedListIterator<CElectronTrack> trackIter(trackList);
CElectronTrack *track;
double x0 = 0;
double y0 = 0;
double z0 = vertex.getZ();
double x1, y1, z1, x2, y2, z2, zFactor;
double x02 = 0;
double y02 = 0;
int trackCounter = 0;
while (track = trackIter()) {
if (track->getClosestSidc1Hit() && track->getClosestSidc2Hit() &&
track->getClosestSidc1Hit()->getAmp() < setup->getMaxDeDxSidc1() &&
track->getClosestSidc2Hit()->getAmp() < setup->getMaxDeDxSidc2() &&
track->getClosestSidc1Hit()->getNPulses() >= setup->getMinPulsesSidc1() &&
track->getClosestSidc2Hit()->getNPulses() >= setup->getMinPulsesSidc2() &&
track->getDeltaRSidc12() < setup->getMaxDeltaRSidc12() &&
track->getDeltaPhiSidc12() < setup->getMaxDeltaPhiSidc12() ) {
trackCounter++;
x1 = track->getClosestSidc1Hit()->getCenterInXYZ().getX();
y1 = track->getClosestSidc1Hit()->getCenterInXYZ().getY();
z1 = track->getClosestSidc1Hit()->getCenterInXYZ().getZ();
x2 = track->getClosestSidc2Hit()->getCenterInXYZ().getX();
y2 = track->getClosestSidc2Hit()->getCenterInXYZ().getY();
z2 = track->getClosestSidc2Hit()->getCenterInXYZ().getZ();
zFactor = (z2-z0)/(z2-z1);
x0 += x2-(x2-x1)*zFactor;
y0 += y2-(y2-y1)*zFactor;
x02 += (x2-(x2-x1)*zFactor)*(x2-(x2-x1)*zFactor);
y02 += (y2-(y2-y1)*zFactor)*(y2-(y2-y1)*zFactor);
}
}
double meanx0 = x0 / trackCounter;
double meany0 = y0 / trackCounter;
double meanx02 = x02 / trackCounter;
double meany02 = y02 / trackCounter;
vertex.setX(meanx0);
vertex.setY(meany0);
vertex.setVarX(meanx02 - meanx0*meanx0);
vertex.setVarY(meany02 - meany0*meany0);
vertex.setNTracks(trackCounter);
}