#include "hades.h"
#include "hevent.h"
#include "heventheader.h"
#include "hiterator.h"
#include "hlinearcategory.h"
#include "hlocation.h"
#include "hmatrixcategory.h"
#include "hparcond.h"
#include "hparset.h"
#include "hrichanalysispar.h"
#include "hrichanalysissim.h"
#include "hrichcalsim.h"
#include "hrichdetector.h"
#include "hrichgeometrypar.h"
#include "hrichhitsim.h"
#include "hrichpadclean.h"
#include "hrichpadlabel.h"
#include "hrichringfindsim.h"
#include "hrichtrack.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "richdef.h"
#include <iomanip>
#include <iostream>
using namespace std;
ClassImp(HRichAnalysisSim)
HRichAnalysisSim::HRichAnalysisSim(void)
{
pRingFindSim = new HRichRingFindSim;
pRichCalSim = NULL;
fIterHitHeader = NULL;
m_pTrackCat = NULL;
pRings = NULL;
}
HRichAnalysisSim::HRichAnalysisSim(const Text_t* name,
const Text_t* title,
Bool_t kSkip)
: HRichAnalysis(name, title, kSkip)
{
pRingFindSim = new HRichRingFindSim;
pRichCalSim = NULL;
fIterHitHeader = NULL;
m_pTrackCat = NULL;
pRings = NULL;
kSkipEvtIfNoRing = kSkip;
}
HRichAnalysisSim::~HRichAnalysisSim(void)
{
if (NULL != pRingFindSim) {
delete pRingFindSim;
pRingFindSim = NULL;
}
if (NULL != fIterHitHeader) {
delete fIterHitHeader;
fIterHitHeader = NULL;
}
}
Bool_t
HRichAnalysisSim::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->buildMatrixCat("HRichCalSim", 1);
if (NULL == m_pCalCat) {
Error("init", "Pointer to HRichCalSim 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->buildLinearCat("HRichHitSim");
if (NULL == m_pHitCat) {
Error("init", "Pointer to HRichHitSim 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;
}
m_pTrackCat = gHades->getCurrentEvent()->getCategory(catRichTrack);
if (NULL == m_pTrackCat) {
m_pTrackCat = pRichDet->buildLinearCat("HRichTrack", 1002);
if (NULL == m_pTrackCat) {
Error("init", "Pointer to HRichTrack category is NULL");
return kFALSE;
} else gHades->getCurrentEvent()
->addCategory(catRichTrack, m_pTrackCat , "Rich");
}
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 HRichAnalysisSim::reinit()
{
UInt_t i, j, k, m, n;
if (kFALSE == initParameters()) {
Error("init", "Parameters not initialized");
return kFALSE;
}
pRingFindSim->init(this);
iPadActive.Set(maxCols * maxRows);
for (i = 0 ; i < (UInt_t)maxCols * maxRows; ++i)
if (getGeometryPar()->getPadsPar()->getPad(i)->getPadActive() > 0)
iPadActive[i] = 1;
else iPadActive[i] = 0;
Int_t iMatrixSurface = 0, iPartSurface = 0;
Int_t iMaskSize = getAnalysisPar()->iRingMaskSize;
Int_t iMaskSizeSquared = iMaskSize * iMaskSize;
for (k = 0; k < (UInt_t)iMaskSizeSquared; k++)
if (getAnalysisPar()->iRingMask[k] == 1) 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 (getAnalysisPar()->iRingMask[k] == 1) iPartSurface++;
}
getGeometryPar()->getPadsPar()->
getPad(i, j)->setAmplitFraction((Float_t)iPartSurface / (Float_t)iMatrixSurface);
}
return kTRUE;
}
Bool_t
HRichAnalysisSim::initParameters()
{
iRingNrTot = 0;
allPairNrTot = 0;
sectorPairNrTot = 0;
maxFiredTotalPads = 3000;
maxCols = getGeometryPar()->getColumns();
maxRows = getGeometryPar()->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 (!getGeometryPar()->getPadsPar()->getPad(col + padOffset)->getPadActive() && col < maxCols) ++col;
if (col == maxCols) {
maxRows = row;
break;
}
pLeftBorder[row] = col;
while (getGeometryPar()->getPadsPar()->getPad(col + padOffset)->getPadActive() && col < maxCols) ++col;
pRightBorder[row] = col - 1;
}
maxPads = maxRows * maxCols;
Int_t i;
Int_t j;
if (pPads) {
for (j = 0; j < 6; j++)
if (pPads[j])
delete [] pPads[j];
delete [] pPads;
}
pPads = new HRichPadSignal * [6];
for (j = 0; j < 6; j++) {
pSectorPads = pPads[j] = new HRichPadSignal[maxPads];
for (i = 0; i < maxPads; pSectorPads[i++].clear());
}
pSectorPads = pPads[iActiveSector];
return kTRUE;
}
Int_t
HRichAnalysisSim::getPadsNTrack1(Int_t padx, Int_t pady, Int_t sec)
{
HLocation loc1;
loc1.set(3, sec, pady, padx);
pRichCalSim = (HRichCalSim*)((HMatrixCategory*)getCalCat())->getObject(loc1);
if (NULL != pRichCalSim)
return pRichCalSim->getTrackId(0);
return 0;
}
Int_t
HRichAnalysisSim::getPadsNTrack2()
{
if (NULL != pRichCalSim)
return pRichCalSim->getTrackId(1);
return 0;
}
Int_t
HRichAnalysisSim::getFlag(Int_t index)
{
HRichTrack *trk = NULL;
trk = (HRichTrack*)((HLinearCategory*)m_pTrackCat)->getObject(index);
if (NULL != trk)
return trk->getFlag();
return -1;
}
Int_t
HRichAnalysisSim::getTrack(Int_t index)
{
HRichTrack *trk = NULL;
trk = ((HRichTrack*)((HLinearCategory*)m_pTrackCat)->getObject(index));
if (NULL != trk)
return trk->getTrack();
return -10;
}
Int_t
HRichAnalysisSim::execute()
{
HRichCal* pCal = NULL;;
Int_t iRingNrTotEvt = 0;
Int_t fPadFiredTot = 0;
Int_t sector = -1;
Int_t padnr = -1;
Int_t j;
Int_t ampl;
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 = static_cast<HRichCal*>(fIter->Next()))) {
if ((ampl = (Int_t)pCal->getCharge()) > 0) {
fPadFired[pCal->getSector()]++;
padnr = pCal->getCol() + pCal->getRow() * maxCols;
if(padnr >= maxPads) {
Warning("execute()","pPads out of range ! signal will be skipped.");
continue;
}
pPads[pCal->getSector()][padnr].setAmplitude(ampl);
}
}
fPadFiredTot = 0;
for (sector = 0; sector < 6; ++sector) {
fPadFiredTot += fPadFired[sector];
if (fPadFired[sector] >= fpAnalysisPar->maxFiredSectorPads) {
Warning("execute", "To many fired pads in sector %i. %i/%i",
sector, fPadFired[sector], fpAnalysisPar->maxFiredSectorPads);
fIter->Reset();
while (NULL != (pCal = static_cast<HRichCal*>(fIter->Next()))) {
if (pCal->getSector() == sector)
pCal->setIsCleanedSector(kTRUE);
}
return 0;
}
}
if (fPadFiredTot > maxFiredTotalPads) {
Warning("execute", "Analysis of event skipped: too many pads fired %i/%i",
fPadFiredTot, maxFiredTotalPads);
fIter->Reset();
while (NULL != (pCal = static_cast<HRichCal*>(fIter->Next()))) {
pCal->setIsCleanedSector(kTRUE);
}
return 0;
}
for (sector = 0; sector < 6; ++sector)
if (getGeometryPar()->getSectorActive(sector) > 0) {
Reset();
SetActiveSector(sector);
iPadFiredNr = fPadFired[sector];
iPadCleanedNr = pPadClean->Execute(this);
iLabelNr = pPadLabel->Execute(this);
if (0 == getAnalysisPar()->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;
pLabelArea[0].iSignature = 0;
}
iRingNr = pRingFindSim->Execute(this);
iRingNrTot += iRingNr;
iRingNrTotEvt += iRingNr;
updateHeaders(sector, gHades->getCurrentEvent()->getHeader()->getEventSeqNumber());
updateHits(sector);
}
fIter->Reset();
while (NULL != (pCal = static_cast<HRichCal*>(fIter->Next()))) {
padnr = pCal->getCol() + pCal->getRow() * maxCols;
if(padnr >= maxPads){
continue;
}
pCal->setIsCleanedSingle(pPads[pCal->getSector()][padnr].getIsCleanedSingle());
pCal->setIsCleanedHigh(pPads[pCal->getSector()][padnr].getIsCleanedHigh());
}
if (kTRUE == kSkipEvtIfNoRing && 0 == iRingNrTotEvt)
return kSkipEvent;
return 0;
}
void
HRichAnalysisSim::updateHits(Int_t nSec)
{
HRichHitSim *hit = NULL;
HLocation loc;
for (Int_t i = 0; i < iRingNr; i++) {
loc.set(1, nSec);
hit = (HRichHitSim *)m_pHitCat->getNewSlot(loc, NULL);
if (hit != NULL) {
hit = new(hit) HRichHitSim;
HRichHitSim & ring = pRings[i];
*hit = ring;
hit->setSector(nSec);
HRichPad * pad = getGeometryPar()->getPadsPar()->
getPad((UInt_t)ring.getRingCenterX(), (UInt_t)ring.getRingCenterY());
if (pad) {
hit->setTheta(pad-> getTheta());
hit->setPhi(pad->getPhi(nSec));
} else {
hit->fPhi = -10;
hit->fTheta = -10;
}
Int_t k = 0;
while (ring.iRingTrack[k]) {
if (hit->track1 == 0) {
hit->track1 = ring.iRingTrack[k];
hit->flag1 = ring.iRingFlag[k];
(hit->weigTrack1)++;
} else {
if (ring.iRingTrack[k] == hit->track1) {
(hit->weigTrack1)++;
} else {
if (hit->track2 == 0) {
hit->track2 = ring.iRingTrack[k];
hit->flag2 = ring.iRingFlag[k];
(hit->weigTrack2)++;
} else {
if (ring.iRingTrack[k] == hit->track2) {
(hit->weigTrack2)++;
} else {
if (hit->track3 == 0) {
hit->track3 = ring.iRingTrack[k];
hit->flag3 = ring.iRingFlag[k];
(hit->weigTrack3)++;
} else (hit->weigTrack3)++;
}
}
}
}
k++;
}
sortTracks(hit);
}
}
}
void
HRichAnalysisSim::sortTracks(HRichHitSim *hit)
{
Int_t tmp[3];
tmp[0] = hit->track1;
tmp[1] = hit->flag1;
tmp[2] = hit->weigTrack1;
if (hit->weigTrack1 < hit->weigTrack2) {
if (hit->weigTrack2 > hit->weigTrack3) {
hit->track1 = hit->track2;
hit->flag1 = hit->flag2;
hit->weigTrack1 = hit->weigTrack2;
if (tmp[2] < hit->weigTrack3) {
hit->track2 = hit->track3;
hit->flag2 = hit->flag3;
hit->weigTrack2 = hit->weigTrack3;
hit->track3 = tmp[0];
hit->flag3 = tmp[1];
hit->weigTrack3 = tmp[2];
} else {
hit->track2 = tmp[0];
hit->flag2 = tmp[1];
hit->weigTrack2 = tmp[2];
}
} else {
hit->track1 = hit->track3;
hit->flag1 = hit->flag3;
hit->weigTrack1 = hit->weigTrack3;
hit->track3 = tmp[0];
hit->flag3 = tmp[1];
hit->weigTrack3 = tmp[2];
}
} else {
if (hit->weigTrack1 < hit->weigTrack3) {
hit->track1 = hit->track3;
hit->flag1 = hit->flag3;
hit->weigTrack1 = hit->weigTrack3;
hit->track3 = hit->track2;
hit->flag3 = hit->flag2;
hit->weigTrack3 = hit->weigTrack2;
hit->track2 = tmp[0];
hit->flag2 = tmp[1];
hit->weigTrack2 = tmp[2];
} else {
if (hit->weigTrack2 < hit->weigTrack3) {
tmp[0] = hit->track2;
tmp[1] = hit->flag2;
tmp[2] = hit->weigTrack2;
hit->track2 = hit->track3;
hit->flag2 = hit->flag3;
hit->weigTrack2 = hit->weigTrack3;
hit->track3 = tmp[0];
hit->flag3 = tmp[1];
hit->weigTrack3 = tmp[2];
}
}
}
if(hit->track2 < 0 && hit->track3 > 0 &&
hit->weigTrack2 > 0 && hit->weigTrack2 == hit->weigTrack3) {
tmp[0] = hit->track2;
tmp[1] = hit->flag2;
tmp[2] = hit->weigTrack2;
hit->track2 = hit->track3;
hit->flag2 = hit->flag3;
hit->weigTrack2 = hit->weigTrack3;
hit->track3 = tmp[0];
hit->flag3 = tmp[1];
hit->weigTrack3 = tmp[2];
}
if(hit->track1 < 0 && hit->track2 > 0 &&
hit->weigTrack2 > 0 && hit->weigTrack1 == hit->weigTrack2) {
tmp[0] = hit->track1;
tmp[1] = hit->flag1;
tmp[2] = hit->weigTrack1;
hit->track1 = hit->track2;
hit->flag1 = hit->flag2;
hit->weigTrack1 = hit->weigTrack2;
hit->track2 = tmp[0];
hit->flag2 = tmp[1];
hit->weigTrack2 = tmp[2];
}
}