CPadChamber.C
//-----------------------------------------------------------------------------
// $Header: /asis/offline/ceres/cool/project/RCS/CPadChamber.C,v 2.5 1997/07/01 10:58:40 lenkeit Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class CPadChamber.
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include <stdio.h>
#include "CPadChamber.h"
#include "CRandom.h"
CPadChamber::CPadChamber(const char* setupFile)
{
//
// Initialize members inherited from CDetector
//
id = PadChamber;
name = "PadChamber";
//
// May be I should comment on this:
// The reason for this massacre is that there exist nothing
// like virtual data member. Each of the defined setup pointers
// has is own scope. With other words there's an ambiguity
// which can only be resolved by either using always the proper
// scope operator or by hacking the code below. I know it looks
// somehow unusual but eases further programming.
//
CDetector::setup = CRichLikeDetector::setup = setup = new CPadChamberSetup;
if (setupFile == 0) {
CString setupName = CString(C_DEFAULT_SETUP_PATH)+
CString(C_SETUPFILE_PADCHAMBER);
setup->read(setupName.data());
}
else
setup->read(setupFile);
//
// Set the proper size of the collections
//
pads.resize(setup->getNPads()+1, setup->getNPads()+1);
lookupItem.resize(setup->getNPads()+1, setup->getNPads()+1);
//
// Read the module map
//
CString modName = CString(C_DEFAULT_SETUP_PATH)+
CString(C_MODULEMAP_PADCHAMBER);
readModulesFromFile(modName.data());
//
// Read the default huffman table
//
CString tabName = CString(C_DEFAULT_SETUP_PATH)+
CString(C_HUFFMAN_TABLE_PADCHAMBER);
readHuffmanTableFromFile(tabName.data());
}
CPadChamber::~CPadChamber() { delete setup; }
CPadChamber::CPadChamber(const CPadChamber &) {}
CPadChamber & CPadChamber::operator = (const CPadChamber &) { return *this; }
void CPadChamber::markSensitiveArea() {}
float CPadChamber::getLeft() const { return 0.5; }
float CPadChamber::getRight() const { return setup->getNPads()+0.5; }
float CPadChamber::getBottom() const { return 0.5; }
float CPadChamber::getTop() const { return setup->getNPads()+0.5; }
// Draw a CPadChamber.
// 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
// v mark sensitive area
// r draw removed pads only ( excludes most of the other pad options )
// i draw cluster ids
//
void CPadChamber::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);
}
}
//
// 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 {
//
// 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.
//
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);
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);
}
}
//
// 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);
}
}
}
}
CBoolean CPadChamber::makeClustersAndDoCleanup()
{
if (!makeClusters(1)) return False;
if (!removeClusters()) return False;
if (!splitClusters()) return False;
return True;
}
CBoolean CPadChamber::removeClusters()
{
if ( clusters.length() <= 0) return False; // nothing to do
int i, j, ioff, iend;
for ( i=0; i<clusters.length(); i++) {
if (!isGoodCluster(*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 CPadChamber::isGoodCluster(CRichLikeCluster& cluster)
{
int ioff = cluster.getFirstPad();
int npads = cluster.getNPads();
int iend = ioff+npads;
//
// collect all properties..
//
int minX = 300; int minY = 300; int maxX = 0; int maxY = 0;
int nsat = 0;
float largestPadAmp = 0;
for (int i=ioff; i<iend; i++) {
if (pads(i)->getAmp() > largestPadAmp) largestPadAmp = pads(i)->getAmp();
if (pads(i)->getAmp() >= setup->getClusterSatPadAmp()) nsat++;
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();
}
int xSize = maxX - minX + 1;
int ySize = maxY - minY + 1;
//
// reject too small clusters (noise)
//
if (npads < setup->getClusterMinPads()) return False;
if (npads == 1 && largestPadAmp < setup->getClusterMinLargestPadAmp()) return False;
//
// reject clusters with too many saturated pads
//
if (nsat > setup->getClusterMaxSat()) return False;
//
// reject clusters with too many pads
//
if (npads >= setup->getClusterMaxPads()) return False;
//
// reject (geometrically) too big clusters
//
if (xSize * ySize > setup->getClusterMaxSize()) return False;
//
// cluster survived all cuts
//
return True;
}
void CPadChamber::makeDigiHits(CMCServer& mcServ)
{
const CDetectorId detId = getDetectorId();
const int NPADS = setup->getNPads();
const double meanAmp = setup->getMCMeanAmplitude();
const double sigHit = setup->getMCHitWidth();
const double sigRad = setup->getMCHitAberration();
const CBoolean useGainMap = setup->getMCUseGainMap();
const CBoolean realNoise = setup->getMCUseRealNoise();
const CBoolean simpleNoise = setup->getMCUseSimpleNoise();
const double meanSigma = setup->getMCMeanSigma();
float minPadAmplitude = (float) setup->getMinPadAmplitude();
float nSigmaCut = (float) setup->getNSigmaCut();
const RWTPtrOrderedVector<CMCDigiHit> *digiHits;
CMCDigiHit *theDigiHit;
float x, y, amp;
double dx, dy;
float hit[5][5];
int i;
CRandom rndmGen;
CPadCPadCoord pad;
CLabXYZCoord xyz;
digiHits = mcServ.getMCDigiHits( detId );
pads.resize( (pads.length() + digiHits->length()*10) ); // assume 10 pads per hit
for (i=0; i<digiHits->length(); i++) {
//
// modify GEANTs x/y according to detector parameters (diffusion, chromatics ... )
// and set an amplitude using the given mean
//
theDigiHit = (*digiHits)(i);
// transform to right z-Position
dx = cos(theDigiHit->getPhi()) * tan(theDigiHit->getTheta()) * setup->getMCZOffset();
dy = sin(theDigiHit->getPhi()) * tan(theDigiHit->getTheta()) * setup->getMCZOffset();
xyz.setX(theDigiHit->getGeantHit().getX() - dx);
xyz.setY(theDigiHit->getGeantHit().getY() - dy);
xyz.setZ(theDigiHit->getGeantHit().getZ());
// transform from cm to pads...
transform(pad,xyz);
theDigiHit->getGeantHit().setX(pad.getX());
theDigiHit->getGeantHit().setY(pad.getY());
theDigiHit->getGeantHit().setZ(setup->getZPosition());
dx = rndmGen.gauss(); // get normal distributed random numbers
dy = rndmGen.gauss();
x = theDigiHit->getGeantHit().getX() + dx*sigRad;
y = theDigiHit->getGeantHit().getY() + dy*sigRad;
if (meanAmp < 1000){
amp = float(-log( drand48()*drand48() ) * meanAmp); // to simulate the funny '95 distribution
}
else{
amp = theDigiHit->getGeantHit().getAmp() * meanAmp; // normal landau distribution
}
theDigiHit->getDigiHit().setX(x); // store info for subsequent use
theDigiHit->getDigiHit().setY(y);
theDigiHit->getDigiHit().setAmp(amp);
} // loop mchits
}
// The following are coordinate transformation routines.
//
// In the transformation from CPadCPadCoord into CLabXYZCoord
// all necessary calibration corrections are performed.
// The related calibration constants are taken from the
// setup object constructed according to the referring setup
// file.
//
void CPadChamber::transform(CPadCPadCoord& pad, const CLabXYZCoord& xyz) const
{
//
// z does not exist in pad coordinates
//
double xOriginInPadUnits = double(setup->getNPads())/2.+0.5 + setup->getXOffset();
double yOriginInPadUnits = double(setup->getNPads())/2.+0.5 + setup->getYOffset();
pad.setX(xyz.getX() / setup->getPadSize() + xOriginInPadUnits);
pad.setY(xyz.getY() / setup->getPadSize() + yOriginInPadUnits);
}
void CPadChamber::transform(CLabXYZCoord& xyz, const CPadCPadCoord& pad) const
{
//
// The z-value gives the absolute
// distance from SiDC-1 since the first SiDC defines
// the origin (0,0,0).
//
double xOriginInPadUnits = double(setup->getNPads())/2.+0.5 + setup->getXOffset();
double yOriginInPadUnits = double(setup->getNPads())/2.+0.5 + setup->getYOffset();
xyz.setX((pad.getX() - xOriginInPadUnits) * setup->getPadSize());
xyz.setY((pad.getY() - yOriginInPadUnits) * setup->getPadSize());
xyz.setZ(setup->getZPosition());
}
void CPadChamber::transform(CPadCPadCoord& pad, const CLabCylinderCoord& cyl) const
{
CLabXYZCoord xyz;
::transform(xyz, cyl);
transform(pad, xyz);
}
void CPadChamber::transform(CLabCylinderCoord& cyl, const CPadCPadCoord& pad) const
{
CLabXYZCoord xyz;
transform(xyz, pad);
::transform(cyl, xyz);
}
CBoolean CPadChamber::rotateHits(CAngle rotationAngle)
{
int i;
RWTPtrOrderedVector<CRichLikeHit> rotatedHits;
//
// Loop over all hits to change the coordinate
//
CRichLikeHit *hit;
CLabCylinderCoord hitPolar;
for ( i = 0; i < hits.length(); i++) {
hit = hits(i);
CPadCPadCoord hitPad(hit->getX(),hit->getY());
transform(hitPolar, hitPad);
hitPolar.setPhi(hitPolar.getPhi() + rotationAngle);
transform(hitPad,hitPolar);
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;
}