CFittedRing.C
//-----------------------------------------------------------------------------
// $Header: /asis/offline/ceres/cool/project/RCS/CFittedRing.C,v 2.4 1997/05/30 12:14:02 messer Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of CFittedRing class.
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include <iomanip.h>
#include <math.h>
#include "CFittedRing.h"
#include "CRichLikeHit.h"
#include "rw/tvsrtvec.h"
CFittedRing::CFittedRing() : CRing()
{
numberOfHits = 0;
variance = 0;
activeCircumference = 0;
kolmogorovTestValue = 0;
chi2 = 0;
numberOfPads = 0;
analogRingSum = 0;
ringCandidate = 0;
}
CFittedRing::~CFittedRing()
{
hits.clear(); // not clearAndDestroy
}
void CFittedRing::setCenter(CRichPadCoord rpad, CRichPolarCoord rpol)
{
polarCenter = rpol;
center = rpad;
}
void CFittedRing::clearHitInformation()
{
hits.clear();
numberOfHits = 0;
}
void CFittedRing::addHit(CRichLikeHit* hit)
{
hits.insert(hit);
numberOfHits++;
}
void CFittedRing::printProperties(ostream& ost) const
{
const ListingNameWidth = 24;
ost.setf(ios::left);
ost << setw(ListingNameWidth) << "ring type:";
switch (type) {
case Electron:
ost << "Electron" << endl;
break;
case Pion:
ost << "Pion" << endl;
break;
case Fake:
ost << "Fake" << endl;
break;
default:
ost << "Unknown" << endl;
break;
}
ost << setw(ListingNameWidth) << "center:" << dec << center << endl;
ost << setw(ListingNameWidth) << " " << dec << polarCenter << endl;
ost << setw(ListingNameWidth) << "radius:" << dec << radius << " pads" << endl;
ost << setw(ListingNameWidth) << "chi2:" << dec << chi2 << endl;
ost << setw(ListingNameWidth) << "fit variance:" << dec << variance << endl;
ost << setw(ListingNameWidth) << "number of hits:" << dec << numberOfHits << endl;
ost << setw(ListingNameWidth) << "active circumference:" << dec << activeCircumference << endl;
ost << setw(ListingNameWidth) << "kolmogorov test value:" << dec << kolmogorovTestValue << endl;
ost << setw(ListingNameWidth) << "number of pads:" << dec << numberOfPads << endl;
ost << setw(ListingNameWidth) << "analog ring sum:" << dec << analogRingSum << endl;
ost.unsetf(ios::left);
}
double CFittedRing::evalActiveSegment(const CCollection<CRichLikeLookupItem>& lookupList,
double minPhi, double maxPhi)
{
double dphi = 0.5/radius;
int totalPads = 0;
int activePads = 0;
int ix, iy;
int ixOld = 0;
int iyOld = 0;
CRichLikeLookupItem *item;
for (double phi = minPhi; phi < maxPhi; phi += dphi) {
ix = int(center.getX()+radius*cos(phi));
iy = int(center.getY()+radius*sin(phi));
if (ix == ixOld && iy == iyOld) continue; // already checked
totalPads++;
item = lookupList(ix, iy);
if (item != 0 && item->getIsSensitive()) activePads++;
ixOld = ix;
iyOld = iy;
}
return (double) activePads/totalPads;
}
void CFittedRing::evalActiveCircumference(const CCollection<CRichLikeLookupItem>& list)
{
activeCircumference = evalActiveSegment(list, 0, TwoPi);
}
void CFittedRing::evalKolmogorovTestValue(const CCollection<CRichLikeLookupItem>& lookupList)
{
//
// Whatever is written below, this function gives segmentation violations.
// The bug is probably some off-by-one-problem in the function
// CProbKolmogorov. Until somebody had time and nerves to look into
// this, the function will remain dummy.
//
cerr << "CFittedRing::evalKolmogorovTestValue()" << endl;
cerr << "\tWARNING" << endl;
cerr << "\tThis function is at present dummy." << endl;
//
//
// This function returns the probability that the distribution
// of the given fit points along the circumference of a ring
// is compatible with a uniform distribution. The return value
// is in the range 0-1 and gives the confidence level for the
// null hypothesis. Values of 'cl' close to zero are taken as
// indicating a small probability of compatibility.
//
// This function uses the CERNLIB routine PROBKL (G102) to
// calculate the CL. In PROBKL the Kolmogorov distribution
// function is calculated which is only valid for large number
// of points. However, a comparison with the former used values
// taken from a table for N<10 shows no significant difference.
//
//if (hits.length() < 2) return; // nothing to check
//
// Get phi of each point and sort them
//
//double nextPhi;
//RWTValSortedVector<double> phi(hits.length());
//for (int i=0; i<hits.length(); i++) {
// nextPhi = atan2(hits(i)->getY()-center.getY(), hits(i)->getX()-center.getX());
// if (nextPhi < 0) nextPhi += TwoPi;
// phi.insert(nextPhi);
//}
//
// We need the active circumference, if not already calculated go and get it.
//
//if (activeCircumference <= 0)
// evalActiveCircumference(lookupList);
//
// Remove inactive area
//
//double lastPhi = 0.;
//double gap = 0.;
//for (i=0; i<phi.length(); i++) {
// gap += (phi(i)-lastPhi)*(1-evalActiveSegment(lookupList, lastPhi, phi(i)));
// lastPhi = phi(i);
// phi(i) -= gap;
//}
//
// Determine Dn
//
//double dn = 0.;
//double d1, d2, d;
//double activeRadians = TwoPi*activeCircumference;
//for (i=0; i<phi.length(); i++) {
// d1 = fabs(phi(i)/activeRadians - double(i)/double(phi.length()));
// d2 = fabs(phi(i)/activeRadians - double(i+1)/double(phi.length()));
// d = (d1 > d2 ? d1 : d2);
// dn = dn > d ? dn : d;
//}
//
// Calculate confidence level (CERNLIB PROBKL G102)
//
//kolmogorovTestValue = CProbKolmogorov(dn*sqrt(phi.length()));
}
void CFittedRing::evalChi2(double r)
{
//
// Get chi2 according to given radius.
// Not divided by degrees-of-freedom.
//
chi2 = 0;
double dist;
for (int i=0; i<hits.entries(); i++) {
dist = sqrt((center.getX()-hits(i)->getX())*(center.getX()-hits(i)->getX()) +
(center.getY()-hits(i)->getY())*(center.getY()-hits(i)->getY()));
chi2 += (dist-r)*(dist-r);
}
}
void CFittedRing::evalRingSum(const CCollection<CRichLikeLookupItem>& list,
const CCollection<CPad>& pads,
double halfMaskWidth)
{
numberOfPads = 0;
analogRingSum = 0;
double maxRadiusSquared = (radius+halfMaskWidth)*(radius+halfMaskWidth);
double minRadiusSquared = (radius-halfMaskWidth)*(radius-halfMaskWidth);
int ix, iy;
double distSquared;
CPad *pad;
double range = radius+halfMaskWidth+1;
int ixmin = int(center.getX()-range);
int ixmax = int(center.getX()+range);
int iymin = int(center.getY()-range);
int iymax = int(center.getY()+range);
for (ix = ixmin; ix <= ixmax; ix++) {
for (iy = iymin; iy <= iymax; iy++) {
distSquared = (ix-center.getX())*(ix-center.getX()) +
(iy-center.getY())*(iy-center.getY());
if (distSquared > minRadiusSquared && distSquared < maxRadiusSquared) {
pad = pads(ix,iy);
if (pad && pad->getAmp() > 0) {
numberOfPads++;
analogRingSum += pad->getAmp() * list(ix,iy)->getGainCalibrationFactor();
}
}
}
}
}
void CFittedRing::getMCGeantTracks(CRich& rich,RWTPtrOrderedVector<CMCGeantTrack>& theTracks)
{
int i, idigi;
CMCGeantTrack *theTrack;
RWTPtrOrderedVector<CMCGeantTrack> tracks;
RWTPtrOrderedVector<CMCDigiHit> theDigiHits;
tracks.clear();
for (i=0; i<getHits().length(); i++ ) {
theDigiHits.clear();
getHits()(i)->getMCDigiHits(rich, theDigiHits);
for (idigi=0; idigi<theDigiHits.length(); idigi++) {
theTrack = theDigiHits(idigi)->getGeantTrack();
tracks.insert( theTrack );
}
}
if (tracks.length() == 0) return; // nothing found... probably an overlay with data
if (tracks.length() > 1) { // more than one track for this ring, find the one with most hits...
int max=0, maxIndex=-1;
for (i=0; i<tracks.length(); i++) {
if (tracks.occurrencesOf( tracks(i) ) > max) {
maxIndex = i;
}
}
theTracks.insert(tracks(maxIndex));
}
for (i=0; i<tracks.length(); i++) { // copy remaining tracks only once into the output list
if ( !theTracks.contains( tracks(i) ) ) {
theTracks.insert(tracks(i));
}
}
tracks.clear();
theDigiHits.clear();
}