CRich.C
//-----------------------------------------------------------------------------
// $Header: /asis/offline/ceres/cool/project/RCS/CRich.C,v 2.18 1997/10/16 09:48:26 messer Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class CRich.
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include <stdio.h>
#include <math.h>
#include <limits.h>
#include "CRich.h"
#include "CHoughWorkArray.h"
#include "CRingFitter.h"
#include "CListIterator.h"
#include "CSortedListIterator.h"
#include "CTrack.h"
#include "CRandom.h"
CRich::CRich()
{
ringList = 0;
candidateList = 0;
houghMask = 0;
}
CRich::~CRich() { delete houghMask; }
float CRich::getLeft() const { return 0.5; }
float CRich::getRight() const { return setup->getNPads()+0.5; };
float CRich::getBottom() const { return 0.5; }
float CRich::getTop() const { return setup->getNPads()+0.5; };
CBoolean CRich::makeClustersAndDoCleanup()
{
if (!makeClusters(5)) return False;
if (!removeBigClusters()) return False;
if (!makeClusters(1)) return False;
removeNoiseClusters();
if (!splitClusters()) return False;
if (!removeSmallClusters()) return False;
return True;
}
void CRich::removeNoiseClusters()
{
int i, j, ioff, iend;
for (i=0; i<clusters.length(); i++) {
if (clusters(i)->getNPads() < setup->getClusterMinPads()) {
ioff = clusters(i)->getFirstPad();
iend = ioff+clusters(i)->getNPads();
for ( j=ioff;j<iend;j++ ) {
pads(j)->setAmp(-fabs(pads(j)->getAmp())); // mark pads
}
clusters(i)->setNPads(-abs(clusters(i)->getNPads()));
}
}
}
CBoolean CRich::removeSmallClusters()
{
if ( clusters.length() <= 0) return False;
int i, j, ioff, iend;
for ( i=0; i<clusters.length(); i++) {
if (!isGoodSmallCluster(*clusters(i))) { // mark clusters found by algorithm
ioff = clusters(i)->getFirstPad();
iend = ioff+clusters(i)->getNPads();
for ( j=ioff;j<iend;j++ ) {
pads(j)->setAmp(-fabs(pads(j)->getAmp())); // mark pads
}
clusters(i)->setNPads(-abs(clusters(i)->getNPads()));
}
}
return True;
}
CBoolean CRich::isGoodSmallCluster(CRichLikeCluster& cluster)
{
if ( cluster.getNPads() < setup->getClusterMinPads() ) return False; // too small
int ioff = cluster.getFirstPad();
int npads = cluster.getNPads();
int iend = ioff+npads;
int nhigh = 0;
int nsat = 0;
//
// collect all properties
//
for (int i=ioff; i<iend; i++) {
if ( pads(i)->getAmp() >= setup->getClusterHighThreshold() ) nhigh++;
if ( pads(i)->getAmp() >= setup->getClusterSatPadAmp() ) nsat++;
}
//
// apply all cuts
//
if ( cluster.getNPads() > setup->getClusterMaxPads() ||
nsat > setup->getClusterMaxSatPads() ||
nhigh > setup->getClusterMaxPadsAboveThreshold() ||
cluster.getNLocMax()==0 ) {
return False;
} else {
return True;
}
}
CBoolean CRich::removeBigClusters()
{
if ( clusters.length() <= 0) return False; // nothing to do
int i, j, ioff, iend;
for ( i=0; i<clusters.length(); i++) {
if (!isGoodBigCluster(*clusters(i))) { // mark clusters found by algorithm
ioff = clusters(i)->getFirstPad();
iend = ioff+clusters(i)->getNPads();
for ( j=ioff;j<iend;j++ ) {
pads(j)->setAmp(-fabs(pads(j)->getAmp())); // mark pads
}
clusters(i)->setNPads(-abs(clusters(i)->getNPads()));
}
}
return True;
}
CBoolean CRich::isGoodBigCluster(CRichLikeCluster& cluster)
{
//
// reject large circular ionization cluster
//
int ioff = cluster.getFirstPad();
int npads = cluster.getNPads();
int iend = ioff+npads;
int nsat = 0;
int lnpads = 1; // to get the same results as in the old lib...
int lminX = 300; int lminY = 300; int lmaxX = 0; int lmaxY = 0;
int minX = 300; int minY = 300; int maxX = 0; int maxY = 0;
float xav = 0; float yav = 0; float xsq = 0; float ysq = 0;
float lAmp = 0;
// collect all properties..
for (int i=ioff; i<iend; i++) {
if ( pads(i)->getAmp()>= setup->getClusterSatPadAmp() ) nsat++;
if ( pads(i)->getAmp() >= setup->getClusterLowThreshold() ) {
lnpads++;
lAmp+=pads(i)->getAmp();
if ( pads(i)->getX()<lminX ) lminX = pads(i)->getX();
if ( pads(i)->getX()>lmaxX ) lmaxX = pads(i)->getX();
if ( pads(i)->getY()<lminY ) lminY = pads(i)->getY();
if ( pads(i)->getY()>lmaxY ) lmaxY = pads(i)->getY();
}
if ( pads(i)->getX()<minX ) minX = pads(i)->getX();
if ( pads(i)->getX()>maxX ) maxX = pads(i)->getX();
if ( pads(i)->getY()<minY ) minY = pads(i)->getY();
if ( pads(i)->getY()>maxY ) maxY = pads(i)->getY();
xav += pads(i)->getX()*pads(i)->getAmp();
yav += pads(i)->getY()*pads(i)->getAmp();
xsq += pads(i)->getX()*pads(i)->getX()*pads(i)->getAmp();
ysq += pads(i)->getY()*pads(i)->getY()*pads(i)->getAmp();
}
int load = int( (float(lAmp) / float(lnpads))+0.5 ); // nint!!!!
int xSize = maxX - minX + 1;
int ySize = maxY - minY + 1;
int clusterWidth = 0; int clusterLength = 0;
int lClusterWidth = 0; int lClusterLength = 0;
if (xSize <= ySize) {
clusterWidth = xSize;
clusterLength = ySize;
} else {
clusterWidth = ySize;
clusterLength = xSize;
}
int lxSize = lmaxX - lminX + 1;
int lySize = lmaxY - lminY + 1;
if (lxSize <= lySize) {
lClusterWidth = lxSize;
lClusterLength = lySize;
} else {
lClusterWidth = lySize;
lClusterLength = lxSize;
}
float ratio = 0; float occup = 0;
if ( nsat>setup->getClusterMaxSat2() ) {
if ( lClusterWidth>setup->getClusterMinSmallSize() ) {
ratio = float(lClusterWidth) / float(lClusterLength); // circular shape
occup = float(lnpads) / float(lxSize) / float(lySize); // large box occup
if (ratio > setup->getClusterMaxRatio() &&
occup > setup->getClusterMaxOccup() ||
load > setup->getClusterMaxLoad2() ) return False;
}
}
//
// reject large saturated ionization tracks
//
float rmsSize = xsq/cluster.getAmpSum() + ysq/cluster.getAmpSum()
- (xav/cluster.getAmpSum())*(xav/cluster.getAmpSum())
- (yav/cluster.getAmpSum())*(yav/cluster.getAmpSum());
if ( clusterWidth>setup->getClusterSatTrackMaxWidth() &&
nsat>setup->getClusterSatTrackMaxSatPads() &&
rmsSize>setup->getClusterMaxRmsSize() &&
load>setup->getClusterMaxLoad() ) return False;
//
// reject very big clusters
//
if ( clusterWidth >= setup->getClusterBigClusterWidth() &&
( nsat >= setup->getClusterMaxSat() ||
npads >= setup->getClusterBigClusterMaxPads() ) ) return False;
//
// reject a line of pads
//
switch (clusterWidth) {
case 1:
if ( clusterLength >= 4 &&
cluster.getRms2() <= setup->getClusterPadLineRms2()) return False;
case 2:
if ( clusterLength >= 8 &&
cluster.getRms2() <= setup->getClusterPadLineRms2()) return False;
case 3:
if ( clusterLength >= 12 &&
cluster.getRms2() <= setup->getClusterPadLineRms2()) return False;
}
//
// reject bad small clusters where all pads are red
//
if ( (cluster.getAmpSum()/float(npads)) > setup->getClusterMeanAmp() ) return False;
if ( (float(nsat)/float(npads)) > setup->getClusterSatPadRatio() ) return False;
//
// reject flat large zones
//
if ( clusterWidth >= 4 &&
cluster.getRms2() <= setup->getClusterMinRms2()) return False;
//
// cluster survived all cuts
//
return True;
}
//
// Draw a Rich like detector.
// The following options are recognized:
// m draw module borders
// s draw spokes
// a draw acceptance
// o show pads from clusters only ('o' and 'c' are mutualy exclusive)
// c draw clusters (colored according to cluster id)
// h draw hits
// n print numeric values of pad amplitudes
// k draw candidates
// v mark sensitive area
// e draw electron rings
// t draw tracks from attached track list
// x draw extended track information (same as 't' but also Sidc segments)
// r draw removed pads only ( excludes most of the other pad options )
// i draw cluster ids
// 0 don't draw pads (except for r option)
// p pion tracks (Green - Rich1 , Yellow - Rich2)
//
void CRich::draw(const char* copt)
{
int i;
CString options(copt);
char number[16];
color(Gray);
clear();
//
// Define beam axis
//
double xcenter = setup->getNPads()/2.+0.5;
double ycenter = xcenter;
//
// First draw modules in black.
// Remember (x,y) in CModule is upper left pad.
//
color(Black);
polyfill(True);
for (i=0; i<modules.length(); i++) {
rect(modules(i)->getX()-0.5,
modules(i)->getY()-modules(i)->getHeight()+0.5,
modules(i)->getX()+modules(i)->getWidth()-0.5,
modules(i)->getY()+0.5);
}
//
// Draw module borders in blue. In addition the index
// of the module is plotted in the upper left corner.
//
if (options.first('m') != RW_NPOS) {
color(Blue);
polyfill(False);
font("futura.l");
leftjustify();
topjustify();
textsize(1, 1);
for (i=0; i<modules.length(); i++) {
rect(modules(i)->getX()-0.5,
modules(i)->getY()-modules(i)->getHeight()+0.5,
modules(i)->getX()+modules(i)->getWidth()-0.5,
modules(i)->getY()+0.5);
move2(modules(i)->getX(), modules(i)->getY());
sprintf(number, "%d\0", i);
drawstr(number);
}
}
//
// Draw spokes
//
float points[5][2];
if (options.first('s') != RW_NPOS) {
if (psMode)
color(White);
else
color(Gray);
polyfill(True);
double halfWidth = setup->getSpokeWidth()/2.;
double phi = setup->getFirstSpokeAngle()*Pi/180.;
double dphi = TwoPi/setup->getNSpokes();
double outerRadius = setup->getNPads()/2.+1;
double innerRadius = outerRadius/4.;
for (i=0; i<setup->getNSpokes(); i++) {
points[0][0] = outerRadius * cos((phi + atan(halfWidth/outerRadius))) + xcenter;
points[0][1] = outerRadius * sin((phi + atan(halfWidth/outerRadius))) + ycenter;
points[1][0] = innerRadius * cos((phi + atan(halfWidth/innerRadius))) + xcenter;
points[1][1] = innerRadius * sin((phi + atan(halfWidth/innerRadius))) + ycenter;
points[2][0] = innerRadius * cos((phi - atan(halfWidth/innerRadius))) + xcenter;
points[2][1] = innerRadius * sin((phi - atan(halfWidth/innerRadius))) + ycenter;
points[3][0] = outerRadius * cos((phi - atan(halfWidth/outerRadius))) + xcenter;
points[3][1] = outerRadius * sin((phi - atan(halfWidth/outerRadius))) + ycenter;
points[4][0] = points[0][0]; points[4][1] = points[0][1];
phi += dphi;
poly2(5, points);
}
}
//
// Show acceptance
//
if (options.first('a') != RW_NPOS) {
if (psMode)
color(Gray);
else
color(White);
polyfill(False);
circle(xcenter, ycenter, setup->getInnerAcceptance());
circle(xcenter, ycenter, setup->getOuterAcceptance());
}
//
// Draw pads and/or clusters
//
CPad *pad;
CBoolean showNumbers = False;
CBoolean showClusterIds = False;
int first, j;
if (options.first('n') != RW_NPOS || options.first('i') != RW_NPOS) {
if (options.first('n') != RW_NPOS ) {
showNumbers = True;
showClusterIds = False;
}
if (options.first('i') != RW_NPOS) {
showClusterIds = True;
showNumbers =False;
}
font("futura.l");
centertext(True);
textsize(0.5, 0.5);
polyfill(False);
}
else
polyfill(True);
if (options.first('o') != RW_NPOS || options.first('c') != RW_NPOS) {
//
// Pads in cluster only
//
CBoolean uniqCluster = options.first('c') != RW_NPOS;
for (i=0; i<clusters.length(); i++) {
if ( clusters(i)->getNPads() < 0 ) continue; // use valid clusters only
first = clusters(i)->getFirstPad();
for (j=first; j<first+clusters(i)->getNPads(); j++) {
pad = pads(j);
if (pad->getAmp() > 0) {
if (uniqCluster) {
//
// Color according to cluster id
//
switch (i%5) {
case 0:
color(Blue);
break;
case 1:
color(Green);
break;
case 2:
color(Yellow);
break;
case 3:
color(Cyan);
break;
case 4:
color(White);
break;
}
if (pad->getIsLocalMax()) // mark local maximum (red)
color(Red);
}
else {
//
// Color according to amplitude
//
color(colors[(int) pad->getAmp()]);
}
rect(pad->getX()-0.5, pad->getY()-0.5,
pad->getX()+0.5, pad->getY()+0.5);
if (showClusterIds) {
sprintf(number, "%d\0", i);
move2(pad->getX(), pad->getY());
drawstr(number);
} else if (showNumbers) {
sprintf(number, "%d\0", (int)pad->getAmp());
move2(pad->getX(), pad->getY());
drawstr(number);
}
}
}
}
}
else if (options.first('r') != RW_NPOS) {
//
// All removed pads (optional as numbers only)
//
for (i=0; i<pads.length(); i++) {
pad = pads(i);
if (pad->getAmp() < 0) {
color(colors[(int) fabs(pad->getAmp())]);
rect(pad->getX()-0.5, pad->getY()-0.5,
pad->getX()+0.5, pad->getY()+0.5);
if (showNumbers) {
sprintf(number, "%d\0", (int)fabs(pad->getAmp()));
move2(pad->getX(), pad->getY());
drawstr(number);
}
}
}
}
else if (options.first('0') == RW_NPOS) { // omit drawing pads for '0' option
//
// All pads (optional as numbers only)
//
for (i=0; i<pads.length(); i++) {
pad = pads(i);
if (pad->getAmp() > 0) {
color(colors[(int) pad->getAmp()]);
rect(pad->getX()-0.5, pad->getY()-0.5,
pad->getX()+0.5, pad->getY()+0.5);
if (showNumbers) {
sprintf(number, "%d\0", (int)pad->getAmp());
move2(pad->getX(), pad->getY());
drawstr(number);
}
}
}
}
//
// Draw hits as white crosses with a small rectangle in the middle.
// If the hits is part of a ring, the hit is marked with a small circle.
//
CBoolean markIfPartOfRing = (options.first('e') != RW_NPOS && ringList);
const float crossWidth = 3;
const float markerWidth = 0.1;
if (options.first('h') != RW_NPOS) {
polyfill(False);
CRichLikeHit *hit;
for (i=0; i<hits.length(); i++) {
hit = hits(i);
if (markIfPartOfRing && hit->getRing())
circle(hit->getX(), hit->getY(), crossWidth/2);
color(White);
move2(hit->getX()-crossWidth/2, hit->getY()-crossWidth/2);
rdraw2(crossWidth, crossWidth);
rmove2(-crossWidth, 0);
rdraw2(crossWidth, -crossWidth);
rect(hit->getX()-markerWidth, hit->getY()-markerWidth,
hit->getX()+markerWidth, hit->getY()+markerWidth);
}
}
//
// Draw candidates as yellow dashed circles
//
if (options.first('k') != RW_NPOS && candidateList != 0) {
CListIterator<CRingCandidate> citer(*candidateList);
polyfill(False);
color(Yellow);
linestyle("10");
setdash(0.015);
while (citer())
circle(citer.key()->getCenter().getX(),
citer.key()->getCenter().getY(),
citer.key()->getRadius());
linestyle("1");
}
//
// Draw electron rings as red circles, mark the center with a red cross
//
CFittedRing *ring;
if (options.first('e') != RW_NPOS && ringList) {
CListIterator<CFittedRing> riter(*ringList);
polyfill(False);
color(Red);
while (ring = riter()) {
circle(ring->getPadCenter().getX(), ring->getPadCenter().getY(), ring->getRadius());
move2(ring->getPadCenter().getX()-crossWidth/2, ring->getPadCenter().getY()-crossWidth/2);
rdraw2(crossWidth, crossWidth);
rmove2(-crossWidth, 0);
rdraw2(crossWidth, -crossWidth);
}
}
//
// Mark sensitive area
//
if (options.first('v') != RW_NPOS) {
polyfill(False);
color(White);
CRichLikeLookupItem *item;
for (i=0; i<lookupItem.length(); i++) {
item = lookupItem(i);
if (item->getIsSensitive()) {
rect(item->getX()-0.5, item->getY()-0.5,
item->getX()+0.5, item->getY()+0.5);
}
}
}
//
// Draw tracks from attached track list.
// Rings of tracks are marked by a number referring to the index
// of the track in the list. Rings belonging to the same track
// are connected by a line.
//
CTrack *track;
CRichPadCoord center, center1, center2;
CRichPadCoord center1n, center2n;
if ((options.first('t') != RW_NPOS || options.first('x') != RW_NPOS) && trackList) {
font("futura.l");
centertext(True);
polyfill(False);
textsize(2*crossWidth, 2*crossWidth);
CSortedListIterator<CTrack> titer(*trackList);
while (track = titer()) {
if(track->getTrackMask() == (CTrack::hasRich1Match|CTrack::hasRich2Match|
CTrack::hasSidc1Match|CTrack::hasSidc2Match|
CTrack::hasVertex))
color(Cyan);
else
color(Blue);
sprintf(number, "%d\0", titer.pos());
//
// Text for Rich1 match if any
// Recolor referring ring according to track mask
//
if (track->getTrackMask()&CTrack::hasRich1Match) {
transform(center1, track->getClosestRich1Ring()->getPolarCenter());
if (id == Rich1)
circle(center1.getX(), center1.getY(),
track->getClosestRich1Ring()->getRadius());
move2(center1.getX(), center1.getY());
drawstr(number);
}
//
// Text for Rich2 match if any
// Recolor referring ring according to track mask
//
if (track->getTrackMask()&CTrack::hasRich2Match) {
transform(center2, track->getClosestRich2Ring()->getPolarCenter());
if (id == Rich2)
circle(center2.getX(), center2.getY(),
track->getClosestRich2Ring()->getRadius());
move2(center2.getX(), center2.getY());
drawstr(number);
}
//
// Visualize match by straight line
//
if (track->getTrackMask()&CTrack::hasRich1Match &&
track->getTrackMask()&CTrack::hasRich2Match) {
move2(center1.getX(), center1.getY());
draw2(center2.getX(), center2.getY());
}
}
//
// Visualize Sidc segment
// only Sidc1 match: triangle
// only Sidc2 match: cross
// Sidc1 & Sidc2 match: star
//
if (options.first('x') != RW_NPOS) {
font("markers");
color(White);
textsize(crossWidth, crossWidth);
titer.reset();
while (track = titer()) {
if (track->getTrackMask()&CTrack::hasVertex) {
transform(center, CRichPolarCoord(track->getThetaSidc(),
track->getPhiSidc()));
move2(center.getX(), center.getY());
if (track->getTrackMask()&CTrack::hasSidc1Match &&
track->getTrackMask()&CTrack::hasSidc2Match)
drawchar('e');
else if (track->getTrackMask()&CTrack::hasSidc1Match)
drawchar('c');
else if (track->getTrackMask()&CTrack::hasSidc2Match)
drawchar('f');
}
}
}
}
//
// Draw the pionTracks on the Riches.
// For having some facilities during the cut procedure, on every Riches are drawed
// two rings, the Green-Rich1 ring and the Yellow_rich2 ring).
//
if ((options.first('p') != RW_NPOS) && trackList) {
font("futura.l");
centertext(True);
polyfill(False);
textsize(2*crossWidth, 2*crossWidth);
CSortedListIterator<CTrack> titer(*trackList);
cout << " number of Pion tracks :" << trackList->length() << endl;
while (track = titer()) {
sprintf(number, "%d\0", i);
//
// Transform from Polar Coordinate to Rich Coordinates in Pad Unit
//
transform(center1, track->getClosestRich1Ring()->getPolarCenter());
transform(center2, track->getClosestRich2Ring()->getPolarCenter());
color(Green); // Rich1
circle(center1.getX(), center1.getY(),
track->getClosestRich1Ring()->getRadius()*setup->getRadiusFactor1());
move2(center1.getX()-crossWidth/2, center1.getY()-crossWidth/2);
rdraw2(crossWidth, crossWidth);
rmove2(-crossWidth, 0);
rdraw2(crossWidth, -crossWidth);
color(Yellow); // Rich2
circle(center2.getX(), center2.getY(),
track->getClosestRich2Ring()->getRadius()*setup->getRadiusFactor2());
move2(center2.getX()-crossWidth/2, center2.getY()-crossWidth/2);
rdraw2(crossWidth, crossWidth);
rmove2(-crossWidth, 0);
rdraw2(crossWidth, -crossWidth);
}
}
}
CBoolean CRich::makeHoughCandidates(int minHough)
{
//
// temporary data types
//
class houghPad {
public:
houghPad(const int size)
{
x = new short[size];
y = new short[size];
length = size;
}
~houghPad() { delete [] x; delete [] y; }
public:
int length;
short *x, *y;
};
//
// anything to do?
//
if (hits.length() <= 0) return False;
//
// make new candidate list for this event
//
delete candidateList;
candidateList = new CList<CRingCandidate>;
//
// parameters for cuts
//
int minAmp;
if (minHough > -1)
minAmp = minHough;
else
minAmp = setup->getHoughMinAmp();
double radius = setup->getAsymptoticRadius(); // asymptotic radius for electrons
//
// Create the hough working arrays and the corresponding list
// (with maximum possible size).
//
static CHoughWorkArray hough;
static CHoughWorkArray smooth;
int maxsize = (hough.maxIndex-hough.minIndex)*(hough.maxIndex-hough.minIndex);
int thisSize = houghMask->length * hits.length();
int listSize = (thisSize > maxsize ? thisSize : maxsize);
houghPad *houghList = new houghPad(listSize);
//
// first Hough transformation
//
int nhough = 0;
int xhit, yhit, ix, iy, i, j;
hough.zeroXY();
for(i=0; i<hits.length(); i++ ) {
xhit = int(hits(i)->getX() + .5) - hough.minIndex; // transform into hough
yhit = int(hits(i)->getY() + .5) - hough.minIndex; // -array-coordinates
for(j=0; j<houghMask->length; j++ ) {
ix = xhit + houghMask->x[j];
iy = yhit + houghMask->y[j];
if (hough.xy[ix][iy] == 0) {
houghList->x[nhough] = ix;
houghList->y[nhough] = iy;
nhough++;
}
hough.xy[ix][iy] += short(hits(i)->getAmp()+0.5);
}
}
//
// sum over 3x3 array i.e. smooth hough array
//
int ixx, iyy;
smooth.zeroXY();
for (i=0; i<nhough; i++) {
ix = houghList->x[i];
iy = houghList->y[i];
for (ixx=ix-1; ixx<=ix+1; ixx++)
for (iyy=iy-1; iyy<=iy+1; iyy++)
smooth.xy[ix][iy] += hough.xy[ixx][iyy];
smooth.xy[ix][iy] /= 10;
}
//
// (fast) zero hough-working-array
//
for (i=0; i<nhough; i++)
hough.xy[houghList->x[i]][houghList->y[i]] = 0;
//
// second Hough transformation
// fill houghlist with points > minAmp
//
nhough = 0;
int iringsum, newAmplitude;
for( i=0; i<hits.length(); i++ ) {
xhit = int(hits(i)->getX() + .5) - hough.minIndex; // transform into hough
yhit = int(hits(i)->getY() + .5) - hough.minIndex; // -array-coordinates
iringsum = 0;
for( j=0; j<houghMask->length; j++ )
iringsum += smooth.xy[xhit+houghMask->x[j]][yhit+houghMask->y[j]];
if (iringsum > 0) {
for( j=0; j<houghMask->length; j++ ) {
ix = xhit + houghMask->x[j];
iy = yhit + houghMask->y[j];
newAmplitude = hough.xy[ix][iy] + (smooth.xy[ix][iy]*1000)/iringsum;
if (newAmplitude >= minAmp && hough.xy[ix][iy] < minAmp) {
houghList->x[nhough] = ix;
houghList->y[nhough] = iy;
nhough++;
}
hough.xy[ix][iy] = newAmplitude;
}
}
}
//
// Fill candidate list looking for local maxima above minAmp in hough array
//
CRichPadCoord candCenter;
float hough1Amplitude, hough2Amplitude;
RWTPtrOrderedVector<CRingCandidate> newCandidates;
for (i=0; i<nhough; i++) {
ix = houghList->x[i];
iy = houghList->y[i];
if (hough.isLocalMax(ix, iy)) {
candCenter.setX(double(ix + hough.minIndex));
candCenter.setY(double(iy + hough.minIndex));
hough1Amplitude = float(smooth.xy[ix][iy]);
hough2Amplitude = float(hough.xy[ix][iy]);
newCandidates.insert(new CRingCandidate(candCenter, radius, CRing::Electron,
hough1Amplitude, hough2Amplitude, 0));
}
}
//
// Fit rings (on hits) using the fast Chernov algorithm.
// Remove unreasonable looking rings.
//
CRingFitter fit;
double minDist2 = setup->getHoughMinDistance() * setup->getHoughMinDistance();
int nhits;
double xCand, yCand, xHit, yHit, hitRadius;
double halfChernovWidth = setup->getHoughChernovWidth()/2;
double fitRadius;
CRichPadCoord deltaCenter, fitCenter;
double deltaCenter2;
for (i=0; i<newCandidates.length(); i++) {
fit.resetAll(); // prepare for next fit
xCand = newCandidates(i)->getCenter().getX();
yCand = newCandidates(i)->getCenter().getY();
nhits = 0;
for (j=0; j<hits.length(); j++) {
xHit = hits(j)->getX();
yHit = hits(j)->getY();
hitRadius = sqrt((xHit-xCand)*(xHit-xCand) + (yHit-yCand)*(yHit-yCand));
if (fabs(hitRadius - radius) <= halfChernovWidth) {
nhits++;
fit.addX(xHit);
fit.addY(yHit);
}
}
if (nhits < setup->getHoughMinFitHit()) continue; // enough hits for later fit?
newCandidates(i)->setNumberOfHits(nhits);
if (!fit.chernov()) continue; // chernov fit was ok?
//
// radius cut
//
fitRadius = fit.getRadius();
if (fabs(fit.getRadius() - radius) > setup->getHoughMaxDeltaRadius()) continue;
newCandidates(i)->setRadius(fit.getRadius());
//
// Too far from expected origin?
//
fitCenter = CRichPadCoord(fit.getXCenter(), fit.getYCenter());
deltaCenter = newCandidates(i)->getCenter() - fitCenter;
deltaCenter2 = deltaCenter.getX()*deltaCenter.getX()+deltaCenter.getY()*deltaCenter.getY();
if (deltaCenter2 > setup->getHoughMaxFitDis()) continue;
newCandidates(i)->setCenter(fitCenter);
//
// If too close to an already existing candidate:
// replace this if hough2Amplitude is larger, otherwise just forget it
//
j=0;
while (j<candidateList->entries()) {
if ((*candidateList)(j)->getHough2Amplitude()) {
deltaCenter = (*candidateList)(j)->getCenter() - newCandidates(i)->getCenter();
deltaCenter2 = deltaCenter.getX()*deltaCenter.getX()+deltaCenter.getY()*deltaCenter.getY();
if (deltaCenter2 < minDist2) {
if (newCandidates(i)->getHough2Amplitude() > (*candidateList)(j)->getHough2Amplitude()) {
delete candidateList->removeAt(j);
}
else {
newCandidates(i)->setHough2Amplitude(0.0);
j++;
}
}
else
j++;
}
}
//
// insert this candidate, if it has survived
//
if (newCandidates(i)->getHough2Amplitude()) {
candidateList->append(new CRingCandidate(*(newCandidates(i))));
}
}
delete houghList;
newCandidates.clearAndDestroy();
return True;
}
CRingMask* CRich::ringMaskWidth1(double ringMaskRadius)
{
double r, rx, ry;
const double rBullshit = 1.e10;
const int maxN = 1000;
CRingMask *tmpMask = new CRingMask(maxN);
//
// start at 0 deg and go in math.positive direction
//
int quadrant = 0;
int n = 0;
int ix = int(ringMaskRadius+0.5);
int iy = 0;
int xDir = -1;
int yDir = 1;
int Xstep = 0;
int Ystep = yDir;
int newXstep, newYstep;
//
// store the start point
//
tmpMask->x[n] = ix;
tmpMask->y[n] = iy;
n++;
for (;;) {
//
// switch directions at 0, 90, 180, 270 deg
//
if (iy == 0) {
xDir = -yDir;
quadrant++;
}
if (ix == 0) {
yDir = xDir;
quadrant++;
}
if (quadrant > 4) break;
//
// find best new approximation of point on ring
//
newXstep = (abs(Xstep) - 1) * (-xDir); // switch between 0 and xDir
newYstep = (abs(Ystep) - 1) * (-yDir); // switch between 0 and yDir
r = sqrt(double((ix+Xstep)*(ix+Xstep) + (iy+Ystep)*(iy+Ystep)));
if (newXstep!=0 || Ystep!=0)
rx = sqrt(double((ix+newXstep)*(ix+newXstep) + (iy+Ystep)*(iy+Ystep)));
else
rx = rBullshit;
if (Xstep!=0 || newYstep!=0)
ry = sqrt(double((ix+Xstep)*(ix+Xstep) + (iy+newYstep)*(iy+newYstep)));
else
ry = rBullshit;
if (fabs(ringMaskRadius - rx) < fabs(ringMaskRadius - ry)) {
if (fabs(ringMaskRadius - rx) < fabs(ringMaskRadius - r)) Xstep = newXstep;
}
else {
if (fabs(ringMaskRadius - ry) < fabs(ringMaskRadius - r)) Ystep = newYstep;
}
ix = ix + Xstep;
iy = iy + Ystep;
//
// store this point
//
if (n < maxN) {
tmpMask->x[n] = ix;
tmpMask->y[n] = iy;
n++;
}
else {
cerr << "CRich::ringMaskWidth1():\n";
cerr << "\tERROR\n";
cerr << "\tSoR masksize exceeds decleration n = " << n << endl;
delete tmpMask;
return 0;
}
}
//
// truncate to correct size
//
CRingMask *newMask = new CRingMask(n);
for (int i=0; i<n; i++) {
newMask->x[i] = tmpMask->x[i];
newMask->y[i] = tmpMask->y[i];
}
delete tmpMask;
return newMask;
}
CRingMask* CRich::ringMask(double ringMaskRadius, double ringMaskWidth)
{
double rmin2 = (ringMaskRadius-ringMaskWidth/2.)*(ringMaskRadius-ringMaskWidth/2.);
double rmax2 = (ringMaskRadius+ringMaskWidth/2.)*(ringMaskRadius+ringMaskWidth/2.);
int mrange = int(ringMaskRadius+ringMaskWidth+0.5);
//
// First pass to get size of mask
//
int n = 0, ix, iy;
double r2;
for (iy=-mrange; iy<=mrange; iy++) {
for (ix=-mrange; ix<=mrange; ix++) {
r2 = ix*ix+iy*iy;
if (r2 >= rmin2 && r2 <= rmax2) n++;
}
}
CRingMask *ringMask = new CRingMask(n); // allocate mask
//
// Fill mask
//
n = 0;
for (iy=-mrange; iy<=mrange; iy++) {
for (ix=-mrange; ix<=mrange; ix++) {
r2 = ix*ix+iy*iy;
if (r2 >= rmin2 && r2 <= rmax2) {
ringMask->x[n] = ix;
ringMask->y[n] = iy;
n++;
}
}
}
return ringMask;
}
void CRich::markSensitiveArea()
{
//
// Sensitive are all items in the lookup table
// unless they are:
// (i) lying under a spoke
// (ii) are not within the acceptance
//
CRichLikeLookupItem *item;
double radius, radius2, phi, x, y, phiSpoke;
double center = setup->getNPads()/2.+0.5;
double rMin2 = setup->getInnerAcceptance()*setup->getInnerAcceptance();
double rMax2 = setup->getOuterAcceptance()*setup->getOuterAcceptance();
int nSpokes = setup->getNSpokes();
double halfWidth = setup->getSpokeWidth()/2.;
double dphi = TwoPi/nSpokes;
double firstPhi = setup->getFirstSpokeAngle()*ToRadian;
//
// Loop over all lookup items
//
for (int k=0; k<lookupItem.length(); k++) {
item = lookupItem(k);
item->setIsSensitive(True);
if (item->getX() > 0 && item->getY() > 0) {
//
// Check if item is within acceptance
//
x = item->getX()-center;
y = item->getY()-center;
radius2 = x*x+y*y;
if (radius2 >= rMin2 && radius2 <= rMax2) {
//
// Check if under spoke
//
radius = sqrt(radius2);
phi = fabs(atan2(y,x)); // phi of item in range 0-Pi
for (phiSpoke = firstPhi; phiSpoke< Pi; phiSpoke += dphi) {
if (fabs(radius*sin(phiSpoke-phi)) < halfWidth) {
item->setIsSensitive(False);
break;
}
}
}
else
item->setIsSensitive(False);
}
else
item->setIsSensitive(False);
}
}
void CRich::printRingProperties(ostream& ost) const
{
ost << "\nList of fitted rings in " << name << endl;
CString line('=',46);
ost << line << endl;
if (ringList) {
for(int i=0; i<ringList->length(); i++) {
ost << "ring " << i << ":\n";
(*ringList)[i]->printProperties(ost);
ost << endl;
}
}
}
CBoolean CRich::applyRingQualityCuts()
{
//
// No cuts defined yet ...
//
return True;
}
int CRich::assignHitsAndRings( CSortedList<CRichLikeHit> &theHitList, CList<CFittedRing> &theRingList)
{
//
// Calculates the distances of each hit to all fits. The hits
// are assigned to the closest fit if the distance is smaller
// than 'minDeviation'. Returns the index of the ring with the
// fewest points or -1 if the ring list is empty. Note that the
// hits are 'added' to the hit list of the rings. Hits already
// in the list remain there.
// If several rings with 'fewest points' exist the index of the
// one with the highest variance is returned.
//
const double minDeviation = setup->getMaxHitRingDistance(); // in pads
//
// Search closest fit for each hit and assign
// the hit to this ring.
//
CRichLikeHit *hit;
CFittedRing *ring, *closestRing;
float closest, distance;
CRichPadCoord center;
CListIterator<CFittedRing> riter(theRingList);
for (int i=0; i<theHitList.length(); i++) {
hit = theHitList(i);
closest = minDeviation;
closestRing = 0;
riter.reset();
while (ring = riter()) {
center = ring->getPadCenter();
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 < closest) {
closest = distance;
closestRing = ring;
}
}
if (closestRing != 0) { // present hit belongs to ring
hit->setRing(closestRing); // hit refers to its ring
closestRing->addHit(hit); // add hit to ring
}
}
//
// Find lowest number of hits per ring
//
riter.reset();
int fewestHits = INT_MAX;
while (ring = riter()) {
if (ring->getNumberOfHits() < fewestHits) {
fewestHits = ring->getNumberOfHits();
}
}
//
// Find - out of the rings with the lowest number of hits -
// the one with the highest variance.
//
riter.reset();
int worstRing = -1;
double highestVariance = 0;
while (ring = riter()) {
if (ring->getNumberOfHits() == fewestHits &&
ring->getVariance() > highestVariance) {
highestVariance = ring->getVariance();
worstRing = riter.pos();
}
}
return worstRing;
}
CBoolean CRich::makeRingFits()
{
if (candidateList == 0 || hits.isEmpty()) return False;
delete ringList; // change for persistency
ringList = new CList<CFittedRing>;
int j;
for (j=0; j<hits.length(); j++) // clear references to rings
hits(j)->setRing(0);
//
// Loop Rich ring candidates
//
int k;
CRichLikeHit *hit;
CFittedRing *ring;
CRingCandidate *candidate;
CRingFitter fitter;
CBoolean isTooClose;
double dx, dy;
double range = 1.25*setup->getAsymptoticRadius();
CRichPolarCoord polarCoord;
CRichPadCoord padCoord;
CListIterator<CRingCandidate> citer(*candidateList);
CListIterator<CFittedRing> riter(*ringList);
while (candidate = citer()) {
fitter.reset();
//
// Accumulate all hits around ring candidate,
// feed them into the fitter.
//
for (k=0; k<hits.length(); k++) {
hit = hits(k);
if (fabs(hit->getX()-candidate->getCenter().getX()) < range &&
fabs(hit->getY()-candidate->getCenter().getY()) < range) {
fitter.addX(hit->getX());
fitter.addY(hit->getY());
fitter.addWeight(1.);
}
}
if (fitter.getNFitPoints() < setup->getMinHitsPerRing()) continue;
//
// Pass start values for center and radius
// and number of iterations.
//
fitter.setMaxIter(100);
fitter.setXCenter(candidate->getCenter().getX());
fitter.setYCenter(candidate->getCenter().getY());
fitter.setRadius(setup->getAsymptoticRadius());
//
// Fit the candidate, perform some simple checks and
// store the new ring in the referring list.
//
if (fitter.robustFixedRadius() && fitter.getError() == 0) {
//
// Check if center is too close to an already existing one
//
isTooClose = False;
riter.reset();
while (ring = riter()) {
dx = ring->getPadCenter().getX()-fitter.getXCenter();
dy = ring->getPadCenter().getY()-fitter.getYCenter();
if (dx*dx + dy*dy < 1) {
isTooClose = True;
break;
}
}
//
// Create and store the new ring
//
if (!isTooClose) {
ring = new CFittedRing;
ring->setType(CRing::Electron);
padCoord = CRichPadCoord(fitter.getXCenter(), fitter.getYCenter());
transform(polarCoord, padCoord);
ring->setCenter(padCoord, polarCoord);
ring->setRadius(fitter.getRadius());
ring->setVariance(fitter.getVariance());
ring->evalActiveCircumference(lookupItem);
ring->setRingCandidate(candidate);
ringList->append(ring);
}
}
}
if (ringList->entries() == 0) return True; // nothing found
//
// Assign hits to fits and vice versa and eliminate fits
// with too few points left. The destructor of CRing will
// clear the hit list (not destroy it).
// Remember that assignHitsAndRings returnd the index of
// the ring with the fewest hits or -1 if the list is empty.
//
int index;
CSortedList<CRichLikeHit> remainingHits;
index = assignHitsAndRings(hits, *ringList);
while (index != -1 &&
(*ringList)(index)->getNumberOfHits() < setup->getMinHitsPerRing()) {
remainingHits = (*ringList)(index)->getHits(); // shallow copy
for (j=0; j<remainingHits.length(); j++) // clear references to rings
remainingHits(j)->setRing(0);
ring = ringList->removeAt(index);
delete ring;
index = assignHitsAndRings(remainingHits, *ringList);
}
remainingHits.clear();
//
// Kolmogorov-Smirnov test, evaluation of chi2
// and analog ring sum.
//
riter.reset();
while (ring = riter()) {
//ring->evalKolmogorovTestValue(lookupItem);
ring->evalRingSum(lookupItem, pads, setup->getHalfRingSumMask());
}
//
// Quality check all rings. If a ring is kicked out at this
// stage, the hits have to be redistributed.
//
CBoolean stillRingsToCheck = True;
int i;
while (stillRingsToCheck) {
stillRingsToCheck = False;
for (i=0; i<ringList->length(); i++) {
ring = (*ringList)(i);
ring->evalChi2(setup->getAsymptoticRadius());
if (ring->getChi2()/(ring->getNumberOfHits()-3) > setup->getMaxChi2PerDofForRing()) {
remainingHits = ring->getHits(); // shallow copy
for (j=0; j<remainingHits.length(); j++) // clear references to rings
remainingHits(j)->setRing(0);
ring = ringList->removeAt(i);
delete ring;
index = assignHitsAndRings(remainingHits, *ringList);
remainingHits.clear();
stillRingsToCheck = True;
break;
}
}
}
return True;
}
//
// The following are coordinate transformation routines.
//
// In transformation from CRichPadCoord into CRichPolarCoord
// all necessary calibrations corrections are performed.
// The related calibrations constants are taken from the
// setup object constructed according to the referring setup
// file.
//
void CRich::transform(CRichXYCoord& xy, const CRichPolarCoord& pol) const
{
xy.setX(pol.getTheta()*cos(pol.getPhi()));
xy.setY(pol.getTheta()*sin(pol.getPhi()));
}
void CRich::transform(CRichPolarCoord& pol, const CRichXYCoord& xy) const
{
pol.setTheta(sqrt(xy.getX()*xy.getX()+xy.getY()*xy.getY()));
pol.setPhi(atan2(xy.getY(), xy.getX()));
}
void CRich::transform(CRichPolarCoord& gp, const CRichPadCoord& rxy) const
{
if (setup->getUseVariableFocalLength())
transformWithVariableFocalLength(gp, rxy);
else
transformWithConstantFocalLength(gp, rxy);
}
void CRich::transform(CRichPadCoord& rxy, const CRichPolarCoord& gp) const
{
if (setup->getUseVariableFocalLength())
transformWithVariableFocalLength(rxy, gp);
else
transformWithConstantFocalLength(rxy, gp);
}
void CRich::transform(CRichXYCoord& gxy, const CRichPadCoord& rxy) const
{
if (setup->getUseVariableFocalLength())
transformWithVariableFocalLength(gxy, rxy);
else
transformWithConstantFocalLength(gxy, rxy);
}
void CRich::transform(CRichPadCoord& rxy, const CRichXYCoord& gxy) const
{
if (setup->getUseVariableFocalLength())
transformWithVariableFocalLength(rxy, gxy);
else
transformWithConstantFocalLength(rxy, gxy);
}
// here follow the protected ones which are used....
void CRich::transformWithConstantFocalLength(CRichPolarCoord& gp, const CRichPadCoord& rxy) const
{
//
// (x [pads], y [pads]) -> (theta [mrad], phi [rad])
// A analytical inverse calculations for CRichPolarCoord -> CRichPadCoord
// does not exist. This is a pure numerical ansatz to solve the problem.
// Note the order in which the different calibration corrections are applied.
//
const double Range = 0.01;
const int MaxIterations = 100;
double center = setup->getNPads()/2.+0.5;
double x = rxy.getX() - center;
double y = rxy.getY() - center;
//
// Apply corrections:
// 1. Correct for xy offset
//
x -= setup->getXOffset();
y -= setup->getYOffset();
double rcm = setup->getPadSize()*sqrt(x*x+y*y);
double twoFocal = 2.*setup->getFocalLength();
double zt = twoFocal - setup->getMirrorDistance();
//
// Regula Falsi
//
double a = rcm/setup->getFudgeFocalLength()-Range;
double b = rcm/setup->getFudgeFocalLength()+Range;
double delta, result_a, result_b;
double tantheta2, zm, rm , rb;
for (int i = 1; i < MaxIterations; i++) {
//
// Get result_b = f(b)
//
tantheta2 = tan(b)*tan(b);
zm = (zt*tantheta2+sqrt((zt*tantheta2)*(zt*tantheta2)-(tantheta2+1.)*
(zt*zt*tantheta2-(twoFocal*twoFocal))))/(tantheta2+1.);
rm = (zm-zt)*tan(b);
rb = rm+(zm-(twoFocal-setup->getDetectorDistance()))*
tan(b-2.*asin(rm/twoFocal));
result_b = rb - rcm;
//
// Get result_a = f(a)
//
tantheta2 = tan(a)*tan(a);
zm = (zt*tantheta2+sqrt((zt*tantheta2)*(zt*tantheta2)-(tantheta2+1.)*
(zt*zt*tantheta2-(twoFocal*twoFocal))))/(tantheta2+1.);
rm = (zm-zt)*tan(a);
rb = rm+(zm-(twoFocal-setup->getDetectorDistance()))*
tan(a-2.*asin(rm/twoFocal));
result_a = rb - rcm;
//
// Prepare next (a,b) until done
//
delta = (result_b-result_a)/(b-a);
a = b;
b -= result_b/delta;
if (fabs(result_b) < FLT_EPSILON) break;
}
b *= 1000; // in mrad
//
// Apply corrections:
// 2. Correction of theta scale
//
b += setup->getDeltaThetaCorrections(0) + setup->getDeltaThetaCorrections(1)*b;
gp.setTheta(b);
gp.setPhi(atan2(y, x));
}
void CRich::transformWithConstantFocalLength(CRichPadCoord& rxy, const CRichPolarCoord& gp) const
{
//
// (theta [mrad], phi [rad]) -> (x [pads], y [pads])
// The referring calculations were done by P. Glaessel.
// Note the order in which the different calibration corrections are applied.
//
double theta = gp.getTheta(); // in mrad
double phi = gp.getPhi(); // in rad
//
// Apply (undo) corrections:
// 2. Correction of theta scale
//
theta = (theta - setup->getDeltaThetaCorrections(0))/
(1 + setup->getDeltaThetaCorrections(1));
theta /= 1000.; // mrad -> rad
double center = setup->getNPads()/2.+0.5; // nominal center in pads
double twofocal = 2*setup->getFocalLength();
double tantheta2 = tan(theta)*tan(theta);
double zt = twofocal - setup->getMirrorDistance();
double zm = (zt*tantheta2 + sqrt((zt*tantheta2)*(zt*tantheta2) - (tantheta2+1.)*
(zt*zt*tantheta2 - twofocal*twofocal)))/(tantheta2 + 1.);
double rm = (zm - zt)*tan(theta);
double rb = rm+(zm - (twofocal - setup->getDetectorDistance()))*
tan(theta - 2.*asin(rm/twofocal));
//
// Apply (undo) corrections:
// 2. Correct for xy offset
//
rxy.setX(rb/setup->getPadSize()*cos(phi)+center+setup->getXOffset());
rxy.setY(rb/setup->getPadSize()*sin(phi)+center+setup->getYOffset());
}
void CRich::transformWithConstantFocalLength(CRichXYCoord& gxy, const CRichPadCoord& rxy) const
{
CRichPolarCoord gp;
transform(gp, rxy);
transform(gxy, gp);
}
void CRich::transformWithConstantFocalLength(CRichPadCoord& rxy, const CRichXYCoord& gxy) const
{
CRichPolarCoord gp;
transform(gp, gxy);
transform(rxy, gp);
}
CList<CRingCandidate>* CRich::detachCandidateList()
{
CList<CRingCandidate>* copyToReturn = candidateList;
candidateList = 0;
return copyToReturn;
}
CList<CFittedRing>* CRich::detachRingList()
{
CList<CFittedRing>* copyToReturn = ringList;
ringList = 0;
return copyToReturn;
}
CBoolean CRich::getSegmentFocalLengthDescription(double phi,
double& focalLength,
double& focalLengthCorrection,
double& focalLengthCorrectionOffset) const {
//
// Function to determine the description of the theta dependend
// focal length from the quantities given in the setup.
// This function performs if necessary the mapping between
// phi angles and an arbitrarily choosen mirror number in the setup
// counting the mirror segments of the Rich detectors.
// It is choosen such that mirror #0 has its upper edge
// at nine o'clock ( parallel to the x-axis )
//
int mirror;
double lowPhi, upperPhi, mirrorPhiRange=2*Pi/setup->getNMirrorSegments();
double mirrorOffset = setup->getMirrorSegmentOffset();
for ( mirror=0; mirror<setup->getNMirrorSegments(); mirror++ ) {
// calculate mirror edges
lowPhi = mirror*mirrorPhiRange - Pi + mirrorOffset;
upperPhi = lowPhi+mirrorPhiRange;
if (phi > Pi + mirrorOffset) phi = phi - 2*Pi;
if (phi < - Pi + mirrorOffset) phi = phi + 2*Pi;
if (lowPhi<phi && phi<=upperPhi) {
focalLength = setup->getSegmentFocalLength(mirror);
focalLengthCorrection = setup->getSegmentFocalLengthCorrection(mirror);
focalLengthCorrectionOffset = setup->getSegmentFocalLengthCorrectionOffset();
return True;
}
}
return False;
}
//
// The following are coordinate transformation routines taking into
// account the mirror quality, namely the change of the focal length.
//
// In transformation from CRichPadCoord into CRichPolarCoord
// all necessary calibrations corrections are performed.
// The related calibrations constants are taken from the
// setup object constructed according to the referring setup
// file.
//
void CRich::transformWithVariableFocalLength(CRichPolarCoord& gp, const CRichPadCoord& rxy) const
{
//
// (x [pads], y [pads]) -> (theta [mrad], phi [rad])
// A analytical inverse calculations for CRichPolarCoord -> CRichPadCoord
// does not exist. This is a pure numerical ansatz to solve the problem.
// Note the order in which the different calibration corrections are applied.
//
const double Range = 0.02;
const int MaxIterations = 100;
double center = setup->getNPads()/2.+0.5;
double x = rxy.getX() - center;
double y = rxy.getY() - center;
//
// Apply corrections:
// 1. Correct for xy offset
//
x -= setup->getXOffset();
y -= setup->getYOffset();
//
// determine phi angle including correction
//
gp.setPhi(atan2(y, x) + setup->getPhiOffset());
//
// determine correction for the mirror at this angle
//
double focalLength, focalLengthCorrection, focalLengthCorrectionOffset;
getSegmentFocalLengthDescription(gp.getPhi(), focalLength, focalLengthCorrection, focalLengthCorrectionOffset);
double rcm = setup->getPadSize()*sqrt(x*x+y*y);
//
// Regula Falsi
//
double a = rcm/setup->getFudgeFocalLength()-Range;
double b = rcm/setup->getFudgeFocalLength()+Range;
double delta, result_a, result_b;
double zMirror, rMirror, rPlane, rZero, rTheta, zTarget;
rZero = 2.* (focalLength + focalLengthCorrection*(0-focalLengthCorrectionOffset));
for (int i = 1; i < MaxIterations; i++) {
//
// Get result_b = f(b)
//
zTarget = rZero - setup->getMirrorDistance();
rTheta = 2.* (focalLength+
focalLengthCorrection*(b-focalLengthCorrectionOffset));
zMirror = rTheta*cos(b-asin(zTarget*sin(b)/rTheta));
rMirror = (zMirror - zTarget)*tan(b);
rPlane = rMirror +
(zMirror - (rZero - setup->getDetectorDistance()))
*tan(b - 2.*asin(rMirror/rTheta));
result_b = rPlane - rcm;
//
// Get result_a = f(a)
//
zTarget = rZero - setup->getMirrorDistance();
rTheta = 2.*(focalLength+
focalLengthCorrection*(a-focalLengthCorrectionOffset));
zMirror = rTheta*cos(a-asin(zTarget*sin(a)/rTheta));
rMirror = (zMirror - zTarget)*tan(a);
rPlane = rMirror +
(zMirror - (rZero - setup->getDetectorDistance()))
*tan(a - 2.*asin(rMirror/rTheta));
result_a = rPlane - rcm;
//
// Prepare next (a,b) until done
//
delta = (result_b-result_a)/(b-a);
a = b;
b -= result_b/delta;
if (fabs(result_b) < FLT_EPSILON) break;
if ( i == MaxIterations-1 )
cerr << "CRich2::transformWithVariableFocalLength: Warning unable to invert function"
<< endl;
}
gp.setTheta(b*1000); // to mrad
}
void CRich::transformWithVariableFocalLength(CRichPadCoord& rxy, const CRichPolarCoord& gp) const
{
//
// (theta [mrad], phi [rad]) -> (x [pads], y [pads])
// The referring calculations were done by C. Voigt.
// Note the order in which the different calibration corrections are applied.
//
double center = setup->getNPads()/2.+0.5; // nominal center in pads
double theta = gp.getTheta(); // in mrad
double phi = gp.getPhi(); // in rad
theta /= 1000.; // mrad -> rad
//
// determine correction for the mirror at this angle
//
double focalLength, focalLengthCorrection, focalLengthCorrectionOffset;
getSegmentFocalLengthDescription(phi, focalLength, focalLengthCorrection, focalLengthCorrectionOffset);
//
// calculate positions on mirror and plane
//
double rZero, zTarget, rTheta, zMirror, rMirror, rPlane;
rZero = 2.* (focalLength+
focalLengthCorrection*(0-focalLengthCorrectionOffset));
zTarget = rZero - setup->getMirrorDistance();
rTheta = 2.* (focalLength+
focalLengthCorrection*(theta-focalLengthCorrectionOffset));
zMirror = rTheta*cos(theta-asin(zTarget*sin(theta)/rTheta));
rMirror = (zMirror - zTarget)*tan(theta);
rPlane = rMirror +
(zMirror - (rZero - setup->getDetectorDistance()))
*tan(theta - 2.*asin(rMirror/rTheta));
//
// Apply (undo) corrections:
// 1. phi correction
// 2. Correct for xy offset
//
phi -= setup->getPhiOffset();
rxy.setX(rPlane/setup->getPadSize()*cos(phi)+center+setup->getXOffset());
rxy.setY(rPlane/setup->getPadSize()*sin(phi)+center+setup->getYOffset());
}
void CRich::transformWithVariableFocalLength(CRichXYCoord& gxy, const CRichPadCoord& rxy) const
{
CRichPolarCoord gp;
transformWithVariableFocalLength(gp, rxy);
transform(gxy, gp);
}
void CRich::transformWithVariableFocalLength(CRichPadCoord& rxy, const CRichXYCoord& gxy) const
{
CRichPolarCoord gp;
transform(gp, gxy);
transformWithVariableFocalLength(rxy, gp);
}
CBoolean CRich::rotateHits(CAngle rotationAngle)
{
int i;
RWTPtrOrderedVector<CRichLikeHit> rotatedHits;
//
// Loop over all hits to change the coordinate
//
CRichLikeHit *hit;
CRichXYCoord hitXY;
CRichPolarCoord hitPolar;
for ( i = 0; i<hits.length(); i++) {
hit = hits(i);
CRichPadCoord hitPad(hit->getX(),hit->getY());
transform(hitXY,hitPad);
transform(hitPolar,hitXY);
hitPolar.setPhi((CAngle)hitPolar.getPhi() + rotationAngle);
transform(hitXY, hitPolar);
transform(hitPad, hitXY);
hit->setX(hitPad.getX());
hit->setY(hitPad.getY());
rotatedHits.insert(hit);
}
//
// now fill the hits in sorted vector
//
hits.clear();
for ( i = 0; i<rotatedHits.length(); i++) {
hit = rotatedHits(i);
hits.insert(hit);
}
rotatedHits.clear();
return True;
}
void CRich::makeDigiHits(CMCServer& mcServ)
{
const CDetectorId detId = getDetectorId();
const int NPADS = setup->getNPads();
const double meanAmp = setup->getMCMeanAmplitude();
const double sigHit = setup->getMCHitWidth();
const double sigRadDiff = setup->getMCHitDiffusion();
const double sigRadAber = setup->getMCHitAberration();
const CBoolean useGainMap = setup->getMCUseGainMap();
const CBoolean realNoise = setup->getMCUseRealNoise();
const CBoolean simpleNoise = setup->getMCUseSimpleNoise();
const double meanSigma = setup->getMCMeanSigma();
const double ringSigma = setup->getMCRingSigma();
float minPadAmplitude = (float) setup->getMinPadAmplitude();
float nSigmaCut = (float) setup->getNSigmaCut();
CRichPadCoord padCoordinate;
CRichPolarCoord polarCoordinate;
double xoffset,yoffset;
double thetaCorrection[2];
const RWTPtrOrderedVector<CMCDigiHit> *digiHits;
CMCDigiHit *theDigiHit;
float amp;
double dxRing, dyRing, dx1, dy1, dx2, dy2, phi, dr2, drr2, range;
int i;
CRandom rndmGen;
digiHits = mcServ.getMCDigiHits( detId );
pads.resize( (pads.length() + digiHits->length()*10) ); // assume 10 pads per hit
//
// simulation of not understandable bad ring center resolution
//
dxRing = rndmGen.gauss() * ringSigma;
dyRing = rndmGen.gauss() * ringSigma;
for (i=0; i<digiHits->length(); i++) {
//
// simulation of transvers diffusion
//
theDigiHit = (*digiHits)(i);
dx1 = rndmGen.gauss() * sigRadDiff;
dy1 = rndmGen.gauss() * sigRadDiff;
//
// simulation of chromatic Aberration
//
range = sigRadAber*sqrt(12)*0.5; // calculate range for rectangle distribution
drr2 = rndmGen.flat(-range,range);
phi = rndmGen.flat(0.,2*Pi);
dr2 = sqrt(range*range-drr2*drr2);
dx2 = dr2*sin(phi);
dy2 = dr2*cos(phi);
//
// modify GEANTs x/y according to detector parameters (diffusion, chromatics ... )
// and set an amplitude using the given mean
//
padCoordinate.setX(theDigiHit->getGeantHit().getX() + dxRing + dx1 + dx2);
padCoordinate.setY(theDigiHit->getGeantHit().getY() + dyRing + dy1 + dy2);
amp = float(-log( drand48() ) * meanAmp);
//
// modify GEANTs x/y according to global detector mismatches
//
xoffset = setup->getXOffset();
yoffset = setup->getYOffset();
thetaCorrection[0] = setup->getDeltaThetaCorrections(0);
thetaCorrection[1] = setup->getDeltaThetaCorrections(1);
setup->setXOffset(0.);
setup->setYOffset(0.);
setup->setDeltaThetaCorrections(0.,0.);
setup->setUseVariableFocalLength(0);
transform(polarCoordinate,padCoordinate);
setup->setXOffset(xoffset);
setup->setYOffset(yoffset);
setup->setDeltaThetaCorrections(thetaCorrection[0],thetaCorrection[1]);
setup->setUseVariableFocalLength(1);
transform(padCoordinate,polarCoordinate);
theDigiHit->getDigiHit().setX(padCoordinate.getX()); // store info for subsequent use
theDigiHit->getDigiHit().setY(padCoordinate.getY());
theDigiHit->getDigiHit().setAmp(amp);
} // loop mchits
}