CSidc.C
//-----------------------------------------------------------------------------
// $Header: /asis/offline/ceres/cool/project/RCS/CSidc.C,v 3.28 1997/10/15 16:33:24 messer Exp $
//
// COOL Program Library
// Copyright (C) CERES collaboration, 1996
//
// Implementation of class CSidc.
//
//-----------------------------------------------------------------------------
#include <iostream.h>
#include <strstream.h>
#include <fstream.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <limits.h>
#include "CDataStream.h"
#include "CSidc.h"
#include "CTrack.h"
#include "CRandom.h"
#include "CSortedListIterator.h"
#include "CEventServer.h"
#include "CCoordinate.h"
#include "CSidcDeltaR.h"
#include "CPulseFitter.h"
#include "CMCServer.h"
#include "CMCDigiHit.h"
#include "CMCParticle.h"
CSidc::CSidc()
{
lookupSidcAnode = 0;
readoutThreshold = 0;
stopPulsePhaseCorrection = 0;
for (int i=0; i<4; i++) {
stopPulseCoG[i] = 0;
stopPulseCoGOffset[i] = 0;
stopPulsePedestal[i] = 0;
}
}
CSidc::~CSidc()
{
delete [] lookupSidcAnode;
lookupItems.clearAndDestroy();
cells.clearAndDestroy();
nims.clearAndDestroy();
}
CBoolean CSidc::unpack(CEventServer &server)
{
int nAnodes = setup->getNAnodes();
int nTimeBins = setup->getNTimeBins();
int nNim = setup->getNNim();
cells.clearAndDestroy(); // clear old collection
nims.clearAndDestroy();
errorCollection.newEvent();
//
// Create the anode lookup table if it doesn't already exist
//
int j, k, nch;
if (lookupSidcAnode == 0) {
lookupSidcAnode = new int[nAnodes+nNim];
for(k=1; k<=4; k++) {
nch = 0;
for (j=1; j<=96; j++) {
if (j%16 == 0) {
lookupSidcAnode[(k-1)*96+j-1] = nAnodes+nch*4+k;
nch++;
} else
lookupSidcAnode[(k-1)*96+j-1] = j+(k-1)*90-nch;
}
}
}
//
// Unpack pedestals if not already loaded.
// The last SoR must be present.
//
const CRawSoR *SoR = server.getSoR();
if (SoR == 0) {
cerr << "CSidc::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tSoR label not found. No further action.\n";
return False;
}
if (lookupItems.length() == 0 || (lastRunNumber != SoR->getRunNumber())) {
if (!unpackPedestals(server)) {
cerr << "CSidc::unpack():\n";
cerr << "\tWARNING\n";
cerr << "\tFailed to unpack pedstals for "
<< name << ". Continue despite warnings.\n";
}
lastRunNumber = SoR->getRunNumber();
}
//
// Get raw data
//
CLabel *label = (CLabel *) server.getLabel(setup->getDataLabel());
if (label == 0) {
cerr << "CSidc::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tData label for " << name << " not found.\n";
return False;
}
//
// No data, nothing to do
//
if (label->getSize() <= 0) {
cerr << "CSidc::unpack():\n";
cerr << "\tWARNING\n";
cerr << "\tNo data to unpack for " << name << endl;
return False;
}
//
// Set default sizes for collection classes
//
cells.resize((size_t) 4*label->getSize()); // ~ 4 cells/long
//
// Short words are swapped within one data word.
// 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 first sync block
//
int index = 0;
index += (unsigned short)evt[index];
//
// Get the hit counter for each chain.
// Note that the SiDCs contains only 8 valid receiver values,
// although 16 are written.
//
const int MaxReceivers = 16;
const int MaxUsedReceivers = 8;
short *chainLength = (short*) &evt[index];
index += MaxReceivers*sizeof(short)/sizeof(long);
//
// Skip second sync block
//
index += (unsigned short) evt[index];
//
// Check the receiver lengths
//
int len = 0;
int i;
for (i=0; i<MaxUsedReceivers; i++) {
len += chainLength[i];
}
if (len*sizeof(short) > (label->getSize()-index)*sizeof(long)) {
cerr << "CSidc::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tSum of all hit counters exceeds actual data length by %d bytes.\n";
errorCollection.remember("dataLengthExceeded");
return False;
}
//
// Loop receivers and unpack data
//
unsigned short* data = (unsigned short*)&evt[index];
index = 0;
CBoolean preSample = False;
CBoolean markerWord;
int first, tbin, wire, lastTbin, lastWire;
float leftAmp, rightAmp;
int leftAnode, rightAnode;
CCell *cell;
for (k=0; k<MaxUsedReceivers; k++) {
if (k >= setup->getFirstCrate() && k < setup->getFirstCrate()+4) {
first = index;
//
// Get first wire & time bin
//
while (!MarkerSet(data[first]) && first < index+chainLength[k]) first++;
if (first != index) {
if (setup->getWarningLevel() < 5) {
cerr << "CSidc::unpack():\n";
cerr << "\tWARNING\n";
cerr << "\tFirst start-of-pulse marker missing in receiver " << k << ".\n";
}
errorCollection.remember("firstStartOfPulseMarkerMissing", k);
}
if (first >= index+chainLength[k]) {
if (setup->getWarningLevel() < 10) {
cerr << "CSidc::unpack():\n";
cerr << "\tWARNING\n";
cerr << "\tNo start-of-pulse marker found in receiver " << k << ".\n";
}
errorCollection.remember("noStartOfPulseMarkerFound", k);
continue;
}
//
// Loop over data of receiver k
//
lastTbin = lastWire = 0;
for (i=first; i<index+chainLength[k]; i++) {
if (MarkerSet(data[i])) {
//
// Get wire & time bins, check for lots of errors
//
tbin = data[i] & ~(~0 << 8); // bit 0-7
wire = (data[i] >> 8) & ~(~0 << 6); // bit 8-13
if (wire < 0 || wire > 47) {
//
// This means in most cases that all bits are set (wire=63).
// In this case we have probably misinterpreted a corrupted
// amplitude information as a start-of-pulse-word.
//
if (setup->getWarningLevel() < 5) {
cerr << "CSidc::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tWire out of range in receiver " << k << " at word " << i-index;
cerr << " of " << chainLength[k] << "." << endl;
cerr << "\tFound " << wire << " but range is limited to [0-47].\n";
}
errorCollection.remember("wireOutOfRange", k);
wire = lastWire;
tbin = lastTbin+1;
markerWord = False;
}
else
markerWord = True;
}
else
markerWord = False;
if (markerWord) {
if (wire < lastWire) {
if (setup->getWarningLevel() < 10) {
cerr << "CSidc::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tWire number not monotonously growing: "
<< lastWire << "->" << wire << endl;
cerr << "\tin receiver " << k << "at word " << i-index;
cerr << " of " << chainLength[k] << "." << endl;
}
errorCollection.remember("wireNumberSmallerThanLast", k);
}
if (wire > lastWire) lastTbin = 0;
lastWire = wire;
if (tbin < lastTbin) {
if (setup->getWarningLevel() < 10) {
cerr << "CSidc::unpack():\n";
cerr << "\tERROR\n";
cerr << "\ttbin not monotonously growing: " << lastTbin << "->" << tbin << endl;
cerr << "\tin receiver " << k << "at word " << i-index;
cerr << " of " << chainLength[k] << "." << endl;
}
errorCollection.remember("timeBinSmallerThanLast", k);
}
lastTbin = tbin;
if (tbin < 0 || tbin >= nTimeBins) {
if (setup->getWarningLevel() < 10) {
cerr << "CSidc::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tStart of pulse out of range in receiver " << k << ". Found " << tbin
<< " but range is limited to [0-" << nTimeBins-1 << "].\n";
}
errorCollection.remember("startOfPulseOutOfRange", k);
return False;
}
leftAnode = lookupSidcAnode[(wire*2+1)+k%4*96-1];
rightAnode = lookupSidcAnode[(wire*2+2)+k%4*96-1];
preSample = True; // next amplitude should be preSample
} // if (markerWord) ...
else {
//
// Get amplitudes of next time bin
//
leftAmp = (float) ((data[i] & ~(~0 << 6))+1);
rightAmp = (float) ((data[i] >> 8 & ~(~0 << 6))+1);
if (leftAmp > 64. || leftAmp < 1.) {
cerr << "CSidc::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tLeft amplitude out of range in receiver " << k << ". Found " << (int)leftAmp
<< " but range is limited to [1-64].\n";
errorCollection.remember("leftAmplitudeOutOfRange", k);
return False;
}
if (rightAmp > 64. || rightAmp < 1.) {
cerr << "CSidc::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tRight amplitude out of range in receiver " << k << ". Found " << (int)rightAmp
<< " but range is limited to [1-64].\n";
errorCollection.remember("rightAmplitudeOutOfRange", k);
return False;
}
if (preSample) {
leftAmp /= 4.; // data compression: 4 presamples are added
rightAmp /= 4.;
preSample = False;
}
//
// Fill into collections
//
if (leftAnode <= nAnodes) {
cell = new CCell(leftAnode-1, tbin, 200.*(leftAmp-1.)/(259.-3.*leftAmp));
if (!cells.insert(cell)) delete cell;
}
else if (leftAnode <= nAnodes + nNim) {
cell = new CCell(leftAnode-nAnodes-1, tbin, 200.*(leftAmp-1.)/(259.-3.*leftAmp));
if (!nims.insert(cell)) delete cell;
}
if (rightAnode <= nAnodes) {
cell = new CCell(rightAnode-1, tbin, 200.*(rightAmp-1.)/(259.-3.*rightAmp));
if (!cells.insert(cell)) delete cell;
}
else if (rightAnode <= nAnodes + nNim) {
cell = new CCell(rightAnode-nAnodes-1, tbin, 200.*(rightAmp-1.)/(259.-3.*rightAmp));
if (!nims.insert(cell)) delete cell;
}
tbin++; // time bin for next cell
}
}
}
index += chainLength[k]; // position of next receiver
}
//
// Get Camac TDC value and evaluate the stop pulse phase correction.
//
stopPulsePhaseCorrection = 0;
const CLabel *ctdc = server.getLabel(CTDC);
if (ctdc == 0) {
cerr << "CSidc::unpack():\n";
cerr << "\tERROR\n";
cerr << "\tData label CTDC not found. Cannot perform stop pulse correction.\n";
errorCollection.remember("noStopPulseCorrection");
return False;
}
const CWord *ctdcData = ctdc->getData();
CWord ctdcTime = ctdcData[0] & 0x0000ffff;
stopPulsePhaseCorrection = ctdcTime/200.;
//
// Now calculate the center of gravity of the stop pulses
//
for (i=0; i<4; i++) {
stopPulseCoG[i] = 0;
stopPulseCoGOffset[i] = 0;
}
double momstop0, momstop1;
int stopTimeBin;
CCell *nimCell;
for (int icrate = 0; icrate<4; icrate++) {
momstop0 = momstop1 = 0;
for (stopTimeBin = 0; stopTimeBin < nTimeBins; stopTimeBin++) {
nimCell = nims(icrate, stopTimeBin);
if ( (nimCell !=0) && (nimCell->getAmp() >= setup->getStopPulseThreshold()) ) {
momstop0 += nimCell->getAmp();
momstop1 += nimCell->getAmp()*stopTimeBin;
}
}
if (momstop0 > 0)
stopPulseCoG[icrate] = momstop1/momstop0;
else
stopPulseCoG[icrate] = 0;
stopPulseCoG[icrate] += stopPulsePhaseCorrection;
if (stopPulseCoG[icrate] <= setup->getLowerStopPulseCoGCut(lastRunNumber, icrate)) // if cog in lower peak
stopPulseCoGOffset[icrate] += 1;
stopPulseCoGOffset[icrate] += stopPulsePhaseCorrection - // TDC correction
setup->getTimeOffsetBTDC();
}
return True;
}
CBoolean CSidc::unpackPedestals(CEventServer &server)
{
//
// Get raw data for pedestals
//
CLabel *label = (CLabel *) server.getLabel(setup->getPedestalLabel());
if (label == 0) {
cerr << "CSidc::unpackPedestals():\n";
cerr << "\tERROR\n";
cerr << "\tPedestal label for " << name << " not found.\n";
return False;
}
//
// Something fishy here. We will make them anew.
//
if (lookupItems.length() != setup->getNAnodes()) lookupItems.clearAndDestroy();
#ifdef ARCH_IS_LITTLE_ENDIAN
label->revertAsciiData();
#endif
//
// Since (for some unknown reason) pedestal label is written in
// ascii we have to use an strstream here.
//
istrstream ist((char*) label->getData(), label->getSize()*sizeof(CWord));
//
// First get the readout threshold
//
int threshold[2];
ist >> threshold[0] >> threshold[1];
if (id == SiDC1)
readoutThreshold = threshold[0];
else
readoutThreshold = threshold[1];
//
// Read the pedestals
//
int icrate = 0;
int anode = -1;
int crate, channel;
float pedestal, sigma;
int minCrate = setup->getFirstCrate(); // first crate for referring data
for (int i=0; i<8*96; i++) {
ist >> crate >> channel >> pedestal >> sigma;
//
// Use only the information related to this SiDC
//
if (crate >= minCrate && crate < minCrate+4) {
//
// Store pedestal ignoring every 16th (nim) channel.
// If old lookupItems exist, simply overwrite them,
// otherwise create new ones.
//
if ((i+1)%16 != 0) {
anode++;
if (lookupItems.length() < anode+1) {
lookupItems.append(new CSidcLookupItem(crate, channel, pedestal, sigma, 0., 1.0));
}
else {
lookupItems(anode)->setCrate(crate);
lookupItems(anode)->setChannel(channel);
lookupItems(anode)->setPedestal(pedestal);
lookupItems(anode)->setSigma(sigma);
}
}
if (anode+1 > setup->getNAnodes()) {
cerr << "CSidc::unpackPedestals():\n";
cerr << "\tERROR\n";
cerr << "\tIncremented anode to illegal value " << anode << ". No further action.\n";
return False;
}
//
// The following values are the stop pulse pedestals.
// First modulo 16 channels in each crate
//
if ((i+1-16)%96 == 0)
stopPulsePedestal[icrate++] = pedestal;
}
if (ist.bad() || ist.eof()) {
cerr << "CSidc::unpackPedestals():\n";
cerr << "\tERROR\n";
cerr << "\tFailed to read all pedestals for " << name << ".\n";
cerr << "\tExpected " << 8*96 << " entries in file but found only " << i << ".\n";
return False;
}
}
return True;
}
CBoolean CSidc::readLookupItemsFromFile(const char* filename)
{
//
// Open file.
//
Cifstream ifs(filename);
if (!ifs) {
cerr << "CSidc::readLookupItemsFromFile():\n";
cerr << "\tERROR\n";
cerr << "\tError opening calibration file '" << filename << endl;
return False;
}
//
// Fill calibration constants from file.
//
int anode, crate, channel;
float pedestal, sigma, t0Correction, gain;
CSidcLookupItem *newItem;
lookupItems.clearAndDestroy();
for (int i=0; i<setup->getNAnodes(); i++) {
ifs >> anode >> crate >> channel >> pedestal >> sigma >> t0Correction >> gain;
if (anode != i) {
cerr << "CSidc::readLookupItemsFromFile():\n";
cerr << "\tERROR\n";
cerr << "\tWrong anode number " << anode << "found. No further action." << endl;
return False;
}
newItem = new CSidcLookupItem(crate, channel, pedestal, sigma, t0Correction, gain);
lookupItems.append(newItem);
}
ifs.close();
return True;
}
void CSidc::setLookupItemsToDefault()
{
//
// Set calibration constants to default values.
//
int crate = 0;
int channel = 0;
CSidcLookupItem *newItem;
lookupItems.clearAndDestroy();
for (int i=0; i<setup->getNAnodes(); i++) {
newItem = new CSidcLookupItem(crate, channel, 5.0, 1.0, 0.0, 1.0);
lookupItems.append(newItem);
//
// every 16th channel is reserved for stoppulses
//
channel++;
if ((channel+1)%16 == 0) channel++;
//
// one crate has 96 channels
//
if (channel == 96) {
channel = 0;
crate++;
}
}
}
CBoolean CSidc::writeLookupItemsToFile(const char* filename)
{
//
// Open file.
//
Cofstream ofs(filename);
if (!ofs) {
cerr << "CSidc::writeLookupItemsToFile():\n";
cerr << "\tERROR\n";
cerr << "\tError opening calibration file '" << filename << endl;
return False;
}
//
// write calibration constants to file.
//
for (int anode=0; anode<setup->getNAnodes(); anode++) {
ofs << anode << " ";
ofs << lookupItems(anode)->getCrate() << " ";
ofs << lookupItems(anode)->getChannel() << " ";
ofs << lookupItems(anode)->getPedestal() << " ";
ofs << lookupItems(anode)->getSigma() << " ";
ofs << lookupItems(anode)->getT0Correction() << " ";
ofs << lookupItems(anode)->getGain() << " ";
ofs << endl;
}
ofs.close();
return True;
}
//
// Routines required by base class CDrawable to set world coordinates
//
float CSidc::getLeft() const { return -setup->getMaxRadius(); }
float CSidc::getRight() const { return setup->getMaxRadius(); };
float CSidc::getBottom() const { return -setup->getMaxRadius(); }
float CSidc::getTop() const { return setup->getMaxRadius(); };
//
// Draw SiDC detector.
// The following options are recognized:
// n plot the numeric values of the cell amplitudes
// k draw clusters as white boxes
// i plot cluster id of each cell
// h Draw hits as simple crosses
//
void CSidc::draw(const char* copt)
{
CString options(copt);
color(Gray);
clear();
//
// Draw waver in black
//
color(Black);
polyfill(True);
circle(0., 0., setup->getMaxRadius());
if (psMode)
color(White);
else
color(Gray);
circle(0., 0., setup->getMinRadius());
//
// Draw cells and (optional) plot their numeric amplitude.
// Since the color table is setup for RICH pulse heights it
// has to be transformed for the smaller range of the SiDC.
//
CCell *cell;
int newcolor;
float points[5][2];
double dphi = 0.5*ToRadian;
double rmin, rmax, phiLeft, phiRight;
double halfCell = 0.5*setup->getDriftVelocity();
CLabCylinderCoord polarCoord;
CLabXYZCoord xyCoord;
char number[16];
CBoolean showNumbers = options.first('n') != RW_NPOS;
CBoolean showCluster = options.first('i') != RW_NPOS;
if (showNumbers || showCluster) {
font("futura.l");
centertext(True);
textsize(halfCell, halfCell);
polyfill(False);
}
else
polyfill(True);
int i;
for (i = 0; i<cells.length(); i++) {
cell = cells(i);
transform(polarCoord, CSidcCellCoord(cell->getAnode(), cell->getTimeBin()));
newcolor = int(256*cell->getAmp()/200);
if (newcolor < colors.length())
color(colors(newcolor));
else
color(Red);
rmin = polarCoord.getRadius()-halfCell;
rmax = polarCoord.getRadius()+halfCell;
phiLeft = double(polarCoord.getPhi())+dphi;
phiRight = double(polarCoord.getPhi())-dphi;
points[0][0] = rmin*cos(phiLeft);
points[0][1] = rmin*sin(phiLeft);
points[1][0] = rmax*cos(phiLeft);
points[1][1] = rmax*sin(phiLeft);
points[2][0] = rmax*cos(phiRight);
points[2][1] = rmax*sin(phiRight);
points[3][0] = rmin*cos(phiRight);
points[3][1] = rmin*sin(phiRight);
points[4][0] = points[0][0]; points[4][1] = points[0][1];
poly2(5, points);
if (showNumbers) {
sprintf(number, "%d\0", (int)cell->getAmp());
::transform(xyCoord, polarCoord);
move2(xyCoord.getX(), xyCoord.getY());
drawstr(number);
}
else if (showCluster) {
sprintf(number, "%d\0", (int)cell->getCluster());
::transform(xyCoord, polarCoord);
move2(xyCoord.getX(), xyCoord.getY());
drawstr(number);
}
}
//
// Draw hits as simple crosses
//
const float crossWidth = 4 * halfCell;
if (options.first('h') != RW_NPOS) {
polyfill(False);
color(White);
CSidcHit *hit;
for (i=0; i<hits.length(); i++) {
hit = hits(i);
move2(hit->getCenterInXYZ().getX()-crossWidth/2, hit->getCenterInXYZ().getY()-crossWidth/2);
rdraw2(crossWidth, crossWidth);
rmove2(-crossWidth, 0);
rdraw2(crossWidth, -crossWidth);
}
}
//
// draw clusters as white boxes
//
CLabCylinderCoord polarCoordMin, polarCoordMax;
if (options.first('c') != RW_NPOS) {
double almostHalfCell = halfCell * 0.85;
double almostDphi = dphi * 0.85;
polyfill(False);
color(White);
for (i=0; i<clusters.length(); i++) {
transform(polarCoordMin, CSidcCellCoord(clusters(i)->getMinAnode(),
clusters(i)->getMinTBin()));
transform(polarCoordMax, CSidcCellCoord(clusters(i)->getMaxAnode(),
clusters(i)->getMaxTBin()));
rmin = polarCoordMin.getRadius()+almostHalfCell;
rmax = polarCoordMax.getRadius()-almostHalfCell;
if (setup->getIsFlipped()) {
phiLeft = double(polarCoordMax.getPhi())-almostDphi;
phiRight = double(polarCoordMin.getPhi())+almostDphi;
}
else {
phiLeft = double(polarCoordMax.getPhi())+almostDphi;
phiRight = double(polarCoordMin.getPhi())-almostDphi;
}
points[0][0] = rmin*cos(phiLeft);
points[0][1] = rmin*sin(phiLeft);
points[1][0] = rmax*cos(phiLeft);
points[1][1] = rmax*sin(phiLeft);
points[2][0] = rmax*cos(phiRight);
points[2][1] = rmax*sin(phiRight);
points[3][0] = rmin*cos(phiRight);
points[3][1] = rmin*sin(phiRight);
points[4][0] = points[0][0]; points[4][1] = points[0][1];
poly2(5, points);
}
}
//
// Draw tracks from attached track list.
// Hits of tracks are marked by a number referring to the index
// of the track in the list. Hits in the other chamber which belong
// to the same track are referred to by red lines.
//
CTrack *track;
CEventCoord evtCoord;
CLabXYZCoord center, center1, center2, centerNext;
if (options.first('t') != RW_NPOS && trackList) {
font("futura.l");
centertext(True);
polyfill(False);
textsize(4*halfCell, 4*halfCell);
CSortedListIterator<CTrack> titer(*trackList);
while (track = titer()) {
sprintf(number, "%d\0", titer.pos());
//
// Text for match to both Sidcs if any
//
if (track->getTrackMask()&CTrack::hasSidc1Match && track->getTrackMask()&CTrack::hasSidc2Match) {
center1 = track->getClosestSidc1Hit()->getCenterInXYZ();
center2 = track->getClosestSidc2Hit()->getCenterInXYZ();
color(Red);
move2(center1.getX(), center1.getY());
drawstr(number);
move2(center1.getX(), center1.getY());
linestyle("10");
setdash(0.01);
draw2(center2.getX(), center2.getY());
linestyle("1");
drawstr(number);
color(Yellow);
if (track->getNextClosestSidc1Hit()) {
centerNext = track->getNextClosestSidc1Hit()->getCenterInXYZ();
move2(center1.getX(), center1.getY());
draw2(centerNext.getX(), centerNext.getY());
}
if (track->getNextClosestSidc2Hit()) {
centerNext = track->getNextClosestSidc2Hit()->getCenterInXYZ();
move2(center2.getX(), center2.getY());
draw2(centerNext.getX(), centerNext.getY());
}
}
//
// Add number in blue for matching RICH-1 ring center
// and draw line to matching sidc1 and sidc2 hits (if any)
//
if (track->getTrackMask()&CTrack::hasRich1Match && track->getTrackMask()&CTrack::hasVertex) {
evtCoord = CEventCoord(track->getClosestRich1Ring()->getPolarCenter().getTheta(),
track->getClosestRich1Ring()->getPolarCenter().getPhi(),
setup->getZPosition() - track->getVertex()->getZ());
track->getVertex()->transform(center, evtCoord);
color(Blue);
//
// Draw track number
//
move2(center.getX(), center.getY());
drawstr(number);
//
// Draw line to matching SiDC-1 hit
//
if (track->getTrackMask()&CTrack::hasSidc1Match) {
move2(center.getX(), center.getY());
center1 = track->getClosestSidc1Hit()->getCenterInXYZ();
linestyle("10");
setdash(0.01);
draw2(center1.getX(), center1.getY());
linestyle("1");
}
if (track->getTrackMask()&CTrack::hasSidc2Match) {
center2 = track->getClosestSidc2Hit()->getCenterInXYZ();
move2(center.getX(), center.getY());
linestyle("10");
setdash(0.01);
draw2(center2.getX(), center2.getY());
linestyle("1");
}
}
}
}
}
CBoolean CSidc::makeClusters(float minAmp)
{
//
// throw away last event
//
clusters.clearAndDestroy();
//
// anything to do?
//
if (cells.length() == 0) return True;
//
// initialize cluster index in cell list
//
int i;
for (i=0; i<cells.length(); i++)
cells(i)->setCluster(0);
CSidcCluster *newCluster;
for (i=0; i<cells.length(); i++) {
if (cells(i)->getAmp() >= minAmp && cells(i)->getCluster() == 0) {
newCluster = new CSidcCluster();
appendCellToCluster(cells(i), newCluster, minAmp);
getNextNeighbourCells(cells(i), newCluster, minAmp);
clusters.insert(newCluster);
}
}
return True;
}
void CSidc::getNextNeighbourCells(CCell* thisCell, CSidcCluster* cluster, float minAmp)
{
int thisAnode = thisCell->getAnode();
int thisTBin = thisCell->getTimeBin();
//
// get all neighbour cells which do not already belong to a cluster
//
int iTBin;
CCell *cell;
for (iTBin=thisTBin-1; iTBin<=thisTBin+1; iTBin+=2) {
if (iTBin >= 0 && iTBin < setup->getNTimeBins()) {
cell = cells(thisAnode,iTBin);
if (cell != 0)
if (appendCellToCluster(cell, cluster, minAmp))
getNextNeighbourCells(cell, cluster, minAmp);
}
}
int iAnode;
for (iAnode=thisAnode-1; iAnode<=thisAnode+1; iAnode+=2) {
cell = cells(iAnode,thisTBin);
if (cell != 0)
if (appendCellToCluster(cell, cluster, minAmp))
getNextNeighbourCells(cell, cluster, minAmp);
}
}
CBoolean CSidc::appendCellToCluster(CCell* cell, CSidcCluster* cluster, float minAmp)
{
if (cell->getAmp() >= minAmp && cell->getCluster() == 0) {
cell->setCluster(cluster);
//
// update the cluster attributes with this cells data
//
cluster->incNCells();
cluster->addSumAmp(cell->getAmp());
int iAnode = cell->getAnode();
int iTBin = cell->getTimeBin();
if (cell->getAmp() > cluster->getMaxAmp()) {
cluster->setMaxAmp(cell->getAmp());
cluster->setTBinOfMax(iTBin);
cluster->setAnodeOfMax(iAnode);
}
if (iAnode > cluster->getMaxAnode()) cluster->setMaxAnode(iAnode);
if (iAnode < cluster->getMinAnode()) cluster->setMinAnode(iAnode);
if (iTBin > cluster->getMaxTBin() ) cluster->setMaxTBin(iTBin);
if (iTBin < cluster->getMinTBin() ) cluster->setMinTBin(iTBin);
return True;
}
else
return False;
}
void CSidc::printHitSummary(ostream& ost)
{
ost << "Sidc Hit Summary\n";
CString line('=',46);
ost << line << endl;
for (int i=0; i<hits.length(); i++) {
ost << i << "\tr\t" << hits(i)->getCenter().getRadius() << "\tphi\t" << hits(i)->getCenter().getPhi()
<< "\tamp\t" << hits(i)->getAmp() << endl;
}
}
CBoolean CSidc::makeFittedHits()
{
//
// erase old stuff
//
pulses.clearAndDestroy();
hits.clearAndDestroy();
//
// anything to do?
//
if (clusters.length() == 0) return True;
//
// create the fitting tools
//
int i, ipulse, oldLength;
CList<CCoGPulse> tmpPulseList;
CPulseFitter pulseFitter(getSetup()->getSinglePulseSigma(0),
getSetup()->getSinglePulseSigma(1),
getSetup()->getMaxDeviationForSinglePulse());
for(i=0; i<clusters.length(); i++) {
if (clusters(i)->getType() >= 0) { // ignore removed clusters
tmpPulseList.clear();
makePulses(clusters(i), tmpPulseList);
oldLength=tmpPulseList.length();
for (ipulse=0; ipulse<oldLength; ipulse++) {
pulseFitter.fitAndSplit(cells, tmpPulseList, ipulse);
}
for (ipulse=0; ipulse<tmpPulseList.length(); ipulse++) {
tmpPulseList(ipulse)->applyStopPulseCorrection(stopPulseCoGOffset);
pulses.insert(tmpPulseList(ipulse));
}
getHitsFromPulses(clusters(i), tmpPulseList);
}
}
doBallisticDeficitCorrection();
//
// We do not want a clearAndDestroy() to be called
// when tmpPulseList goes out of scope.
//
tmpPulseList.clear();
return True;
}
void CSidc::makePulses(CSidcCluster *cluster, CList<CCoGPulse>& pulsesList )
{
int minAnode = cluster->getMinAnode();
int maxAnode = cluster->getMaxAnode();
int minTime = cluster->getMinTBin();
int maxTime = cluster->getMaxTBin();
float amp, lastAmp;
CBoolean insidePulse;
CCoGPulse *newPulse = 0;
int anode, time;
for (anode=minAnode; anode<=maxAnode; anode++) {
insidePulse = False;
amp = 0;
for (time=minTime; time<=maxTime; time++) {
lastAmp = amp;
if (cells(anode,time))
amp = cells(anode,time)->getAmp();
else
amp = 0;
if (amp >= readoutThreshold && cluster == cells(anode,time)->getCluster()) {
if (!insidePulse) { // new pulse found
if (time == maxTime) continue; // one time bin at the end, skip it
newPulse = new CCoGPulse(anode, time);
insidePulse = True;
}
else {
if (time == maxTime || cells(anode,time+1) == 0) {
newPulse->setTStop(time);
insidePulse = False;
pulsesList.insert(newPulse);
}
}
}
else if (insidePulse) { // amp is below threshold
newPulse->setTStop(time - 1); // last timebin was the end of pulse
insidePulse = False;
pulsesList.insert(newPulse);
}
}
}
}
CBoolean CSidc::makeCoGHits()
{
//
// erase old stuff
//
pulses.clearAndDestroy();
hits.clearAndDestroy();
//
// anything to do?
//
if (clusters.length() == 0) return True;
int i, ipulse;
CList<CCoGPulse> tmpPulseList;
for(i=0; i<clusters.length(); i++) {
if (clusters(i)->getType() >= 0) { // ignore removed clusters
tmpPulseList.clear();
makeSplittedPulses(clusters(i), tmpPulseList);
for (ipulse=0; ipulse<tmpPulseList.length(); ipulse++) {
tmpPulseList(ipulse)->calculateCoG(*this);
tmpPulseList(ipulse)->applyStopPulseCorrection(stopPulseCoGOffset);
}
getHitsFromPulses(clusters(i), tmpPulseList);
for (ipulse=0; ipulse<tmpPulseList.length(); ipulse++) {
pulses.insert(tmpPulseList(ipulse));
}
}
}
doBallisticDeficitCorrection();
//
// We do not want a clearAndDestroy() to be called
// when tmpPulseList goes out of scope.
//
tmpPulseList.clear();
return True;
}
//
// Reconstruct pulses in a cluster
//
// This routine searches for local maxima for a given anode in the cluster
// and fills the pulse list.
// If there are local minima (=> multiple pulses) the timebin of the minimum
// is not in the pulse.
void CSidc::makeSplittedPulses(CSidcCluster *cluster, CList<CCoGPulse>& pulseList)
{
int minAnode = cluster->getMinAnode();
int maxAnode = cluster->getMaxAnode();
int minTime = cluster->getMinTBin();
int maxTime = cluster->getMaxTBin();
int wasMultiPulse; // offset to eliminate the use of local minima
float amp, lastAmp;
CBoolean insidePulse; // just inside a pulse ?
CBoolean localMaximumFound; // already found a local maximum ?
CCoGPulse *newPulse = 0;
CBoolean subtractPedestals = lookupItems.length() >= setup->getNAnodes() ? True : False;
int anode, time;
for (anode=minAnode; anode<=maxAnode; anode++) {
wasMultiPulse = 0;
insidePulse = False;
localMaximumFound = False;
amp = 0;
for (time=minTime; time<=maxTime; time++) {
lastAmp = amp;
if (cells(anode,time - wasMultiPulse) != 0)
amp = cells(anode,time - wasMultiPulse)->getAmp();
else
amp = 0;
if (amp >= readoutThreshold &&
cluster == cells(anode,time - wasMultiPulse)->getCluster()) {
if (!insidePulse) { // new pulse found
if (time == maxTime) continue; // one time bin at the end, skip it
newPulse = new CCoGPulse(anode, time - wasMultiPulse);
if (subtractPedestals)
newPulse->setMaxAmp(amp - lookupItems(anode)->getPedestal());
else
newPulse->setMaxAmp(amp);
insidePulse = True;
}
else {
if (!localMaximumFound && amp >= lastAmp ) { // search for maximum amplitude
if (subtractPedestals)
newPulse->setMaxAmp(amp - lookupItems(anode)->getPedestal());
else
newPulse->setMaxAmp(amp);
newPulse->setMaxTime(time);
if (time < maxTime) continue; // next tbin if not at the end
}
if (amp <= lastAmp && time < maxTime && // there was a local maximum or
cells(anode,time+1-wasMultiPulse) != 0 &&
cells(anode,time+1-wasMultiPulse)->getAmp() < lastAmp) {
localMaximumFound = True; // we are already at the end
continue;
}
if (time == maxTime || cells(anode,time+1-wasMultiPulse) == 0)
newPulse->setTStop(time); // if we are at the end, last tbin is the end
else {
if(localMaximumFound) {
newPulse->setTStop(time - 2); // go back two tbins to get last tbin because
wasMultiPulse = 1; // we are already at the first bin of a new pulse
}
else
continue;
} // wasMultiPulse corrects tstart for next pulse
localMaximumFound = False;
insidePulse = False;
pulseList.insert(newPulse);
}
}
else if (insidePulse) { // amp is below threshold
newPulse->setTStop(time - 1); // last timebin was the end of pulse
localMaximumFound = False;
insidePulse = False;
pulseList.insert(newPulse);
}
}
}
}
//
// Construct hits out of pulses
//
void CSidc::getHitsFromPulses(CSidcCluster *cluster, CList<CCoGPulse>& pulseList)
{
double timeWindowForAdjacentPulses = setup->getTimeWindowForHit();
int maxAnodesInHit = setup->getAnodeWindowForHit();
int minTimeBinsOverThreshold = setup->getMinTimeBinsForHit();
double momAno0, momAno1, momAno2, sumAmp;
double cogTot; // 'width' of hit
double maxAnoAmp; // buffer to store anode with highest amp
int idMaxAnoAmp; // anode with highest amp for 'width' calc.
int firstPulse, lastAnode, nPulses, maxTBins, i, j, jMax;
double sigma2;
CSidcHit *newHit;
CSidcCellCoord cell;
CLabCylinderCoord center;
//
// loop over pulses in cluster
//
for (i=0; i<pulseList.length(); i++) {
if (pulseList(i)->getIsUsedForHit()) continue;
newHit = new CSidcHit();
pulseList(i)->setIsUsedForHit(True);
pulseList(i)->setHit(newHit);
firstPulse = i;
nPulses = 1;
cogTot = pulseList(i)->getCenter() * pulseList(i)->getAmp(); // cog in time is weigthed with amp
sumAmp = pulseList(i)->getSumAmp();
momAno0 = pulseList(i)->getAmp();
momAno1 = pulseList(i)->getAmp() * pulseList(i)->getAnode();
momAno2 = pulseList(i)->getAmp() * pulseList(i)->getAnode() * pulseList(i)->getAnode();
maxAnoAmp = pulseList(i)->getAmp();
idMaxAnoAmp = i;
lastAnode = pulseList(i)->getAnode();
maxTBins = pulseList(i)->getTStop() - pulseList(i)->getTStart();
//
// compare to residuing pulses
//
jMax = i+1;
for (j=i+1; j<pulseList.length(); j++) {
if (!pulseList(j)->getIsUsedForHit() && // not yet used
fabs(pulseList(i)->getCenter()-pulseList(j)->getCenter())<timeWindowForAdjacentPulses &&
pulseList(j)->getAnode() == lastAnode + 1 && // neigbouring anodes
fabs(pulseList(i)->getAnode()-pulseList(j)->getAnode())<maxAnodesInHit) {
jMax = j;
if(pulseList(j)->getTStop()-pulseList(j)->getTStart()>maxTBins)
maxTBins = pulseList(j)->getTStop() - pulseList(j)->getTStart();
pulseList(j)->setIsUsedForHit(True);
pulseList(j)->setHit(newHit);
lastAnode = pulseList(j)->getAnode();
nPulses++;
cogTot += pulseList(j)->getCenter() * pulseList(j)->getAmp();
sumAmp += pulseList(j)->getSumAmp();
momAno0 += pulseList(j)->getAmp(); // cog in anodes
momAno1 += pulseList(j)->getAmp() * pulseList(j)->getAnode();
momAno2 += pulseList(j)->getAmp() * pulseList(j)->getAnode() * pulseList(j)->getAnode();
if (pulseList(j)->getAmp() > maxAnoAmp) { // update position of biggest amp
maxAnoAmp = pulseList(j)->getAmp();
idMaxAnoAmp = j;
}
}
}
if (maxTBins < minTimeBinsOverThreshold) momAno0=0;
if (momAno0 > 0) {
newHit->setCluster(cluster);
newHit->setFirstPulse(firstPulse);
newHit->setNPulses(nPulses);
newHit->setAmp(momAno0);
newHit->setSumAmp(sumAmp);
newHit->setLengthOfLongestPulse(maxTBins);
//
// calculate and store radius and phi
//
cell.setTimeBin(cogTot/momAno0);
cell.setAnode(momAno1/momAno0);
transform(center, cell);
newHit->setCenter(center);
//
// calculate and store sigma of radius and phi
//
cell.setTimeBin(cell.getTimeBin() + pulseList(idMaxAnoAmp)->getRmsTime());
sigma2 = momAno2/momAno0-momAno1/momAno0*momAno1/momAno0;
if (sigma2 >= 0.)
cell.setAnode(cell.getAnode() + sqrt(sigma2));
transform(center, cell);
center.setRadius(fabs(center.getRadius()-newHit->getCenter().getRadius()));
center.setPhi(fabs(center.getPhi()-newHit->getCenter().getPhi()));
newHit->setSigmaCenter(center);
hits.insert(newHit);
}
else {
for (j=i; j<jMax; j++)
if (pulseList(j)->getHit() == newHit) pulseList(j)->setHit(0);
delete newHit;
}
}
}
CBoolean CSidc::doBallisticDeficitCorrection()
{
//
// Transform the amplitude of each hit in the hit list
//
double oldAmp, newAmp, r;
double bd[4];
int i;
for (i=0; i<4; i++)
bd[i] = setup->getBallisticDeficit(i);
for (i=0; i<hits.length(); i++) {
oldAmp = hits(i)->getAmp();
r = hits(i)->getCenter().getRadius();
newAmp = oldAmp * (1+(bd[0]-(bd[1]+bd[2]*r+bd[3]*r*r))/bd[1]);
hits(i)->setAmp(newAmp);
}
return True;
}
void CSidc::makeDigiHits(CMCServer& mcServ)
{
const CDetectorId detId = getDetectorId();
const int nAnodes = setup->getNAnodes();
const int nTimeBins = setup->getNTimeBins();
const double rmax = setup->getMaxRadius();
const double rmin = setup->getMinRadius();
const double tmax = setup->getMaxTimeBin();
const double tmin = setup->getMinTimeBin();
const double minCell = setup->getMinCellAmplitude();
const double phioffset = setup->getPhiOffset();
const CBoolean isflipped = setup->getIsFlipped();
const double x0 = setup->getXOffset();
const double y0 = setup->getYOffset();
const double z0 = setup->getZPosition(); // position of sidc, (0,0,0) = center of sidc1
const double zoffset = setup->getMCZOffset(); // use this to correct for shift in coordinates
const double meanAmp = setup->getMCMeanAmplitude();
const CBoolean useNoise = setup->getMCUseNoise();
const double sigRad = setup->getMCSigmaRadius();
const double sigPhi = setup->getMCSigmaPhi();
int i;
double dr;
CAngle dphi;
const RWTPtrOrderedVector<CMCDigiHit> *digiHit;
double amp;
CLabCylinderCoord polar;
CLabXYZCoord xy;
CSidcCellCoord at;
CRandom rndmGen;
digiHit = mcServ.getMCDigiHits( detId );
cells.resize( (cells.length() + digiHit->length()*8) ); // assume 8 longs per hit
for (i=0; i<digiHit->length(); i++) {
//
// modify GEANTs x/y according to detector parameters
// simple parametrisation used here..
//
xy.setX( (*digiHit)(i)->getGeantHit().getX() );
xy.setY( (*digiHit)(i)->getGeantHit().getY() );
xy.setZ( z0 );
dr = rndmGen.gauss(); // get normal distributed random numbers
dphi = rndmGen.gauss();
::transform(polar, xy); // transform from xy to polar
polar.setRadius( polar.getRadius() + dr *sigRad ); // apply resolution
polar.setPhi ( polar.getPhi () + dphi*sigPhi );
::transform(xy, polar); // transform from polar to xy
amp = (*digiHit)(i)->getGeantHit().getAmp() * 1.e5; // convert from dEdx to FADC counts (so far arbitrary)
(*digiHit)(i)->getDigiHit().setX(xy.getX()); // update info
(*digiHit)(i)->getDigiHit().setY(xy.getY());
(*digiHit)(i)->getDigiHit().setZ(xy.getZ()); // first shot, needs proper function for SiDC-2 later...
(*digiHit)(i)->getDigiHit().setAmp(amp);
} // end loop all digiHits
}
CBoolean CSidc::digitize(CMCServer& mcServ, CMCDigiMode digiMode){
const CDetectorId detId = getDetectorId();
const double diffusion = setup->getMCDiffusionConst();
const double cmtorad = setup->getMCCmToRad();
const double rMax = setup->getMaxRadius();
const double rMin = setup->getMinRadius();
const double GeVtoe = setup->getMCGevToElectrons();
const double eToVolt = setup->getMCElectronsToVolt();
const double noiseForHit = setup->getMCNoise();
const double sigScale = setup->getMCPedestal();
const double sigelec = setup->getMCSigmaElectronic();
const CBoolean useNoise = setup->getMCUseNoise();
const double vdrift = setup->getDriftVelocity();
const int scanThresh = setup->getMCThreshhold(); // in counts
const double sqrt2 = 1.4142135624;
CRandom noisegenerator;
CLabCylinderCoord cylinderHit;
CLabXYZCoord xyHit;
CSidcCellCoord cellHit;
CCell *cell;
double hitRadius, hitTBin, dEdx;
double hitAnode, anodeAmplitude1, anodeAmplitude2, cellAmplitude1, cellAmplitude2;
double cellAmplitudeBefore1, cellAmplitudeBefore2, overlayAmplitude1, overlayAmplitude2;
int anodeNumber, currentAnodeNumber1, currentAnodeNumber2, tBinNumber;
double sigAnode, sigTime;
double anodePedestal1, anodePedestal2, anodeSigma1, anodeSigma2, baseline;
int iAnode, iAnode2, iDigi, iTBin, icrate;
const RWTPtrOrderedVector<CMCDigiHit> *digiHit;
digiHit = mcServ.getMCDigiHits( detId );
RWTValOrderedVector<int> startReadout; // holds the start-
RWTValOrderedVector<int> stopReadout; // and stoppoints for the scanner
RWTValOrderedVector<int> presamples;
if ( digiMode == CMCResetMode ) {
cells.clearAndDestroy();
}
for (iDigi=0; iDigi<digiHit->length(); iDigi++) { //loop over DigiHits
//
// unpack hit information and transform it in useful Coordinates
//
xyHit.setX(((*digiHit)(iDigi)->getDigiHit().getX()));
xyHit.setY(((*digiHit)(iDigi)->getDigiHit().getY()));
xyHit.setZ((*digiHit)(iDigi)->getDigiHit().getZ());
::transform(cylinderHit, xyHit);
transform(cellHit,cylinderHit);
hitRadius = cylinderHit.getRadius(); // in cm
if((hitRadius>rMax) || (hitRadius<rMin)) continue; // skip hits which are not on SIDC
hitAnode = cellHit.getAnode(); // center of hit in Anodes (float)
anodeNumber = int (hitAnode); // number of center anode (int)
hitTBin = cellHit.getTimeBin(); // center of hit in time bins (float)
icrate = anodeNumber/90;
hitTBin = hitTBin - stopPulseCoGOffset[icrate]; // simulate stop pulse correction
tBinNumber = int (hitTBin + 0.5); // number of center time bin (int)
dEdx = (*digiHit)(iDigi)->getGeantHit().getAmp(); // in GeV
dEdx = GeVtoe*dEdx; // in electrons
//
// calculate sigma of electron cloud in anodes and time bins
//
sigAnode = sqrt (2*diffusion * rMax/vdrift * (rMax/hitRadius - 1)) * cmtorad * ToDegree;
sigTime = sqrt (2*diffusion * (rMax - hitRadius)/vdrift) / vdrift;
sigTime = sqrt (sigTime*sigTime + sigelec*sigelec);
//
// distribute the electron cloud over the anodes
// in steps of 2 Anodes at the same time, to simulate the scanner
//
baseline = noisegenerator.gauss() * noiseForHit;
for(iAnode = -2; iAnode <3; iAnode = iAnode + 2){
currentAnodeNumber1 = anodeNumber + iAnode;
if (currentAnodeNumber1 < 0) currentAnodeNumber1 = 360 - currentAnodeNumber1;
if (currentAnodeNumber1 > 359) currentAnodeNumber1 = currentAnodeNumber1 - 360;
if ( (currentAnodeNumber1 % 2) == 0){
currentAnodeNumber2 = currentAnodeNumber1 + 1;
iAnode2 = iAnode + 1;
}
else {
currentAnodeNumber2 = currentAnodeNumber1 - 1;
iAnode2 = iAnode - 1;
}
if (currentAnodeNumber2 < 0) cout << "!!!!!!! attention anode2 out of range :" << currentAnodeNumber2 << endl;
if (currentAnodeNumber2 > 359) cout << "!!!!!!! attention anode2 out of range :" << currentAnodeNumber2 << endl;
//
// calculate the amplitudes for each anode,
// simulate ballistic deficit and transform electrons to mVolt
//
anodeAmplitude1 = calculateAmplitude((double) iAnode, hitAnode, sigAnode, dEdx);
anodeAmplitude1 = anodeAmplitude1 * ballisticDeficitFunction(hitRadius) * eToVolt * 1000.;
anodeAmplitude2 = calculateAmplitude((double) iAnode2, hitAnode, sigAnode, dEdx);
anodeAmplitude2 = anodeAmplitude2 * ballisticDeficitFunction(hitRadius) * eToVolt * 1000.;
//
// check if dead anode
//
CBoolean deadAnode;
int deadFadcChannel, idead;
deadAnode=False;
for (idead = 0; idead<setup->getNumberOfDeadFadcChannel(); idead++){
deadFadcChannel = setup->getDeadFadcChannel(idead);
if(deadFadcChannel == currentAnodeNumber1) anodeAmplitude1 = 0;
if(deadFadcChannel == currentAnodeNumber2) anodeAmplitude2 = 0;
}
if (anodeAmplitude1 == 0 && anodeAmplitude2 == 0) continue;
//
// distribute the electron cloud over time bins and fill into cell array
//
anodePedestal1 = 200.*(lookupItems(currentAnodeNumber1)->getPedestal()-1.) /
(259.-3.*lookupItems(currentAnodeNumber1)->getPedestal()); // transform pedestal to mVolt
anodeSigma1 = 200.*(lookupItems(currentAnodeNumber1)->getSigma()-1.) /
(259.-3.*lookupItems(currentAnodeNumber1)->getSigma()); // transform sigma to mVolt
cellAmplitudeBefore1 = 0;
anodePedestal2 = 200.*(lookupItems(currentAnodeNumber2)->getPedestal()-1.) /
(259.-3.*lookupItems(currentAnodeNumber2)->getPedestal()); // transform pedestal to mVolt
anodeSigma2 = 200.*(lookupItems(currentAnodeNumber2)->getSigma()-1.) /
(259.-3.*lookupItems(currentAnodeNumber2)->getSigma()); // transform sigma to mVolt
cellAmplitudeBefore2 = 0;
for(iTBin = tBinNumber-6; iTBin < tBinNumber+7; iTBin++){
cellAmplitude1 = anodeAmplitude1 * 0.5 * (erf((iTBin + 0.5 - hitTBin)/(sqrt(2)*sigTime)) -
erf((iTBin - 0.5 - hitTBin)/(sqrt(2)*sigTime)));
cellAmplitude2 = anodeAmplitude2 * 0.5 * (erf((iTBin + 0.5 - hitTBin)/(sqrt(2)*sigTime)) -
erf((iTBin - 0.5 - hitTBin)/(sqrt(2)*sigTime)));
if( cells( currentAnodeNumber1, iTBin) !=0 ){ // cell already exists
cellAmplitude1 = cellAmplitude1 + baseline;
overlayAmplitude1 = cells( currentAnodeNumber1, iTBin)->getAmp(); // in mVolt
cellAmplitude1 = cellAmplitude1 + overlayAmplitude1;
if (cellAmplitude1 > 200) cellAmplitude1 = 200; // dynamic region of the Fadc 0 - 200mV
cellAmplitude1 = cellAmplitude1/(0.05+0.00075*cellAmplitude1)*0.064; // transform to channels
cellAmplitudeBefore1 = cellAmplitude1;
cellAmplitude1 = 200.*(cellAmplitude1-1.)/(259.-3.*cellAmplitude1); // retransform to mVolt to simulate binning
cells( currentAnodeNumber1, iTBin)->setAmp( cellAmplitude1);
cells( currentAnodeNumber1, iTBin)->addMCDigiHit((*digiHit)(iDigi));
}
if( cells( currentAnodeNumber2, iTBin) !=0 ){ // cell already exists
cellAmplitude2 = cellAmplitude2 + baseline;
overlayAmplitude2 = cells( currentAnodeNumber2, iTBin)->getAmp(); // in mVolt
cellAmplitude2 = cellAmplitude2 + overlayAmplitude2;
if (cellAmplitude2 > 200) cellAmplitude2 = 200; // dynamic region of the Fadc 0 - 200mV
cellAmplitude2 = cellAmplitude2/(0.05+0.00075*cellAmplitude2)*0.064; // transform to channels
cellAmplitudeBefore2 = cellAmplitude2;
cellAmplitude2 = 200.*(cellAmplitude2-1.)/(259.-3.*cellAmplitude2); // retransform to mVolt to simulate binning
cells( currentAnodeNumber2, iTBin)->setAmp( cellAmplitude2);
cells( currentAnodeNumber2, iTBin)->addMCDigiHit((*digiHit)(iDigi));
}
// for new cells you first have to check if minimum 2 tbins in one of
// the 2 anodes are above the scanner threshold
if (cellAmplitude1 < scanThresh && cellAmplitudeBefore1 < scanThresh &&
cellAmplitude2 < scanThresh && cellAmplitudeBefore2 < scanThresh) continue;
if( cells( currentAnodeNumber1, iTBin) == 0 ){ // new cell
cellAmplitude1 = cellAmplitude1 + anodePedestal1; // add pedestals as stored in SOR
cellAmplitude1 = noisegenerator.gauss(cellAmplitude1,anodeSigma1 * sigScale); // add noise as stored in SOR
cellAmplitude1 = cellAmplitude1 + baseline;
if (cellAmplitude1 > 200) cellAmplitude1 = 200; // dynamic region of the Fadc 0 - 200 mV
cellAmplitude1 = cellAmplitude1/(0.05+0.00075*cellAmplitude1)*0.064; // transform to channels
cellAmplitudeBefore1 = cellAmplitude1;
cellAmplitude1 = 200.*(cellAmplitude1-1.)/(259.-3.*cellAmplitude1); // retransform to mVolt to simulate binning
cell = new CCell(currentAnodeNumber1, iTBin, cellAmplitude1); // creation of the cell,
cell->addMCDigiHit((*digiHit)(iDigi));
if (!cells.insert(cell)) delete cell;
}
if( cells( currentAnodeNumber2, iTBin) == 0 ){ // new cell
cellAmplitude2 = cellAmplitude2 + anodePedestal2; // add pedestals as stored in SOR
cellAmplitude2 = noisegenerator.gauss(cellAmplitude2,anodeSigma2 * sigScale); // add noise as stored in SOR
cellAmplitude2 = cellAmplitude2 + baseline;
if (cellAmplitude2 > 200) cellAmplitude2 = 200; // dynamic region of the Fadc 0 - 200 mV
cellAmplitude2 = cellAmplitude2/(0.05+0.00075*cellAmplitude2)*0.064; // transform to channels
cellAmplitudeBefore2 = cellAmplitude2;
cellAmplitude2 = 200.*(cellAmplitude2-1.)/(259.-3.*cellAmplitude2); // retransform to mVolt to simulate binning
cell = new CCell(currentAnodeNumber2, iTBin, cellAmplitude2); // creation of the cell,
cell->addMCDigiHit((*digiHit)(iDigi));
if (!cells.insert(cell)) delete cell;
}
}
}
} // end loop all DigiHits
return True;
}
double CSidc::ballisticDeficitFunction(double rMax){
const CDetectorId detId = getDetectorId();
double ballistic_deficit_function;
//
// Dabrowski Chip Specification
//
if (detId == SiDC1){
ballistic_deficit_function=0.79098+0.055940*rMax+0.0029307*rMax*rMax;
}
if (detId == SiDC2){
ballistic_deficit_function=0.61446+0.055907*rMax+0.0086416*rMax*rMax;
}
return ballistic_deficit_function;
}
CBoolean CSidc::makeHitsFromMC(CMCServer& mcServ, CMCDigiMode digiMode)
{
const CDetectorId detId = getDetectorId();
const int nAnodes = setup->getNAnodes();
const int nTimeBins = setup->getNTimeBins();
const double rmax = setup->getMaxRadius();
const double rmin = setup->getMinRadius();
const double tmax = setup->getMaxTimeBin();
const double tmin = setup->getMinTimeBin();
const double minCell = setup->getMinCellAmplitude();
const double phioffset = setup->getPhiOffset();
const CBoolean isflipped = setup->getIsFlipped();
const double x0 = setup->getXOffset();
const double y0 = setup->getYOffset();
const double z0 = setup->getZPosition(); // position of sidc, (0,0,0) = center of sidc1
const double zoffset = setup->getMCZOffset(); // use this to correct for shift in coordinates
if ( digiMode == CMCResetMode ) {
hits.clearAndDestroy(); // avoid endless overlays
}
int i;
const RWTPtrOrderedVector<CMCDigiHit> *digiHit;
CSidcHit *hit;
CLabCylinderCoord polar;
CLabXYZCoord xy;
digiHit = mcServ.getMCDigiHits( detId );
for (i=0; i<digiHit->length(); i++) {
xy.setX( (*digiHit)(i)->getDigiHit().getX() ); // for an alternative...
xy.setY( (*digiHit)(i)->getDigiHit().getY() );
xy.setZ( (*digiHit)(i)->getDigiHit().getZ() );
::transform(polar, xy); // transform from xy to polar
hit = new CSidcHit;
hit->setCenterInXYZ(xy);
hit->setAmp( (*digiHit)(i)->getDigiHit().getAmp() );
hits.insert(hit);
} // end loop all digiHits
return True;
}
CBoolean CSidc::readDeltaRadiusFromFile(const char* filename)
{
//
// Open file.
//
Cifstream ifs(filename);
if (!ifs) {
cerr << "CSidc::readDeltaRadiusFromFile():\n";
cerr << "\tERROR\n";
cerr << "\tError opening deltaradius-calibration-file '" << filename << endl;
return False;
}
//
// Fill calibration constants from file.
//
double radius, deltaradius;
CLookupDeltaRadius *newdeltaradius;
deltaRLookupTable.clearAndDestroy();
for (int i=0; i<211; i++) {
ifs >> radius >> deltaradius;
if ((ifs.bad())|| (ifs.eof())) {
cerr << "CSidc::readDeltaRadiusFromFile():\n";
cerr << "\tERROR\n";
cerr << "\tFile too short. No further action." << endl;
return False;
}
newdeltaradius = new CLookupDeltaRadius(radius, deltaradius);
deltaRLookupTable.append(newdeltaradius);
}
ifs.close();
return True;
}
//
// The following are coordinate transformation routines.
//
// In the transformation from CSidcCellCoord into CLabCylinderCoord
// all necessary calibrations corrections are performed.
// The related calibrations constants are taken from the
// setup object constructed according to the referring setup
// file.
//
void CSidc::transform(CSidcCellCoord& sc, const CLabCylinderCoord& spIn) const
{
double maxRadius = setup->getMaxRadius();
double minTimeBin = setup->getMinTimeBin();
double driftVelocity = setup->getDriftVelocity();
double phiOffset = setup->getPhiOffset();
double timeBin, anode;
CLabXYZCoord xyz;
//
// First undo Tilt correction (see below)
//
CAngle phiOfTilt = setup->getTiltCorrection(0) * ToRadian;
double sinOfTiltAngle = setup->getTiltCorrection(1) / maxRadius;
double cosOfTiltAngle = cos(asin(sinOfTiltAngle));
// revert the angle of tilt, that's all to undo the correction ...
sinOfTiltAngle = - sinOfTiltAngle;
CLabCylinderCoord sp = spIn;
sp.setPhi(sp.getPhi() - phiOfTilt);
::transform(xyz, sp);
double xBeforeTilt = xyz.getX();
double zBeforeTilt = xyz.getZ();
xyz.setX(xBeforeTilt*cosOfTiltAngle - zBeforeTilt*sinOfTiltAngle);
xyz.setZ(xBeforeTilt*sinOfTiltAngle + zBeforeTilt*cosOfTiltAngle);
::transform(sp, xyz);
sp.setPhi(sp.getPhi() + phiOfTilt);
//
// Apply (undo) corrections
// 1. Correct for xy offset
// Both corrections are to be applied in xyz coordinates
// so we first have to transform the coordinates, apply
// the corrections and then transform back into cylinder
// coordinates. Since cell coordinates do not have a z
// value the z(xy) corrections can be ignored.
//
CLabCylinderCoord spnew;
::transform(xyz, sp);
xyz.setX(xyz.getX() - setup->getXOffset()); // correct for x offset
xyz.setY(xyz.getY() - setup->getYOffset()); // correct for y offset
::transform(spnew, xyz);
double phi = spnew.getPhi();
if (setup->getIsFlipped())
anode = phiOffset - ToDegree*phi;
else
anode = ToDegree*phi - phiOffset;
if (anode < 1) anode = anode + 360;
anode -= 1; // -1 for [1-360] -> [0-359]
//
// undo anodewise correction for t0 offset
//
int intanode = int(anode);
if (intanode < 0) intanode = 0;
if (intanode > setup->getNAnodes()-1) intanode = setup->getNAnodes()-1;
double r0corr = lookupItems(intanode)->getT0Correction();
double radius = spnew.getRadius() + r0corr;
// correct for local variations, ap-08-02-97
// if at some point the following correction becomes large, we need to implement
// a proper version. for now it's enough to assume deltaR(rSidc) == deltaR(rCorrected).
radius += deltaRcorrection->getDeltaR(radius);
timeBin = minTimeBin+(maxRadius-radius)/driftVelocity;
sc.setAnode(anode);
sc.setTimeBin(timeBin);
}
void CSidc::transform(CLabCylinderCoord& sp, const CSidcCellCoord& sc) const
{
double maxRadius = setup->getMaxRadius();
double minTimeBin = setup->getMinTimeBin();
double driftVelocity = setup->getDriftVelocity();
double phiOffset = setup->getPhiOffset();
double radius, phi;
double timeBin = sc.getTimeBin();
double anode = sc.getAnode();
//
// correct anodewise for t0 offset
//
int intanode = int(anode);
if (intanode < 0) intanode = 0;
if (intanode > setup->getNAnodes()-1) intanode = setup->getNAnodes()-1;
double r0corr = lookupItems(intanode)->getT0Correction();
radius = maxRadius-(timeBin-minTimeBin)*driftVelocity - r0corr;
radius -= deltaRcorrection->getDeltaR(radius); // correct for local variations, ap-08-02-97
if (setup->getIsFlipped())
phi = phiOffset - (anode+1); // +1 for [0-359] -> [1-360]
else
phi = phiOffset + (anode+1);
if (phi > 360) phi = phi-360;
if (phi < 0) phi = phi+360;
phi *= ToRadian;
sp.setPhi(phi);
sp.setRadius(radius);
sp.setZ(setup->getZPosition());
//
// Apply geometry corrections
// 1. Correct z because of waver deformations
// 2. Correct for xy offset
// 3. Correct for tilt
// Both corrections are to be applied in xyz coordinates
// so we first have to transform the coordinates, apply
// the corrections and then transform back into cylinder
// coordinates. The order of the corrections is important.
// to 1.) the obtained z value gives the absolute
// distance from SiDC-1. Since all corrections constants
// for SiDC-1 are 0 the z comes out 0 which is what we
// want since the first SiDC defines the origin (0,0,0).
// The value obtained from setup->getZPosition() is only
// an average which should be used only as a fallback
// solution. The correction function was calculated by
// Peter Glaessel in Jan/Feb 1996.
//
CLabXYZCoord xyz;
::transform(xyz, sp);
double z = (setup->getDeformationCorrection(0) +
setup->getDeformationCorrection(1)*xyz.getX()*xyz.getX() +
setup->getDeformationCorrection(2)*xyz.getY()*xyz.getY() +
setup->getDeformationCorrection(3)*xyz.getY()*xyz.getX());
xyz.setZ(z);
xyz.setX(xyz.getX() + setup->getXOffset()); // correct for x offset
xyz.setY(xyz.getY() + setup->getYOffset()); // correct for y offset
::transform(sp, xyz); // back to cylinder coordinates
//
// Tilt correction:
// In the e-pair-run 1996 the aluminium ring on which both sidcs are
// mounted was inserted into the carbon-(target)tube with a tilt.
// Thus we have to rotate around the origin to correct this.
// To make this an easy two-dim rotation we first transform into
// the tilt-plane (phi->phi+phiOfTilt) and undo this after the correction.
//
CAngle phiOfTilt = setup->getTiltCorrection(0) * ToRadian;
double sinOfTiltAngle = setup->getTiltCorrection(1) / maxRadius;
double cosOfTiltAngle = cos(asin(sinOfTiltAngle));
sp.setPhi(sp.getPhi() - phiOfTilt);
::transform(xyz, sp);
double xBeforeTilt = xyz.getX();
double zBeforeTilt = xyz.getZ();
xyz.setX(xBeforeTilt*cosOfTiltAngle - zBeforeTilt*sinOfTiltAngle);
xyz.setZ(xBeforeTilt*sinOfTiltAngle + zBeforeTilt*cosOfTiltAngle);
::transform(sp, xyz);
sp.setPhi(sp.getPhi() + phiOfTilt);
}