CSidcTrackingStrategy.C
//-----------------------------------------------------------------------------
// $Header: /asis/offline/ceres/cool/project/RCS/CSidcTrackingStrategy.C,v 3.5 1997/06/03 10:52:44 messer Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class CSidcTrackingStrategy.
//
//-----------------------------------------------------------------------------
#include "CSidcTrackingStrategy.h"
#include "CVertex.h"
#include "CSidc.h"
#include "CListIterator.h"
#include "CSortedListIterator.h"
#include "CCyclicIterator.h"
#include "CElectronTrack.h"
CSidcTrackingStrategy::CSidcTrackingStrategy(CVertex* vtx, CSidc* s1,
CSidc* s2, const char* file)
: CTrackingStrategy(file), vertex(vtx), sidc1(s1), sidc2(s2)
{
sidcTracks = new CSortedList<CTrack>;
}
CBoolean CSidcTrackingStrategy::getSidcTracks(CSortedList<CTrack>& trackList, unsigned long mask)
{
//
// Create deep copy of selected tracks
//
if (!sidcTracks) {
cerr << "CSidcTrackingStrategy::getSidcTracks():" << endl;
cerr << "\tERROR:" << endl;
cerr << "\t The internal track list does not exist." << endl;
cerr << "\tNo further action." << endl;
return False;
}
CSortedListIterator<CTrack> sidcTrackIter(*sidcTracks);
CTrack *sidcTrack;
while (sidcTrack = sidcTrackIter()) {
if ((sidcTrack->getTrackMask()&mask) == mask) {
trackList.insert(new CTrack(*(sidcTrack)));
}
}
return True;
}
CSidcTrackingStrategy::~CSidcTrackingStrategy()
{
delete sidcTracks;
}
// Reconstruct SiDC segments
// =========================
//
// Reconstruct SiDC tracks using both SiDCs and the
// found vertex. Therefore it is mandatory to reconstruct
// a vertex prior to this routine.
//
// Method:
// Track segments are reconstructed using the SiDC-2
// hits and the vertex. For each track segment the
// predicted hit position in SiDC-1 is calculated and
// compared with the actual hits in a restricted region
// given by the 'phiMatchWindowSidc12' and 'radiusMatchWindowSidc12'
// parameter of the setup. If a match to sidc1 is found,
// this is flagged in the 'segmentMask' of the segment.
// The phi and theta are then the mean value of both
// contibuting hits.
// Closest and next closest hits of both sidc1 and sidc2
// are stored with the segment.
// For the remaining hits in Sidc1 which have so far not matched
// the closest and next closest hits in Sidc2 stored as well.
//
CBoolean CSidcTrackingStrategy::makeSidcTracks()
{
if (!sidcTracks) {
cerr << "CSidcTrackingStrategy::makeSidcTracks():" << endl;
cerr << "\tERROR:" << endl;
cerr << "\t The internal track list does not exist." << endl;
cerr << "\tNo further action." << endl;
return False;
}
//
// remove old stuff
//
sidcTracks->clearAndDestroy();
//
// anything to do?
//
if (sidc1==0 || sidc2==0 || vertex==0 ||
sidc1->getHits().length() == 0 || sidc2->getHits().length() == 0)
return True;
//
// temporary array to mark hits as used
//
RWTValVector<CBoolean> sidc1HitIsUsed(sidc1->getHits().length(), False);
//
// some variables needed in the loops below...
//
int jClosest = -1;
double dist2Min, drForClosest, dPhiForClosest;
CTrack *newSegment;
CEventCoord hit2Event, hit1Event;
CSidcHit *thisSidc2Hit;
CSidcHit *closestHit = 0;
CSidcHit *nextClosestHit = 0;
//
// Loop sidc2 hits and create the track segments
//
int i;
for (i=0; i<sidc2->getHits().length(); i++) {
thisSidc2Hit = sidc2->getHits()(i);
//
// By default every hit in sidc2 is a track segment.
//
newSegment = new CTrack;
newSegment->setVertex(vertex);
newSegment->setClosestSidc2Hit(thisSidc2Hit);
newSegment->setTrackMask(CTrack::hasVertex | CTrack::hasSidc2Match);
vertex->transform(hit2Event, thisSidc2Hit->getCenter());
getClosestHitInSidc(thisSidc2Hit, hit2Event, sidc2, closestHit, nextClosestHit,
jClosest, drForClosest, dPhiForClosest, dist2Min);
newSegment->setNextClosestSidc2Hit(nextClosestHit);
//
// Find closest and next closest hits to the projection into sidc1.
//
getClosestHitInSidc(thisSidc2Hit, hit2Event, sidc1, closestHit, nextClosestHit,
jClosest, drForClosest, dPhiForClosest, dist2Min);
//
// If the best match is good enough, store sidc1 data in this segment
// If not, theta and phi of this segment are those of the sidc2-hit
//
newSegment->setClosestSidc1Hit(closestHit);
newSegment->setDeltaRSidc12(drForClosest);
newSegment->setDeltaPhiSidc12(dPhiForClosest);
newSegment->setNextClosestSidc1Hit(nextClosestHit);
if (dist2Min < 1) {
sidc1HitIsUsed[jClosest] = True;
newSegment->setTrackMask(newSegment->getTrackMask() | CTrack::hasSidc1Match);
vertex->transform(hit1Event, closestHit->getCenter());
newSegment->setThetaSidc((hit1Event.getTheta() + hit2Event.getTheta()) / 2);
newSegment->setPhiSidc(average(hit2Event.getPhi(), hit1Event.getPhi()));
} else {
newSegment->setThetaSidc(hit2Event.getTheta());
newSegment->setPhiSidc(hit2Event.getPhi());
}
sidcTracks->insert(newSegment);
} // end loop over sidc2 hits
//
// Now look for closest and next closest hit in sidc2 for those hits in
// sidc1 which have so far not been included into a segment.
//
CSidcHit *thisSidc1Hit;
for (i=0; i<sidc1->getHits().length(); i++) {
if (!sidc1HitIsUsed[i]) {
//
// Calculate theta and phi of this sidc1-hit in event-coordinates i.e. relative
// to the reconstructed vertex
//
thisSidc1Hit = sidc1->getHits()(i);
vertex->transform(hit1Event, thisSidc1Hit->getCenter());
//
// By default every remaining hit in sidc1 will now also become a track segment.
// As we cannot expect a match to sidc2 here, theta of the segment is always
// theta of the hit in sidc1
//
newSegment = new CTrack;
newSegment->setVertex(vertex);
newSegment->setClosestSidc1Hit(sidc1->getHits()(i));
newSegment->setTrackMask(CTrack::hasVertex | CTrack::hasSidc1Match);
getClosestHitInSidc(thisSidc1Hit, hit1Event, sidc1, closestHit, nextClosestHit,
jClosest, drForClosest, dPhiForClosest, dist2Min);
newSegment->setNextClosestSidc1Hit(nextClosestHit);
newSegment->setThetaSidc(hit1Event.getTheta());
newSegment->setPhiSidc(hit1Event.getPhi());
//
// Find closest and next closest hits to the projection into sidc2.
//
getClosestHitInSidc(thisSidc1Hit, hit1Event, sidc2, closestHit, nextClosestHit,
jClosest, drForClosest, dPhiForClosest, dist2Min);
newSegment->setClosestSidc2Hit(closestHit);
newSegment->setDeltaRSidc12(drForClosest);
newSegment->setDeltaPhiSidc12(dPhiForClosest);
newSegment->setNextClosestSidc2Hit(nextClosestHit);
sidcTracks->insert(newSegment);
}
}
return True;
}
//
// Look for the closest and next closest hit with respect to the given 'referenceHit'
// which is not too far away in phi.
// The (presumably in phi sorted) hit-list is scanned from the 'left' to the 'right'
// with the best match in phi minus the maximum dphi as starting point.
// The so found closest and next closest
// hits are returned by the according pointer arguments of this function.
//
void CSidcTrackingStrategy::getClosestHitInSidc(const CSidcHit* referenceHit, CEventCoord referenceEvent,
const CSidc* sidc, CSidcHit*& closestHit,
CSidcHit*& nextClosestHit, int& jClosest,
double& drForClosest, double& dPhiForClosest,
double& dist2Min)
{
double PhiMatchFactor = setup->getPhiMatchFactorSidc12();
double sidcMatchR = setup->getRadiusMatchWindowSidc12();
double sidcMatchPhi = setup->getPhiMatchWindowSidc12();
double sidcMatchR2 = sidcMatchR*sidcMatchR;
double sidcMatchPhi2 = sidcMatchPhi*sidcMatchPhi;
double dr, dist2;
double dPhi;
double dist2Next = PhiMatchFactor * PhiMatchFactor;
dist2Min = PhiMatchFactor * PhiMatchFactor;
CAngle maxDPhi = CAngle(sidcMatchPhi) * PhiMatchFactor;
//
// Calculate start and stop index in the hit list using the built-in binary search.
//
CLabCylinderCoord leftCoord(referenceHit->getCenter());
CLabCylinderCoord rightCoord(referenceHit->getCenter());
leftCoord.setPhi(referenceHit->getCenter().getPhi() - maxDPhi);
rightCoord.setPhi(referenceHit->getCenter().getPhi() + maxDPhi);
CSidcHit leftHit(leftCoord, referenceHit->getAmp());
CSidcHit rightHit(rightCoord, referenceHit->getAmp());
int jStart = sidc->getHits().getIndexOfBestMatch(leftHit);
int jStop = sidc->getHits().getIndexOfBestMatch(rightHit);
//
// loop over hits in the sidc given as argument
//
CCyclicIterator<CSidcHit, CSortedList<CSidcHit> > hitIterator(sidc->getHits());
CSidcHit* thisHit;
CLabCylinderCoord referenceLab;
closestHit = nextClosestHit = 0;
hitIterator.setPosition(jStart);
while(hitIterator.getPosition()-1 != jStop) {
//
// Calculate elliptic distance always using the correct z-coordinate.
//
thisHit = hitIterator.next();
referenceEvent.setZ(thisHit->getCenter().getZ() - vertex->getZ());
vertex->transform(referenceLab, referenceEvent);
dPhi = double(referenceLab.getPhi() - thisHit->getCenter().getPhi());
dr = referenceLab.getRadius() - thisHit->getCenter().getRadius();
dist2 = dr*dr/sidcMatchR2 + dPhi*dPhi/sidcMatchPhi2;
if (dist2 < dist2Min) {
jClosest = hitIterator.getPosition()-1;
nextClosestHit = closestHit;
closestHit = thisHit;
dist2Next = dist2Min;
dist2Min = dist2;
drForClosest = dr;
dPhiForClosest = dPhi;
}
else if (dist2 < dist2Next) {
nextClosestHit = thisHit;
dist2Next = dist2;
}
}
}
CSortedList<CTrack>* CSidcTrackingStrategy::detachInternalSidcTrackList()
{
CSortedList<CTrack> *copyToReturn = sidcTracks;
sidcTracks = 0;
return copyToReturn;
}