#include "hades.h"
#include "hcategory.h"
#include "hdebug.h"
#include "hevent.h"
#include "heventheader.h"
#include "hiterator.h"
#include "hlinearcategory.h"
#include "hparset.h"
#include "hrichanalysis.h"
#include "hrichanalysispar.h"
#include "hrichcal.h"
#include "hrichdetector.h"
#include "hrichgeometrypar.h"
#include "hrichhitheader.h"
#include "hrichpadclean.h"
#include "hrichpadlabel.h"
#include "hrichringfind.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "richdef.h"
#include <iomanip>
#include <iostream>
using namespace std;
ClassImp(HRichAnalysis)
HRichAnalysis::HRichAnalysis() : HReconstructor()
{
clear();
pPadClean = new HRichPadClean;
pPadLabel = new HRichPadLabel;
pRingFind = new HRichRingFind;
}
HRichAnalysis::HRichAnalysis(const Text_t* name, const Text_t* title, Bool_t kSkip)
: HReconstructor(name, title)
{
clear();
pPadClean = new HRichPadClean;
pPadLabel = new HRichPadLabel;
pRingFind = new HRichRingFind;
kSkipEvtIfNoRing = kSkip;
}
HRichAnalysis::~HRichAnalysis()
{
if (NULL != pPadClean) {
delete pPadClean;
pPadClean = NULL;
}
if (NULL != pPadLabel) {
delete pPadLabel;
pPadLabel = NULL;
}
if (NULL != pRingFind) {
delete pRingFind;
pRingFind = NULL;
}
if (NULL != pLabelArea) {
delete [] pLabelArea;
pLabelArea = NULL;
}
if (NULL != pPads) {
for (Int_t i = 0; i < 6; i++)
if (pPads[i]) delete [] pPads[i];
delete [] pPads;
pPads = NULL;
}
if (NULL != pLeftBorder) {
delete [] pLeftBorder;
pLeftBorder = NULL;
}
if (NULL != pRightBorder) {
delete [] pRightBorder;
pRightBorder = NULL;
}
if (NULL != fIter) {
delete fIter;
fIter = NULL;
}
if (NULL != fIterHitHeader) {
delete fIterHitHeader;
fIterHitHeader = NULL;
}
}
void HRichAnalysis::clear()
{
pLabelArea = NULL;
pRings = NULL;
pPads = NULL;
pSectorPads = NULL;
fIter = NULL;
fIterHitHeader = NULL;
pLeftBorder = NULL;
pRightBorder = NULL;
m_pCalCat = NULL;
m_pHitCat = NULL;
m_pHitHdrCat = NULL;
fpAnalysisPar = NULL;
fpGeomPar = NULL;
hithdrLoop = NULL;
for (Int_t i = 0; i < 6; ++i)
fPadFired[i] = 0;
iClustersCleaned.Reset();
allPairNrTot = 0;
iActiveSector = 0;
iClusterCleanedNr = 0;
iFakeLocalMax4 = 0;
iFakeLocalMax8 = 0;
iFakePad = 0;
iLabelNr = 0;
iPadCleanedNr = 0;
iPadFiredNr = 0;
iRingNr = 0;
iRingNrTot = 0;
sectorPairNrTot = 0;
}
void
HRichAnalysis::Reset()
{
iClustersCleaned.Reset();
iClusterCleanedNr = 0;
iFakeLocalMax4 = 0;
iFakeLocalMax8 = 0;
iFakePad = 0;
iLabelNr = 0;
iPadCleanedNr = 0;
iPadFiredNr = 0;
iRingNr = 0;
if (pLabelArea) {
delete [] pLabelArea;
pLabelArea = NULL;
}
}
Bool_t
HRichAnalysis::init()
{
HRichDetector* pRichDet = NULL;
HRuntimeDb* rtdb = NULL;
HRichAnalysisPar* pAnalysisPar = NULL;
HRichGeometryPar* pGeomPar = NULL;
pRichDet = static_cast<HRichDetector*>(gHades->getSetup()
->getDetector("Rich"));
if (NULL == pRichDet) {
Error("init", "Pointer to HRichDetector is NULL");
return kFALSE;
}
m_pCalCat = gHades->getCurrentEvent()->getCategory(catRichCal);
if (NULL == m_pCalCat) {
m_pCalCat = pRichDet->buildCategory(catRichCal);
if (NULL == m_pCalCat) {
Error("init", "Pointer to HRichCal category is NULL");
return kFALSE;
} else
gHades->getCurrentEvent()
->addCategory(catRichCal, m_pCalCat, "Rich");
}
fIter = static_cast<HIterator*>(m_pCalCat->MakeIterator());
if (NULL == fIter) {
Error("init", "Pointer to HRichCalSim iterator is NULL");
return kFALSE;
}
m_pHitCat = gHades->getCurrentEvent()->getCategory(catRichHit);
if (NULL == m_pHitCat) {
m_pHitCat = pRichDet->buildCategory(catRichHit);
if (NULL == m_pHitCat) {
Error("init", "Pointer to HRichHit category is NULL");
return kFALSE;
} else gHades->getCurrentEvent()
->addCategory(catRichHit, m_pHitCat, "Rich");
}
m_pHitHdrCat = gHades->getCurrentEvent()->getCategory(catRichHitHdr);
if (NULL == m_pHitHdrCat) {
m_pHitHdrCat = pRichDet->buildCategory(catRichHitHdr);
if (NULL == m_pHitHdrCat) {
Error("init", "Pointer to HRichHitHdr category is NULL");
return kFALSE;
} else gHades->getCurrentEvent()
->addCategory(catRichHitHdr, m_pHitHdrCat, "Rich");
}
fIterHitHeader = static_cast<HIterator*>(m_pHitHdrCat->MakeIterator("native"));
if (NULL == fIterHitHeader) {
Error("init", "Pointer to HRichHitHdr iterator is NULL");
return kFALSE;
}
rtdb = gHades->getRuntimeDb();
pAnalysisPar = static_cast<HRichAnalysisPar*>(rtdb->
getContainer("RichAnalysisParameters"));
if (NULL == pAnalysisPar) {
Error("init", "Pointer to HRichAnalysisPar parameters is NULL");
return kFALSE;
}
setAnalysisPar(pAnalysisPar);
pGeomPar = static_cast<HRichGeometryPar*>(rtdb->
getContainer("RichGeometryParameters"));
if (NULL == pGeomPar) {
Error("init", "Pointer to RichGeometryParameters parameters is NULL");
return kFALSE;
}
setGeometryPar(pGeomPar);
if (kFALSE == initParameters()) {
Error("init", "Parameters not initialized");
return kFALSE;
}
if (kFALSE == pPadClean->init()) {
Error("init", "HRichPadClean did not initialize");
return kFALSE;
}
return kTRUE;
}
Bool_t
HRichAnalysis::reinit()
{
UInt_t i, j, k, m, n;
Int_t iMatrixSurface = 0;
Int_t iPartSurface = 0;
Int_t iMaskSize = 0;
Int_t iMaskSizeSquared = 0;
if (kFALSE == initParameters()) {
Error("reinit", "Parameters not initialized");
return kFALSE;
}
pRingFind->init(this);
iPadActive.Set(maxCols * maxRows);
for (i = 0; i < (UInt_t)maxCols * maxRows; ++i)
if (((HRichGeometryPar*)fpGeomPar)->getPadsPar()->getPad(i)->getPadActive() > 0)
iPadActive[i] = 1;
else
iPadActive[i] = 0;
iMaskSize = ((HRichAnalysisPar*)fpAnalysisPar)->iRingMaskSize;
iMaskSizeSquared = iMaskSize * iMaskSize;
for (k = 0; k < (UInt_t)iMaskSizeSquared; ++k)
if (1 == ((HRichAnalysisPar*)fpAnalysisPar)->iRingMask[k])
iMatrixSurface++;
for (j = 0; j < (UInt_t)maxRows; j++)
for (i = 0; i < (UInt_t)maxCols; i++) {
iPartSurface = 0;
for (k = 0; k < (UInt_t)iMaskSizeSquared; k++) {
m = (k % iMaskSize) - iMaskSize / 2;
n = (k / iMaskSize) - iMaskSize / 2;
if (!IsOut(i, j, m, n))
if (1 == ((HRichAnalysisPar*)fpAnalysisPar)->iRingMask[k])
iPartSurface++;
}
((HRichGeometryPar*)fpGeomPar)->getPadsPar()->
getPad(i, j)->setAmplitFraction((Float_t)iPartSurface / (Float_t)iMatrixSurface);
}
return kTRUE;
}
Bool_t
HRichAnalysis::initParameters()
{
HRichGeometryPar *pGeomPar = getGeometryPar();
iRingNrTot = 0;
allPairNrTot = 0;
sectorPairNrTot = 0;
maxFiredTotalPads = 3000;
maxCols = pGeomPar->getColumns();
maxRows = pGeomPar->getRows();
if (pLeftBorder)
delete [] pLeftBorder;
pLeftBorder = new Short_t[maxRows];
if (pRightBorder)
delete [] pRightBorder;
pRightBorder = new Short_t[maxRows];
for (Int_t row = 0; row < maxRows; ++row) {
Int_t col = 0;
Int_t padOffset = row * maxCols;
while (!pGeomPar->getPadsPar()->getPad(col + padOffset)->getPadActive() && col < maxCols)
++col;
if (col == maxCols) {
maxRows = row;
break;
}
pLeftBorder[row] = col;
while (pGeomPar->getPadsPar()->getPad(col + padOffset)->getPadActive() && col < maxCols)
++col;
pRightBorder[row] = col - 1;
}
maxPads = maxRows * maxCols;
pPads = new HRichPadSignal * [6];
for (Int_t j = 0; j < 6; ++j) {
pSectorPads = pPads[j] = new HRichPadSignal[maxPads];
for (Int_t i = 0; i < maxPads; pSectorPads[i++].clear());
}
pSectorPads = pPads[iActiveSector];
return kTRUE;
}
Bool_t
HRichAnalysis::finalize()
{
return kTRUE;
}
Int_t
HRichAnalysis::execute()
{
HRichCal* pCal = NULL;
Bool_t skipSector[] = {0, 0, 0, 0, 0, 0};
Int_t iRingNrTotEvt = 0;
Int_t fPadFiredTot = 0;
Int_t lastRingNr = 0;
Int_t ampl = -1;
Int_t sector = -1;
Int_t padnr = -1;
const Int_t eventNr = gHades->getCurrentEvent()->getHeader()->getEventSeqNumber();
Int_t j;
secWithCurrent = -1;
Reset();
for (sector = 0; sector < 6; ++sector) {
if (fPadFired[sector] > 0) {
fPadFired[sector] = 0;
HRichPadSignal* pSecPads = pPads[sector];
for (j = 0; j < maxPads; pSecPads[j++].clear());
}
}
fIter->Reset();
while (NULL != (pCal = (HRichCal *)fIter->Next())) {
if ((ampl = (Int_t)pCal->getCharge()) > 0) {
sector = pCal->getSector();
fPadFired[sector]++;
padnr = pCal->getCol() + pCal->getRow() * maxCols;
if(padnr >= maxPads) {
Warning("execute()","pPads out of range ! signal will be skipped.");
continue;
}
pPads[sector][padnr].setAmplitude(ampl);
}
}
for (sector = 6; sector--; fPadFiredTot += fPadFired[sector]) {
if (fPadFiredTot > maxFiredTotalPads) {
Warning("execute", "Analysis of event %d skipped: too many pads fired %i/%i",
eventNr, fPadFiredTot, maxFiredTotalPads);
fIter->Reset();
while (NULL != (pCal = static_cast<HRichCal*>(fIter->Next()))) {
pCal->setIsCleanedSector(kTRUE);
}
return 0;
}
if (fPadFired[sector] >= fpAnalysisPar->maxFiredSectorPads) {
Warning("execute", "Event %d: too many fired pads in sector %i. %i/%i",
eventNr, sector, fPadFired[sector], fpAnalysisPar->maxFiredSectorPads);
fIter->Reset();
while (NULL != (pCal = static_cast<HRichCal*>(fIter->Next()))) {
if (pCal->getSector() == sector)
pCal->setIsCleanedSector(kTRUE);
}
skipSector[sector] = kTRUE;
}
}
lastRingNr = iRingNrTot;
for (sector = 0; sector < 6; ++sector) {
if (kTRUE == skipSector[sector])
continue;
if (((HRichGeometryPar*)fpGeomPar)->getSectorActive(sector) > 0) {
Reset();
SetActiveSector(sector);
iPadFiredNr = fPadFired[sector];
iPadCleanedNr = pPadClean->Execute(this);
iLabelNr = pPadLabel->Execute(this);
if (0 == ((HRichAnalysisPar*)fpAnalysisPar)->isActiveLabelPads) {
iLabelNr = 1;
pLabelArea = new HRichLabel[1];
pLabelArea[0].iLeftX = 0;
pLabelArea[0].iRightX = maxCols;
pLabelArea[0].iLowerY = 0;
pLabelArea[0].iUpperY = maxRows;
pLabelArea[0].iLabeledPadsNr = maxPads;
pLabelArea[0].iFiredPadsNr = iPadFiredNr;
}
iRingNr = pRingFind->Execute(this);
iRingNrTot += iRingNr;
iRingNrTotEvt += iRingNr;
if (iRingNr > 1) {
++sectorPairNrTot;
}
updateHeaders(sector, eventNr);
updateHits(sector);
}
}
fIter->Reset();
while (NULL != (pCal = (HRichCal *)fIter->Next())) {
sector = pCal->getSector();
padnr = pCal->getCol() + pCal->getRow() * maxCols;
if(padnr >= maxPads) {
continue;
}
pCal->setIsCleanedSingle(pPads[sector][padnr].getIsCleanedSingle());
pCal->setIsCleanedHigh(pPads[sector][padnr].getIsCleanedHigh());
}
if (iRingNrTot > lastRingNr + 1) {
++allPairNrTot;
}
if (kSkipEvtIfNoRing && iRingNrTotEvt == 0)
return kSkipEvent;
return 0;
}
void
HRichAnalysis::updateHits(Int_t nSec)
{
HRichHit* hit = NULL;
HLocation loc;
for (Int_t i = 0; i < iRingNr; i++) {
loc.set(1, nSec);
hit = (HRichHit*)m_pHitCat->getNewSlot(loc);
if (NULL != hit) {
hit = new(hit) HRichHit;
*hit = pRings[i];
hit->setEventNr(gHades->getCurrentEvent()->getHeader()->getEventSeqNumber());
hit->setSector(nSec);
if (hit->fX > 1000.)
Info("updateHits", "fX : %f, %i fY : %f, %i",
hit->fX, hit->iRingX, hit->fY, hit->iRingY);
HRichPad* pad = ((HRichGeometryPar*)fpGeomPar)->getPadsPar()->
getPad((UInt_t)hit->iRingX, (UInt_t)hit->iRingY);
hit->fPhi = pad->getPhi(nSec);
hit->fTheta = pad->getTheta();
Float_t delta = hit->fPadX - float(hit->iRingX);
if(delta > 0.01F) {
HRichPad* pad1 = ((HRichGeometryPar*)fpGeomPar)->getPadsPar()->
getPad((UInt_t)hit->iRingX+1, (UInt_t)hit->iRingY);
hit->fPhi += (pad1->getPhi(nSec) - hit->fPhi)*delta;
} else if(delta < -0.01F){
HRichPad* pad1 = ((HRichGeometryPar*)fpGeomPar)->getPadsPar()->
getPad((UInt_t)hit->iRingX-1, (UInt_t)hit->iRingY);
hit->fPhi -= (pad1->getPhi(nSec) - hit->fPhi)*delta;
}
delta = hit->fPadY - float(hit->iRingY);
if(delta > 0.01F) {
HRichPad* pad1 = ((HRichGeometryPar*)fpGeomPar)->getPadsPar()->
getPad((UInt_t)hit->iRingX, (UInt_t)hit->iRingY+1);
hit->fTheta += (pad1->getTheta() - hit->fTheta)*delta;
} else if(delta < -0.01F){
HRichPad* pad1 = ((HRichGeometryPar*)fpGeomPar)->getPadsPar()->
getPad((UInt_t)hit->iRingX, (UInt_t)hit->iRingY-1);
hit->fTheta -= (pad1->getTheta() - hit->fTheta)*delta;
}
}
}
}
void
HRichAnalysis::updateHeaders(const Int_t nSec,
const Int_t nEvNr)
{
HRichHitHeader* hithdr = NULL;
HLocation loc;
fIterHitHeader->Reset();
hithdrLoop = NULL;
while (NULL != (hithdrLoop = (HRichHitHeader *)fIterHitHeader->Next())) {
secWithCurrent = hithdrLoop->getSector();
if (nSec == secWithCurrent) {
hithdr = hithdrLoop;
break;
}
}
if (hithdr == NULL) {
loc.set(1, nSec);
hithdr = (HRichHitHeader*)(m_pHitHdrCat)->getNewSlot(loc);
if (NULL != hithdr)
hithdr = new(hithdr) HRichHitHeader;
}
if (NULL != hithdr) {
hithdr->setSector(nSec);
hithdr->setEventNr(nEvNr);
hithdr->setExecutedAnalysis(1);
hithdr->iPadFiredNr = iPadFiredNr;
hithdr->iPadCleanedNr = iPadCleanedNr;
hithdr->iClusterCleanedNr = iClusterCleanedNr;
hithdr->iClustersCleaned = iClustersCleaned;
hithdr->iRingNr = iRingNr;
hithdr->iFakePad = iFakePad;
hithdr->iFakeLocalMax4 = iFakeLocalMax4;
hithdr->iFakeLocalMax8 = iFakeLocalMax8;
} else {
Warning("updateHeaders", "No Header Object could be created");
}
}
void
HRichAnalysis::SetActiveSector(Int_t sectornr)
{
if (sectornr == iActiveSector)
return;
if (sectornr > 5 || sectornr < 0) {
Error("SetActiveSector", "Sector number %d out of range (0..5)!", sectornr);
return;
} else if (0 == ((HRichGeometryPar*)fpGeomPar)->getSectorActive(sectornr)) {
Error("SetActiveSector", "Sector number %d is not present (and cannot be active)!", sectornr);
return;
} else {
iActiveSector = sectornr;
pSectorPads = pPads[sectornr];
}
}
Int_t
HRichAnalysis::SetNextSector()
{
Int_t i = iActiveSector;
while (i < 6) {
i++;
if (((HRichGeometryPar*)fpGeomPar)->getSectorActive(i) > 0) {
pSectorPads = pPads[i];
return (iActiveSector = i);
}
}
Error("SetNextSector", "No more sectors (last sector %d)", iActiveSector);
return iActiveSector;
}
Bool_t
HRichAnalysis::IsOut(Int_t nowPad, Int_t dx, Int_t dy)
{
if (nowPad <= 0)
return 1;
Int_t x = nowPad % maxCols;
Int_t y = nowPad / maxCols;
if ((x + dx) >= 0 && (x + dx) < maxCols &&
(y + dy) >= 0 && (y + dy) < maxRows &&
iPadActive[nowPad+dx + maxCols*dy]) {
return 0;
} else {
return 1;
}
}
Bool_t
HRichAnalysis::IsOut(Int_t x, Int_t y, Int_t dx, Int_t dy)
{
if ((x + dx) >= 0 && (x + dx) < maxCols &&
(y + dy) >= 0 && (y + dy) < maxRows &&
iPadActive[x + dx + maxCols *(y + dy)]) {
return 0;
} else {
return 1;
}
}
void
HRichAnalysis::Streamer(TBuffer &R__b)
{
if (R__b.IsReading()) {
Version_t R__v = R__b.ReadVersion();
if (R__v) { }
R__b >> iPadFiredNr;
R__b >> iPadCleanedNr;
R__b >> iRingNr;
R__b >> pRings;
R__b >> iClusterCleanedNr;
iClustersCleaned.Streamer(R__b);
R__b >> iFakePad;
R__b >> iFakeLocalMax4;
R__b >> iFakeLocalMax8;
} else {
R__b.WriteVersion(HRichAnalysis::IsA());
R__b << iPadFiredNr;
R__b << iPadCleanedNr;
R__b << iRingNr;
R__b << pRings;
R__b << iClusterCleanedNr;
iClustersCleaned.Streamer(R__b);
R__b << iFakePad;
R__b << iFakeLocalMax4;
R__b << iFakeLocalMax8;
}
}