CExternalElectronTrackingStrategy.C
//-----------------------------------------------------------------------------
// $Header: /asis/offline/ceres/cool/project/RCS/CExternalElectronTrackingStrategy.C,v 1.6 1997/06/03 10:56:10 messer Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class CExternalElectronTrackingStrategy.
//
//-----------------------------------------------------------------------------
#include "rw/tvvector.h"
#include "CExternalElectronTrackingStrategy.h"
#include "CListIterator.h"
#include "CSortedListIterator.h"
#include "CDeflection.h"
CExternalElectronTrackingStrategy::CExternalElectronTrackingStrategy(CVertex* vtx, CSidc* s1,
CSidc* s2, CRich1* r1,
CRich2 *r2, CPadChamber* pc,
const char* file)
: CSidcTrackingStrategy(vtx, s1, s2, file), rich1(r1), rich2(r2), padChamber(pc)
{
richPadcTracks = new CList<CElectronTrack>;
tracks = new CSortedList<CElectronTrack>;
}
CExternalElectronTrackingStrategy::~CExternalElectronTrackingStrategy()
{
delete richPadcTracks;
delete tracks;
}
CBoolean CExternalElectronTrackingStrategy::makeTracks()
{
//
// Check for existence of internal lists. Maybe one of them has been detached.
//
if (!sidcTracks || !richPadcTracks || !tracks) {
cerr << "CExternalElectronTrackingStrategy::makeTracks():" << endl;
cerr << "\tERROR:" << endl;
cerr << "\t One of the internal track lists does not exist" << endl;
cerr << "\tNo further action." << endl;
return False;
}
//
// clear the list of tracks
//
tracks->clearAndDestroy();
//
// make sidc tracks if necessary
//
if (sidcTracks->length() == 0)
if (!makeSidcTracks()) return False;
//
// make rich and padc tracks if necessary
//
if (richPadcTracks->length() == 0)
if (!makeRichPadcTracks()) return False;
//
// link everything
//
if (!linkAllTracks()) return False;
return True;
}
CBoolean CExternalElectronTrackingStrategy::getTracks(CSortedList<CElectronTrack>& list,
unsigned long mask)
{
if (!sidcTracks || !richPadcTracks || !tracks) {
cerr << "CExternalElectronTrackingStrategy::getTracks():" << endl;
cerr << "\tERROR:" << endl;
cerr << "\t One of the internal track lists does not exist" << endl;
cerr << "\tNo further action." << endl;
return False;
}
CSortedListIterator<CTrack> sidcTrackIter(*sidcTracks);
CListIterator<CElectronTrack> richPadcTrackIter(*richPadcTracks);
CSortedListIterator<CElectronTrack> trackIter(*tracks);
CTrack *sidcTrack;
CElectronTrack *richPadcTrack, *track;
//
// First see in which list we have to look for the stuff
//
CBoolean needsSidcTracks = mask&CTrack::hasVertex||mask&CTrack::hasSidc1Match||
mask&CTrack::hasSidc2Match;
CBoolean needsRichPadcTracks = mask&CTrack::hasRich1Match||mask&CTrack::hasRich2Match||
mask&CTrack::hasPadCMatch;
CBoolean needsTracks = needsSidcTracks && needsRichPadcTracks;
CBoolean needsThemAll = !needsSidcTracks && !needsRichPadcTracks;
//
// Scan the internal list and create copies of the matching entries
//
if (needsTracks || needsThemAll) {
//
// Look in track list only
//
trackIter.reset();
while (track = trackIter()) {
if ((track->getTrackMask()&mask) == mask)
list.insert(new CElectronTrack(*track));
}
}
if (needsSidcTracks || needsThemAll) {
//
// Look in list of sidc segments only.
// Merge SiDC track into electron track
//
sidcTrackIter.reset();
while (sidcTrack = sidcTrackIter()) {
if ((sidcTrack->getTrackMask()&mask) == mask) {
track = new CElectronTrack;
*((CTrack*)track) = *(sidcTrack);
list.insert(track);
}
}
}
if (needsRichPadcTracks || needsThemAll) {
//
// Look in list of rich segments only
//
richPadcTrackIter.reset();
while (richPadcTrack = richPadcTrackIter()) {
if ((richPadcTrack->getTrackMask()&mask) == mask)
list.insert(new CElectronTrack(*richPadcTrack));
}
}
return True;
}
CBoolean CExternalElectronTrackingStrategy::makeRichPadcTracks()
{
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 = 0;
CListIterator<CFittedRing> *rich2Iter = 0;
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, indexPadc;
CDeflection padcDeflection;
CPadCPadCoord padcPredictor, padcHitCenter;
CLabXYZCoord padcXYZ;
CEventXYZCoord eventXYZ;
CEventCoord padcPredictorPolar, padcHitCenterPolar;
double bestMatchToPadcTheta, bestMatchToPadcThetaForBest;
CAngle bestMatchToPadcPhi, bestMatchToPadcPhiForBest;
double dxRichPadc, matchToPadc, bestMatchToPadc, secondBestMatchToPadc;
CBoolean matchesToPadc, smallestDeltaThetaMatchesToPadc;
CRichLikeHit padcPredictorHit;
CRichLikeHit *closestPadcHit, *nextClosestPadcHit;
CRichLikeHit *closestPadcHitForBest, *nextClosestPadcHitForBest;
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
smallestDeltaThetaMatchesToPadc = False;
indexOfBestInRich1 = -1; // -1 == no further match
indexOfBestInRich2 = -1;
closestPadcHitForBest = 0;
nextClosestPadcHitForBest = 0;
//
// 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 Rich1 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 then 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;
//
// Calculate prediction point in the padChamber for this ring combination
// and get the best matching hits. We assume that the list of padc-hits
// is sorted in x and one can therefore apply the built-in binary-search.
//
if (padChamber && vertex) {
padcDeflection.calculatePCPrediction(theta1, (double)phi1, theta2, (double)phi2,
(CPadChamber&) *padChamber, (CVertex&) *vertex);
padcPredictor = padcDeflection.getPadcPredictor();
padcPredictorHit.setX(padcPredictor.getX() - setup->getMaxDeltaXYRichPadc());
closestPadcHit = nextClosestPadcHit = 0;
bestMatchToPadc = secondBestMatchToPadc = 1000.;
indexPadc = padChamber->getHits().getIndexOfBestMatch(padcPredictorHit);
dxRichPadc = -setup->getMaxDeltaXYRichPadc();
while (dxRichPadc < setup->getMaxDeltaXYRichPadc() &&
indexPadc < padChamber->getHits().length()) {
dxRichPadc = padChamber->getHits()(indexPadc)->getX() - padcPredictor.getX();
padcHitCenter.setX(padChamber->getHits()(indexPadc)->getX());
padcHitCenter.setY(padChamber->getHits()(indexPadc)->getY());
matchToPadc = abs(padcPredictor - padcHitCenter);
if (matchToPadc < bestMatchToPadc) {
secondBestMatchToPadc = bestMatchToPadc;
bestMatchToPadc = matchToPadc;
nextClosestPadcHit = closestPadcHit;
closestPadcHit = padChamber->getHits()(indexPadc);
padChamber->transform(padcXYZ, padcPredictor);
vertex->transform(eventXYZ, padcXYZ);
::transform(padcPredictorPolar,eventXYZ);
padChamber->transform(padcXYZ, padcHitCenter);
vertex->transform(eventXYZ, padcXYZ);
::transform(padcHitCenterPolar,eventXYZ);
bestMatchToPadcTheta = padcPredictorPolar.getTheta() - padcHitCenterPolar.getTheta();
bestMatchToPadcPhi = padcPredictorPolar.getPhi() - padcHitCenterPolar.getPhi();
}
else if (matchToPadc < secondBestMatchToPadc) {
secondBestMatchToPadc = matchToPadc;
nextClosestPadcHit = padChamber->getHits()(indexPadc);
}
indexPadc++;
}
if (bestMatchToPadc < setup->getMaxDeltaXYRichPadc())
matchesToPadc = True;
else
matchesToPadc = False;
}
else
matchesToPadc = False;
//
// Store indices of best matching rings, favouring those with a
// match to the padChamber
//
if ((matchesToPadc && !smallestDeltaThetaMatchesToPadc) ||
(matchesToPadc == smallestDeltaThetaMatchesToPadc &&
fabs(dtheta) < fabs(smallestDeltaTheta))) {
smallestDeltaThetaMatchesToPadc = matchesToPadc;
smallestDeltaTheta = dtheta;
indexOfBestInRich1 = indexRich1;
indexOfBestInRich2 = indexRich2;
closestPadcHitForBest = closestPadcHit;
nextClosestPadcHitForBest = nextClosestPadcHit;
bestMatchToPadcThetaForBest = bestMatchToPadcTheta;
bestMatchToPadcPhiForBest = bestMatchToPadcPhi;
signOfTrack = deflection < 0 ? -1 : 1;
dphiOfBest = deflection;
pOfBest = p;
}
}
}
//
// Store best track so far, or stop if no more matching
// 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;
if (smallestDeltaThetaMatchesToPadc) {
segment->setTrackMask(CTrack::hasRich1Match|CTrack::hasRich2Match|CTrack::hasPadCMatch);
segment->setClosestPadCHit(closestPadcHitForBest);
segment->setNextClosestPadCHit(nextClosestPadcHitForBest);
segment->setDeltaThetaRich2PadC(bestMatchToPadcThetaForBest);
segment->setDeltaPhiRich2PadC((double) bestMatchToPadcPhiForBest);
}
else
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.)) ));
}
richPadcTracks->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);
richPadcTracks->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);
richPadcTracks->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.
//
CListIterator<CElectronTrack> richPadcTrackIter(*richPadcTracks);
CElectronTrack *richPadcTrack;
CRichXYCoord myXY, otherXY;
int indexOfClosest;
double distance, closestDistance;
const CFittedRing *thisRing;
const CFittedRing *nextRing;
while (richPadcTrack = richPadcTrackIter()) {
//
// Next closest ring in Rich-1
//
if (richPadcTrack->getTrackMask()&CTrack::hasRich1Match) {
thisRing = richPadcTrack->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
richPadcTrack->setNextClosestRich1Ring((*(rich1->getRingList()))(indexOfClosest));
}
//
// Next closest ring in Rich-2
//
if (richPadcTrack->getTrackMask()&CTrack::hasRich2Match) {
thisRing = richPadcTrack->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
richPadcTrack->setNextClosestRich2Ring((*(rich2->getRingList()))(indexOfClosest));
}
}
delete rich1Iter;
delete rich2Iter;
return True;
}
CBoolean CExternalElectronTrackingStrategy::linkAllTracks()
{
CListIterator<CElectronTrack> richPadcTrackIter(*richPadcTracks);
CSortedListIterator<CTrack> sidcTrackIter(*sidcTracks);
double thetaRich;
CAngle phiRich;
CElectronTrack *richSegment;
CTrack *sidcSegment;
CElectronTrack *track;
int indexOfBestMatch;
double distance2, smallestDistance2, scale;
double dphi, dtheta, dphiOfBest, dthetaOfBest;
CRichPolarCoord polarRichCoord;
unsigned long fullSidcTrackMask = CTrack::hasVertex|CTrack::hasSidc1Match|CTrack::hasSidc2Match;
CBoolean matchesToBothSidcs, bestMatchesToBothSidcs;
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 (richPadcTracks->length() == 0 || sidcTracks->length() == 0)
return True;
else {
//
// Loop RICH segments
//
while (richSegment = richPadcTrackIter()) {
//
// 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 << "CExternalElectronTrackingStrategy::linkallTracks():\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
//
// Find the best matching SiDC segment, always favouring those which
// have a match to both SiDCs.
//
indexOfBestMatch = -1;
smallestDistance2 = FLT_MAX;
sidcTrackIter.reset();
bestMatchesToBothSidcs = False;
while (sidcSegment = sidcTrackIter()) {
if (fabs(sidcSegment->getPhiSidc() - phiRich) < phiWindow &&
fabs(sidcSegment->getThetaSidc() - thetaRich) < thetaWindow) {
if ((sidcSegment->getTrackMask()&fullSidcTrackMask) == fullSidcTrackMask)
matchesToBothSidcs = True;
else
matchesToBothSidcs = False;
//
// 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 ((matchesToBothSidcs && !bestMatchesToBothSidcs) ||
(matchesToBothSidcs == bestMatchesToBothSidcs &&
distance2 < smallestDistance2)) {
if (matchesToBothSidcs)
bestMatchesToBothSidcs = True;
else
bestMatchesToBothSidcs = False;
smallestDistance2 = distance2;
indexOfBestMatch = sidcTrackIter.pos();
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->insert(track);
}
}
}
return True;
}