CCandidatory.C
//-----------------------------------------------------------------------------
// $Header: /tmp_mnt/asis/offline/ceres/cool/project/RCS/CCandidatory.C,v 1.7 1997/04/25 15:05:38 messer Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class ...
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include "CCandidatory.h"
#include "CDeflection.h"
#include "CListIterator.h"
#include "CSortedListIterator.h"
#include "CCyclicIterator.h"
CCandidatory::CCandidatory(const char *setupFile)
{
candidateList = 0;
setup = new CCandidatorySetup;
if (setupFile == 0) {
CString setupName = CString(C_DEFAULT_SETUP_PATH)+
CString(C_SETUPFILE_CANDIDATORY);
setup->read(setupName.data());
}
else
setup->read(setupFile);
}
CCandidatory::~CCandidatory()
{
delete setup;
delete candidateList;
}
CCandidatory::CCandidatory(const CCandidatory &)
{
}
CCandidatory & CCandidatory::operator = (const CCandidatory &)
{
return *this;
}
void CCandidatory::listSetup(ostream& ost) const
{
ost << "\nSetup listing of CCandidatory " << endl;
CString line('=', 46);
ost << line << endl;
setup->list(ost);
}
void CCandidatory::copyCandidatesToRings(CRich& rich)
{
CListIterator<CRingCandidate> candIter(*(rich.getCandidateList()));
CRingCandidate *thisCandidate;
CRichPadCoord padCenter;
CRichPolarCoord polarCenter;
CFittedRing *newRing;
CList<CFittedRing> *ringList = new CList<CFittedRing>;
while (thisCandidate = candIter()) {
padCenter = thisCandidate->getCenter();
rich.transform(polarCenter, padCenter);
newRing = new CFittedRing;
newRing->setCenter(padCenter, polarCenter);
newRing->setType(CRing::Electron);
newRing->setRadius(thisCandidate->getRadius());
ringList->append(newRing);
}
//
// hand the control over those rings to rich1
//
delete rich.detachRingList();
rich.attachRingList(ringList);
}
CBoolean CCandidatory::makeHoughCandidates(CSortedList<CElectronTrack>& sidcTracks, CRich1& rich1)
{
//
// anything to do?
//
if (rich1.getCandidateList()->length() == 0) return True;
if (sidcTracks.length() == 0) {
rich1.getCandidateList()->clearAndDestroy();
return True;
}
CCyclicIterator<CElectronTrack, CSortedList<CElectronTrack> > sidcTrackIter(sidcTracks);
CListIterator<CRingCandidate> candIter(*(rich1.getCandidateList()));
CRingCandidate *thisCandidate;
CElectronTrack *thisSidcTrack;
CElectronTrack startTrack, stopTrack;
int jStart, jStop;
CRichPolarCoord candPolar, trackPolar;
CRichPadCoord trackPad;
CAngle maxDPhi = CAngle(setup->getHoughMaxDeltaXYRich1Sidc() * (2.*5./1000.));
CBoolean thisCandidateHasAnSidcMatch;
if (!candidateList) candidateList = new CList<CRingCandidate>;
candidateList->clearAndDestroy();
//
// Copy all candidates which have a match to an sidc track segment
// to the internal candidate list. We assume that the sidc tracks are
// sorted in phi and one apply a binary search. This is guarantied by
// the way in which these segments are made. See CSidcTrackingStrategy for details.
//
while (thisCandidate = candIter()) {
rich1.transform(candPolar, thisCandidate->getCenter());
thisCandidateHasAnSidcMatch = False;
sidcTrackIter.reset();
startTrack.setPhiSidc(candPolar.getPhi() - maxDPhi);
stopTrack.setPhiSidc(candPolar.getPhi() + maxDPhi);
jStart = sidcTracks.getIndexOfBestMatch(startTrack);
jStop = sidcTracks.getIndexOfBestMatch(stopTrack);
sidcTrackIter.setPosition(jStart);
while (sidcTrackIter.getPosition()-1 != jStop) {
thisSidcTrack = sidcTrackIter.next();
//
// transform into required coordinates
//
trackPolar.setTheta(thisSidcTrack->getThetaSidc());
trackPolar.setPhi(thisSidcTrack->getPhiSidc());
rich1.transform(trackPad, trackPolar);
//
// If a match was found, leave the track-loop.
//
if (abs(trackPad - thisCandidate->getCenter()) < setup->getHoughMaxDeltaXYRich1Sidc()) {
thisCandidateHasAnSidcMatch = True;
break;
}
}
if (thisCandidateHasAnSidcMatch)
candidateList->append(new CRingCandidate(*thisCandidate));
}
//
// hand the control over the so found candidates to rich1
//
delete rich1.detachCandidateList();
rich1.attachCandidateList(candidateList);
candidateList = 0;
return True;
}
CBoolean CCandidatory::makeHoughCandidates(CSortedList<CElectronTrack>& externalTracks, CRich2& rich2,
CPadChamber& padChamber)
{
//
// anything to do?
//
if (rich2.getCandidateList()->length() == 0) return True;
if (externalTracks.length() == 0) {
rich2.getCandidateList()->clearAndDestroy();
return True;
}
//
// variables needed in the loops below
//
CSortedListIterator<CElectronTrack> trackIter(externalTracks);
CListIterator<CRingCandidate> candIter(*(rich2.getCandidateList()));
CRingCandidate *thisCandidate;
CElectronTrack *thisTrack;
CRichPolarCoord candPolar;
double thetaBefor, phiBefor, thetaAfter, phiAfter;
CPadCPadCoord padcPredictor, padcHitCenter;
CDeflection deflection;
CBoolean thisCandidateHasAPadCMatch;
if (!candidateList) candidateList = new CList<CRingCandidate>;
candidateList->clearAndDestroy();
//
// copy all candidates which have a match to
// an external track segment to the internal candidate list
//
while (thisCandidate = candIter()) {
thisCandidateHasAPadCMatch = False;
rich2.transform(candPolar, thisCandidate->getCenter());
trackIter.reset();
while (thisTrack = trackIter()) {
thetaBefor = thisTrack->getClosestRich1Ring()->getPolarCenter().getTheta();
phiBefor = (double) thisTrack->getClosestRich1Ring()->getPolarCenter().getPhi();
thetaAfter = candPolar.getTheta();
phiAfter = (double) candPolar.getPhi();
deflection.calculatePCPrediction(thetaBefor, phiBefor, thetaAfter, phiAfter,
padChamber, (CVertex&)*(thisTrack->getVertex()));
padcPredictor = deflection.getPadcPredictor();
//
// If a match was found, leave the track-loop.
//
padcHitCenter.setX(thisTrack->getClosestPadCHit()->getX());
padcHitCenter.setY(thisTrack->getClosestPadCHit()->getY());
if (abs(padcPredictor - padcHitCenter) <
setup->getHoughMaxDeltaXYRich2Padc()) {
thisCandidateHasAPadCMatch = True;
break;
}
}
if (thisCandidateHasAPadCMatch)
candidateList->append(new CRingCandidate(*thisCandidate));
}
//
// hand the control over the so found candidates to rich2
//
delete rich2.detachCandidateList();
rich2.attachCandidateList(candidateList);
candidateList = 0;
return True;
}
CBoolean CCandidatory::makeHoughCandidates(CVertex& vertex, CRich1& rich1, CRich2& rich2,
CPadChamber& padChamber)
{
//
// anything to do?
//
if (rich1.getRingList()->length() == 0) return True;
if (rich2.getCandidateList()->length() == 0) return True;
if (padChamber.getHits().length() == 0) return True;
CListIterator<CFittedRing> rich1RingIter(*(rich1.getRingList()));
CListIterator<CRingCandidate> rich2CandIter(*(rich2.getCandidateList()));
CFittedRing *thisRich1Ring;
double theta1, phi1, theta2, phi2;
CRingCandidate *thisRich2Candidate;
CRichPolarCoord rich2Polar;
CDeflection padcDeflection;
CPadCPadCoord padcPredictor, padcHitCoord;
CRichLikeHit padcPredictorHit;
double padcMatchWindow = setup->getHoughMaxDeltaXYRich2Padc();
CBoolean thisCandidateHasAPadcMatch;
if (!candidateList) candidateList = new CList<CRingCandidate>;
candidateList->clearAndDestroy();
int iHit;
//
// Loop over all previously found rich2 candidates.
//
while (thisRich2Candidate = rich2CandIter()) {
thisCandidateHasAPadcMatch = False;
rich2.transform(rich2Polar, thisRich2Candidate->getCenter());
theta2 = rich2Polar.getTheta();
phi2 = rich2Polar.getPhi();
rich1RingIter.reset();
//
// Loop over all rich1 rings and see if the combination of the present
// rich2 candidate with a rich1 ring matches to the padchamber.
// The theta match is set to a fixed value.
// Leave the loop as soon as a match is found.
//
while (thisRich1Ring = rich1RingIter()) {
theta1 = thisRich1Ring->getPolarCenter().getTheta();
if (fabs(theta1 - theta2) > setup->getRich12ThetaMatch()) continue;
phi1 = thisRich1Ring->getPolarCenter().getPhi();
padcDeflection.calculatePCPrediction(theta1, phi1, theta2, phi2, padChamber, vertex);
padcPredictor = padcDeflection.getPadcPredictor();
//
// Loop over the in x sorted hit-list in the padchamber,
// beginning at the furthermost possible index (predictor - matchwindow).
// Leave the loop as soon as a match is found.
//
padcPredictorHit.setX(padcPredictor.getX() - padcMatchWindow);
iHit = padChamber.getHits().getIndexOfBestMatch(padcPredictorHit);
while (iHit < padChamber.getHits().length() &&
padChamber.getHits()(iHit)->getX()-padcPredictor.getX() < padcMatchWindow) {
padcHitCoord.setX(padChamber.getHits()(iHit)->getX());
padcHitCoord.setY(padChamber.getHits()(iHit)->getY());
if (abs(padcPredictor - padcHitCoord) < padcMatchWindow) {
thisCandidateHasAPadcMatch = True;
break;
}
iHit++;
}
if (thisCandidateHasAPadcMatch) break;
}
if (thisCandidateHasAPadcMatch)
candidateList->append(new CRingCandidate(*thisRich2Candidate));
}
//
// hand the control over the so found candidates to rich2
//
delete rich2.detachCandidateList();
rich2.attachCandidateList(candidateList);
candidateList = 0;
return True;
}
CBoolean CCandidatory::makeHoughCandidates(CVertex& vertex, CRich2& rich2, CPadChamber& padChamber)
{
//
// anything to do?
//
if (rich2.getCandidateList()->length() == 0) return True;
if (padChamber.getHits().length() == 0) return True;
CListIterator<CRingCandidate> rich2CandIter(*(rich2.getCandidateList()));
CRingCandidate *thisRich2Candidate;
CRichPolarCoord rich2Polar;
double theta2, phi2, thetaPadc, phiPadc, deflection, p, dTheta;
double maxDeltaTheta, maxDeflection;
CLabXYZCoord padcXYZCoord;
CEventCoord padcEventCoord;
CBoolean thisCandidateHasAPadcMatch;
if (!candidateList) candidateList = new CList<CRingCandidate>;
candidateList->clearAndDestroy();
int iHit;
//
// Loop over all previously found rich2 candidates.
//
while (thisRich2Candidate = rich2CandIter()) {
thisCandidateHasAPadcMatch = False;
rich2.transform(rich2Polar, thisRich2Candidate->getCenter());
theta2 = rich2Polar.getTheta();
phi2 = rich2Polar.getPhi();
//
// For each padChamber hit check whether it matches within a
// butterly to this candidate. As soon as a match is found we
// leave the loop.
//
for (iHit=0; iHit<padChamber.getHits().length(); iHit++) {
padcXYZCoord.setX(padChamber.getHits()(iHit)->getX());
padcXYZCoord.setY(padChamber.getHits()(iHit)->getY());
padcXYZCoord.setZ(padChamber.getSetup()->getZPosition());
vertex.transform(padcEventCoord, padcXYZCoord);
thetaPadc = padcEventCoord.getTheta();
phiPadc = padcEventCoord.getPhi();
deflection = double(phiPadc - phi2);
//
// Get momentum of of this combination.
// For no-field data the lowest momentum according to the
// given pt-cut is assumed (later needed for butterfly).
//
if (setup->getFieldFlag())
p = rich2PadcPhiDeflectionAt1GeV(theta2)/fabs(deflection);
else
p = setup->getPtCutForButterfly()/sin(theta2/1000.);
//
// Check if padchamber hit 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.
//
dTheta = thetaPadc - theta2;
maxDeltaTheta = setup->getButterflySigmas() * padChamberButterfly(p);
if (setup->getFieldFlag())
maxDeflection = rich2PadcPhiDeflectionAt1GeV(theta2)/(setup->getPtCutForButterfly()/
sin(theta2/1000.));
else
maxDeflection = maxDeltaTheta/sin(theta2/1000.);
if (fabs(deflection) < maxDeflection || fabs(dTheta) < maxDeltaTheta) {
thisCandidateHasAPadcMatch = True;
break;
}
}
if (thisCandidateHasAPadcMatch)
candidateList->append(new CRingCandidate(*thisRich2Candidate));
}
//
// hand the control over the so found candidates to rich2
//
delete rich2.detachCandidateList();
rich2.attachCandidateList(candidateList);
candidateList = 0;
return True;
}
double CCandidatory::padChamberButterfly(double p)
{
//
// Returns the one sigma fluctuation in theta due
// to hit resolution in the tracking detectors and multiple
// scattering. Units are mrad.
//
double ringResolution2 = setup->getRich2CenterResolution() *
setup->getRich2CenterResolution() +
setup->getPadcHitResolution() *
setup->getPadcHitResolution();
double multScatter2 = setup->getMultipleScattering() *
setup->getMultipleScattering();
return sqrt(ringResolution2 + multScatter2 / (p*p));
}
double CCandidatory::rich2PadcPhiDeflectionAt1GeV(double theta)
{
//
// Returns phi deflection for 1 GeV vs. theta (mrad)
//
return ((setup->getPhiDeflectionParameters(0) +
setup->getPhiDeflectionParameters(1) * theta +
setup->getPhiDeflectionParameters(2) * theta * theta)*
setup->getScaleFactorForPadcPhiDeflection());
}