CRichLikeDetector.C
//-----------------------------------------------------------------------------
// $Header: /asis/offline/ceres/cool/project/RCS/CRichLikeDetector.C,v 3.6 1997/07/16 15:11:46 messer Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class CRichLikeDetector.
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include <strstream.h>
#include <fstream.h>
#include <math.h>
#include <stdlib.h>
#include "CDataStream.h"
#include "CEventServer.h"
#include "CBitStream16.h"
#include "CRichLikeDetector.h"
#include "CRichLikeSetup.h"
#include "CRichLikeCluster.h"
#include "CCoordinate.h"
#include "rw/tstack.h"
#include "rw/tpordvec.h"
#include "rw/tvordvec.h"
#include "CCollection.h"
#include "CPad.h"
#include "CMCServer.h"
#include "CMCDigiHit.h"
#include "CMCParticle.h"
#include "CRandom.h"
CHuffmanCode::CHuffmanCode()
{
lengthInBits = 0;
code = 0;
}
CHuffmanCode::CHuffmanCode(unsigned int lb, unsigned short dc)
{
lengthInBits = lb;
code = dc;
}
CRichLikeDetector::CRichLikeDetector()
{
huffmanTable = 0;
bigHuffmanTable = 0;
}
CRichLikeDetector::~CRichLikeDetector()
{
delete [] huffmanTable;
delete [] bigHuffmanTable;
pads.clearAndDestroy();
lookupItem.clearAndDestroy();
clusters.clearAndDestroy();
modules.clearAndDestroy();
hits.clearAndDestroy();
}
CBoolean CRichLikeDetector::unpack(CEventServer& server)
{
//
// General parameters and initialisations
//
const int maxChains = 16;
const int bitsPerPad = 9;
const int maxLookup = setup->getMaxLookupPerChain() * setup->getMaxUsedChains();
const int maxNumberOfPads = maxLookup;
pads.clearAndDestroy(); // clear old collections
errorCollection.newEvent();
//
// Check existence of XY lookup table, huffman table, pedestals and sigma
// Unpack those if necessary. The last SoR must be present.
//
const CRawSoR *SoR = server.getSoR();
if (SoR == 0) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tSoR label not found. No further action.\n";
return False;
}
if (lookupItem.length() == 0 || (lastRunNumber != SoR->getRunNumber())) {
if (!unpackXYLookupTable(server)) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tFailed to unpack xy lookup table for "
<< name << ". No further action.\n";
return False;
}
if (!unpackPedestals(server)) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tWARNING\n";
cerr << "\tFailed to unpack pedstals for "
<< name << ". Continue despite warnings.\n";
}
if (!unpackHuffmanTable(server)) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tWARNING\n";
cerr << "\tFailed to unpack huffman table for "
<< name << ". Use default table.\n";
}
lastRunNumber = SoR->getRunNumber();
}
//
// Get raw data
//
CLabel *label = (CLabel *) server.getLabel(setup->getDataLabel());
if (label == 0) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tData label for " << name << " not found.\n";
return False;
}
//
// Nothing to do
//
if (label->getSize() <= 0) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tNo data to unpack for " << name << endl;
return False;
}
//
// Set default sizes for collection classes
//
pads.resize(label->getSize()*sizeof(CWord)*2);
//
// Get label.
// Short words are swapped within one data word. (2,1),(4,3),...
// To ease further proccessing we swap all words in the label
// already here in one go.
//
if (!label->isSwapped()) label->swap();
#ifdef ARCH_IS_LITTLE_ENDIAN
if (label->isSwapped()) label->swap(); // swap back for little-endians
#endif
const CWord *evt = label->getData();
//
// Skip sync data
//
int index = 0;
index += (unsigned short) evt[index];
//
// Get the local hit counter for each chain.
// Hit counters give measure in short words.
// Note that the pad chamber contains only 8
// valid receiver values. To keep the handling
// common for the further unpacking procedure
// we explicitely zero the last 8 receivers.
//
short *chainLength = (short*) &evt[index];
index += maxChains*sizeof(short)/sizeof(CWord);
//
// Skip sync data
//
index += (unsigned short) evt[index];
//
// Check the hit counter.
//
int len = 0;
for (int i=0; i<setup->getMaxUsedChains(); i++) {
len += chainLength[i];
if (chainLength[i] > setup->getMaxPadsPerChain()*bitsPerPad) {
if (setup->getWarningLevel() < 10) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tLength of chain (receiver) #" << i << " in " << name << endl;
cerr << "\tExceeds max allowed size of " << setup->getMaxPadsPerChain() << " pads\n";
cerr << "\tsize = " << chainLength[i]*sizeof(short) << endl;
}
errorCollection.remember("chainTooLong", i);
return False;
}
}
if (len*sizeof(short) > (label->getSize()-index)*sizeof(CWord)) {
if (setup->getWarningLevel() < 10) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tSum of all " << name << " chain lengths exceeds actual data";
cerr << " length by " << len*sizeof(short)-(label->getSize()-index)*sizeof(long) << " bytes\n";
}
errorCollection.remember("chainLengthExceedsData");
return False;
}
//
// Unpack all chains/receivers, decode them, get x,y and subtract pedestals.
//
CBoolean nextPatternMustBeAmplitude;
int chainLengthInBytes;
unsigned int lookupPosition = 0;
unsigned int validBits;
//unsigned int maxCodeLength = 9;
unsigned int maxCodeLength = 16;
unsigned int next16Bit;
unsigned long padsPerChain;
CBitStream16 *bstream;
unsigned long pattern, nZero;
float amplitude;
float minPadAmplitude = (float) setup->getMinPadAmplitude();
float nSigmaCut = (float) setup->getNSigmaCut();
unsigned int window;
const unsigned short *chain = (unsigned short*) &evt[index]; // pointer to current chain/receiver
for (int ichain = 0; ichain < setup->getMaxUsedChains(); ichain++) {
chainLengthInBytes = chainLength[ichain]*sizeof(short);
padsPerChain = 0; // reset counter for this chain
nextPatternMustBeAmplitude = False;
validBits = bitsPerPad*((chainLengthInBytes*8)/bitsPerPad); // align buffer for n*BitsPerPad
bstream = new CBitStream16(chain, validBits); // make bit stream object out of chain
while (!bstream->eos() && !bstream->bad() && padsPerChain < setup->getMaxPadsPerChain()) {
//
// Be careful at end of bitstream! Reduce size of window to avoid bstream->bad()
//
if (validBits - bstream->count() >= 16)
window = maxCodeLength;
else
window = (unsigned int)(validBits - bstream->count());
//
// Look (peek) at next 16 bits from stream and get encoded value
//
next16Bit = (unsigned int)(bstream->peek(window));
pattern = (unsigned long) bigHuffmanTable[next16Bit] & 0x0000ffff;
//pattern = (unsigned int)(bstream->next(window));
//
// check whether this is a legal huffman code
//
if (pattern == 0x0000ffff) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tERROR\n";
cerr << '\t' << name << " chain " << ichain << " : illegal huffman code." << endl;
break;
}
//
// move pointer by length of codeword
//
bstream->next(huffmanTable[pattern].lengthInBits);
//
// pattern contains # of zero pads
//
if (pattern & 0x100) {
//
// The checks for correct alternation of zero supression words and pad
// amplitudes as well as for chain overflow are reliable methods to
// secure overall data integrity. However, if the problem occurs in the
// last pad of a chain, it can be ignored.
//
if (nextPatternMustBeAmplitude) {
if (padsPerChain < setup->getMaxPadsPerChain()-10) {
if (setup->getWarningLevel() < 5) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tERROR\n";
cerr << '\t' << name << " chain " << ichain << " expected amplitude";
cerr << " but got counter at pad no. " << padsPerChain << ".\n";
}
errorCollection.remember("counterNotExpected", ichain);
}
padsPerChain = setup->getMaxPadsPerChain();
continue;
}
nZero = (pattern == 0x100 ? 256 : pattern & 0xff); // 0 means 256 pads with amp=0
if (nZero < 256) nextPatternMustBeAmplitude = True;
if (padsPerChain+nZero > setup->getMaxPadsPerChain()) {
if (padsPerChain < setup->getMaxPadsPerChain()-10) {
if (setup->getWarningLevel() < 5) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tOverflow at pad " << padsPerChain << " in chain " << ichain;
cerr << " of " << name << endl;
cerr << "\tApplying next decoded counter value of " << nZero << endl;
cerr << "\twould exceed max number of pads per chain by ";
cerr << padsPerChain+nZero-setup->getMaxPadsPerChain() << " pads.\n";
}
errorCollection.remember("chainOverflow", ichain);
}
padsPerChain = setup->getMaxPadsPerChain();
continue;
}
padsPerChain += nZero;
lookupPosition += (int)nZero;
}
else { // pattern contains valid amplitude
nZero = 0;
nextPatternMustBeAmplitude = False;
if (lookupPosition >= maxLookup) {
if (setup->getWarningLevel() < 10) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tCan't lookup xy for pad in " << name << ".\n";
cerr << "\tRequest exceeds size of lookup table.\n";
}
errorCollection.remember("sizeOfLookupTableExceeded");
delete bstream;
return False;
}
//
// Subtract pedestals and check for thresholds.
// Two different kinds of thresholds are applied:
// 1. constant thresholds as defined by minPadAmplitude
// in the setup file
// 2. threshold as defined by nSigmaCut (setup file) times
// the referring sigma of the pad as stored in the data.
// Note the >= below.
//
amplitude = pattern-lookupItem[lookupPosition]->getPedestal();
if (amplitude >= minPadAmplitude &&
amplitude >= nSigmaCut*lookupItem[lookupPosition]->getSigma()) {
if (lookupItem[lookupPosition]->getX() > 0 &&
lookupItem[lookupPosition]->getY() > 0 ) {
CPad *newpad = new CPad(lookupItem[lookupPosition]->getX(),
lookupItem[lookupPosition]->getY(),
amplitude);
if (!pads.insert(newpad)) delete newpad; // cannot be inserted
}
}
padsPerChain++;
lookupPosition++;
}
}
delete bstream;
if (int(validBits)-setup->getMaxPadsPerChain()*bitsPerPad > bitsPerPad) {
if (setup->getWarningLevel() < 5) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tWARNING\n";
cerr << '\t' << "Reached limit of %d pads in chain " << ichain << endl;
cerr << '\t' << "of " << name << " but " << validBits-setup->getMaxPadsPerChain()*bitsPerPad;
cerr << " aligned bits are left in buffer.\n";
}
errorCollection.remember("bitsLeftInChain", ichain);
}
if (padsPerChain < setup->getMaxPadsPerChain()-10 && chainLength[ichain] > 0) {
if (setup->getWarningLevel() < 5) {
cerr << "CRichLikeDetector::unpack():\n";
cerr << "\tWARNING\n";
cerr << "\tExpected " << setup->getMaxPadsPerChain() << " pads in chain " << ichain;
cerr << " of " << name << " but extracted only " << padsPerChain << " pads.\n";
if (nZero > 0)
cerr << "\tLast sequence was ZeroSup with nZero = " << nZero << ".\n";
else
cerr << "\tLast decoded sequence was an amplitude.\n";
}
errorCollection.remember("notEnoughPadsInChain", ichain);
}
lookupPosition = (ichain+1)*setup->getMaxLookupPerChain();
chain += chainLength[ichain];
}
return True;
}
CBoolean CRichLikeDetector::unpackPedestals(CEventServer& server)
{
//
// Get raw data for pedestals
//
CLabel *label = (CLabel *) server.getLabel(setup->getPedestalLabel());
if (label == 0) {
cerr << "CRichLikeDetector::unpackPedestals():\n";
cerr << "\tERROR\n";
cerr << "\tPedestal label for " << name << " not found.\n";
return False;
}
//
// Check if all data are available.
// Minimum minDataLength is measured in short integers (4 extra words).
//
int minDataLength = setup->getMaxLookupPerChain()*setup->getMaxUsedChains()+4;
if (label->getSize()*2 < minDataLength) {
cerr << "CRichLikeDetector::unpackPedestals():\n";
cerr << "\tERROR\n";
cerr << "\tNot enough data for " << name << ": length is " << label->getSize()*2
<< " and should be >= " << minDataLength << endl;
return False;
}
#ifdef ARCH_IS_LITTLE_ENDIAN
// raw data are a 'stream-of-shorts'. swap them for little-endian architectures
if (!label->isSwapped()) label->swap();
#endif
//
// Data is stored in short words
//
short *data = (short *) label->getData();
//
// Check for marker words
//
int marker1 = data[0]; // these are for checking of the (long) <-> (short) matching
int marker2 = data[1];
if (marker1 != 1 || marker2 != 2) {
cerr << "CRichLikeDetector::unpackPedestals():\n";
cerr << "\tERROR\n";
cerr << "\tGot wrong marker words for pedestals of " << name
<< " (" << marker1 << ',' << marker2 << ")"
<< " instead of (1,2).\n";
return False;
}
//
// Read and store data.
// Lower and upper threshold (data[2] and data[3]) are skipped.
//
int i;
for (i=4; i<minDataLength; i++)
lookupItem[i-4]->setPedestal(float(data[i])/10.);
//
// Get raw data for sigma
//
label = (CLabel *) server.getLabel(setup->getSigmaLabel());
if (label == 0) {
cerr << "CRichLikeDetector::unpackPedestals():\n";
cerr << "\tERROR\n";
cerr << "\tSigmaLabel for " << name << " not found.\n";
return False;
}
//
// Check if all data are available.
// Minimum minDataLength is measured in short integers (6 extra words).
//
minDataLength = setup->getMaxLookupPerChain()*setup->getMaxUsedChains()+6;
if (label->getSize()*2 < minDataLength) {
cerr << "CRichLikeDetector::unpackXYLookupTable():\n";
cerr << "\tERROR\n";
cerr << "\tNot enough data for " << name << ": length is " << label->getSize()*2
<< " and should be >= " << minDataLength << endl;
return False;
}
#ifdef ARCH_IS_LITTLE_ENDIAN
// raw data are a 'stream-of-shorts'. swap them for little-endian architectures
if (!label->isSwapped()) label->swap();
#endif
//
// Data is stored in short words
//
data = (short *) label->getData();
//
// Check for marker words
//
marker1 = data[0]; // these are for checking of the (long) <-> (short) matching
marker2 = data[1];
if (marker1 != 1 || marker2 != 2) {
cerr << "\tCRichLikeDetector::unpackPedestals():\n";
cerr << "\tERROR\n";
cerr << "\tGot wrong sigma marker words for " << name
<< " (" << marker1 << "," << marker2 << ")"
<< " instead of (1,2)\n";
return False;
}
//
// Read and store data.
// Skip data[2] - data[5] (threshold mode, used sigmas etc.)
//
for (i=6 ; i<minDataLength ; i++)
lookupItem[i-6]->setSigma(float(data[i])/10.);
return True;
}
CBoolean CRichLikeDetector::unpackXYLookupTable(CEventServer& server)
{
//
// Get raw data
//
CLabel *label = (CLabel *) server.getLabel(setup->getXYLookupLabel());
if (label == 0) {
cerr << "CRichLikeDetector::unpackXYLookupTable():\n";
cerr << "\tERROR\n";
cerr << "\tXYLookupLabel for " << name << " not found.\n";
return False;
}
//
// Check if enough data is available.
// minDataLength is measured in long integers (short x + short y).
//
int minDataLength = setup->getMaxLookupPerChain()*setup->getMaxUsedChains();
if (label->getSize() < minDataLength) {
cerr << "CRichLikeDetector::unpackXYLookupTable():\n";
cerr << "\tERROR\n";
cerr << "\tNot enough data for " << name << ": length is " << label->getSize()
<< " and should be >= " << minDataLength << endl;
return False;
}
#ifdef ARCH_IS_LITTLE_ENDIAN
// raw data are a 'stream-of-shorts'. swap them for little-endian architectures
if (!label->isSwapped()) label->swap();
#endif
//
// Read and store xy lookup items, clear old collection
//
short *data = (short *) label->getData();
lookupItem.clearAndDestroy();
for (int i=0; i<minDataLength*2; i+=2) {
if (!lookupItem.insert(new CRichLikeLookupItem(data[i], data[i+1]))) {
cerr << "CRichLikeDetector::unpackXYLookupTable():\n";
cerr << "\tFATAL ERROR\n";
cerr << "\tCouldn't insert xy lookup item in list.\n";
cerr << "\tRejected item: (" << data[i] << ',' << data[i+1] << ")\n";
cerr << "\tLookup list cannot be used for decoding. Aborting session.\n";
abort();
return False;
}
}
//
// Mark the sensitive items. This is detector dependent
// and is influenced by spokes borders etc..
//
markSensitiveArea();
return True;
}
CBoolean CRichLikeDetector::unpackHuffmanTable(CEventServer& server)
{
//
// Get raw data
//
CLabel *label = (CLabel *) server.getLabel(setup->getHuffmanLabel());
if (label == 0) {
cerr << "CRichLikeDetector::unpackHuffmanTable():\n";
cerr << "\tERROR\n";
cerr << "\thuffmanLabel for " << name << " not found.\n";
return False;
}
#ifdef ARCH_IS_LITTLE_ENDIAN
label->revertAsciiData();
#endif
//
// Since the huffman label is written in
// ascii we have to use an strstream here.
//
istrstream ist((char*) label->getData(), label->getSize()*sizeof(CWord));
//
// Fill Huffmann table from label data.
//
int dummy;
if (!huffmanTable) huffmanTable = new CHuffmanCode[512];
for (int i=0; i<512; i++)
ist >> dummy >> huffmanTable[i].lengthInBits >> huffmanTable[i].code >> hex >> dummy >> dec;
//
// Create unambigous lookup table.
//
if (!expandHuffmanTable())
return False;
return True;
}
//-----------------------------------------------------------------------------
// Read Huffman table from file
// Format of the (ascii) file must be:
// ADC value no. of bits Huffman code(dec) Huffman code(hex)
//-----------------------------------------------------------------------------
CBoolean CRichLikeDetector::readHuffmanTableFromFile(const char* filename)
{
//
// Open file.
//
Cifstream ifs(filename);
if (!ifs) {
cerr << "CRichLikeDetector::readHuffmanTableFromFile():\n";
cerr << "\tERROR\n";
cerr << "\tError opening huffman table file '" << filename << endl;
return False;
}
//
// Fill Huffmann table from file.
//
int dummy;
if (!huffmanTable) huffmanTable = new CHuffmanCode[512];
for (int i=0; i<512; i++)
ifs >> dummy >> huffmanTable[i].lengthInBits >> huffmanTable[i].code >> hex >> dummy >> dec;
ifs.close();
//
// Create unambigous lookup table.
//
if (!expandHuffmanTable())
return False;
return True;
}
CBoolean CRichLikeDetector::expandHuffmanTable()
{
//
// check existence of huffman table
//
if (!huffmanTable) return False;
//
// Create unambigous lookup table.
//
if (!bigHuffmanTable) bigHuffmanTable = new unsigned short[0x10000];
int i;
for(i=0;i<0x10000L;i++) {
bigHuffmanTable[i] = (unsigned short) (0xffff); // clear for redundancy check
}
unsigned long adr;
CBoolean errors = False;
for(i=0;i<512;i++){
if (huffmanTable[i].lengthInBits == 0 || huffmanTable[i].lengthInBits > 16) {
cerr << "<E> in huffman table: lengthInBits = " << huffmanTable[i].lengthInBits
<< " for i= " << i << endl;
errors = True;
}
else {
//
//------------------------------------------------------------------
// fill lookup table with all possible completitions of the code word:
// code : 4bit 0010 -> xxxxxxxxxxxx0010
// . |
// . +-> xxx...xxx all possible codes
//
for(adr=huffmanTable[i].code; adr<0x10000L; adr+=1L<<huffmanTable[i].lengthInBits) {
if (bigHuffmanTable[adr] == (unsigned short) (0xffff)) // if lookuptable entry is free
bigHuffmanTable[adr] = (unsigned short)(i & 0x01ff); // write the value to the addresses
else {
cerr << "<E> in Huffman table entry occupied adr= "
<< hex << adr << " value " << bigHuffmanTable[adr] << dec << endl;
errors = True;
}
}
}
}
if (errors)
return False;
return True;
}
CBoolean CRichLikeDetector::readModulesFromFile(const char* filename)
{
Cifstream ifs(filename);
if (!ifs) {
cerr << "CRichLikeDetector::readModulesFromFile():\n";
cerr << "\tERROR\n";
cerr << "\tError opening module map file '" << filename << ".\n";
return False;
}
int x, y, w, h;
while (ifs.good() && !ifs.eof()) {
ifs >> x >> y >> w >> h;
modules.insert(new CModule(x, y, w, h));
}
return True;
}
CBoolean CRichLikeDetector::readGainmapFromFile(const char* filename)
{
//
// The gainmap contains the following entries (binary format):
// (short x) (short y) (float factor)
// Only entries with factors != 1 are listed
// The calibration factors get stored in the CRichLikeLookupItem list.
//
if (lookupItem.length() == 0) {
cerr << "CRichLikeDetector::readGainmapFromFile():\n";
cerr << "\tERROR\n";
cerr << "\tThe xy lookup table does not exist. Cannot store gainmap." << endl;
return False;
}
Cifstream ifs(filename);
if (!ifs) {
cerr << "CRichLikeDetector::readGainmapFromFile():\n";
cerr << "\tERROR\n";
cerr << "\tError opening gainmap file '" << filename << ".\n";
return False;
}
short x, y;
float factor;
CRichLikeLookupItem *item;
while (ifs.good() && !ifs.eof()) {
ifs.read((char *)&x, sizeof(short));
ifs.read((char *)&y, sizeof(short));
ifs.read((char *)&factor, sizeof(float));
if (item = lookupItem(x, y))
item->setGainCalibrationFactor(factor);
}
return True;
}
CBoolean CRichLikeDetector::readModuleGainmapFromFile(const char* filename)
{
Cifstream ifs(filename);
if (!ifs) {
cerr << "CRichLikeDetector::readGainFromFile():\n";
cerr << "\tERROR\n";
cerr << "\tError opening module gain file '" << filename << ".\n";
return False;
}
float gain;
int i=0;
while (ifs.good() && !ifs.eof()) {
ifs >> gain;
if(i < modules.length() ) {
modules(i)->setGain(gain);
i++;
}
else {
cerr << "CRichLikeDetector::readGainFromFile():\n";
cerr << "\tERROR\n";
cerr << "\tMore values in gain file than existing modules" << ".\n";
return False;
}
}
return True;
}
int CRichLikeDetector::countNeighbours(CPad& pad, float threshold)
{
int padsAboveThreshold = 0;
int ix = pad.getX();
int iy = pad.getY();
int lx,ly;
for ( lx=ix-1; lx<ix+2; lx++ ) {
for ( ly=iy-1; ly<iy+2; ly++ ) {
if ( pads(lx,ly) == 0 ) continue;
if ( (pads(lx,ly)->getAmp()>threshold) ) padsAboveThreshold++;
}
}
return padsAboveThreshold;
}
void CRichLikeDetector::printHitSummary(ostream& ost)
{
ost << "\nRichLikeDetector Hit Summary\n";
CString line('=',46);
ost << line << endl;
for (int i=0; i<hits.length(); i++) {
ost << i << "\tx\t" << hits(i)->getX() << "\ty\t" << hits(i)->getY()
<< "\tamp\t" << hits(i)->getAmp() << endl;
}
}
CBoolean CRichLikeDetector::makeClusters(float minAmp)
{
enum padstatus {UNUSED, USED};
RWTValOrderedVector<CPad*> newpads; // this vector contains the sorted pads afterwards
RWTStack<CPad*, RWTValOrderedVector<CPad*> > padStack; // stack for amp>=minAmp
RWTStack<CPad*, RWTValOrderedVector<CPad*> > lowPadStack; // stack for amp<minAmp
const float startAmp = 300;
if ( clusters.length()>0 ) clusters.clearAndDestroy(); // remove all old cluster
int i;
for ( i=0; i<pads.length(); i++ ) {
pads(i)->setTmp(UNUSED);
}
for ( i=0; i<pads.length(); i++ ) {
if ( pads(i)->getAmp()>=minAmp && pads(i)->getTmp()==UNUSED ) {
clusters.append(new CRichLikeCluster);
clusters.last()->setFirstPad(newpads.length()); // get entry to new sorted list
clusters.last()->setMinAmpLocMax(startAmp);
padStack.push(pads(i)); // put pad onto the stack
pads(i)->setTmp(USED);
while ( !padStack.isEmpty() ) {
newpads.insert(padStack.pop()); // append ??? change CCollection....
newpads.last()->setIsLocalMax(True); // default
newpads.last()->setClusterNumber(clusters.length()-1);
int ix = newpads.last()->getX();
int iy = newpads.last()->getY();
float amp = pads(ix,iy)->getAmp();
int ly = iy;
int lx;
for ( lx=ix-1; lx<ix+2; lx+=2 ) { // loop neighbours in x
if ( pads(lx,ly) == 0 ) continue; // pad must exist !!!!
if ( (pads(lx,ly)->getAmp()>amp) ||
(pads(lx,ly+1) != 0 && // check diagonal
pads(lx,ly+1)->getAmp()>amp ) || // connected neighbours
(pads(lx,ly-1) != 0 &&
pads(lx,ly-1)->getAmp()>amp ) ) {
newpads.last()->setIsLocalMax(False);
}
if ( pads(lx,ly)->getAmp() >= minAmp &&
pads(lx,ly)->getTmp() == UNUSED ) {
padStack.push(pads(lx,ly));
pads(lx,ly)->setTmp(USED);
}
} // end loop neighbours in x
lx = ix;
for ( ly=iy-1; ly<iy+2; ly+=2 ) { // loop neighbours in y
if ( pads(lx,ly) == 0 ) continue; // pad must exist !!!!
if ( (pads(lx,ly)->getAmp()>amp) ||
(pads(lx+1,ly) != 0 && // check diagonal
pads(lx+1,ly)->getAmp()>amp ) || // connected neighbours
(pads(lx-1,ly) != 0 &&
pads(lx-1,ly)->getAmp()>amp ) ) {
newpads.last()->setIsLocalMax(False);
}
if ( pads(lx,ly)->getAmp() >= minAmp &&
pads(lx,ly)->getTmp() == UNUSED ) {
padStack.push(pads(lx,ly));
pads(lx,ly)->setTmp(USED);
}
} // end loop neighbours in y
//
// now update all cluster properties
//
clusters.last()->incrNPads();
clusters.last()->updateAmp(newpads.last()->getAmp());
if ( newpads.last()->getIsLocalMax() )
clusters.last()->incrNLocMax();
if ( newpads.last()->getIsLocalMax() &&
newpads.last()->getAmp()<clusters.last()->getMinAmpLocMax() )
clusters.last()->setMinAmpLocMax((int)newpads.last()->getAmp());
} // end while !padstack.isEmpty()
clusters.last()->rms2Calc(); // since all pads are collected, update rms2 here....
} else if ( pads(i)->getAmp() < minAmp ) {
lowPadStack.push(pads(i)); // collect lowPads here
} // end if ( pads(i)->getAmp()>=minAmp && ...
} // end loop all pads
//
// clean up programming utilities
//
//concatenate lists
//
while ( !lowPadStack.isEmpty() ) {
newpads.insert(lowPadStack.pop());
newpads.last()->setClusterNumber(-1);
}
//
// finally update the pointer list in the collection class
//
if ( pads.length()!= newpads.length() ) {
cerr << "CRichLikeDetector::makeClusters():\n";
cerr << "\t FATAL ERROR\n";
cerr << "length of pointer lists mixed up!\n";
cerr << "Aborting program now\n";
abort();
} else {
for (i=0;i<newpads.length();i++) {
pads(i) = newpads(i);
}
}
return True;
}
void CRichLikeDetector::printClusterSummary()
{
if ( clusters.length()==0 ) {
cerr << "CRichLikeDetector::printClusterSummary():\n";
cerr << "\tWARNING\n";
cerr << "\tCluster list is empty.\n";
} else {
cout << "Cluster Properties" << endl;
CString line('=', 46);
cout << line << endl;
cout << "cluster npads firstPad nlocmax minAmpLocMax " << endl;
for ( int i=0; i<clusters.length(); i++ ) {
cout << i << "\t"
<< clusters(i)->getNPads() << "\t"
<< clusters(i)->getFirstPad() << "\t"
<< clusters(i)->getNLocMax() << "\t"
<< clusters(i)->getMinAmpLocMax()<< "\t" << endl;
}
}
}
CBoolean CRichLikeDetector::splitClusters(int minPad )
{
CBoolean oneClusterDidSplit = True;
while (oneClusterDidSplit) {
oneClusterDidSplit = False;
for ( int i=0; i<clusters.length(); i++ ) {
if ( clusters(i)->getNPads()>minPad && clusters(i)->getNLocMax()>1 ) {
if ( clusters(i)->split(clusters, pads, i)>1 ) oneClusterDidSplit = True;
}
}
}
return True;
}
CBoolean CRichLikeDetector::makeHitsFromClusters()
{
int ioff,iend;
CRichLikeHit *newHit;
CBoolean makeHitFromCluster;
hits.clearAndDestroy();
for ( int i=0; i<clusters.length(); i++) {
if ( clusters(i)->getNPads()<=0 ) continue;
makeHitFromCluster = False;
ioff = clusters(i)->getFirstPad();
iend = ioff+clusters(i)->getNPads();
if (clusters(i)->getNLocMax()==1) { // only one locmax
makeHitFromCluster = True;
} else if (clusters(i)->getNLocMax()==2) { // check if just neighbored
if (clusters(i)->getNPads()<=3) {
makeHitFromCluster = True;
}else{
float dist2 = 0;
int j,k;
for ( j=ioff; j<iend; j++) {
if ( !pads(j)->getIsLocalMax() ) continue;
for (k=j; k<iend; k++) {
if ( !pads(k)->getIsLocalMax() ) continue;
dist2 = ((pads(j)->getX()-pads(k)->getX())*
(pads(j)->getX()-pads(k)->getX()) +
(pads(j)->getY()-pads(k)->getY())*
(pads(j)->getY()-pads(k)->getY()) );
}
if (dist2>0) break;
}
if ( dist2 <= 2 ) { // treat as one
makeHitFromCluster = True;
}else{ // check if one is at
if ( (countNeighbours(*pads(j),0)==6 && // cluster border
countNeighbours(*pads(k),0)<=5 ) ||
(countNeighbours(*pads(j),0)<=5 &&
countNeighbours(*pads(k),0)==6 ) ) {
makeHitFromCluster = True;
}
}
}
}
else {
//
// other cases we will examine later, since they are at percent level...
//
}
if (makeHitFromCluster) {
float sum=0, x=0, y=0, x2=0, y2=0;
for (int j=ioff; j<iend; j++) {
sum = sum+pads(j)->getAmp();
x = x+pads(j)->getAmp()*pads(j)->getX();
y = y+pads(j)->getAmp()*pads(j)->getY();
x2 = x2+pads(j)->getAmp()*pads(j)->getX()*pads(j)->getX();
y2 = y2+pads(j)->getAmp()*pads(j)->getY()*pads(j)->getY();
}
newHit = new CRichLikeHit;
newHit->setCluster(clusters(i));
newHit->setX(x/sum);
newHit->setY(y/sum);
newHit->setAmp(sum);
newHit->setRms2((x2+y2-x*x-y*y)/sum);
hits.insert(newHit);
}
} //enddo cluster_loop
return True;
}
void CRichLikeDetector::killModule(int i)
{
//
// Mark all pads of module killed by setting
// the referring pad amplitudes negative
//
int x, y;
CPad *pad;
if (i < 0 || i > modules.length()) {
cerr << "CRichLikeDetector::killModule():\n";
cerr << "\tERROR\n";
cerr << "\tModule id " << i << " is out of range. Nothing killed.\n";
return;
}
//
// Remember x,y of module is upper left pad
//
for (x = modules(i)->getX() ; x < modules(i)->getX()+modules(i)->getWidth() ; x++) {
for (y = modules(i)->getY()-modules(i)->getHeight()+1; y <= modules(i)->getY(); y++) {
pad = pads(x, y);
if (pad) pad->setAmp( -fabs(pad->getAmp()) );
}
}
}
void CRichLikeDetector::killPads(int x1, int y1, int x2, int y2)
{
//
// Mark pad killed by setting the pad amplitude negative
//
int x, y;
CPad *pad;
for (x = x1; x <= x2; x++) {
for (y = y1; y <= y2; y++) {
pad = pads(x, y);
if (pad) pad->setAmp( -fabs(pad->getAmp()) );
}
}
}
CBoolean CRichLikeDetector::makeCogAndClusterHits()
{
const int minPadNumberForCogAlgorithm = 10;
enum { useThis=99, alreadyUsed=101 };
int indexFirst, indexLast, indexOfLocalMax;
int k, j;
float sum, x, y, x2, y2;
CRichLikeHit *newHit, *lastHit;
hits.clearAndDestroy();
for (int i=0; i<clusters.length(); i++) {
if ( clusters(i)->getNPads() <= 0 ) continue; // cluster was removed
//
// Get index range of pads referring to the cluster
//
indexFirst = clusters(i)->getFirstPad();
indexLast = indexFirst+clusters(i)->getNPads();
//
// Create the hits.
// Small and large clusters are treated differently.
//
if (clusters(i)->getNLocMax() == 1 &&
clusters(i)->getNPads() < minPadNumberForCogAlgorithm ) {
//
// Case I:
// Small cluster with less than minPadNumberForCogAlgorithm
// can be easily treated without the CoG algorithm to speed
// up the hit finding
//
x = y = x2 = y2 = sum = 0;
for (j=indexFirst; j<indexLast; j++) {
sum = sum+pads(j)->getAmp();
x = x+pads(j)->getAmp()*pads(j)->getX();
y = y+pads(j)->getAmp()*pads(j)->getY();
x2 = x2+pads(j)->getAmp()*pads(j)->getX()*pads(j)->getX();
y2 = y2+pads(j)->getAmp()*pads(j)->getY()*pads(j)->getY();
}
if (sum > 0) {
x /=sum; y /=sum; x2 /=sum; y2 /=sum;
}
else
x = y = x2 = y2 = 0;
newHit = new CRichLikeHit;
newHit->setCluster(clusters(i));
newHit->setX(x);
newHit->setY(y);
newHit->setAmp(sum);
newHit->setRms2(x2+y2-x*x-y*y);
hits.insert(newHit);
}
else {
//
// Case II:
// The good old Cog algorithm for larger clusters and those
// with several local maxima.
//
indexOfLocalMax = indexFirst;
for (k=indexFirst; k<indexLast; k++) {
pads(k)->setTmp(useThis);
if (pads(k)->getIsLocalMax()) indexOfLocalMax = k;
}
//
// Local maximum comes first
//
if (pads(indexOfLocalMax)->getTmp() == useThis) {
if (lastHit = findOneHit(pads(indexOfLocalMax), i)) {
//
// Mark pads in 3x3 around center pad as 'used'
//
for (j = indexOfLocalMax+1; j < indexLast; j++ ) {
if (abs(pads(j)->getX()-lastHit->getX()) <= 1 &&
abs(pads(j)->getY()-lastHit->getY()) <= 1)
pads(j)->setTmp(alreadyUsed);
}
}
}
//
// Now get the rest
//
for (k = indexFirst; k < indexLast; k++) {
if (k == indexOfLocalMax) continue; // already done
if (pads(k)->getTmp() == useThis) {
if (lastHit = findOneHit(pads(k), i)) {
//
// Mark pads in 3x3 around center pad as 'used'
//
for (j = k+1; j < indexLast; j++ ) {
if (abs(pads(j)->getX()-lastHit->getX()) <= 1 &&
abs(pads(j)->getY()-lastHit->getY()) <= 1)
pads(j)->setTmp(alreadyUsed);
}
}
}
}
}
//
// Avoid second usage by neighboured cluster
//
for (j=indexFirst; j<indexLast; j++)
pads(j)->setTmp(alreadyUsed);
}
return True;
}
void CRichLikeDetector::findCenterOfGravity(int ix, int iy, float& x, float& y, float& ampSum,
float& x2, float& y2, CPad*& ampMaxPad)
{
int lx, ly;
float amp = 0;
float ampMax = 0;
x = y = x2 = y2 = ampSum = 0;
ampMaxPad = 0;
for ( lx=ix-1; lx<=ix+1; lx++ ) {
for ( ly=iy-1; ly<=iy+1; ly++ ) {
if ( pads(lx,ly) == 0 || pads(lx,ly)->getAmp() <= 0 )
continue;
else {
amp = pads(lx,ly)->getAmp();
ampSum += amp;
x += amp*lx;
y += amp*ly;
x2 += amp*lx*lx;
y2 += amp*ly*ly;
if ( amp > ampMax ) {
ampMaxPad = pads(lx,ly);
ampMax = amp;
}
}
}
}
if (ampSum>0) {
x /=ampSum;
y /=ampSum;
x2 /=ampSum;
y2 /=ampSum;
}
else
x = y = x2 = y2 = 0;
}
CRichLikeHit* CRichLikeDetector::findOneHit(const CPad* startPad, int clusterNumber)
{
const int MaxIterations = 7;
CPad* ampMaxPad = 0;
CBoolean hitTooNear;
CBoolean tryAgain=True;
float x, y, x2, y2, ampSum;
CRichLikeHit *newHit = 0;
x = y = x2 = y2 = ampSum = 0;
int ix = startPad->getX();
int iy = startPad->getY();
float amp = startPad->getAmp();
for (int count=0; (tryAgain && count < MaxIterations); count++ ) {
tryAgain = False;
findCenterOfGravity(ix, iy, x, y, ampSum, x2, y2, ampMaxPad);
if (int(x+0.5) != ix) {
ix = int(x+0.5);
tryAgain = True;
}
if (int(y+0.5) != iy) {
iy = int(y+0.5);
tryAgain = True;
}
}
//
// Check if the hit make sense
// If the center pad is not a local maximum the pad must
// have at least 'hitMinCenterPadPercentage' of the total
// amplitude of the hit (~10-20%).
//
if (!ampMaxPad) return False;
if (!pads(ix,iy)) return False;
if (!pads(ix,iy)->getIsLocalMax())
if (setup->getHitMinCenterPadPercentage()*ampSum > pads(ix,iy)->getAmp()) return False;
//
// Check if too close to already existing one
//
hitTooNear = False;
for (int j=0; j<hits.length(); j++) {
if (abs(hits(j)->getX()-x) < 1 && abs(hits(j)->getY()-y) < 1) {
hitTooNear = True;
break;
}
}
//
// Create a new entry in the hit list
//
if (!hitTooNear) {
newHit = new CRichLikeHit;
newHit->setCluster(clusters(clusterNumber));
newHit->setX(x);
newHit->setY(y);
newHit->setAmp(ampSum);
newHit->setRms2(x2+y2-x*x-y*y);
hits.insert(newHit);
}
return newHit;
}
double CRichLikeDetector::getAverageHitAmplitude()
{
double sum = 0;
for (int i=0; i<hits.length(); i++) sum += hits(i)->getAmp();
return sum/hits.length();
}
//----------------------------------------------------------------------------------------
//
// Dummy function (temporary only)
//
//----------------------------------------------------------------------------------------
CBoolean CRichLikeDetector::makeHitsFromMC(CMCServer&)
{
cerr << "CRichLikeDetector::makeHitsFromMC():\n";
cerr << "\tWARNING\n";
cerr << "\tDummy version of makeHitsFromMC() called. No action taken.\n";
return False;
}
double erf2(float x0, float y0, int ix, int iy, float sighit)
{
/**********************************************************
* *
* Returns value of integrated 2-dim Gaussian with *
* means XMEAN,YMEAN around (ix/iy) +/- 0.5 and *
* sigma SIGMA. *
* *
* original version TU *
* *
**********************************************************/
const float sqrt2 = 1.4142135624;
double x1 = double(ix - 0.5);
double x2 = double(ix + 0.5);
double y1 = double(iy - 0.5);
double y2 = double(iy + 0.5);
return ( ( erf((x2-x0)/(sqrt2*sighit))-
erf((x1-x0)/(sqrt2*sighit)) ) *
( erf((y2-y0)/(sqrt2*sighit))-
erf((y1-y0)/(sqrt2*sighit)) )/4. );
}
double invers_cosh(float x0, float y0, int ix, int iy, float fwhm)
{
/**********************************************************
* *
* Returns value of integrated 2-dim 1/cosh with *
* means XMEAN,YMEAN around (ix,iy) +/- 0.5 and *
* full width half maximum of fwhm. *
* *
**********************************************************/
float L = fwhm/1.68; // dipl thesis S.Tapprogge S.69
float fac = Pi/(2*L);
double x1 = double(ix - 0.5);
double x2 = double(ix + 0.5);
double y1 = double(iy - 0.5);
double y2 = double(iy + 0.5);
return ( (atan(exp(fac*(x2-x0)))-
atan(exp(fac*(x1-x0)))) *
(atan(exp(fac*(y2-y0)))-
atan(exp(fac*(y1-y0)))) / (Pi/2) / (Pi/2));
}
CBoolean CRichLikeDetector::digitize(CMCServer& mcServ, CMCDigiMode digiMode)
{
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();
const double noiseFactor = setup->getMCNoiseFactor();
const double pedShift = setup->getMCPedestalShift();
float minPadAmplitude = (float) setup->getMinPadAmplitude();
float nSigmaCut = (float) setup->getNSigmaCut();
if ( digiMode == CMCResetMode ) {
pads.clearAndDestroy();
}
const RWTPtrOrderedVector<CMCDigiHit> *digiHits;
float x, y, amp;
int ix, iy;
float hit[7][7];
double sig, camexError;
float gaincal, noise;
CRandom rndmGen;
CPad *pad;
int i, ii, jj;
int id_module;
digiHits = mcServ.getMCDigiHits( detId );
pads.resize( (pads.length() + digiHits->length()*10) ); // assume 10 pads per hit
for (i=0; i<digiHits->length(); i++) {
x = (*digiHits)(i)->getDigiHit().getX();
y = (*digiHits)(i)->getDigiHit().getY();
amp = (*digiHits)(i)->getDigiHit().getAmp();
//
// get gainfactor
//
gaincal = 1.;
for (id_module=0; id_module<modules.length(); id_module++){
if(modules(id_module)->isInModule((int) x,(int) y)) {
gaincal = modules(id_module)->getGain();
break;
}
}
amp = amp/gaincal; // local gain calibration
//
// next is to distribute the hit over a 7x7 array of pads:
//
for (ii=0; ii<7; ii++) {
ix = int( x + 0.5 ) + ii - 3;
for (jj=0; jj<7; jj++) {
iy = int( y + 0.5 ) + jj - 3;
hit[ii][jj] = (amp * invers_cosh( x, y, ix, iy, sigHit)); // use 2-d 1/cosh-function
camexError = hit[ii][jj] * rndmGen.gauss() * noiseFactor;
hit[ii][jj] = hit[ii][jj] + camexError; // error caused by camex
}
}
//
// and here we add the 7x7 array to the pad-array:
//
for (ii=0; ii<7; ii++) {
ix = int( x + 0.5 ) + ii - 3;
for (jj=0; jj<7; jj++) {
iy = int( y + 0.5 ) + jj - 3;
if ( ix >= 1 && ix <= NPADS &&
iy >= 1 && iy <= NPADS ) { // within array ...
if ( lookupItem(ix,iy) == 0 ) continue; // no entry for this x/y
if ( !lookupItem(ix,iy)->getIsSensitive() ) continue; // ignore pads under spokes and outside acceptance
noise = 0; // pre-set
if ( simpleNoise ) { // simple noise mode requested ...
sig = rndmGen.gauss(); // get normal distributed random number
noise = sig * meanSigma; // ... set noise
} else if ( realNoise ) {
sig = rndmGen.gauss(); // get normal distributed random number
noise = sig * lookupItem(ix, iy)->getSigma();
// set noise according to SOR info
}
if ( pads(ix,iy) != 0 ) { // pad already existing (i.e. above threshold !!)
amp = int( hit[ii][jj] + pads(ix,iy)->getAmp());
amp = (amp > 255) ? 255 : amp;
pads(ix,iy)->setAmp(amp); // update amplitude of pad
} else { // pad does not exist
amp = int( hit[ii][jj] + noise + pedShift);
amp = (amp > 255) ? 255 : amp;
if (amp >= minPadAmplitude && // check if amplitude is above threshold
amp >= nSigmaCut*lookupItem(ix,iy)->getSigma()) { // as in unpack() ...
pad = new CPad(ix, iy, amp);
if (pad == 0) {
cerr << "CRichLikeDetector::digitize():\n"
<< "\tFATAL\n"
<< "\tcannot allocate new pad for x/y/a = " << ix << "/" << iy << "/" << amp << endl
<< "\taborting." << endl;
abort();
}
pad->addMCDigiHit( (*digiHits)(i) );
if ( !pads.insert(pad) ) delete pad;
}
} // if pad exists
} // x/y is within limits of array
}
}
} // loop mchits
//
// Next is to add noise to all pads (if requested). This is a very costly part of the simulation,
// since we do this for about 65000 pads (ignoring the 10 % which have been 'hit' before).
//
int nnois = 0;
float lupSig;
if ( (digiMode == CMCResetMode) && // don't make any noise in overlay mode
(realNoise || simpleNoise) ) { // fill the other pads according to their sigma
sig = 0;
for (i=0; i<lookupItem.length(); i++) {
ix = int(lookupItem(i)->getX()); if (ix <= 0) continue; // can't be > NPADS
iy = int(lookupItem(i)->getY()); if (iy <= 0) continue;
if (pads(ix, iy) == 0) { // fill noise only for empty pads
if (realNoise)
lupSig = lookupItem(ix, iy)->getSigma();
else
lupSig = meanSigma;
sig = rndmGen.gauss(); // get normal distributed random number
amp = float(sig) * lupSig; // set noise according to SOR info
if (amp >= minPadAmplitude && // check if amplitude is above threshold
amp >= nSigmaCut*lupSig) { // as in unpack() ...
nnois++;
pad = new CPad(ix, iy, amp);
if (pad == 0) {
cerr << "CRichLikeDetector::digitize():\n"
<< "\tFATAL\n"
<< "\tcannot allocate new pad for x/y/a = "
<< ix << "/" << iy << "/" << amp << endl
<< "\taborting." << endl;
abort();
}
if ( !pads.insert(pad) ) delete pad;
}
}
} // end for i<lookupItem.length()
}
return True;
} // end CRichLikeDetector::digitize