using namespace std;
#include "TRandom.h"
#include <time.h>
#include "hwalldigitizer.h"
#include "walldef.h"
#include "hdebug.h"
#include "hades.h"
#include "hiterator.h"
#include "hruntimedb.h"
#include "hspectrometer.h"
#include "hwalldetector.h"
#include "hwalldigipar.h"
#include "hwallgeompar.h"
#include "hgeantwall.h"
#include "hwallrawsim.h"
#include "hevent.h"
#include "hcategory.h"
#include "hlocation.h"
#include "hwallrawsimfilter.h"
#include <iostream>
#include <iomanip>
void HWallDigitizer::initParContainer() {
pWallDigiPar=(HWallDigiPar *)gHades->getRuntimeDb()->getContainer("WallDigiPar");
}
HWallDigitizer::HWallDigitizer(void) {
fGeantCat = 0;
fRawCat = 0;
pWallDigiPar = 0;
pWallGeomPar = 0;
fLoc.set(1,0);
iterGeant = 0;
iterWallRaw = 0;
}
HWallDigitizer::HWallDigitizer(const Text_t *name,const Text_t *title) :
HReconstructor(name,title) {
fGeantCat = 0;
fRawCat = 0;
pWallDigiPar = 0;
pWallGeomPar = 0;
fLoc.set(1,0);
iterGeant = 0;
iterWallRaw = 0;
}
HWallDigitizer::~HWallDigitizer(void) {
}
Bool_t HWallDigitizer::init(void) {
time_t curtime;
fGeantCat = gHades->getCurrentEvent()->getCategory(catWallGeantRaw);
if (!fGeantCat) {
Warning("init()","No category catWallGeantRaw found!") ;
} else iterGeant = (HIterator *)fGeantCat->MakeIterator("native");
fRawCat = gHades->getCurrentEvent()->getCategory(catWallRaw);
if (!fRawCat) {
HWallDetector* wall=(HWallDetector*)(gHades->getSetup()->getDetector("Wall"));
if(wall){
fRawCat=wall->buildMatrixCategory("HWallRawSim",0.5F);
if (fRawCat) gHades->getCurrentEvent()->addCategory(catWallRaw,fRawCat,"Wall");
}
}
if(fRawCat) iterWallRaw = (HIterator *)fRawCat ->MakeIterator("native");
else Warning("init()","No category catWallRaw !") ;
if(fGeantCat&&fRawCat) initParContainer();
time(&curtime);
return kTRUE;
}
Int_t HWallDigitizer::execute(void) {
if(!iterGeant || !iterWallRaw) return 0;
const Float_t quantEff = 0.24;
const Float_t photYield = 11100.0;
const Float_t timeResWall[3] = {0.42, 0.52, 0.67};
const Float_t relAmplResol = 0.08;
const Float_t minEnerRelease = 1.8;
HGeantWall* geant = 0;
HWallRawSim* raw = 0;
HWallRawSimFilter fRawFilter;
TString cellName,str;
Float_t slopeTDC;
Int_t thrCFD, thrADC;
Int_t cell = -1;
fLoc.set(1,-1);
Int_t numTrack, numTrack1 = -1, numTrack2 = -1;
Float_t trackLen;
Float_t time, charge;
Int_t timeCh, chargeCh;
Float_t prevTime, prevCharge;
Float_t timeResol, amplResol, chargeRef;
Float_t geaTof = 0.;
Float_t geaTof1[383];
Float_t geaTof2[383];
Float_t geaEner = 0.;
Float_t geaX = 0.;
Float_t geaY = 0.;
Float_t geaMom = 0.;
iterGeant->Reset();
while ((geant = (HGeantWall *)iterGeant->Next()) != 0) {
fLoc[0] = geant->getCell();
raw = (HWallRawSim*) fRawCat->getObject(fLoc);
if(raw) {
raw->incNHit();
numTrack1 = raw->getNTrack1();
numTrack2 = raw->getNTrack2();
prevTime = raw->getTime();
prevCharge = raw->getCharge();
}
else {
prevTime = 100000.;
prevCharge = 0.;
raw = (HWallRawSim*) fRawCat->getNewSlot(fLoc);
if(!raw) continue;
raw = new(raw) HWallRawSim;
raw->setNHit(1);
}
cell = fLoc[0];
slopeTDC = pWallDigiPar->getTDC_Slope(cell);
geant->getHit(geaEner,geaX,geaY,geaTof,geaMom,trackLen);
numTrack = geant->getTrack();
cellName = pWallGeomPar->getCellName(fLoc[0]);
str=cellName.Replace(0,3,"",0);
str.Replace(1,3,"",0);
time = geaTof;
timeResol = timeResWall[atoi(str)-1];
time = gRandom->Gaus(time,timeResol);
timeCh = (Int_t) (time/slopeTDC);
if (timeCh < 0) timeCh = 0;
if (timeCh > 4095) timeCh = 4095;
time = ( Float_t ) (timeCh);
if(raw->getNHit() > 1){
if(geaTof < geaTof1[fLoc[0]]){
numTrack2 = numTrack1;
numTrack1 = numTrack;
geaTof2[fLoc[0]] = geaTof1[fLoc[0]];
geaTof1[fLoc[0]] = geaTof;
} else {
if(geaTof < geaTof2[fLoc[0]]){
numTrack2 = numTrack;
geaTof2[fLoc[0]] = geaTof;
}
}
}
if(time > prevTime) time = prevTime;
charge = geaEner*photYield*quantEff;
amplResol = charge*relAmplResol;
charge = gRandom->Gaus(charge,amplResol);
chargeRef = photYield*quantEff*minEnerRelease;
chargeCh = (Int_t) ((charge/chargeRef)*256. + prevCharge);
if (chargeCh < 0) chargeCh = 0;
if (chargeCh > 4095) chargeCh = 4095;
charge = (Float_t)chargeCh;
raw->setTime(time);
raw->setCharge(charge);
raw->setCell((Int_t) fLoc[0]);
if(raw->getNHit() > 1){
raw->setNTrack1(numTrack1);
raw->setNTrack2(numTrack2);
} else {
raw->setNTrack1(numTrack);
raw->setNTrack2(-1);
geaTof1[fLoc[0]] = geaTof;
geaTof2[fLoc[0]] = 100000.;
}
}
iterWallRaw->Reset();
while ( (raw = (HWallRawSim *)iterWallRaw->Next()) != NULL) {
fLoc[0] = raw->getCell();
time = raw->getTime();
charge = raw->getCharge();
if(((Int_t)time) >= 4095) raw->setTime(0.0);
cell=fLoc[0];
thrCFD = (Int_t)pWallDigiPar->getCFD_Threshold(cell);
thrADC = (Int_t)pWallDigiPar->getADC_Threshold(cell);
if(((Int_t)charge) < thrCFD){
raw->setTime(0.0);
if(((Int_t)charge) < thrADC){
raw->setCharge(0.0);
}
}
}
fRawCat->filter(fRawFilter);
return 0;
}
ClassImp(HWallDigitizer)