CElectronTrackingStrategy.C
//-----------------------------------------------------------------------------
// $Header: /cool/project/RCS/CElectronTrackingStrategy.C,v 2.3 1997/02/24 11:38:19 lenkeit Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class CElectronTrackingStrategy.
//
//-----------------------------------------------------------------------------
#include "rw/tvvector.h"
#include "CElectronTrackingStrategy.h"
#include "CListIterator.h"
CElectronTrackingStrategy::CElectronTrackingStrategy(CVertex* vtx, CSidc* s1,
CSidc* s2, CRich1* r1,
CRich2 *r2, const char* file)
: CSidcTrackingStrategy(vtx, s1, s2, file), rich1(r1), rich2(r2)
{ /* nop */ }
CElectronTrackingStrategy::~CElectronTrackingStrategy()
{
richTracks.clearAndDestroy();
tracks.clearAndDestroy();
}
CBoolean CElectronTrackingStrategy::makeTracks()
{
richTracks.clearAndDestroy();
tracks.clearAndDestroy();
if (!makeSidcTracks()) return False;
if (!makeRichTracks()) return False;
if (!linkRichAndSidcTracks()) return False;
return True;
}
CBoolean CElectronTrackingStrategy::getTracks(CList<CElectronTrack>& list,
unsigned long mask)
{
int i;
CElectronTrack *track;
//
// First see in which list we have to look for the stuff
//
CBoolean needsSidcSegments = mask&CTrack::hasVertex||mask&CTrack::hasSidc1Match||
mask&CTrack::hasSidc2Match;
CBoolean needsRichSegments = mask&CTrack::hasRich1Match||mask&CTrack::hasRich2Match;
CBoolean needsTracks = needsSidcSegments && needsRichSegments;
CBoolean needsThemAll = !needsSidcSegments && !needsRichSegments;
//
// Scan the internal list and create copies of the matching entries
//
if (needsTracks || needsThemAll) {
//
// Look in track list only
//
for (i=0; i<tracks.length(); i++) {
if ((tracks(i)->getTrackMask()&mask) == mask)
list.append(new CElectronTrack(*tracks(i)));
}
}
if (needsSidcSegments || needsThemAll) {
//
// Look in list of sidc segments only.
// Merge SiDC track into electron track
//
for (i=0; i<sidcTracks->length(); i++) {
if (((*sidcTracks)(i)->getTrackMask()&mask) == mask) {
track = new CElectronTrack;
*((CTrack*)track) = *((*sidcTracks)(i));
list.append(track);
}
}
}
if (needsRichSegments || needsThemAll) {
//
// Look in list of rich segments only
//
for (i=0; i<richTracks.length(); i++) {
if ((richTracks(i)->getTrackMask()&mask) == mask)
list.append(new CElectronTrack(*richTracks(i)));
}
}
return True;
}
CBoolean CElectronTrackingStrategy::makeRichTracks()
{
enum { notUsed = 0,
usedForNegativeDeflectedTrack,
usedForPositiveDeflectedTrack,
used };
//
// Check if there's anything to track
//
if (!rich1 && !rich2) return True;
//
// Create some internal lists, initialized to enum 'notUsed'.
// If no rings exist it's just an empty list, i.e.,
// usedRich1Rings.length() == 0
//
size_t length = (rich1 && rich1->getRingList()) ? rich1->getRingList()->length() : 0;
RWTValVector<int> usedRich1Rings(length, notUsed);
length = (rich2 && rich2->getRingList()) ? rich2->getRingList()->length() : 0;
RWTValVector<int> usedRich2Rings(length, notUsed);
//
// Create ring lists iterators
//
CListIterator<CFittedRing> *rich1Iter;
CListIterator<CFittedRing> *rich2Iter;
if (usedRich1Rings.length())
rich1Iter = new CListIterator<CFittedRing>(*(rich1->getRingList()));
if (usedRich2Rings.length())
rich2Iter = new CListIterator<CFittedRing>(*(rich2->getRingList()));
//
// First match Rich1 & Rich2 rings
// The phi of the ring centers are treated as CAngles, which
// saves a lot of checks.
//
CAngle phi1, phi2;
double theta1, theta2;
CFittedRing *rich1Ring, *rich2Ring;
int indexRich1, indexRich2;
double deflection; // phi deflection in rad
double p; // momentum of track
double dtheta;
double maxDeltaTheta;
double maxDeflection;
double smallestDeltaTheta;
double dphiOfBest;
double pOfBest;
int indexOfBestInRich1;
int indexOfBestInRich2;
int signOfTrack;
CElectronTrack *segment;
//
// Iterate to find successive the best matching rings
//
while (1) {
smallestDeltaTheta = FLT_MAX; // init to find best match
indexOfBestInRich1 = -1; // -1 == no further match
indexOfBestInRich2 = -1;
//
// Loop Rich2 rings
//
for (indexRich2 = 0; indexRich2 < usedRich2Rings.length(); indexRich2++) {
if (usedRich2Rings(indexRich2) == used) continue; // already used, skip
//
// Get center of Rich2 rings
//
rich2Ring = (*(rich2->getRingList()))(indexRich2);
phi2 = rich2Ring->getPolarCenter().getPhi();
theta2 = rich2Ring->getPolarCenter().getTheta();
//
// Loop Rich1 rings
//
for (indexRich1 = 0; indexRich1 < usedRich1Rings.length(); indexRich1++) {
//
// Get and transform ring center of Rich2 rings
//
rich1Ring = (*(rich1->getRingList()))(indexRich1); // pointer to ring
phi1 = rich1Ring->getPolarCenter().getPhi();
theta1 = rich1Ring->getPolarCenter().getTheta();
//
// Calculate deflection.
// For nominal field direction: electrons have deflection < 0
//
deflection = (double)(phi2-phi1);
//
// Check if already used with same sign.
// Note that a Rich1 ring can be used only twice (conversion) if
// the second Rich2 fit has opposite sign than the previous stored
// one.
//
if (usedRich1Rings(indexRich1) == used) continue;
if (usedRich1Rings(indexRich1) == usedForNegativeDeflectedTrack && deflection < 0) continue;
if (usedRich1Rings(indexRich1) == usedForPositiveDeflectedTrack && deflection > 0) continue;
//
// Get momentum of track.
// For no-field data the lowest momentum according to the
// given pt-cut is assumed (later needed for butterfly).
//
if (setup->getFieldFlag())
p = phiDeflectionAt1GeV(theta1)/fabs(deflection);
else
p = setup->getPtCut()/sin(theta1/1000.);
//
// Evaluate theta deflection corrected for 2-order field effect
//
dtheta = theta2-theta1;
if (setup->getFieldFlag()) dtheta += secondOrderFieldEffect(theta1, deflection);
//
// Check if Rich1 ring is within butterfly:
// Note that for no field the range is determined by the spread
// given by the lowest momentum according to the pt-cut.
//
maxDeltaTheta = setup->getButterflySigmas() * richButterfly(p);
if (setup->getFieldFlag())
maxDeflection = phiDeflectionAt1GeV(theta1)/(setup->getPtCut()/sin(theta1/1000.));
else
maxDeflection = maxDeltaTheta/sin(theta2/1000.);
if (fabs(deflection) > maxDeflection || fabs(dtheta) > maxDeltaTheta) continue;
//
// Store indices of best matching rings
//
if (fabs(dtheta) < fabs(smallestDeltaTheta)) {
smallestDeltaTheta = dtheta;
indexOfBestInRich1 = indexRich1;
indexOfBestInRich2 = indexRich2;
signOfTrack = deflection < 0 ? -1 : 1;
dphiOfBest = deflection;
pOfBest = p;
}
}
}
//
// Store best track so far, or stop if no more matchings
// rings were found. We do not perform any calculations
// here but fill the track list with all the info we have.
//
if (indexOfBestInRich1 == -1 && indexOfBestInRich2 == -1) break;
segment = new CElectronTrack;
segment->setTrackMask(CTrack::hasRich1Match | CTrack::hasRich2Match);
segment->setClosestRich1Ring((*(rich1->getRingList()))(indexOfBestInRich1));
segment->setClosestRich2Ring((*(rich2->getRingList()))(indexOfBestInRich2));
segment->setDeltaThetaRich12(smallestDeltaTheta);
segment->setDeltaPhiRich12(dphiOfBest);
//
// Some member make only sense with field
//
if (setup->getFieldFlag()) {
segment->setType(signOfTrack > 0 ? CTrack::Positron : CTrack::Electron);
theta1 = segment->getClosestRich1Ring()->getPolarCenter().getTheta();
phi1 = segment->getClosestRich1Ring()->getPolarCenter().getPhi();
segment->setMomentum(C4Momentum(ElectronMass,
C3Momentum(pOfBest*sin(theta1/1000.)*cos(phi1),
pOfBest*sin(theta1/1000.)*sin(phi1),
pOfBest*cos(theta1/1000.)) ));
}
richTracks.append(segment);
//
// Mark as used. Rich1 rings may be used twice (Vs) but only
// for tracks with opposite sign.
//
if (usedRich1Rings(indexOfBestInRich1) == notUsed)
usedRich1Rings(indexOfBestInRich1) =
signOfTrack > 0 ? usedForPositiveDeflectedTrack : usedForNegativeDeflectedTrack;
else
usedRich1Rings(indexOfBestInRich1) = used;
usedRich2Rings(indexOfBestInRich2) = used;
}
//
// Now all matching rings are put together.
// Check if there are any unresolved rings left and add
// them to the track list.
//
for (indexRich1 = 0; indexRich1 < usedRich1Rings.length(); indexRich1++) {
if (usedRich1Rings(indexRich1) == notUsed) {
segment = new CElectronTrack;
segment->setClosestRich1Ring((*(rich1->getRingList()))(indexRich1));
segment->setTrackMask(CTrack::hasRich1Match);
richTracks.append(segment);
}
}
for (indexRich2 = 0; indexRich2 < usedRich2Rings.length(); indexRich2++) {
if (usedRich2Rings(indexRich2) == notUsed) {
segment = new CElectronTrack;
segment->setClosestRich2Ring((*(rich2->getRingList()))(indexRich2));
segment->setTrackMask(CTrack::hasRich2Match);
richTracks.append(segment);
}
}
//
// Go through the completed list of Rich tracks and add
// the missing information about the next closest rings.
// The next closest is looked up only if the closest exist.
//
CRichXYCoord myXY, otherXY;
int indexOfClosest;
double distance, closestDistance;
const CFittedRing *thisRing;
const CFittedRing *nextRing;
for (int i = 0; i<richTracks.length(); i++) {
//
// Next closest ring in Rich-1
//
if (richTracks(i)->getTrackMask()&CTrack::hasRich1Match) {
thisRing = richTracks(i)->getClosestRich1Ring();
rich1->transform(myXY, thisRing->getPolarCenter());
closestDistance = FLT_MAX;
indexOfClosest = -1;
rich1Iter->reset();
while (nextRing = (*rich1Iter)()) {
if (nextRing == thisRing) continue;
rich1->transform(otherXY, nextRing->getPolarCenter());
distance = abs(otherXY-myXY);
if (distance < closestDistance) {
closestDistance = distance;
indexOfClosest = rich1Iter->pos();
}
}
if (indexOfClosest != -1) // add next closest R1 ring if any
richTracks(i)->setNextClosestRich1Ring((*(rich1->getRingList()))(indexOfClosest));
}
//
// Next closest ring in Rich-2
//
if (richTracks(i)->getTrackMask()&CTrack::hasRich2Match) {
thisRing = richTracks(i)->getClosestRich2Ring();
rich2->transform(myXY, thisRing->getPolarCenter());
closestDistance = FLT_MAX;
indexOfClosest = -1;
rich2Iter->reset();
while (nextRing = (*rich2Iter)()) {
if (nextRing == thisRing) continue;
rich2->transform(otherXY, nextRing->getPolarCenter());
distance = abs(otherXY-myXY);
if (distance < closestDistance) {
closestDistance = distance;
indexOfClosest = rich2Iter->pos();
}
}
if (indexOfClosest != -1) // add next closest R2 ring if any
richTracks(i)->setNextClosestRich2Ring((*(rich2->getRingList()))(indexOfClosest));
}
}
return True;
}
CBoolean CElectronTrackingStrategy::linkRichAndSidcTracks()
{
int i, k;
double thetaRich;
CAngle phiRich;
CElectronTrack *richSegment;
CTrack *sidcSegment;
CElectronTrack *track;
int indexOfBestMatch;
double distance2, smallestDistance2, scale;
double dphi, dtheta, dphiOfBest, dthetaOfBest;
CRichPolarCoord polarRichCoord;
double phiWindow = setup->getPhiMatchWindowRichSidc(); // in rad
double thetaWindow = setup->getThetaMatchWindowRichSidc(); // in mrad
//
// First the easy cases:
// No RICH- or no SiDC-segments or no segments at all.
// Then the hard case:
// Merge matching SiDC and RICH segments together.
//
if (richTracks.length() == 0 || sidcTracks->length() == 0)
return True;
else {
//
// Loop RICH segments
//
for (i=0; i<richTracks.length(); i++) {
richSegment = richTracks(i);
//
// Get theta and phi of track.
// Prefer Rich1 coordinates (if available)
//
if (richSegment->getTrackMask()&CTrack::hasRich1Match)
polarRichCoord = richSegment->getClosestRich1Ring()->getPolarCenter();
else if (richSegment->getTrackMask()&CTrack::hasRich2Match)
polarRichCoord = richSegment->getClosestRich2Ring()->getPolarCenter();
else {
cerr << "CTrackManager::linkRichAndSidcTracks():\n";
cerr << "\tERROR:\n";
cerr << "\tFound RICH segment without any ring. This is fatal.\n";
cerr << "\tNo further action.\n";
return False;
}
thetaRich = polarRichCoord.getTheta(); // in mrad
phiRich = polarRichCoord.getPhi(); // in rad
//
// Loop SiDC segments
//
indexOfBestMatch = -1;
smallestDistance2 = FLT_MAX;
for (k=0; k<sidcTracks->length(); k++) {
sidcSegment = (*sidcTracks)(k);
if (fabs(sidcSegment->getPhiSidc() - phiRich) < phiWindow &&
fabs(sidcSegment->getThetaSidc() - thetaRich) < thetaWindow) {
//
// Calculate the squared distance in theta units
//
scale = sin(thetaRich/1000.) * 1000.; // phi rad -> theta mrad
dtheta = sidcSegment->getThetaSidc() - thetaRich;
dphi = (double) (sidcSegment->getPhiSidc() - phiRich);
distance2 = dtheta*dtheta + dphi*dphi*scale*scale;
if (distance2 < smallestDistance2) {
smallestDistance2 = distance2;
indexOfBestMatch = k;
dphiOfBest = dphi;
dthetaOfBest = dtheta;
}
}
}
//
// If matching segment is found merge the two segments and put
// them into the track list.
//
if (indexOfBestMatch != -1) {
//
// Create as copy of RICH track
//
track = new CElectronTrack(*richSegment);
//
// Merge data from SiDC track
//
sidcSegment = (*sidcTracks)(indexOfBestMatch);
track->setVertex(sidcSegment->getVertex());
track->setClosestSidc1Hit(sidcSegment->getClosestSidc1Hit());
track->setClosestSidc2Hit(sidcSegment->getClosestSidc2Hit());
track->setNextClosestSidc1Hit(sidcSegment->getNextClosestSidc1Hit());
track->setNextClosestSidc2Hit(sidcSegment->getNextClosestSidc2Hit());
track->setThetaSidc(sidcSegment->getThetaSidc());
track->setPhiSidc(sidcSegment->getPhiSidc());
//
// Add info concerning RICH-SiDC match
//
track->setDeltaRSidc12(sidcSegment->getDeltaRSidc12());
track->setDeltaPhiSidc12(sidcSegment->getDeltaPhiSidc12());
track->setTrackMask(track->getTrackMask()|sidcSegment->getTrackMask());
track->setDeltaThetaSidcRich(dthetaOfBest);
track->setDeltaPhiSidcRich(dphiOfBest);
tracks.append(track);
}
}
}
return True;
}