#include "hmdcdrifttimepar.h"
#include "hades.h"
#include "hspectrometer.h"
#include "hruntimedb.h"
#include "hmdccal2parsim.h"
#include "hmdccal2par.h"
#include "hmdcgetcontainers.h"
ClassImp(HMdcDriftTimeParBin)
ClassImp(HMdcDriftTimeParAngle)
ClassImp(HMdcDriftTimeParMod)
ClassImp(HMdcDriftTimeParSec)
ClassImp(HMdcDriftTimePar)
HMdcDriftTimeParBin::HMdcDriftTimeParBin(void) {
c0 = 0.;
c1 = 0.;
c2 = 0.;
c3 = 0.;
d0 = 0.;
d1 = 0.;
d2 = 0.;
d3 = 0.;
}
void HMdcDriftTimeParBin::print(void) const {
printf("DriftTime = %.5g + %.5g*alpha + %.5g*dist + %.5g*alpha*dist\n",
c0,c1,c2,c3);
printf("ErrDriftTime = %.5g + %.5g*alpha + %.5g*dist + %.5g*alpha*dist\n",
c0err,c1err,c2err,c3err);
}
void HMdcDriftTimeParBin::printTable(void) const {
printf("| %10.4g | %10.4g | %10.4g | %10.4g || %10.4g | %10.4g | %10.4g | %10.4g |\n",
c0,c1,c2,c3,c0err,c1err,c2err,c3err);
}
HMdcDriftTimeParMod::HMdcDriftTimeParMod(Int_t nAnBins, Int_t nDsBins) {
nAngleBins = nAnBins;
lAngleBin = nAngleBins-1;
angleBinSize = 90./nAngleBins;
lastAngleBin = lAngleBin*angleBinSize;
invAngleBinSize = 1./angleBinSize;
nDistBins = nDsBins;
lDistBin = nDistBins-1;
setDistBinSize(0.1);
}
void HMdcDriftTimeParMod::setDistBinSize(Double_t dBSize) {
distBinSize = dBSize;
lastBinDist = lDistBin*distBinSize;
invDistBinSize = 1./distBinSize;
}
void HMdcDriftTimeParMod::setTimeBinSize(Double_t tBSize) {
timeBinSize = tBSize;
lastBinTime = lDistBin*timeBinSize;
invTimeBinSize = 1./timeBinSize;
}
Bool_t HMdcDriftTimeParMod::initContainer(HMdcCal2ParModSim& fC2PModS,HMdcCal2ParMod& fC2PMod,
Double_t slOut,Double_t scaleErr) {
slopeOutside = slOut;
for(Int_t anBin=0; anBin<nAngleBins; anBin++) {
for(Int_t dsBin=0; dsBin<nDistBins; dsBin++) {
HMdcDriftTimeParBin* bin=at(anBin)->at(dsBin);
if(bin==0) return kFALSE;
HMdcCal2ParAngleSim& rAng1 = fC2PModS[anBin];
HMdcCal2ParAngleSim& rAng2 = anBin<lAngleBin ? fC2PModS[anBin+1] : rAng1;
fillDriftTime(rAng1,rAng2,bin,anBin,dsBin);
fillDriftTimeErr(rAng1,rAng2,bin,anBin,dsBin,scaleErr);
HMdcCal2ParAngle& tAng1 = fC2PMod[anBin];
HMdcCal2ParAngle& tAng2 = anBin<lAngleBin ? fC2PMod[anBin+1] : tAng1;
fillDriftDist(tAng1,tAng2,bin,anBin,dsBin);
}
}
return kTRUE;
}
void HMdcDriftTimeParMod::print(void) {
for(Int_t anBin=0; anBin<nAngleBins; anBin++) {
for(Int_t dsBin=0; dsBin<nDistBins; dsBin++) {
HMdcDriftTimeParBin* bin=at(anBin)->at(dsBin);
if(bin==0) continue;
printf("%3i %3i ",anBin,dsBin);
bin->printTable();
}
}
}
void HMdcDriftTimeParMod::fillDriftTime(HMdcCal2ParAngleSim& rAng1,
HMdcCal2ParAngleSim& rAng2, HMdcDriftTimeParBin* bin,
Int_t anBin, Int_t dsBin) {
Double_t k0=90./angleBinSize-anBin;
if(dsBin==lDistBin) {
Double_t y1=rAng1.getDriftTime1(lDistBin);
Double_t y2=rAng2.getDriftTime1(lDistBin);
Double_t k1=y1-y2;
bin->setC0(y1-lastBinDist*slopeOutside - k0*k1);
bin->setC1(k1/angleBinSize);
bin->setC2(slopeOutside);
bin->setC3(0.);
} else {
Double_t y1=rAng1.getDriftTime1(dsBin);
Double_t y2=rAng2.getDriftTime1(dsBin);
Double_t y3=rAng2.getDriftTime1(dsBin+1);
Double_t y4=rAng1.getDriftTime1(dsBin+1);
Double_t k1=y1-y2;
Double_t k2=y1-y4;
Double_t k3=k1+y3-y4;
k1+=dsBin*k3;
bin->setC0(y1 + dsBin*k2 - k0*k1);
bin->setC1(k1/angleBinSize);
bin->setC2((k0*k3-k2)/distBinSize);
bin->setC3(-k3/distBinSize/angleBinSize);
}
}
void HMdcDriftTimeParMod::fillDriftTimeErr(HMdcCal2ParAngleSim& rAng1,
HMdcCal2ParAngleSim& rAng2, HMdcDriftTimeParBin* bin,Int_t anBin, Int_t dsBin,Double_t scaleErr) {
Double_t k0=90./angleBinSize-anBin;
if(dsBin==lDistBin) {
Double_t y1=rAng1.getDriftTime1Error(lDistBin) * scaleErr;
Double_t y2=rAng2.getDriftTime1Error(lDistBin) * scaleErr;
Double_t k1=y1-y2;
bin->setC0Err(y1 - k0*k1);
bin->setC1Err(k1/angleBinSize);
bin->setC2Err(0.);
bin->setC3Err(0.);
} else {
Double_t y1=rAng1.getDriftTime1Error(dsBin) * scaleErr;
Double_t y2=rAng2.getDriftTime1Error(dsBin) * scaleErr;
Double_t y3=rAng2.getDriftTime1Error(dsBin+1) * scaleErr;
Double_t y4=rAng1.getDriftTime1Error(dsBin+1) * scaleErr;
Double_t k1=y1-y2;
Double_t k2=y1-y4;
Double_t k3=k1+y3-y4;
k1+=dsBin*k3;
bin->setC0Err(y1 + dsBin*k2 - k0*k1);
bin->setC1Err(k1/angleBinSize);
bin->setC2Err((k0*k3-k2)/distBinSize);
bin->setC3Err(-k3/distBinSize/angleBinSize);
}
}
void HMdcDriftTimeParMod::fillDriftDist(HMdcCal2ParAngle& tAng1,HMdcCal2ParAngle& tAng2,
HMdcDriftTimeParBin* bin,Int_t anBin, Int_t tmBin) {
Double_t k0=90./angleBinSize-anBin;
if(tmBin==lDistBin) {
Double_t y1 = tAng1.getDistance(lDistBin/10,lDistBin%10);
Double_t y2 = tAng2.getDistance(lDistBin/10,lDistBin%10);
Double_t k1 = y1-y2;
bin->setD0(y1-lastBinTime/slopeOutside - k0*k1);
bin->setD1(k1/angleBinSize);
bin->setD2(1./slopeOutside);
bin->setD3(0.);
} else {
Double_t y1 = tAng1.getDistance(tmBin/10,tmBin%10);
Double_t y2 = tAng2.getDistance(tmBin/10,tmBin%10);
Double_t y3 = tAng2.getDistance((tmBin+1)/10,(tmBin+1)%10);
Double_t y4 = tAng1.getDistance((tmBin+1)/10,(tmBin+1)%10);
Double_t k1 = y1-y2;
Double_t k2 = y1-y4;
Double_t k3 = k1+y3-y4;
k1 += tmBin*k3;
bin->setD0(y1 + tmBin*k2 - k0*k1);
bin->setD1(k1/angleBinSize);
bin->setD2((k0*k3-k2)/timeBinSize);
bin->setD3(-k3/timeBinSize/angleBinSize);
}
}
Bool_t HMdcDriftTimeParMod::testSizes(Int_t anSz, Int_t dsSz) const {
return anSz==nAngleBins && dsSz==nDistBins;
}
HMdcDriftTimeParSec::HMdcDriftTimeParSec(void) : TObjArray(4) {
}
Bool_t HMdcDriftTimeParSec::initContainer(HMdcCal2ParSecSim& fC2PSecS,HMdcCal2ParSec& fC2PSec,
Int_t sec, Double_t slOut,Double_t* scaleError) {
sector= sec;
for(Int_t m=0; m<4; m++) if(HMdcGetContainers::getObject()->isModActive(sec,m)) {
HMdcCal2ParModSim& fC2PModS = fC2PSecS[m];
if(&fC2PModS == 0) return kFALSE;
HMdcCal2ParMod& fC2PMod = fC2PSec[m];
if(&fC2PModS == 0) return kFALSE;
Int_t nAngleBins = fC2PModS.getSize();
Int_t nDistBins = 100;
Double_t distBinSize = 0.1;
Double_t timeBinSize = 4.0;
if(At(m)==0) AddAt(new HMdcDriftTimeParMod(nAngleBins,nDistBins),m);
else if(!at(m)->testSizes(nAngleBins,nDistBins)) {
delete at(m);
AddAt(new HMdcDriftTimeParMod(nAngleBins,nDistBins),m);
}
at(m)->setDistBinSize(distBinSize);
at(m)->setTimeBinSize(timeBinSize);
if(!at(m)->initContainer(fC2PModS,fC2PMod,slOut,scaleError[m])) return kFALSE;
}
return kTRUE;
}
void HMdcDriftTimeParSec::print(void) {
for(Int_t m=0; m<4; m++) if(HMdcGetContainers::getObject()->isModActive(sector,m)) {
printf("\n HMdcDriftTimePar: %i sector %i module \n",sector+1,m+1);
printf("---bin---------------------------------------------------------");
printf( "-----------------------------------------------------\n");
printf("ang. dis. | c0 | c1 | c2 | c3 |");
printf( "| c0err | c1err | c2err | c3err |\n");
printf("----------|------------|------------|------------|------------|");
printf( "|------------|------------|------------|------------|\n");
at(m)->print();
}
}
HMdcDriftTimeParSec::~HMdcDriftTimeParSec(void) {
Delete();
}
HMdcDriftTimePar* HMdcDriftTimePar::fMdcDriftTimePar = NULL;
HMdcDriftTimePar::HMdcDriftTimePar(void) : TObjArray(6) {
fMdcDriftTimePar = this;
pCal2ParSim = (HMdcCal2ParSim*)gHades->getRuntimeDb()->getContainer("MdcCal2ParSim");
pCal2Par = (HMdcCal2Par*)gHades->getRuntimeDb()->getContainer("MdcCal2Par");
isInited = kFALSE;
setScalerTime1Err(1.,1.,1.,1.);
}
HMdcDriftTimePar* HMdcDriftTimePar::getObject(void) {
if(fMdcDriftTimePar == NULL) fMdcDriftTimePar = new HMdcDriftTimePar();
return fMdcDriftTimePar;
}
void HMdcDriftTimePar::deleteCont(void) {
if(fMdcDriftTimePar == NULL) return;
delete fMdcDriftTimePar;
fMdcDriftTimePar = NULL;
}
void HMdcDriftTimePar::setScalerTime1Err(Double_t m0,Double_t m1,Double_t m2,Double_t m3) {
scaleError[0] = m0;
scaleError[1] = m1;
scaleError[2] = m2;
scaleError[3] = m3;
}
Bool_t HMdcDriftTimePar::initContainer(void) {
if( isInited ) return kTRUE;
if(pCal2ParSim == NULL || !HMdcGetContainers::isInited(pCal2ParSim)) return kFALSE;
if(pCal2Par == NULL || !HMdcGetContainers::isInited(pCal2Par)) return kFALSE;
Double_t slopeOutside = pCal2ParSim->getSlopeOutside();
for(Int_t s=0; s<6; s++) if( HMdcGetContainers::getObject()->isSecActive(s) ) {
HMdcCal2ParSecSim& fC2PSecS = (*pCal2ParSim)[s];
if(&fC2PSecS == NULL) return kFALSE;
HMdcCal2ParSec& fC2PSec = (*pCal2Par)[s];
if(&fC2PSec == NULL) return kFALSE;
if(At(s)==0) AddAt(new HMdcDriftTimeParSec(),s);
if( !at(s)->initContainer(fC2PSecS,fC2PSec,s,slopeOutside,scaleError) ) return kFALSE;
}
isInited = kTRUE;
return kTRUE;
}
void HMdcDriftTimePar::print(void) {
for(Int_t s=0; s<6; s++) if(HMdcGetContainers::getObject()->isSecActive(s)) at(s)->print();
}
HMdcDriftTimePar::~HMdcDriftTimePar(void) {
Delete();
}