CPionTrackingStrategy.C
//-----------------------------------------------------------------------------
// $Header: CPionTrackingStrategy.C,v 2.1 1996/08/21 21:59:49 ceretto
//
// COOL Program Library // Copyright (C) CERES collaboration, 1996
//
// Implementation of class CPionTrackingStrategy.
//
//-----------------------------------------------------------------------------
#include "rw/tvvector.h"
#include "CPionTrackingStrategy.h"
#include "CSidc1.h"
#include "CSidc2.h"
#include "CPadChamber.h"
#include "CRich1.h"
#include "CRich2.h"
#include "CDeflection.h"
#include "CList.h"
#include "CSortedList.h"
#include "CListIterator.h"
#include "CSortedListIterator.h"
#include "CCyclicIterator.h"
#include "CElectronTrack.h"
CPionTrackingStrategy::CPionTrackingStrategy(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)
{
radiusVector.reshape(setup->getVectorSize());
xCoord.reshape(1000);
yCoord.reshape(1000);
binPos.reshape(1000);
richTracks = new CList<CTrack>; // analog to RichSegment
sidcPadTracks = new CSortedList<CTrack>; // sidc-pad tracks
tracks = new CSortedList<CTrack>; // final tracks
}
CPionTrackingStrategy::~CPionTrackingStrategy()
{
delete richTracks;
delete sidcPadTracks;
delete tracks;
}
CBoolean CPionTrackingStrategy::makeTracks()
{
if (!sidcTracks || !richTracks || !sidcPadTracks || !tracks) {
cerr << "CPionTrackingStrategy::makeTracks():" << endl;
cerr << "\tERROR:" << endl;
cerr << "\t One of the internal track lists does not exist" << endl;
cerr << "\tNo further action." << endl;
return False;
}
richTracks->clearAndDestroy();
sidcPadTracks->clearAndDestroy();
tracks->clearAndDestroy();
if (sidcTracks->length() == 0)
if (!makeSidcTracks()) return False;
if (!makePadChamberSidcTracks()) return False;
if (!makeRichTracksAndLinkToSidcPadTracks()) return False;
return True;
}
CBoolean CPionTrackingStrategy::getTracks(CSortedList<CTrack>& list, unsigned long mask)
{
if (!sidcTracks || !richTracks || !sidcPadTracks || !tracks) {
cerr << "CPionTrackingStrategy::getTracks():" << endl;
cerr << "\tERROR:" << endl;
cerr << "\t One of the internal track lists does not exist" << endl;
cerr << "\tNo further action." << endl;
return False;
}
CTrack *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
//
CSortedListIterator<CTrack> titer(*tracks);
CTrack *thisTrack;
while (thisTrack = titer()) {
if ((thisTrack->getTrackMask()&mask) == mask)
list.append(new CTrack(*thisTrack));
}
}
if (needsSidcSegments || needsThemAll) {
//
// Look in list of sidc segments only.
// Merge SiDC track into electron track
//
CSortedListIterator<CTrack> siter(*sidcTracks);
CTrack *sidcTrack;
while (sidcTrack = siter()) {
if ((sidcTrack->getTrackMask()&mask) == mask) {
track = new CTrack;
*((CTrack*)track) = *(sidcTrack);
list.append(track);
}
}
}
if (needsRichSegments || needsThemAll) {
//
// Look in list of rich segments only
//
CListIterator<CTrack> riter(*richTracks);
CTrack *richTrack;
while (richTrack = riter()) {
if ((richTrack->getTrackMask()&mask) == mask)
list.append(new CTrack(*richTrack));
}
}
return True;
}
CTrack* CPionTrackingStrategy::findBestPadCMatch(CTrack* thisTrack, CPadCPadCoord predictor,double rich2Theta, CAngle rich2Phi )
{
double thetaPadC;
CAngle phiPadC;
double distance;
int bestMatch = -1;
int next = -1;
double closestDistance = setup->getPCFensterMatch();
double nextClosestDistance = closestDistance;
CLabXYZCoord hitXYZ;
CEventCoord hitCyl;
CPadCPadCoord hitPosition;
CRichLikeHit *hit;
//
// loop over all Pad Chamber Hits
//
for (int k=0; k < padChamber->getHits().length(); k++) {
hit = padChamber->getHits()(k);
hitPosition.setX(hit->getX());
hitPosition.setY(hit->getY());
padChamber->transform(hitXYZ, hitPosition);
vertex->transform(hitCyl, hitXYZ);
distance = abs(predictor - hitPosition);
if( distance < closestDistance ) {
nextClosestDistance = closestDistance;
closestDistance = distance;
next = bestMatch;
bestMatch = k;
thetaPadC = hitCyl.getTheta();
phiPadC = hitCyl.getPhi();
}
else if(distance < nextClosestDistance) {
next = k;
nextClosestDistance = distance;
}
}
if(bestMatch > -1) {
thisTrack->setTrackMask(thisTrack->getTrackMask()|CTrack::hasPadCMatch);
thisTrack->setClosestPadCHit(padChamber->getHits()(bestMatch));
thisTrack->setThetaPadC(thetaPadC);
thisTrack->setPhiPadC(phiPadC);
thisTrack->setDeltaThetaSidcPadC(thisTrack->getThetaSidc() - thetaPadC);
thisTrack->setDeltaPhiSidcPadC(thisTrack->getPhiSidc() - phiPadC);
thisTrack->setDeltaThetaRich2PadC(rich2Theta - thetaPadC);
thisTrack->setDeltaPhiRich2PadC(rich2Phi - phiPadC);
if (next > -1)
thisTrack->setNextClosestPadCHit(padChamber->getHits()(next));
}
else {
thisTrack->setTrackMask(thisTrack->getTrackMask() & ~(CTrack::hasPadCMatch));
}
return thisTrack;
}
CBoolean CPionTrackingStrategy::makePadChamberSidcTracks()
{
//
// anything to do?
//
if (!sidcTracks || !sidcPadTracks) {
cerr << "CPionTrackingStrategy::makePadChamberSidcTracks():" << endl;
cerr << "\tERROR:" << endl;
cerr << "\t One of the internal track lists does not exist" << endl;
cerr << "\tNo further action." << endl;
return False;
}
if (sidcTracks->length() == 0 || padChamber->getHits().length() == 0) return True;
//
// make Pad chamber segment for every Sidc Segment
// the window is : (50 mrad (phi) * 6 mrad (theta))*4
//
int jStart, jStop;
double thetaSidc, thetaPad;
CAngle phiSidc, phiPad;
CRichLikeHit *hit;
CTrack *sidcSegment;
CTrack *track;
CTrack startTrack, stopTrack;
double dphi,dtheta;
CPadCPadCoord thisPadPad;
CLabXYZCoord thisPadXYZ;
CEventCoord thisPadEvent;
CCyclicIterator<CTrack, CSortedList<CTrack> > sidcSegmentsIterator(*sidcTracks);
//
// Loop over all PAD hits
//
for (int k = 0; k < padChamber->getHits().length(); k++) {
hit = padChamber->getHits()(k);
thisPadPad.setX(hit->getX());
thisPadPad.setY(hit->getY());
padChamber->transform(thisPadXYZ, thisPadPad);
vertex->transform(thisPadEvent, thisPadXYZ); // transform in Event Coord
thetaPad = thisPadEvent.getTheta(); // PadC hit coordinate
phiPad = thisPadEvent.getPhi();
//
// Loop over the sidcSegments to find all those which match to this padc hit.
// We apply a binary search in phi, assuming that the list of sidcSegments
// is sorted in phi.
//
startTrack.setPhiSidc(phiPad - (CAngle)setup->getPhiWindow());
stopTrack.setPhiSidc(phiPad + (CAngle)setup->getPhiWindow());
jStart = sidcTracks->getIndexOfBestMatch(startTrack);
jStop = sidcTracks->getIndexOfBestMatch(stopTrack);
sidcSegmentsIterator.setPosition(jStart);
while(sidcSegmentsIterator.getPosition()-1 != jStop) {
sidcSegment = sidcSegmentsIterator.next();
if(sidcSegment->getTrackMask() == (CTrack::hasSidc1Match | CTrack::hasSidc2Match |
CTrack::hasVertex)) {
thetaSidc = sidcSegment->getThetaSidc(); // in mrad
phiSidc = sidcSegment->getPhiSidc(); // in rad
dphi = (double) (phiSidc - phiPad) ; // dphi (sidc - pad)
dtheta = thetaSidc - thetaPad; // dtheta (sidc - pad)
if (fabs(dphi) < setup->getPhiWindow() &&
fabs(dtheta) < setup->getThetaWindow()) {
track= new CTrack(*sidcSegment);
track->setClosestPadCHit(hit);
track->setThetaPadC(thetaPad);
track->setPhiPadC(phiPad);
track->setDeltaThetaSidcPadC(dtheta);
track->setDeltaPhiSidcPadC(dphi);
track->setTrackMask(sidcSegment->getTrackMask()|CTrack::hasPadCMatch);
sidcPadTracks->append(track);
}
}
}
}
return True;
}
int CPionTrackingStrategy::returnIndexOfMaxVector()
{
int maxIndex = 0;
int maxContent = 0;
int vectorSize = setup->getVectorSize();
for ( int k = 0; k < vectorSize ; k++) {
if(radiusVector(k) > maxContent){
maxContent = radiusVector(k);
maxIndex = k;
}
}
return maxIndex;
}
CRingCandidate* CPionTrackingStrategy::searchForCandidate(CRich* rich, CTrack* track)
{
if (rich->getHits().isEmpty()) return 0;
int nhitsForChernov;
int indexOfMax, nbin ;
int intVectorIndex;
double fitRadius,relRad;
double factor = 1.5;
double maxRange = 2.*rich->getSetup()->getAsymptoticRadius();
CRichLikeHit *hit;
CRingCandidate *candidate;
CRichPolarCoord thisRichPolarCoord;
CRichPadCoord candCenter,thisRichPadCoord;
radiusVector = 0; // reset the Radius Vector
xCoord = 0; // reset for every Rich and tracks
yCoord = 0;
binPos = 0;
numberOfHits = 0;
CAngle dphiRich = (track->getDeltaPhiSidcPadC()*rich->getSetup()->getDeltaPhiFactor())*
setup->getFieldFlag();
CAngle phiRich = track->getPhiSidc() - dphiRich;
double thetaRich = track->getThetaPadC(); // prefer the theta of pad Chamber
thisRichPolarCoord.setTheta(thetaRich);
thisRichPolarCoord.setPhi(phiRich);
rich->transform(thisRichPadCoord,thisRichPolarCoord); // Sollposition in Rich
double x0 = thisRichPadCoord.getX(); // predictor
double y0 = thisRichPadCoord.getY();
CRichLikeHit richPredictorHit(x0 - maxRange, 0, 0);
int k = rich->getHits().getIndexOfBestMatch(richPredictorHit);
double dx = -maxRange;
double dy, distance;
while (dx < maxRange && k < rich->getHits().length()) {
hit = rich->getHits()(k);
dx = hit->getX()-x0;
dy = hit->getY()-y0;
distance =sqrt(dx*dx + dy*dy);
if(distance < maxRange) {
relRad = distance/rich->getSetup()->getAsymptoticRadius();
intVectorIndex = (int)(10*(relRad/factor)*(relRad/factor));
intVectorIndex += 1; // correct position
xCoord(numberOfHits) = hit->getX();
yCoord(numberOfHits) = hit->getY();
binPos(numberOfHits) = intVectorIndex;
if (intVectorIndex < setup->getVectorSize() - 1 ) {
radiusVector(intVectorIndex) += 1; // increase by one
numberOfHits++;
}
}
k++;
}
//
// use the hit that are in the max bin (+ o - 1) for Chernov fit
//
indexOfMax = returnIndexOfMaxVector();
nhitsForChernov = 0;
fit.resetAll();
double x1, y1;
for (int d1 = 0; d1 < numberOfHits; d1++) {
x1 = xCoord(d1);
y1 = yCoord(d1);
nbin = binPos(d1);
if( nbin >= indexOfMax - 1 && nbin <= indexOfMax + 1 ) {
nhitsForChernov += 1;
fit.addX(x1);
fit.addY(y1);
}
}
if (nhitsForChernov < setup->getMinHits()) return 0;
if (!fit.chernov()) return 0; // chernov fit was ok?
candCenter.setX(fit.getXCenter());
candCenter.setY(fit.getYCenter());
fitRadius = fit.getRadius();
candidate = new CRingCandidate(candCenter, fitRadius, CRing::Pion, 0., 0., numberOfHits);
return candidate;
}
CBoolean CPionTrackingStrategy::makeRichTracksAndLinkToSidcPadTracks()
{
//
// anything to do?
//
if (!sidcTracks || !richTracks || !sidcPadTracks || !tracks) {
cerr << "CPionTrackingStrategy::makeRichTracksAndLinkToSidcPadTracks():" << endl;
cerr << "\tERROR:" << endl;
cerr << "\t One of the internal track lists does not exist" << endl;
cerr << "\tNo further action." << endl;
return False;
}
if (sidcPadTracks->length() == 0) return True;
CTrack *sidcPadTrack;
CTrack *completeTrack;
CSortedListIterator<CTrack> spiter(*sidcPadTracks);
CRingCandidate *candidate1, *candidate2;
CFittedRing *ring1, *ring2;
while (sidcPadTrack = spiter()) {
candidate1 = candidate2 = 0;
ring1 = ring2 = 0;
completeTrack = 0;
if(candidate1 = searchForCandidate((CRich1*)rich1,sidcPadTrack)) {
if(ring1 = makePionFreeRadiusFit((CRich1*)rich1, candidate1)) {
if(candidate2 = searchForCandidate((CRich2*)rich2, sidcPadTrack)){
if(ring2 = makePionFreeRadiusFit((CRich2*)rich2, candidate2)) {
if(completeTrack = isTrueTrack(ring1, ring2, sidcPadTrack)) {
tracks->append(completeTrack);
}
}
}
}
}
delete candidate1;
delete candidate2;
if (!completeTrack) {
delete ring1;
delete ring2;
}
}
return True;
}
CFittedRing* CPionTrackingStrategy::makePionFreeRadiusFit(CRich* rich, CRingCandidate* candidate)
{
int hitDensity;
CFittedRing *ring;
CRichPadCoord padCoord;
CRichPolarCoord polarCoord;
double distance;
double range = 1.25*rich->getSetup()->getAsymptoticRadius();
//
// Accumulate all hits around ring candidate,
// feed them into the fitter.
//
hitDensity = numberOfHits;
fit.resetAll();
for (int k = 0; k < numberOfHits; k++) {
distance= sqrt((xCoord(k) - candidate->getCenter().getX())*
(xCoord(k) - candidate->getCenter().getX()) +
(yCoord(k) - candidate->getCenter().getY())*
(yCoord(k) - candidate->getCenter().getY()));
if (distance < range) {
fit.addX(xCoord(k));
fit.addY(yCoord(k));
fit.addWeight(1.);
}
}
if (fit.getNFitPoints() <= setup->getMinHits()) return 0;
//
// Pass start values for center and radius
// and number of iterations.
//
fit.setXCenter(candidate->getCenter().getX());
fit.setYCenter(candidate->getCenter().getY());
fit.setRadius(candidate->getRadius());
//
// Fit the candidate, perform some simple checks and
// store the new ring in the referring list.
//
if (fit.robustFreeRadius()) {
padCoord = CRichPadCoord(fit.getXCenter(), fit.getYCenter());
rich->transform(polarCoord, padCoord);
//
// store new ring and return it if it is good else return 0 pointer
//
ring = new CFittedRing();
ring->setType(CRing::Pion);
ring->setCenter(padCoord, polarCoord);
ring->setRadius(fit.getRadius());
ring->setVariance(fit.getVariance());
ring->setHitDensity(hitDensity);
return ring;
}
else {
return 0;
}
}
CTrack* CPionTrackingStrategy::isTrueTrack(CFittedRing* ring1,CFittedRing* ring2, CTrack* SPTrack)
{
int k, signOfTrack;
int nhit1,nhit2;
CAngle phi1,phi2;
double distance;
double theta1,theta2;
double dtheta12, dphi12, waitedDeltaPhi;
double RelRadius1,RelRadius2;
CBoolean isTooClose = False;
CTrack *isThisTrack;
CTrack *newTrack;
CTrack *newSegment;
CDeflection deflection;
CPadCPadCoord predictor;
RelRadius1 = ring1->getRadius()/setup->getRich1AsymptoticRadius();
RelRadius2 = ring2->getRadius()/setup->getRich2AsymptoticRadius();
phi1 = ring1->getPolarCenter().getPhi();
phi2 = ring2->getPolarCenter().getPhi();
theta1 = ring1->getPolarCenter().getTheta();
theta2 = ring2->getPolarCenter().getTheta();
double thetaMeanNoField = (theta1+theta2)/2.;
CAngle phiMeanNoField = average(phi1,phi2);
if(!setup->getFieldFlag()) { // if No Field => take the mean values !!!!
phi1 = phi2 = phiMeanNoField;
theta1 = theta2 = thetaMeanNoField;
}
deflection.calculatePCPrediction(theta1, phi1, theta2, phi2, (CPadChamber&)*padChamber,
(CVertex&)*vertex);
predictor = deflection.getPadcPredictor(); // in pad Coord
dphi12 = (double)(phi2 - phi1);
dtheta12 = theta2 - theta1;
signOfTrack = (dphi12 < 0) ? -1 : 1 ;
double tmpRadius = (RelRadius2 + RelRadius1)/2.;
if (tmpRadius > 1.) tmpRadius = 1;
waitedDeltaPhi = (0.03214*sqrt(1 - tmpRadius*tmpRadius))*setup->getFieldFlag(); // theoric value for deflection
//
// decide if the two ring can be a pion track => apply cut
//
if(fabs(fabs(dphi12)- waitedDeltaPhi) < setup->getWaitedPhiWindow() &&
fabs((RelRadius2/RelRadius1) - 0.9656) < setup->getRelRadiusAllowedDifference() &&
fabs(dtheta12)< setup->getThetaWindow() &&
RelRadius1 < setup->getCutRelRadius1() &&
RelRadius2 < setup->getCutRelRadius2()) {
calculateNumberOfHits((CRich1*)rich1,ring1);
calculateNumberOfHits((CRich2*)rich2,ring2);
ring1->evalChi2(ring1->getRadius());
ring2->evalChi2(ring2->getRadius());
nhit1 = ring1->getNumberOfHits();
nhit2 = ring2->getNumberOfHits();
if(nhit1 < setup->getMinHits() || nhit2 < setup->getMinHits()) return 0;
//
// check if already in the list or it is too close to another one
//
k = 0;
while (k < tracks->length()) {
isThisTrack = (*tracks)(k);
distance = abs(ring1->getPadCenter() - isThisTrack->getClosestRich1Ring()->getPadCenter());
if(distance < (setup->getRich1AsymptoticRadius())/2.) {
if(dtheta12 < isThisTrack->getDeltaThetaRich12()){
delete tracks->removeAt(k);
}
else{
isTooClose = True;
k++;
}
}
else
k++;
}
if(isTooClose) return 0;
SPTrack = findBestPadCMatch(SPTrack,predictor,theta2, phi2); // best PadChamber Match
if((SPTrack->getTrackMask()&CTrack::hasPadCMatch) != CTrack::hasPadCMatch) return 0;
newSegment = new CTrack;
newSegment->setTrackMask(CTrack::hasRich1Match | CTrack::hasRich2Match);
newSegment->setClosestRich1Ring(ring1);
newSegment->setClosestRich2Ring(ring2);
newSegment->setDeltaThetaRich12(dtheta12);
newSegment->setDeltaPhiRich12(dphi12);
if(setup->getFieldFlag()) // it makes sense only with field
newSegment->setType(signOfTrack > 0 ? CTrack::PiPlus : CTrack::PiMinus);
richTracks->append(newSegment);
newTrack = new CTrack; // place for a new track
*newTrack = *SPTrack; // fill with old track value (Sidc - Pad)
newTrack->setTrackMask(SPTrack->getTrackMask()|newSegment->getTrackMask());
newTrack->setClosestRich1Ring(ring1);
newTrack->setClosestRich2Ring(ring2);
newTrack->setDeltaThetaRich12(dtheta12);
newTrack->setDeltaPhiRich12(dphi12);
newTrack->setDeltaThetaSidcRich(SPTrack->getThetaSidc() -
ring1->getPolarCenter().getTheta());
newTrack->setDeltaPhiSidcRich((CAngle) SPTrack->getPhiSidc() -
ring1->getPolarCenter().getPhi());
if(setup->getFieldFlag()) // it makes sense only with field
newTrack->setType(signOfTrack > 0 ? CTrack::PiPlus : CTrack::PiMinus) ;
return newTrack;
}
else
return 0;
}
CBoolean CPionTrackingStrategy::calculateNumberOfHits(CRich* rich, CFittedRing* ring)
{
float distance;
CRichLikeHit *hit;
CRichPadCoord center;
center = ring->getPadCenter();
for (int i=0; i < rich->getHits().length(); i++) {
hit = rich->getHits()(i);
distance = sqrt((hit->getX() - center.getX())*(hit->getX() - center.getX()) +
(hit->getY() - center.getY())*(hit->getY() - center.getY())) -
ring->getRadius();
distance = fabs(distance);
if(distance < rich->getSetup()->getMaxHitRingDistance()){
ring->addHit(hit); // add hit to ring
}
}
return True;
}