//*-- Authors: Leszek Kidon & Jacek Otwinowski
//*-- Last Modified: 08/07/2001 (Marcin Jaskula)
//*-- Last Modified: 17/11/2005 (Jacek Otwinowski)
//*-- Last Modified: 05/12/2006 (Jacek Otwinowski)
//_HADES_CLASS_DESCRIPTION
/////////////////////////////////////////////////////////////////////////////
//
// HShowerPadDigitizer
//
// This class digitizes the shower pads. For each fired wire it calculates all
// the pads which the charge couples to. For each pad the track numbers of the
// particles that fired the pad are stored in a linear category (HShowerTrack).
// This category is sortable, in particular the track numbers are sorted by
// the respective pad address.
// All the fired pads are stored in a matrix category (catShowerRawMatr,
// this matrix category is used only for the simulated data).
//
// The Shower digitization is split into several tasks as shown
// in the flow diagram below.
//
// ---------------------- //
// | HShowerUnpacker | //
// | (embedding mode) | //
// | | ------------------ //
// ---------------------- | | HGeantShower | //
// | ------------------ //
// | //
// | ------------------ -------------> ---------------------- //
// | | HGeantWire | <------------ | HShowerHitDigitizer | //
// | ------------------ ---------------------- //
// | //
// ------------- ------------------ -------------> ----------------------- //
// -- | HShowerRaw | | HShowerRawMatr | <------------ | HShowerPadDigitizer | //
// | ------------- ------------------ |( creates track objects| //
// | | for real tracks in | //
// ---------------------- ------------------ | embedding mode too) | //
// | HShowerCalibrater | | HShowerTrack | <------------ ----------------------- //
// | (embedding mode) | ------------------ //
// ---------------------- ----------------------- //
// | ------------------ ----------> | HShowerCopy | //
// -------------------> | HShowerCal | <------------ |(add charge of real hit| //
// ------------------ | in embedding too ) | //
// ----------------------- //
// ------------------ -------------> ----------------------- //
// | HShowerHitHdr | <----------- | HShowerHitFinder | //
// ------------------ ----------------------- //
// ------------------ | //
// | HShowerPID | <------------| //
// ------------------ | //
// ------------------ | //
// | HShowerHit | <------------| //
// ------------------ //
// //
// ------------------ -------------> ------------------------ //
// | HShowerHitTrack | <------------ | HShowerHitTrackMatcher | //
// ------------------ ------------------------ //
//
//
//
//
// In the case of TRACK EMBEDDING of simulated tracks into
// experimental data the real data are written by the HShowerUnpacker into
// HShowerRaw category. The real hits are taken into
// account by the digitizer (adding of charges). The embedding mode is recognized
// automatically by analyzing the
// gHades->getEmbeddingMode() flag.
// Mode ==0 means no embedding
// ==1 realistic embedding (first real or sim hit makes the game)
// ==2 keep GEANT tracks (as 1, but GEANT track numbers will always
// win against real data. besides the tracknumber the output will
// be the same as in 1)
//
//////////////////////////////////////////////////////////////////////////////
#include "hshowerpaddigitizer.h"
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "hruntimedb.h"
#include "hevent.h"
#include "hspectrometer.h"
#include "hdetector.h"
#include "hshowerdetector.h"
#include "hcategory.h"
#include "hmatrixcategory.h"
#include "hlinearcategory.h"
#include "hlocation.h"
#include "hshowergeantwire.h"
#include "hshowerwire.h"
#include "hshowercal.h"
#include "hshowertrack.h"
#include "hshowerpad.h"
#include "hshowerframe.h"
#include "hshowerdigidetpar.h"
#include "hshowergeometry.h"
#include "hdebug.h"
#include "hades.h"
#include "showerdef.h"
#define MIN(A,B) (((A) <= (B)) ? (A) : (B))
#define MAX(A,B) (((A) >= (B)) ? (A) : (B))
#define TABLE_SIZE 10
ClassImp(HShowerPadDigitizer)
HShowerPadDigitizer::HShowerPadDigitizer(Text_t *name,Text_t *title) :
HShowerDigitizer(name,title)
{
fTrackIter = NULL;
fChannelCoeff = 256.0 / 60.0; // QDC calibration factor
m_pSimulPar= NULL;
fShowerRawMatrIter = NULL;
}
HShowerPadDigitizer::HShowerPadDigitizer()
{
fTrackIter = NULL;
m_pSimulPar= NULL;
fShowerRawMatrIter = NULL;
}
HShowerPadDigitizer::~HShowerPadDigitizer(void) {
if (fTrackIter) delete fTrackIter;
if (fShowerRawMatrIter) delete fShowerRawMatrIter;
}
Bool_t HShowerPadDigitizer::init()
{
// creates ShowerGeantWire(input), ShowerRawMatr(output) and ShowerTrack
// categories and adds them to current event
// creates an iterator for ShowerTrack category
printf("initialization of shower pad digitizer n");
HCategory *pCat;
HShowerDetector *pShowerDet = (HShowerDetector*)gHades->getSetup()->getDetector("Shower");
pCat=gHades->getCurrentEvent()->getCategory(catShowerGeantWire);
if (!pCat) {
pCat=pShowerDet->buildCategory(catShowerGeantWire);
if (!pCat) return kFALSE;
else gHades->getCurrentEvent()->addCategory(catShowerGeantWire, pCat, "Shower");
}
setInCat(pCat);
pCat=gHades->getCurrentEvent()->getCategory(catShowerRawMatr);
if (!pCat) {
pCat=pShowerDet->buildCategory(catShowerRawMatr);
if (!pCat) return kFALSE;
else {
gHades->getCurrentEvent()->addCategory(catShowerRawMatr, pCat, "Shower");
}
}
setOutCat(pCat);
if(pCat) fShowerRawMatrIter=(HIterator*)pCat->MakeIterator("native");
if(gHades->getEmbeddingMode()>0)
{ // we need to get the HShowerCal catgeory
// in addition
pCat=gHades->getCurrentEvent()->getCategory(catShowerCal);
if (!pCat) {
Error("init()","No catShowerCal in input ... needed for embedding!");
return kFALSE;
}
if(pCat) fShowerCalIter= (HIterator*)pCat->MakeIterator("native");
}
pCat=gHades->getCurrentEvent()->getCategory(catShowerTrack);
if (!pCat) {
pCat=pShowerDet->buildCategory(catShowerTrack);
if (!pCat) return kFALSE;
else gHades->getCurrentEvent()->addCategory(catShowerTrack, pCat, "Shower");
}
setTrackCat(pCat);
fTrackIter = (HIterator*)getTrackCat()->MakeIterator("native");
int res = HShowerDigitizer::init();
//init ShowerPar container
m_pSimulPar = (HShowerSimulPar*)gHades->getRuntimeDb()->getContainer("ShowerSimulPar");
if (!m_pSimulPar)
{
::Error("HShowerPadDigitizer::init()","ShowerSimulPar not initialised");
return kFALSE;
}
return res;
}
HShowerPadDigitizer& HShowerPadDigitizer::operator=(HShowerPadDigitizer &c) {
// It should have been done
return c;
}
Int_t HShowerPadDigitizer::execute()
{
if(gHades->getEmbeddingMode()>0)
{
//-------------------------------------------------
// loop over ShowerCal category and create the
// corresponding HShowerTrack for the real data objects.
HLocation loc;
loc.set(4,0,0,0);
HShowerCal *pCal;
fShowerCalIter->Reset();
HShowerRaw pRaw;
while ((pCal=(HShowerCal *)fShowerCalIter->Next())!=0)
{
loc[0]=pCal->getSector();
loc[1]=pCal->getModule();
loc[2]=pCal->getRow();
loc[3]=pCal->getCol();
pRaw.setSector(loc[0]);
pRaw.setModule(loc[1]);
pRaw.setRow(loc[2]);
pRaw.setCol(loc[3]);
updateTrack(&pRaw,loc,gHades->getEmbeddingRealTrackId());
}
}
//-------------------------------------------------
TObject *pHit;
HShowerRawMatr *pRawMatr;
lNrEvent++;
fIter->Reset();
Int_t n = 0;
Double_t fQ;
while((pHit = fIter->Next()))
{
digitize(pHit);
n++;
}
// gain and efficiency
fShowerRawMatrIter->Reset();
while((pRawMatr = (HShowerRawMatr*)fShowerRawMatrIter->Next()))
{
// calculate gain
fQ = gainCharge(pRawMatr);
pRawMatr->setCharge(fQ);
// effciency model good for post1 and post2
if(!checkEfficiency(pRawMatr)) {
pRawMatr->setCharge(0.0);
}
}
//printf("chamber %d, gain %4.2f n",nChamber,m_pSimulPar->getGain(nChamber));
sort(); // this is used to call the track sort() in pad digitizer
return 0;
}
Bool_t HShowerPadDigitizer::digitize(TObject *pHit) {
HLocation fLoc;
#if DEBUG_LEVEL>2
gDebuger->enterFunc("HShowerPadDigitizer::execute");
gDebuger->message("Hit cat points to %p",pHit);
#endif
HShowerGeantWire *pGeantWire = (HShowerGeantWire*) pHit;
HShowerDigiDetPar *pDetPar = (HShowerDigiDetPar*) getDigiParSet();
if (!pDetPar)
{
Error("HShowerPadDigitizer::digitize",
"Digitization parameters not initialised!");
return 0;
}
if (pGeantWire) {
//////////////////////////////////////////////////////////
digiPads(pGeantWire);
/////////////////////////////////////////////////////////
}
#if DEBUG_LEVEL>2
gDebuger->leaveFunc("HShowerPadDigitizer::execute");
#endif
return (pGeantWire) ? kTRUE : kFALSE;
}
Int_t HShowerPadDigitizer::sort() {
// if (((HLinearCategory*)getTrackCat())->IsSortable())
((HLinearCategory*)getTrackCat())->sort();
return 0;
}
void HShowerPadDigitizer::updateTrack(HShowerRaw *pShowerRaw, HLocation& loc, Int_t nTrack) {
// this function eliminates double track for the same address
// it means that with one hit should be connected no more than one pad address and one track
HShowerTrack *pShowerTrack=NULL;
Int_t nAddress = pShowerRaw->getAddress();
fTrackIter->Reset();
while((pShowerTrack = (HShowerTrack*)fTrackIter->Next()))
if ((pShowerTrack->getAddress()==nAddress) &&
(pShowerTrack->getTrack()==nTrack))
return;
pShowerTrack = (HShowerTrack*)((HLinearCategory*)getTrackCat())
->getNewSlot(loc);
if (pShowerTrack) {
pShowerTrack = new(pShowerTrack) HShowerTrack;
pShowerTrack->setAddress(pShowerRaw->getAddress());
pShowerTrack->setTrack(nTrack);
}
}
void HShowerPadDigitizer::updatePad(HShowerPad *pPad, Float_t fIndQ, Int_t nSect, Int_t nMod, Int_t nTrack)
{
Int_t nRow, nCol;
HShowerRaw *pShowerRaw=NULL;
HLocation loc;
Float_t fUpdateThreshold = ((HShowerDigiDetPar*)getDigiParSet())
->getUpdatePadThreshold();
pPad->getPadPos(&nRow, &nCol);
loc.set(4, nSect, nMod, nRow, nCol);
pShowerRaw = (HShowerRaw*)((HMatrixCategory*)getOutCat())->getObject(loc);
if (!pShowerRaw)
{
pShowerRaw = (HShowerRaw*)((HMatrixCategory*)getOutCat())->getSlot(loc);
if (pShowerRaw != NULL)
{
pShowerRaw = new(pShowerRaw) HShowerRaw;
pShowerRaw->setSector(loc[0]);
pShowerRaw->setModule(loc[1]);
pShowerRaw->setRow(loc[2]);
pShowerRaw->setCol(loc[3]);
}
}
if(pShowerRaw)
{
pShowerRaw->addCharge(fIndQ);
if(pShowerRaw->getCharge() >= fUpdateThreshold)
updateTrack(pShowerRaw, loc, nTrack);
}
}
Float_t HShowerPadDigitizer::calcCharge(Float_t charge, Float_t dist,
Float_t Xd, Float_t Yd, Float_t Xu, Float_t Yu) {
//charge density matrix is calculated
const float twoPI = 6.28318530718;
return ((charge/twoPI)*(
atan(Xd*Yd/(dist*sqrt(dist*dist+Xd*Xd+Yd*Yd))) -
atan(Xd*Yu/(dist*sqrt(dist*dist+Xd*Xd+Yu*Yu))) +
atan(Xu*Yu/(dist*sqrt(dist*dist+Xu*Xu+Yu*Yu))) -
atan(Xu*Yd/(dist*sqrt(dist*dist+Xu*Xu+Yd*Yd)))));
}
/* ************************************************************************** */
void HShowerPadDigitizer::digiPads(HShowerGeantWire* pWireHit)
{
Float_t fX, fY, fQ;
Float_t factor;
Float_t sum;
Int_t nPadX, nPadY, nPadRange;
Float_t pChargeTable[MAX_PADS_DIST];
int iMaxInArray, m;
HShowerPad *pTmpPad;
HShowerPad *pCPad;
Float_t fPadThreshold = ((HShowerDigiDetPar *)getDigiParSet())
->getPadThreshold();
Float_t fUpdateThreshold = ((HShowerDigiDetPar *)getDigiParSet())
->getUpdatePadThreshold();
HShowerPadTab *pPadParam = ((HShowerGeometry *)getGeometry())
->getPadParam(pWireHit->getModule());
// don't calculate if charge belowe threshold
if((fQ = pWireHit->getCharge()) < fUpdateThreshold)
return;
pWireHit->getXY(&fX, &fY);
//row = Y,col = X
if((pCPad = pPadParam->getPad(fX, fY)) == NULL)
return;
pCPad->getPadPos(&nPadY, &nPadX); //row = Y; col = X
iMaxInArray = MAX_PADS_DIST;
nPadRange = (MAX_PADS_DIST - 1) / 2;
if(numericalCalc(pCPad, pWireHit, pChargeTable) != 0)
return;
// calculate sum for the first renormalization
sum = 0.0f;
for(m = 0; m < iMaxInArray; m++)
sum += pChargeTable[m];
if(sum <= 0.0f)
return;
factor = fQ / sum;
for(m = 0; m < iMaxInArray; m++)
{
pChargeTable[m] *= factor;
}
// calculate sum for a renormalization
sum = 0.0f;
for(m = 0; m < iMaxInArray; m++)
{
if(pChargeTable[m] >= fPadThreshold)
sum += pChargeTable[m];
else
pChargeTable[m] = 0.0f;
}
// if sum is below threshold put all charge to the central pad
if(sum < fPadThreshold)
{
for(m = 0; m < iMaxInArray; m++)
pChargeTable[m] = 0.0f;
pChargeTable[nPadRange] = fQ;
}
else
{
// renormalization: sum of the charge on the pads
// must be equal to the charge on the wire
factor = fQ / sum;
for(m = 0; m < iMaxInArray; m++)
{
if(pChargeTable[m] > 0.0f)
pChargeTable[m] *= factor;
}
}
for(m = 0; m < iMaxInArray; m++)
{
pTmpPad = pPadParam->getPad(nPadY, nPadX + m - nPadRange);
if(((m == nPadRange) || (pChargeTable[m] >= fPadThreshold))
&& (pTmpPad != NULL) && (pTmpPad->getPadFlag()))
{
updatePad(pTmpPad, pChargeTable[m], pWireHit->getSector(),
pWireHit->getModule(), pWireHit->getTrack());
}
}
}
//---------------------------------------------------------
Double_t HShowerPadDigitizer::gainCharge(HShowerRawMatr *pRawMatr)
{
// gain charge in simulation
Int_t chamber = pRawMatr->getSector()+pRawMatr->getModule()*6;
Double_t charge = pRawMatr->getCharge()*m_pSimulPar->getGain(chamber);
return charge;
}
//---------------------------------------------------------
Bool_t HShowerPadDigitizer::checkEfficiency(HShowerRawMatr *pRawMatr)
{
// check efficiency
// model good for post converters
// if pre-converter
if(pRawMatr->getModule() == 0) return kTRUE;
Int_t chamber = pRawMatr->getSector()+pRawMatr->getModule()*6;
Double_t threshold = m_pSimulPar->getThreshold(chamber);
//printf("chamber %d, threshold %4.2f, charge %4.2f n",chamber,threshold,charge);
if(pRawMatr->getCharge() > threshold) return kTRUE;
else return kFALSE;
}
/* ************************************************************************** */
Int_t HShowerPadDigitizer::calcLimit(Float_t fCor, Int_t nMatrixRange,
Float_t fBoxSize)
{
//charge density matrix is delimited to the sensitive area which covers fired pads
Int_t k, nLimit;
k = (Int_t)((fabs(fCor) + 0.5*fBoxSize)/fBoxSize);
if (fCor < 0 )
nLimit = (k>nMatrixRange) ? 0 : nMatrixRange - k;
else
nLimit = (k>nMatrixRange) ? 2*nMatrixRange : nMatrixRange + k;
return nLimit;
}
void HShowerPadDigitizer::moveCoord(HShowerPad *pPad,Float_t distWire, Float_t fDx, Float_t fDy,
Float_t *corXld, Float_t *corYld, Float_t *corXlu, Float_t *corYlu,
Float_t *corXrd, Float_t *corYrd, Float_t *corXru, Float_t *corYru) {
// moves charge density matrix to the place where is placed induced charge on the sense wire
*corXld = pPad->fXld - fDx;
// *corYld = pPad->fYld - fDy;
*corYld = fDy - fDy - distWire;
*corXlu = pPad->fXlu - fDx;
// *corYlu = pPad->fYlu - fDy;
*corYlu = fDy - fDy + distWire;
*corXrd = pPad->fXrd - fDx;
// *corYrd = pPad->fYrd - fDy;
*corYrd = fDy - fDy - distWire;
*corXru = pPad->fXru - fDx;
// *corYru = pPad->fYru - fDy;
*corYru = fDy - fDy + distWire;
}
void HShowerPadDigitizer::analyticCalc(HShowerPad *pPad, HShowerGeantWire* pWireHit)
{
// this function calculates charge induced on pad analytically
// it is useful only for square or rectangular pads
Float_t fX, fY;
Float_t fIndQ;
Float_t corXld, corYld, corXlu, corYlu;
Float_t corXrd, corYrd, corXru, corYru;
pWireHit->getXY(&fX, &fY);
Int_t nModule = pWireHit->getModule();
HShowerWireTab *pfWire = ((HShowerGeometry*)getGeometry())->
getWireTab(nModule);
Float_t distWire = pfWire->getDistWire();
// cout << " distWire " << distWire << endl;
moveCoord(pPad, distWire, fX, fY, &corXld, &corYld, &corXlu, &corYlu,
&corXrd, &corYrd, &corXru, &corYru);
fIndQ = calcCharge(pWireHit->getCharge(),
((HShowerDigiDetPar*)getDigiParSet())->getPlaneDist(),
corXld, corYld, corXru, corYru);
updatePad(pPad, fIndQ, pWireHit->getSector(),
pWireHit->getModule(), pWireHit->getTrack());
}
//------------------------------------------------------------------------------
int HShowerPadDigitizer::numericalCalc(HShowerPad *pPad,
HShowerGeantWire* pWireHit, Float_t pfValues[MAX_PADS_DIST])
{
Int_t iModule;
HShowerFrame *pFrame;
Int_t iMatrixRange;
Float_t fBoxSize;
Float_t fDistWire;
Float_t fAl, fAr, fBl, fBr;
Float_t fXd, fXu, fYd, fYu;
Float_t fXmin, fXmax;
Float_t fYmin, fYmax;
Float_t fWireX, fWireY;
Float_t fY, fX, fQ;
Int_t i;
HShowerFrameCorner *pCorner;
if((pPad == NULL) || (pPad->getPadFlag() == 0))
return -1;
iModule = pWireHit->getModule();
pFrame = ((HShowerGeometry *)getGeometry())->getFrame(iModule);
fDistWire = (((HShowerGeometry *)getGeometry())->getWireTab(iModule))
->getDistWire();
iMatrixRange = ((HShowerDigiDetPar *)getDigiParSet())->getMatrixRange();
fBoxSize = ((HShowerDigiDetPar *)getDigiParSet())->getBoxSize();
const Float_t* pfChargeMatrix = ((HShowerDigiDetPar *)getDigiParSet())->getChargeMatrix();
// left border
if((pPad->getPadFlag() == 2) && (pPad->fXld < 0.0))
{
pCorner = pFrame->getCorner(0);
fXd = pCorner->getX();
fYd = pCorner->getY();
pCorner = pFrame->getCorner(1);
fXu = pCorner->getX();
fYu = pCorner->getY();
}
else
{
fXd = pPad->fXld;
fYd = pPad->fYld;
fXu = pPad->fXlu;
fYu = pPad->fYlu;
}
fAl = (fXu - fXd) / (fYu - fYd);
fBl = fXu - fAl * fYu;
// right border
if((pPad->getPadFlag() == 2) && (pPad->fXld >= 0.0))
{
pCorner = pFrame->getCorner(3);
fXd = pCorner->getX();
fYd = pCorner->getY();
pCorner = pFrame->getCorner(2);
fXu = pCorner->getX();
fYu = pCorner->getY();
}
else
{
fXd = pPad->fXrd;
fYd = pPad->fYrd;
fXu = pPad->fXru;
fYu = pPad->fYru;
}
fAr = (fXu - fXd) / (fYu - fYd);
fBr = fXu - fAr * fYu;
// range of the charge matrix
pWireHit->getXY(&fWireX, &fWireY);
fXmin = fWireX - iMatrixRange * fBoxSize;
fXmax = fWireX + iMatrixRange * fBoxSize;
fYmin = fWireY - fDistWire + 0.5 * fBoxSize;
fYmax = fWireY + fDistWire - 0.5 * fBoxSize;
memset(pfValues, 0, MAX_PADS_DIST * sizeof(Float_t));
for(fY = fYmin; fY <= fYmax; fY += fBoxSize)
{
// left border
if((fX = fAl * fY + fBl) >= fXmin && (fWireX-fX) >= 0.0)
{
i = (int)((fX - fXmin) / fBoxSize);
if(i>2*iMatrixRange||i<0)
{
cout << " left border fWireX-fX " << fWireX-fX << endl;
cout << " fXd,fYd,fXu,fYu " << pPad->fXld <<","<< pPad->fYld << ","<< pPad->fXlu << ","<< pPad->fYlu << endl;
Error("numericalCalc()","left border matrix : max %i, ind %i, fX %7.4f, fXmin %7.4f, fWireX %7.4f,fBoxSize %5.3f, PadNr %i",2*iMatrixRange,i,fX,fXmin,fWireX,fBoxSize,pPad->getPadNr());
continue;
}
pfValues[0] += pfChargeMatrix[i];
fQ = pfChargeMatrix[i];
}
else
fQ = 0.0f;
// right border
if((fX = fAr * fY + fBr) <= fXmax && (fWireX-fX) <= 0.0)
{
i = (int)((fX - fXmin) / fBoxSize);
if(i>2*iMatrixRange||i<0)
{
cout << " right border fWireX-fX " << fWireX-fX << endl;
cout << " fXd,fYd,fXu,fYu " << pPad->fXrd <<","<< pPad->fYrd << ","<< pPad->fXru << ","<< pPad->fYru << endl;
Error("numericalCalc()","right border matrix : max %i, ind %i, fX %7.4f, fXmin %7.4f, fWireX %7.4f, fBoxSize %5.3f, PadNr %i",2*iMatrixRange,i,fX,fXmin,fWireX,fBoxSize,pPad->getPadNr());
continue;
}
pfValues[1] += pfChargeMatrix[i] - fQ;
pfValues[2] += pfChargeMatrix[2 * iMatrixRange] - pfChargeMatrix[i];
}
else
pfValues[1] += pfChargeMatrix[2 * iMatrixRange] - fQ;
}
return 0;
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.