#include "hemccalibrater.h"
#include "emcdef.h"
#include "hemcraw.h"
#include "hemccalqa.h"
#include "hemccalsim.h"
#include "hemcdetector.h"
#include "hemccalpar.h"
#include "hemccellstatuspar.h"
#include "hemccalibraterpar.h"
#include "hstartdef.h"
#include "hstart2hit.h"
#include "hades.h"
#include "hcategory.h"
#include "hdebug.h"
#include "hevent.h"
#include "hiterator.h"
#include "hruntimedb.h"
#include "hrun.h"
#include "hspectrometer.h"
#include <TMath.h>
#include <iomanip>
#include <iostream>
#include <stdlib.h>
#include <string.h>
using namespace std;
#define SELECT_FIRST_OR_LAST // Select the first or the last pair, if commented, the last pair will be stored
ClassImp(HEmcCalibrater)
HEmcCalibrater::HEmcCalibrater(void) {
pRawCat = NULL;
pCalCat = NULL;
iter = NULL;
pCalpar = NULL;
embedding = 0;
}
HEmcCalibrater::HEmcCalibrater(const Text_t *name, const Text_t *title)
: HReconstructor(name, title) {
pRawCat = NULL;
pCalCat = NULL;
iter = NULL;
pCalpar = NULL;
embedding = 0;
}
HEmcCalibrater::~HEmcCalibrater(void) {
if (NULL != iter) {
delete iter;
iter = NULL;
}
}
Bool_t HEmcCalibrater::init(void) {
HEmcDetector* pDet = (HEmcDetector*)gHades->getSetup()->getDetector("Emc");
if (!pDet) {
Error("init", "No Emc found.");
return kFALSE;
}
pCalpar = (HEmcCalPar*)gHades->getRuntimeDb()->getContainer("EmcCalPar");
if (!pCalpar) return kFALSE;
pStatuspar = (HEmcCellStatusPar*)gHades->getRuntimeDb()->getContainer("EmcCellStatusPar");
if (!pStatuspar) return kFALSE;
pCalibpar = (HEmcCalibraterPar*)gHades->getRuntimeDb()->getContainer("EmcCalibraterPar");
if (!pCalibpar) return kFALSE;
pRawCat = gHades->getCurrentEvent()->getCategory(catEmcRaw);
if (!pRawCat) {
Warning("init()", "HEmcRaw category not available!");
}
if(gHades->getEmbeddingMode()>0) embedding = 1;
pCalCat = (HCategory*)(((HEvent*)(gHades->getCurrentEvent()))->getCategory(catEmcCal));
if(pCalCat == NULL) {
if (embedding == 0) pCalCat = pDet->buildMatrixCategory("HEmcCal",0.5);
else pCalCat = pDet->buildMatrixCategory("HEmcCalSim",0.5);
if(pCalCat == NULL) return kFALSE;
gHades->getCurrentEvent()->addCategory(catEmcCal,pCalCat,"Emc");
}
pStartHitCat = (HCategory*)(((HEvent*)(gHades->getCurrentEvent()))->getCategory(catStart2Hit));
if(pStartHitCat == NULL) {
Info("HEmcCalibrater::init()", "No Start2Hit category, no Start time correction\n");
}
#ifdef WITHCALQA
pQACat = (HCategory*)(((HEvent*)(gHades->getCurrentEvent()))->getCategory(catEmcCalQA));
if(pQACat == NULL) {
pQACat = pDet->buildMatrixCategory("HEmcCalQA",0.5);
if(pQACat == NULL) return kFALSE;
gHades->getCurrentEvent()->addCategory(catEmcCalQA, pQACat,"Emc");
}
#endif
if(pRawCat) iter = (HIterator*)pRawCat->MakeIterator();
else iter = 0;
loc.set(2, 0, 0);
fActive = kTRUE;
return kTRUE;
}
Float_t twcCorrection(Float_t t, Float_t * twcPars)
{
if (t > twcPars[2])
return twcPars[0] + twcPars[1]/TMath::Sqrt(t - twcPars[2]);
return 0.0;
}
Int_t HEmcCalibrater::execute(void) {
if(!pRawCat) return 0;
HEmcRaw *pRaw = NULL;
HEmcCal *pCal = NULL;
Int_t sec = -1;
Int_t cell = 0;
Int_t pos = 0;
Char_t row = -1;
Char_t col = -1;
UInt_t nfasthits = 0, nslowhits=0;
Float_t rawTimeLeading = 0.;
Float_t rawTimeTrailing = 0.;
Float_t twcTimeLeading = 0.;
Float_t rawTimeLeadingS = 0.;
Float_t rawTimeTrailingS = 0.;
Float_t rawWidth = 0.;
Float_t calTimeLeading = 0.;
Float_t parCellData[8] = {1.F, 0.F, 0.F, 1.F, 0.F, 0.F, 0.F, 0.F};
Float_t start_time = 0.0;
if (pStartHitCat && pStartHitCat->getEntries()>0) {
HStart2Hit *pStartHit = (HStart2Hit *) pStartHitCat->getObject(0);
if (pStartHit && pStartHit->getFlag())
start_time = pStartHit->getTime();
}
Float_t diff_min = pCalibpar->getMatchWindowMin();
Float_t diff_max = pCalibpar->getMatchWindowMax();
iter->Reset();
while ((pRaw = (HEmcRaw *)iter->Next()) != 0) {
Float_t rawLead = 0.F;
Float_t rawTot = 0.F;
Float_t calTime = 0.F;
Float_t calEnergy = 0.F;
nfasthits = pRaw->getFastMultiplicity();
nslowhits = pRaw->getSlowMultiplicity();
if (nfasthits == 0 && nslowhits == 0)
continue;
pRaw->getAddress(sec, cell, row, col);
if (sec < 0)
continue;
pos = HEmcDetector::getPositionFromCell(cell);
if (pStatuspar->getCellStatus(sec, pos) == 0)
continue;
loc[0] = sec;
loc[1] = cell;
HEmcCalParCell &pPar = (*pCalpar)[sec][cell];
pPar.getData(parCellData);
Int_t hits_cnt = 0;
#ifdef WITHCALQA
HEmcCalQA *pQA = (HEmcCalQA*)pQACat->getObject(loc);
if (!pQA) {
pQA = (HEmcCalQA *)pQACat->getSlot(loc);
if (pQA) {
pQA = new(pQA) HEmcCalQA;
pQA->setAddress((Char_t)sec, (UChar_t)cell, row, col);
} else {
Error("execute()", "Can't get slot for CalQA sector=%i, cell=%i", sec, cell);
return -1;
}
} else {
Error("execute()", "Slot already exists for CalQA sector=%i, cell=%i", sec, cell);
return -1;
}
#endif
UInt_t j = 0;
UInt_t k = 0;
UInt_t l = 0;
for (UInt_t i = 0; i < nfasthits; ++i)
{
rawTimeLeading = pRaw->getFastTimeLeading(i);
if (0.0 == rawTimeLeading) continue;
#ifdef WITHCALQA
for (; j < nfasthits; ++j)
{
rawTimeTrailing = pRaw->getFastTimeTrailing(j);
if (0.0 == rawTimeTrailing) continue;
if (rawTimeTrailing < rawTimeLeading) continue;
pQA->addFastHit(rawTimeLeading * parCellData[0] + parCellData[1], (rawTimeTrailing - rawTimeLeading) * parCellData[0]);
break;
}
#endif
for (; k < nslowhits; ++k)
{
rawTimeLeadingS = pRaw->getSlowTimeLeading(k);
if (0.0 == rawTimeLeadingS) continue;
Float_t diff = rawTimeLeadingS - rawTimeLeading;
if (diff < diff_min) continue;
if (diff > diff_max) break;
#ifdef WITHCALQA
for (UInt_t _l = l; _l < nslowhits; ++_l)
{
rawTimeTrailingS = pRaw->getSlowTimeTrailing(_l);
if (0.0 == rawTimeTrailingS) continue;
if (rawTimeTrailingS < rawTimeLeadingS) continue;
pQA->addSlowHit(rawTimeLeadingS * parCellData[0] + parCellData[1], (rawTimeTrailingS - rawTimeLeadingS) * parCellData[0]);
break;
}
#endif
for (; l < nslowhits; ++l)
{
rawTimeTrailingS = pRaw->getSlowTimeTrailing(l);
if (0.0 == rawTimeTrailingS) continue;
if (rawTimeTrailingS < rawTimeLeadingS) continue;
rawWidth = rawTimeTrailingS - rawTimeLeading;
Float_t twc = twcCorrection(rawWidth, (parCellData+5));
twcTimeLeading = rawTimeLeading - twc;
calTimeLeading = parCellData[0] * twcTimeLeading + parCellData[1];
Float_t calWidth = rawWidth * parCellData[0];
#ifdef WITHCALQA
pQA->addMatchedHit(parCellData[0] * rawTimeLeading + parCellData[1], calWidth);
#endif
++hits_cnt;
#ifdef SELECT_FIRST_OR_LAST
if (hits_cnt > 1) break;
#endif
Float_t corrTime = calTimeLeading - start_time;
rawLead = rawTimeLeading;
rawTot = rawWidth;
calTime = corrTime;
calEnergy = parCellData[2]+exp( parCellData[3]+ calWidth*parCellData[4]);
break;
}
break;
}
}
if (hits_cnt)
{
pCal = (HEmcCal*)pCalCat->getObject(loc);
if (!pCal) {
pCal = (HEmcCal *)pCalCat->getSlot(loc);
if (pCal) {
if(embedding == 0) pCal = new(pCal) HEmcCal;
else pCal = new(pCal) HEmcCalSim;
pCal->setAddress((Char_t)sec, (UChar_t)cell, row, col);
} else {
Error("execute()", "Can't get slot sector=%i, cell=%i", sec, cell);
return -1;
}
} else {
Error("execute()", "Slot already exists for sector=%i, cell=%i", sec, cell);
return -1;
}
pCal->setLead(rawLead);
pCal->setWidth(rawTot);
pCal->setTime(calTime);
pCal->setEnergy(calEnergy);
pCal->setNHits(hits_cnt);
}
}
return 0;
}