//_HADES_CLASS_DESCRIPTION ///////////////////////////////////////////////////////// // // HShowerDigiPar - parameters of digitalisation // ///////////////////////////////////////////////////////// #include "hshowerdigipar.h" #include "hhistconverter.h" #include "hparamlist.h" #include "hpario.h" #include "hdetpario.h" #include "TMath.h" #include "TRandom.h" #include #include #include #include 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); // 6*3*32*32 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) { // removes all histograms 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; // intitializes the container from an input HDetParIo* input=inp->getDetParIo("HCondParIo"); if (input) rc=input->init(this,set); if (rc && changed) { rc=recreateHistograms(); } return rc; } //------------------------------------------------------------------------------ void HShowerDigiPar::putParams(HParamList* l) { // Puts all params to the parameter list of type HParamList // (which is used by the ASCII and Oracle io); 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; //printf("%d: %f\n", i, fM); 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); // if(iBin >= nGlobalEffBins) iBin = nGlobalEffBins - 1; 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()); }