#include "hshowerdigipar.h"
#include "hhistconverter.h"
#include "hparamlist.h"
#include "hpario.h"
#include "hdetpario.h"
#include "TMath.h"
#include "TRandom.h"
#include <iostream>
#include <iomanip>
#include <stdlib.h>
#include <cmath>
using namespace std;
ClassImp(HShowerDigiPar)
HShowerDigiPar::HShowerDigiPar(const Char_t* name,const Char_t* title,
const Char_t* context)
: HParCond(name,title,context) {
gain.Set(18);
threshold.Set(18);
effScaleMap.Set(18432);
clear();
pfChargeMatrix = NULL;
phEff = NULL;
for(Int_t s = 0; s < 6; s ++) {
for(Int_t m = 0; m < 3; m ++) {
ph2QvB[s][m] = NULL;
pArrayQvB[s][m] = NULL;
nQvBxbins[s][m] = 0;
}
}
nMatrixRange = 0;
chargeMatrixSize = 0;
nGlobalEffBins = 0;
}
HShowerDigiPar::~HShowerDigiPar() {
removeHistograms();
}
void HShowerDigiPar::clear(void) {
gain.Reset(0.f);
threshold.Reset(0.f);
effScaleMap.Reset(1.f);
globalEff.Set(0);
for(Int_t s = 0; s < 6; s ++) {
for(Int_t m = 0; m < 3; m ++) {
chargeVsBeta[s][m].Set(0);
}
}
fThickDet = 1.0f;
nMatrixRange = 0;
fBoxSize = 1.0f;
fChargeSlope = 1.0f;
fPlaneDist = 1.0f;
fPadThreshold = 0.0f;
fUpdatePadThreshold = 0.0f;
status=kFALSE;
resetInputVersions();
changed=kFALSE;
}
void HShowerDigiPar::removeHistograms(void) {
removeChargeHistograms();
removeEfficiencyHistograms();
removeQvBHistograms();
}
void HShowerDigiPar::removeChargeHistograms(void) {
if (pfChargeMatrix != NULL) {
delete [] pfChargeMatrix;
pfChargeMatrix = NULL;
}
chargeMatrixSize = 0;
}
void HShowerDigiPar::removeEfficiencyHistograms(void) {
if (phEff != NULL) {
delete phEff;
phEff = NULL;
}
nGlobalEffBins = 0;
}
void HShowerDigiPar::removeQvBHistograms(void) {
nGlobalEffBins = 0;
for(Int_t s = 0; s < 6; s ++) {
for(Int_t m = 0; m < 3; m ++) {
if (ph2QvB[s][m] != NULL) {
delete ph2QvB[s][m];
ph2QvB[s][m] = NULL;
}
if (pArrayQvB[s][m] != NULL) {
pArrayQvB[s][m]->Delete();
delete pArrayQvB[s][m];
pArrayQvB[s][m] = NULL;
}
nQvBxbins[s][m]= 0;
}
}
}
Bool_t HShowerDigiPar::init(HParIo* inp,Int_t* set) {
Bool_t rc=kFALSE;
HDetParIo* input=inp->getDetParIo("HCondParIo");
if (input) rc=input->init(this,set);
if (rc && changed) {
rc=recreateHistograms();
}
return rc;
}
void HShowerDigiPar::putParams(HParamList* l) {
if (!l) return;
l->add("gain" ,gain);
l->add("threshold" ,threshold);
l->add("fThickDet" ,fThickDet);
l->add("nMatrixRange" ,nMatrixRange);
l->add("fBoxSize" ,fBoxSize);
l->add("fChargeSlope" ,fChargeSlope);
l->add("fPlaneDist" ,fPlaneDist);
l->add("fPadThreshold" ,fPadThreshold);
l->add("fUpdatePadThreshold" ,fUpdatePadThreshold);
l->add("effScaleMap" ,effScaleMap);
l->add("globalEff" ,globalEff);
for(Int_t s = 0; s < 6; s ++) {
for(Int_t m = 0; m < 3; m ++) {
l->add(Form("chargeVsBeta%i%i",s,m),chargeVsBeta[s][m]);
}
}
}
Bool_t HShowerDigiPar::getParams(HParamList* l) {
if (!l) return kFALSE;
if( !(l->fill("gain" ,&gain))) return kFALSE;
if( !(l->fill("threshold" ,&threshold))) return kFALSE;
if (!(l->fill("fThickDet" ,&fThickDet))) return kFALSE;
if (!(l->fill("nMatrixRange" ,&nMatrixRange))) return kFALSE;
if (!(l->fill("fBoxSize" ,&fBoxSize))) return kFALSE;
if (!(l->fill("fChargeSlope" ,&fChargeSlope))) return kFALSE;
if (!(l->fill("fPlaneDist" ,&fPlaneDist))) return kFALSE;
if (!(l->fill("fPadThreshold" ,&fPadThreshold))) return kFALSE;
if (!(l->fill("fUpdatePadThreshold",&fUpdatePadThreshold))) return kFALSE;
if (!(l->fill("effScaleMap" ,&effScaleMap))) return kFALSE;
if (!(l->fill("globalEff" ,&globalEff))) return kFALSE;
for(Int_t s = 0; s < 6; s ++) {
for(Int_t m = 0; m < 3; m ++) {
if (!(l->fill(Form("chargeVsBeta%i%i",s,m),&chargeVsBeta[s][m]))) return kFALSE;
}
}
return kTRUE;
}
Bool_t HShowerDigiPar::recreateHistograms() {
removeHistograms();
Bool_t rc=initChargeMatrix();
if (rc) rc=initEffHistogram();
if (rc) rc=initSumVersBetaHistograms();
return rc;
}
void HShowerDigiPar::setChargeMatrix(Int_t nRange, const Float_t *pMatrix) {
if(pfChargeMatrix != NULL) delete [] pfChargeMatrix;
nMatrixRange = nRange;
if (nMatrixRange == 0) pfChargeMatrix = NULL;
else {
chargeMatrixSize = 2 * nMatrixRange + 1;
pfChargeMatrix = new Float_t[chargeMatrixSize];
if (pMatrix != NULL) {
memcpy(pfChargeMatrix, pMatrix, sizeof(Float_t) * chargeMatrixSize);
} else
memset(pfChargeMatrix, 0, sizeof(Float_t) * chargeMatrixSize);
}
}
Float_t HShowerDigiPar::calcCharge(Float_t fCharge, Float_t fDist,
Float_t fXd, Float_t fYd, Float_t fXu, Float_t fYu) {
return ((fCharge / TMath::TwoPi()) *
(atan(fXd * fYd / (fDist * sqrt(fDist * fDist + fXd * fXd + fYd * fYd)))
- atan(fXd * fYu / (fDist * sqrt(fDist * fDist + fXd * fXd + fYu * fYu)))
+ atan(fXu * fYu / (fDist * sqrt(fDist * fDist + fXu * fXu + fYu * fYu)))
- atan(fXu * fYd / (fDist * sqrt(fDist * fDist + fXu * fXu + fYd * fYd)))
));
}
Bool_t HShowerDigiPar::initChargeMatrix(void) {
Float_t cornerX = -0.5 * fBoxSize;
Float_t cornerY = cornerX;
setChargeMatrix(nMatrixRange);
if (pfChargeMatrix==NULL) {
Error("initChargeMatrix()", "No charge matrix");
return kFALSE;
}
for(Int_t i = 0; i <= nMatrixRange; i++) {
pfChargeMatrix[nMatrixRange - i] = calcCharge(1., fPlaneDist,
cornerX + i * fBoxSize,
cornerY,
cornerX + (i + 1) * fBoxSize,
cornerY + fBoxSize);
Float_t fM = (Float_t)(1.0 - fChargeSlope * i / nMatrixRange);
if(fM < 0.0f) fM = 0.0f;
pfChargeMatrix[nMatrixRange - i] *= fM;
pfChargeMatrix[nMatrixRange + i] = pfChargeMatrix[nMatrixRange - i];
}
for(Int_t i = 1; i < chargeMatrixSize; i++) {
pfChargeMatrix[i] += pfChargeMatrix[i - 1];
}
return kTRUE;
}
Bool_t HShowerDigiPar::initEffHistogram(void) {
phEff = (TH1D*)HHistConverter::createHist("phEff",globalEff);
if(phEff == NULL) {
Error("initEffHistogram()", "No efficiency histogram");
return kFALSE;
}
nGlobalEffBins = phEff->GetNbinsX();
return kTRUE;
}
Float_t HShowerDigiPar::getEfficiency(Int_t sec, Int_t mod, Int_t row, Int_t col,
Float_t fBeta) {
Int_t iBin = phEff->GetXaxis()->FindFixBin((Double_t)fBeta);
return (Float_t)(phEff->GetBinContent(iBin) * effScaleMap[padIndex(sec,mod,row,col)]);
}
Bool_t HShowerDigiPar::checkEfficiency(Int_t sec, Int_t mod, Int_t row,Int_t col,
Float_t fBeta) {
Float_t f = getEfficiency(sec,mod,row,col,fBeta);
return ((f >= 1.0f) || ((f > 0.0f) && (f > gRandom->Rndm())));
}
Bool_t HShowerDigiPar::initSumVersBetaHistograms(void) {
for(Int_t s = 0; s < 6; s ++) {
for(Int_t m = 0; m < 3; m ++) {
TH2D *ph2 = (TH2D*)HHistConverter::createHist(Form("ph2QvB%i%i",s,m),
chargeVsBeta[s][m]);
if(ph2 == NULL) {
Error("%s", "No QvB histogram for sector %i, module %i",s,m);
return kFALSE;
}
Int_t nBins = ph2->GetNbinsX();
TObjArray *pArr = new TObjArray(nBins);
pArr->SetOwner();
for(Int_t iBeta = 0; iBeta < nBins; iBeta++) {
TH1D *ph1 = ph2->ProjectionY(Form("sumVersBeta%i%i_%i",s,m,iBeta),
iBeta, iBeta+1);
ph1->SetDirectory(NULL);
ph1->SetTitle(Form("sum_%03.0f_%03.0f",
ph2->GetXaxis()->GetBinLowEdge(iBeta),
ph2->GetXaxis()->GetBinUpEdge(iBeta)));
pArr->AddAt(ph1, iBeta);
}
ph2QvB[s][m] = ph2;
nQvBxbins[s][m] = nBins;
pArrayQvB[s][m] = pArr;
}
}
return kTRUE;
}
Float_t HShowerDigiPar::getCharge(Int_t sec,Int_t mod, Float_t fBeta) {
TAxis* axis = ph2QvB[sec][mod]->GetXaxis();
Int_t iBin = axis->FindFixBin((Double_t)fBeta);
if(iBin > nQvBxbins[sec][mod]) iBin = nQvBxbins[sec][mod];
if(iBin<0) iBin = 1;
Double_t width = axis->GetBinWidth(iBin);
TH1D *ph1 = (TH1D*)(pArrayQvB[sec][mod]->At(--iBin));
if(ph1 == NULL) return 0.0f;
Double_t fQ = ph1->GetRandom() + (2 * gRandom->Rndm() - 1) * width;
return ((fQ < 0.0) ? 0.0f : (Float_t)fQ);
}
void HShowerDigiPar::setEffScaleMap(TArrayF& arr) {
Int_t n=arr.GetSize();
if (n==18432) {
effScaleMap.Set(n,arr.GetArray());
} else {
Error("setEffScaleMap(TArrayF&)","Array size %i differs from 6*3*32*32=18432",n);
}
}
void HShowerDigiPar::setGlobalEff(TArrayD& arr) {
removeEfficiencyHistograms();
globalEff.Set(arr.GetSize(),arr.GetArray());
}
void HShowerDigiPar::setChargeVsBeta(Int_t sec, Int_t mod, TArrayD& arr) {
removeQvBHistograms();
chargeVsBeta[sec][mod].Set(arr.GetSize(),arr.GetArray());
}