#include "hrich700digipar.h"
#include "hparamlist.h"
#include "richdef.h"
#include "hhistconverter.h"
#include "TRandom.h"
#include "TMath.h"
#include <sstream>
#include <iostream>
ClassImp(HRich700DigiPar)
using namespace std;
HRich700DigiPar::HRich700DigiPar(const Char_t* name,const Char_t* title,
const Char_t* context)
: HParCond(name,title,context)
{
clear();
}
HRich700DigiPar::~HRich700DigiPar()
{
}
void HRich700DigiPar::clear()
{
fArrayPmt.Reset();
fArrayZVertices.Reset();
for (UInt_t i = 0; i < fArrayThetaMean.size(); i++) {
fArrayThetaMean[i].Reset();
}
fArrayThetaMean.clear();
fhxyThetaMean.clear();
fArrayThetaTransParamsPoly.Reset();
fArrayThetaTransParamsGeo.Reset();
fArrayPhiAlign.Reset();
fArrayThetaAlign.Reset();
fhPhiAlign = nullptr;
fhThetaAlign = nullptr;
fArrayInvPhiAlign.Reset();
fArrayInvThetaAlign.Reset();
fhInvPhiAlign = nullptr;
fhInvThetaAlign = nullptr;
resetInputVersions();
}
void HRich700DigiPar::printParam(void)
{
printf("----------------------------------------------------------------------\n");
printf("HRich700DigiPar: \n");
printf("fNofPixelsInRow = %i \n" , fNofPixelsInRow);
printf("fPmtSize = %f \n" , fPmtSize);
printf("fPmtGap = %f \n" , fPmtGap);
printf("fPmtSensSize = %f \n" , fPmtSensSize);
printf("fCollectionEfficiency = %f \n" , fCollectionEfficiency);
printf("fCrossTalkProbability = %f \n" , fCrossTalkProbability);
printf("fNofNoiseHits = %i \n" , fNofNoiseHits);
printf("----------------------------------------------------------------------\n");
printf("Sec PmtId IndX IndY x y z pmtType theta phi \n");
for(Int_t i = 0; i < fArrayPmt.GetSize()/NPARPMT; i ++){
Int_t id =getPMTId ( (Float_t)fArrayPmt[i*NPARPMT+3], (Float_t) fArrayPmt[i*NPARPMT+4]);
if(id!= (Int_t)fArrayPmt[i*NPARPMT+0] ) cout <<"missmatch : "<<id<<" "<<(Int_t)fArrayPmt[i*NPARPMT+0]<<endl;
printf("%3i %5i %4i %4i %10.3f %10.3f %10.3f %5i %10.3f %10.3f \n",
getSector((Float_t)fArrayPmt[i*NPARPMT+3],(Float_t)fArrayPmt[i*NPARPMT+4]),
(Int_t)fArrayPmt[i*NPARPMT+0],
(Int_t)fArrayPmt[i*NPARPMT+1],
(Int_t)fArrayPmt[i*NPARPMT+2],
fArrayPmt[i*NPARPMT+3],
fArrayPmt[i*NPARPMT+4],
fArrayPmt[i*NPARPMT+5],
(Int_t)fArrayPmt[i*NPARPMT+6],
fArrayPmt[i*NPARPMT+7],
fArrayPmt[i*NPARPMT+8]
);
}
for(Int_t i = 0; i < fArrayZVertices.GetSize(); i ++){
printf("fArrayZVertices[%i] = %f \n" , i, (Float_t)fArrayZVertices[i]);
}
for(UInt_t i = 0; i < fArrayThetaMean.size(); i ++){
TArrayD arTmp;
TH2F*htmp = (TH2F*)HHistConverter::createHist("htmp" ,fArrayThetaMean[i]);
HHistConverter::fillArray(htmp,arTmp,Form("fArrayThetaMean[%i]",i),10,10,kTRUE);
delete htmp;
}
for(Int_t i = 0; i < fArrayThetaTransParamsPoly.GetSize(); i ++){
printf("fArrayThetaTransParamsPoly[%i] = %f \n" , i, (Float_t)fArrayThetaTransParamsPoly[i]);
}
for(Int_t i = 0; i < fArrayThetaTransParamsGeo.GetSize(); i ++){
printf("fArrayThetaTransParamsGeo[%i] = %f \n" , i, (Float_t)fArrayThetaTransParamsGeo[i]);
}
TArrayD tmpArr;
TH2F* tmpHist = (TH2F*)HHistConverter::createHist("tmpHist", fArrayPhiAlign);
HHistConverter::fillArray(tmpHist,tmpArr, "fArrayPhiAlign", 10, 10, kTRUE);
delete tmpHist;
tmpHist = (TH2F*)HHistConverter::createHist("tmpHist", fArrayThetaAlign);
HHistConverter::fillArray(tmpHist,tmpArr, "fArrayThetaAlign", 10, 10, kTRUE);
delete tmpHist;
tmpHist = (TH2F*)HHistConverter::createHist("tmpHist", fArrayInvPhiAlign);
HHistConverter::fillArray(tmpHist,tmpArr, "fArrayInvPhiAlign", 10, 10, kTRUE);
delete tmpHist;
tmpHist = (TH2F*)HHistConverter::createHist("tmpHist", fArrayInvThetaAlign);
HHistConverter::fillArray(tmpHist,tmpArr, "fArrayInvThetaAlign", 10, 10, kTRUE);
delete tmpHist;
printf("----------------------------------------------------------------------\n");
}
string HRich700DigiPar::getStringForParTxtFile()
{
stringstream ss;
ss << "##############################################################################" << endl;
ss << "# Class: HRich700DigiPar" << endl;
ss << "# Context: Rich700DigiParProduction" << endl;
ss << "##############################################################################" << endl;
ss << "[Rich700DigiPar]" << endl;
ss << "//----------------------------------------------------------------------------" << endl;
ss << "fNofPixelsInRow: Int_t " << fNofPixelsInRow << endl;
ss << "fPmtSize: Double_t " << fPmtSize << endl;
ss << "fPmtGap: Double_t " << fPmtGap << endl;
ss << "fPmtSensSize: Double_t " << fPmtSensSize << endl;
ss << "fCollectionEfficiency: Double_t " << fCollectionEfficiency << endl;
ss << "fCrossTalkProbability: Double_t " << fCrossTalkProbability << endl;
ss << "fNofNoiseHits: Int_t " << fNofNoiseHits << endl;
return ss.str();
}
Bool_t HRich700DigiPar::init(HParIo* inp,Int_t* set)
{
Bool_t rc = HParCond::init(inp,set);
if (rc && changed){
fillMaps();
}
return rc;
}
void HRich700DigiPar::putParams(HParamList* l)
{
if (!l) return;
l->add("fNofPixelsInRow" , fNofPixelsInRow);
l->add("fPmtSize" , fPmtSize);
l->add("fPmtGap" , fPmtGap);
l->add("fPmtSensSize" , fPmtSensSize);
l->add("fCollectionEfficiency", fCollectionEfficiency);
l->add("fCrossTalkProbability", fCrossTalkProbability);
l->add("fNofNoiseHits", fNofNoiseHits);
l->add("fArrayPmt" , fArrayPmt);
l->add("fArrayZVertices" , fArrayZVertices);
for (UInt_t i = 0; i < fArrayThetaMean.size(); i++) {
l->add(Form("fArrayThetaMean[%i]",i) , fArrayThetaMean[i]);
}
l->add("fArrayThetaTransParamsPoly", fArrayThetaTransParamsPoly);
l->add("fArrayThetaTransParamsGeo", fArrayThetaTransParamsGeo);
l->add("fArrayPhiAlign", fArrayPhiAlign);
l->add("fArrayThetaAlign", fArrayThetaAlign);
l->add("fArrayInvPhiAlign", fArrayInvPhiAlign);
l->add("fArrayInvThetaAlign", fArrayInvThetaAlign);
}
Bool_t HRich700DigiPar::getParams(HParamList* l)
{
if (!l) return kFALSE;
if(!( l->fill("fNofPixelsInRow", &fNofPixelsInRow))) return kFALSE;
if(!( l->fill("fPmtSize" , &fPmtSize))) return kFALSE;
if(!( l->fill("fPmtGap" , &fPmtGap))) return kFALSE;
if(!( l->fill("fPmtSensSize" , &fPmtSensSize))) return kFALSE;
if(!( l->fill("fCollectionEfficiency", &fCollectionEfficiency))) return kFALSE;
if(!( l->fill("fCrossTalkProbability", &fCrossTalkProbability))) return kFALSE;
if(!( l->fill("fNofNoiseHits", &fNofNoiseHits))) return kFALSE;
if(!( l->fill("fArrayPmt" , &fArrayPmt))) return kFALSE;
if(!( l->fill("fArrayZVertices" , &fArrayZVertices))) return kFALSE;
fArrayThetaMean.resize(fArrayZVertices.GetSize());
for (Int_t i = 0; i < fArrayZVertices.GetSize(); i++) {
if(!( l->fill(Form("fArrayThetaMean[%i]",i), &fArrayThetaMean[i]))) return kFALSE;
}
if(!( l->fill("fArrayThetaTransParamsPoly", &fArrayThetaTransParamsPoly))) return kFALSE;
if(!( l->fill("fArrayThetaTransParamsGeo", &fArrayThetaTransParamsGeo))) return kFALSE;
if(!( l->fill("fArrayPhiAlign", &fArrayPhiAlign))) return kFALSE;
if(!( l->fill("fArrayThetaAlign", &fArrayThetaAlign))) return kFALSE;
if(!( l->fill("fArrayInvPhiAlign", &fArrayInvPhiAlign))) return kFALSE;
if(!( l->fill("fArrayInvThetaAlign", &fArrayInvThetaAlign))) return kFALSE;
return kTRUE;
}
void HRich700DigiPar::fillMaps(void)
{
for (UInt_t i = 0; i < fhxyThetaMean.size(); i++) {
if(fhxyThetaMean[i]) delete fhxyThetaMean[i];
}
fhxyThetaMean.clear();
for (UInt_t i = 0; i < fArrayThetaMean.size(); i++) {
fhxyThetaMean.push_back((TH2F*)HHistConverter::createHist(Form("fhxyThetaMean[%i]",i),fArrayThetaMean[i]) );
}
if(fhPhiAlign) delete fhPhiAlign;
fhPhiAlign = (TH2D*)HHistConverter::createHist("fhPhiAlign", fArrayPhiAlign);
if(fhThetaAlign) delete fhThetaAlign;
fhThetaAlign = (TH2D*)HHistConverter::createHist("fhThetaAlign", fArrayThetaAlign);
if(fhInvPhiAlign) delete fhInvPhiAlign;
fhInvPhiAlign = (TH2D*)HHistConverter::createHist("fhInvPhiAlign", fArrayInvPhiAlign);
if(fhInvThetaAlign) delete fhInvThetaAlign;
fhInvThetaAlign = (TH2D*)HHistConverter::createHist("fhInvThetaAlign", fArrayInvThetaAlign);
fPmtDataMapPmtId.clear();
fPmtDataMapXY .clear();
fMaxX = -1000;
fMaxY = -1000;
for(Int_t i = 0; i < fArrayPmt.GetSize()/NPARPMT; i++){
Int_t j = i*NPARPMT;
HRich700PmtData pmtData;
pmtData.fPmtId = (Int_t)fArrayPmt[j+0];
pmtData.fIndX = (Int_t)fArrayPmt[j+1];
pmtData.fIndY = (Int_t)fArrayPmt[j+2];
pmtData.fX = fArrayPmt[j+3];
pmtData.fY = fArrayPmt[j+4];
pmtData.fZ = fArrayPmt[j+5];
pmtData.fPmtType = static_cast<HRich700PmtTypeEnum>(fArrayPmt[j+6]);
pmtData.fTheta = fArrayPmt[j+7];
pmtData.fPhi = fArrayPmt[j+8];
if(pmtData.fX > fMaxX) fMaxX = pmtData.fX ;
if(pmtData.fY > fMaxY) fMaxY = pmtData.fY ;
pmtData.fSector = getSector( (Float_t)pmtData.fX, (Float_t)pmtData.fY);
fPmtDataMapPmtId[pmtData.fPmtId] = pmtData;
pair<Int_t,Int_t> xyInd(pmtData.fIndX, pmtData.fIndY);
fPmtDataMapXY[xyInd] = pmtData;
}
}
Int_t HRich700DigiPar::getSector(Float_t x, Float_t y)
{
if(x==0&&y==0) {
Warning("getSector()","x and y equal zero, return sector 0");
return 0;
}
Float_t phi = TMath::ACos(x/sqrt(x*x+y*y));
if (y<0) phi = 2.*TMath::Pi()-phi;
Int_t sectorPhi = (Int_t)(phi/1.0471975)-1;
if (sectorPhi == -1) sectorPhi = 5;
return sectorPhi;
}
void HRich700DigiPar::getLocation(Int_t pmtId, Float_t x, Float_t y, Int_t *loc, Bool_t silent)
{
UInt_t ufNofPixelsInRow = (UInt_t)fNofPixelsInRow;
loc[0]=-1; loc[1]=-1; loc[2]=-1;
map<Int_t, HRich700PmtData>::iterator it = fPmtDataMapPmtId.find(pmtId);
if (it == fPmtDataMapPmtId.end()){
if(!silent) Warning("getLocation","No PMT with id: %i!", pmtId);
return;
}
HRich700PmtData pmtData = it->second;
Double_t dx = fPmtSensSize / (Double_t)fNofPixelsInRow;
Double_t halfPmtSensSize = fPmtSensSize / 2.;
if (x < -halfPmtSensSize || x > halfPmtSensSize || y < -halfPmtSensSize || y > halfPmtSensSize) {
return;
}
UInt_t indX = (UInt_t)((x + halfPmtSensSize) / dx);
UInt_t indY = (UInt_t)((y + halfPmtSensSize) / dx);
if (indX >= ufNofPixelsInRow){
indX = ufNofPixelsInRow - 1;
}
if (indY >= ufNofPixelsInRow){
indY = ufNofPixelsInRow - 1;
}
Int_t sectorPhi = getSector(x,y);
loc[0] = sectorPhi;
loc[1] = indX + pmtData.fIndX * ufNofPixelsInRow;
loc[2] = indY + pmtData.fIndY * ufNofPixelsInRow;
}
void HRich700DigiPar::pmtIdPixelToColRowSec(Int_t pmtId,Int_t pixel,Int_t& sec,Int_t& col,Int_t& row, Bool_t silent)
{
UInt_t ufNofPixelsInRow = (UInt_t)fNofPixelsInRow;
col =-1;
row =-1;
sec =-1;
if(pixel<1){
Warning("pmtIdPixelToColRowSec","pixel < 1 !");
}
map<Int_t, HRich700PmtData>::iterator it = fPmtDataMapPmtId.find(pmtId);
if (it == fPmtDataMapPmtId.end()){
if(!silent) Warning("pmtIdPixelToColRowSec","No PMT with id: %i!", pmtId);
return;
}
pixel-=1;
HRich700PmtData pmtData = it->second;
UInt_t pix_x = pixel % fNofPixelsInRow;
UInt_t pix_y = pixel / fNofPixelsInRow;
col = pix_x + pmtData.fIndX * ufNofPixelsInRow;
row = pix_y + pmtData.fIndY * ufNofPixelsInRow;
Double_t dx = fPmtSensSize / ufNofPixelsInRow;
Double_t halfPmtSensSize = fPmtSensSize / 2.;
Double_t xLoc = (pix_x + 0.5) * dx - halfPmtSensSize;
Double_t yLoc = (pix_y + 0.5) * dx - halfPmtSensSize;
Float_t x = xLoc + pmtData.fX;
Float_t y = yLoc + pmtData.fY;
sec = getSector(x,y);
}
pair<Double_t, Double_t> HRich700DigiPar::getXY(Int_t* loc,Bool_t silent)
{
UInt_t ufNofPixelsInRow = (UInt_t)fNofPixelsInRow;
Int_t indX = loc[1];
Int_t indY = loc[2];
Int_t pmtIndX = (Int_t) indX / ufNofPixelsInRow;
Int_t pmtIndY = (Int_t) indY / ufNofPixelsInRow;
map<pair<Int_t, Int_t>, HRich700PmtData>::iterator it = fPmtDataMapXY.find(make_pair(pmtIndX, pmtIndY));
if (it == fPmtDataMapXY.end()){
if(!silent)Warning("getXY()","No PMT with XY: %i %i!",pmtIndX,pmtIndY);
return pair<Double_t, Double_t>(0., 0.);
}
HRich700PmtData pmtData = it->second;
Int_t locIndX = indX % ufNofPixelsInRow;
Int_t locIndY = indY % ufNofPixelsInRow;
Double_t dx = fPmtSensSize / ufNofPixelsInRow;
Double_t halfPmtSensSize = fPmtSensSize / 2.;
Double_t xLoc = (locIndX + 0.5) * dx - halfPmtSensSize;
Double_t yLoc = (locIndY + 0.5) * dx - halfPmtSensSize;
pair<Double_t, Double_t> xy(xLoc + pmtData.fX, yLoc + pmtData.fY);
return xy;
}
pair<Double_t, Double_t> HRich700DigiPar::getPmtCenter(Int_t pmtId)
{
return pair<Double_t, Double_t>(fPmtDataMapPmtId[pmtId].fX, fPmtDataMapPmtId[pmtId].fY);
}
vector<pair<Double_t, Double_t> > HRich700DigiPar::getPmtCenters()
{
vector<pair<Double_t, Double_t> > v;
v.reserve(fPmtDataMapPmtId.size());
for(map<Int_t,HRich700PmtData>::iterator it = fPmtDataMapPmtId.begin(); it != fPmtDataMapPmtId.end(); it++) {
v.push_back(pair<Double_t,Double_t>(it->second.fX, it->second.fY));
}
return v;
}
vector<pair<Int_t, Int_t> > HRich700DigiPar::getDirectNeighbourPixels(Int_t col, Int_t row)
{
UInt_t ufNofPixelsInRow = (UInt_t)fNofPixelsInRow;
std::vector<pair<Int_t, Int_t> > v;
UInt_t locCol = (UInt_t)(col % ufNofPixelsInRow);
UInt_t locRow = (UInt_t)(row % ufNofPixelsInRow);
if (locCol < ufNofPixelsInRow - 1) {
v.push_back(make_pair(col + 1, row));
}
if (locCol > 0) {
v.push_back(make_pair(col - 1, row));
}
if (locRow < ufNofPixelsInRow - 1) {
v.push_back(make_pair(col, row + 1));
}
if (locRow > 0) {
v.push_back(make_pair(col, row - 1));
}
return v;
}
vector<pair<Int_t, Int_t> > HRich700DigiPar::getDiagonalNeighbourPixels(Int_t col, Int_t row)
{
UInt_t ufNofPixelsInRow = (UInt_t)fNofPixelsInRow;
std::vector<pair<Int_t, Int_t> > v;
UInt_t locCol = (UInt_t)(col % ufNofPixelsInRow);
UInt_t locRow = (UInt_t)(row % ufNofPixelsInRow);
if (locCol < ufNofPixelsInRow - 1 && locRow < ufNofPixelsInRow - 1) {
v.push_back(make_pair(col + 1, row + 1));
}
if (locCol > 0 && locRow > 0) {
v.push_back(make_pair(col - 1, row - 1));
}
if (locCol < ufNofPixelsInRow - 1 && locRow > 0) {
v.push_back(make_pair(col + 1, row - 1));
}
if (locCol > 0 && locRow < ufNofPixelsInRow - 1) {
v.push_back(make_pair(col - 1, row + 1));
}
return v;
}
vector<pair<Int_t, Int_t> > HRich700DigiPar::getNoisePixels(UInt_t nofNoisePixels)
{
UInt_t ufNofPixelsInRow = (UInt_t)fNofPixelsInRow;
std::vector<pair<Int_t, Int_t> > v;
while(kTRUE) {
Int_t col = (Int_t)(gRandom->Rndm() * RICH700_MAX_COLS);
Int_t row = (Int_t)(gRandom->Rndm() * RICH700_MAX_ROWS);
Int_t pmtIndX = (Int_t) col / ufNofPixelsInRow;
Int_t pmtIndY = (Int_t) row / ufNofPixelsInRow;
map<pair<Int_t, Int_t>, HRich700PmtData>::iterator it = fPmtDataMapXY.find(make_pair(pmtIndX, pmtIndY));
if (it != fPmtDataMapXY.end()){
v.push_back(make_pair(col, row));
if (v.size() >= nofNoisePixels) break;
}
}
return v;
}
Int_t HRich700DigiPar::getPMTId(Float_t x, Float_t y)
{
if(x>fMaxX) x = fMaxX;
if(y>fMaxY) y = fMaxY;
if(x<-fMaxX) x = -fMaxX;
if(y<-fMaxY) y = -fMaxY;
Int_t pmtIndX = ( (fMaxX+fPmtSize/2.) + x)/(fPmtSize+fPmtGap);
Int_t pmtIndY = ( (fMaxY+fPmtSize/2.) + y)/(fPmtSize+fPmtGap);
map<pair<Int_t, Int_t>, HRich700PmtData>::iterator it = fPmtDataMapXY.find(make_pair(pmtIndX, pmtIndY));
if (it != fPmtDataMapXY.end()){
return it->second.fPmtId;
} else {
return -1;
}
}
Int_t HRich700DigiPar::getPMTId(Int_t col, Int_t row)
{
Int_t pmtIndX = col/fNofPixelsInRow;
Int_t pmtIndY = row/fNofPixelsInRow;
map<pair<Int_t, Int_t>, HRich700PmtData>::iterator it = fPmtDataMapXY.find(make_pair(pmtIndX, pmtIndY));
if (it != fPmtDataMapXY.end()){
return it->second.fPmtId;
} else {
Warning("getPMTId()","No PMT found for (col,row) = (%i,%i).",col,row);
return -1;
}
}
HRich700PmtData* HRich700DigiPar::getPMTData(Int_t pmtid)
{
map<Int_t,HRich700PmtData>::iterator it = fPmtDataMapPmtId.find(pmtid);
if(it == fPmtDataMapPmtId.end()) return NULL;
return &fPmtDataMapPmtId[pmtid];
}
Int_t HRich700DigiPar::getInterpolatedSectorThetaPhi(Float_t x, Float_t y, Float_t& theta,Float_t& phi)
{
theta = 0;
phi = 0;
Int_t sector = getSector(x,y);
phi = atan2(y, x)*TMath::RadToDeg();
phi = phi>0 ? phi : 360+phi;
theta = fhxyThetaMean[0] ->Interpolate(x,y) ;
return sector;
}
Int_t HRich700DigiPar::getInterpolatedSectorThetaPhi(Float_t x, Float_t y, Float_t zv, Float_t& theta,Float_t& phi)
{
theta = 0;
phi = 0;
Int_t sector = getSector(x,y);
phi = atan2(y, x)*TMath::RadToDeg();
phi = phi>0 ? phi : 360+phi;
if (zv >= fArrayZVertices[0]) {
theta = fhxyThetaMean[0] ->Interpolate(x,y);
} else if (zv <= fArrayZVertices[fArrayZVertices.GetSize() - 1]) {
theta = fhxyThetaMean[fhxyThetaMean.size() - 1] ->Interpolate(x,y);
} else {
Int_t iZ0 = 0;
Int_t iZ1 = 0;
for (Int_t i = 0; fArrayZVertices.GetSize(); i++){
if (zv > fArrayZVertices[i]){
iZ1 = i;
iZ0 = iZ1 - 1;
break;
}
}
Double_t z0 = fArrayZVertices[iZ0];
Double_t z1 = fArrayZVertices[iZ1];
Double_t vTheta0 = fhxyThetaMean[iZ0]->Interpolate(x, y);
Double_t vTheta1 = fhxyThetaMean[iZ1]->Interpolate(x, y);
theta = vTheta0 + (zv - z0) * (vTheta1 - vTheta0) / (z1 - z0);
}
return sector;
}
Int_t HRich700DigiPar::getInterpolatedSectorThetaPhiAnalytical(Float_t x, Float_t y, Float_t zv, Float_t& theta,Float_t& phi)
{
theta = 0;
phi = 0;
Int_t sector = getSector(x,y);
phi = atan2(y, x)*TMath::RadToDeg();
phi = phi>0 ? phi : 360+phi;
Double_t innerPlaneWidth = fArrayThetaTransParamsGeo[0];
Double_t meanTargetWRTMirror = fArrayThetaTransParamsGeo[1];
Double_t TargetCenter = fArrayThetaTransParamsGeo[2];
Double_t P[16][5];
Int_t i = 0;
for (Int_t i1 = 0; i1 < 16; i1++) {
for (Int_t i2 = 0; i2 < 5; i2++) {
P[i1][i2] = fArrayThetaTransParamsPoly[i];
i++;
}
}
Bool_t isInnerPlane = (std::abs(x)<innerPlaneWidth && std::abs(y)<innerPlaneWidth);
Double_t z = (zv < -120 || zv > 20)?meanTargetWRTMirror + TargetCenter : meanTargetWRTMirror + zv;
Double_t z2 = z * z;
Double_t z3 = z2 * z;
Double_t z4 = z3 * z;
Double_t P0I = P[0][0] + P[0][1] * z + P[0][2] * z2 + P[0][3] * z3 + P[0][4] * z4;
Double_t P1I = P[1][0] + P[1][1] * z + P[1][2] * z2 + P[1][3] * z3 + P[1][4] * z4;
Double_t P2I = P[2][0] + P[2][1] * z + P[2][2] * z2 + P[2][3] * z3 + P[2][4] * z4;
Double_t P3I = P[3][0] + P[3][1] * z + P[3][2] * z2 + P[3][3] * z3 + P[3][4] * z4;
Double_t P4I = P[4][0] + P[4][1] * z + P[4][2] * z2 + P[4][3] * z3 + P[4][4] * z4;
Double_t P5I = P[5][0] + P[5][1] * z + P[5][2] * z2 + P[5][3] * z3 + P[5][4] * z4;
Double_t P6I = P[6][0] + P[6][1] * z + P[6][2] * z2 + P[6][3] * z3 + P[6][4] * z4;
Double_t P7I = P[7][0] + P[7][1] * z + P[7][2] * z2 + P[7][3] * z3 + P[7][4] * z4;
Double_t P0O = P[8][0] + P[8][1] * z + P[8][2] * z2 + P[8][3] * z3 + P[8][4] * z4;
Double_t P1O = P[9][0] + P[9][1] * z + P[9][2] * z2 + P[9][3] * z3 + P[9][4] * z4;
Double_t P2O = P[10][0] + P[10][1] * z + P[10][2] * z2 + P[10][3] * z3 + P[10][4] * z4;
Double_t P3O = P[11][0] + P[11][1] * z + P[11][2] * z2 + P[11][3] * z3 + P[11][4] * z4;
Double_t P4O = P[12][0] + P[12][1] * z + P[12][2] * z2 + P[12][3] * z3 + P[12][4] * z4;
Double_t P5O = P[13][0] + P[13][1] * z + P[13][2] * z2 + P[13][3] * z3 + P[13][4] * z4;
Double_t P6O = P[14][0] + P[14][1] * z + P[14][2] * z2 + P[14][3] * z3 + P[14][4] * z4;
Double_t P7O = P[15][0] + P[15][1] * z + P[15][2] * z2 + P[15][3] * z3 + P[15][4] * z4;
Double_t r = sqrt(x*x + y*y);
Double_t r2 = r * r;
Double_t r3 = r2 * r;
Double_t r4 = r3 * r;
Double_t r5 = r4 * r;
Double_t r6 = r5 * r;
Double_t r7 = r6 * r;
if(isInnerPlane ) {
theta = P0I + P1I * r + P2I * r2 + P3I * r3 + P4I * r4 + P5I * r5 + P6I * r6 + P7I * r7;
} else {
theta = P0O + P1O * r + P2O * r2 + P3O * r3 + P4O * r4 + P5O * r5 + P6O * r6 + P7O * r7;
}
return sector;
}
void HRich700DigiPar::getRingCenterXY(Float_t theta, Float_t phi, Float_t zv, Float_t& x,Float_t& y)
{
Double_t innerPlaneWidth = fArrayThetaTransParamsGeo[0];
Double_t meanTargetWRTMirror = fArrayThetaTransParamsGeo[1];
Double_t TargetCenter = fArrayThetaTransParamsGeo[2];
Double_t mirrorRadius = fArrayThetaTransParamsGeo[3];
Double_t PMTPlane1Z = fArrayThetaTransParamsGeo[4];
Double_t PMTPlane2Z = fArrayThetaTransParamsGeo[5];
Double_t z = (zv < -120 || zv > 20)?meanTargetWRTMirror + TargetCenter : meanTargetWRTMirror + zv;
Double_t phiModNeg45 = fmod(phi, 90);
if(phiModNeg45 > 45) {
phiModNeg45 = fmod(phi, 90) - 90;
}
Double_t innerPlaneWidthAtPhi = innerPlaneWidth / cos(phiModNeg45 * TMath::DegToRad() );
Double_t phiRad = phi*TMath::DegToRad();
Double_t ct = cos(theta*TMath::DegToRad() );
Double_t st = sin(theta*TMath::DegToRad() );
Double_t zz = z*z;
Double_t rr = mirrorRadius*mirrorRadius;
Double_t xInner = (2*PMTPlane1Z*zz*st*st*st + 2*sqrt(-zz*st*st + rr)*PMTPlane1Z*z*ct*st - PMTPlane1Z*rr*st - rr*z*st)
/(2*zz*ct*st*st - 2*sqrt(-zz*st*st + rr)*z*st*st - rr*ct);
if (xInner < innerPlaneWidthAtPhi){
x = xInner*cos(phiRad);
y = xInner*sin(phiRad);
return;
}
double xOuter = (2*PMTPlane2Z*zz*st*st*st + 2*sqrt(-zz*st*st + rr)*PMTPlane2Z*z*ct*st - PMTPlane2Z*rr*st - rr*z*st)
/(2*zz*ct*st*st - 2*sqrt(-zz*st*st + rr)*z*st*st - rr*ct);
if (xOuter > innerPlaneWidthAtPhi){
x = xOuter*cos(phiRad);
y = xOuter*sin(phiRad);
return;
}
x = innerPlaneWidthAtPhi*cos(phiRad);
y = innerPlaneWidthAtPhi*sin(phiRad);
}
void HRich700DigiPar::getAlignedThetaPhi(const Float_t theta, const Float_t phi
, Float_t& thetaCor, Float_t& phiCor)
{
thetaCor = theta;
phiCor = phi;
if(theta > 90 || theta < 10) {
return;
}
else {
Float_t mod_phi = fmod(phi,360);
phiCor -= fhPhiAlign->Interpolate(theta, mod_phi);
thetaCor -= fhThetaAlign->Interpolate(theta, mod_phi);
return;
}
}
void HRich700DigiPar::getAlignedThetaPhiInv(const Float_t theta, const Float_t phi
, Float_t & thetaCor, Float_t & phiCor)
{
thetaCor = theta;
phiCor = phi;
if(theta > 90 || theta < 10) {
return;
}
else {
Float_t mod_phi = fmod(phi,360);
phiCor -= fhInvPhiAlign->Interpolate(theta, mod_phi);
thetaCor -= fhInvThetaAlign->Interpolate(theta, mod_phi);
return;
}
}
Bool_t HRich700DigiPar::getInterpolatedThetaPhiPMT(Float_t x, Float_t y, Float_t& theta,Float_t& phi)
{
Int_t pmtid = getPMTId(x,y);
if(pmtid < 0) return kFALSE;
HRich700PmtData& data = fPmtDataMapPmtId[pmtid];
Int_t iX = data.fIndX;
Int_t iY = data.fIndY;
Double_t x0 = data.fX;
Double_t y0 = data.fY;
Double_t phi0 = data.fPhi;
Double_t theta0 = data.fTheta;
theta = theta0;
phi = phi0;
map<pair<Int_t, Int_t>, HRich700PmtData>::iterator itX, itY, itXY;
Double_t x1,x2,y1,y2;
Double_t p11,p12,p21,p22;
Double_t t11,t12,t21,t22;
if(x <= x0 && y >= y0){
itX = fPmtDataMapXY.find(make_pair(iX-1, iY));
itY = fPmtDataMapXY.find(make_pair(iX, iY+1));
itXY = fPmtDataMapXY.find(make_pair(iX-1, iY+1));
} else if(x >= x0 && y >= y0){
itX = fPmtDataMapXY.find(make_pair(iX+1, iY));
itY = fPmtDataMapXY.find(make_pair(iX, iY+1));
itXY = fPmtDataMapXY.find(make_pair(iX+1, iY+1));
} else if(x >= x0 && y <= y0){
itX = fPmtDataMapXY.find(make_pair(iX+1, iY));
itY = fPmtDataMapXY.find(make_pair(iX, iY-1));
itXY = fPmtDataMapXY.find(make_pair(iX+1, iY-1));
} else if(x <= x0 && y <= y0){
itX = fPmtDataMapXY.find(make_pair(iX-1, iY));
itY = fPmtDataMapXY.find(make_pair(iX, iY-1));
itXY = fPmtDataMapXY.find(make_pair(iX-1, iY-1));
}
x1 = x0;
y1 = y0;
p11 = phi0;
t11 = theta0;
if (itX != fPmtDataMapXY.end()){
x2 = itX->second.fX;
p21 = itX->second.fPhi;
t21 = itX->second.fTheta;
} else {
x2 = x0;
p21 = phi0;
t21 = theta0;
}
if (itY != fPmtDataMapXY.end()){
y2 = itY->second.fY;
p12 = itY->second.fPhi;
t12 = itY->second.fTheta;
} else {
y2 = y0;
p12 = phi0;
t12 = theta0;
}
if (itXY != fPmtDataMapXY.end()){
p22 = itXY->second.fPhi;
t22 = itXY->second.fTheta;
} else {
p22 = phi0;
t22 = theta0;
}
Double_t d = 1.0*(x2-x1)*(y2-y1);
if (d == 0) {
theta = theta0;
phi = phi0;
return kTRUE;
}
theta = (Float_t)(1.0*t11/d*(x2-x)*(y2-y)+1.0*t21/d*(x-x1)*(y2-y)+1.0*t12/d*(x2-x)*(y-y1)+1.0*t22/d*(x-x1)*(y-y1));
phi = (Float_t)(1.0*p11/d*(x2-x)*(y2-y)+1.0*p21/d*(x-x1)*(y2-y)+1.0*p12/d*(x2-x)*(y-y1)+1.0*p22/d*(x-x1)*(y-y1));
return kTRUE;
}
Int_t HRich700DigiPar::getSectorPixels(Int_t col,Int_t row)
{
Int_t pmtIndX = (Int_t) col / fNofPixelsInRow;
Int_t pmtIndY = (Int_t) row / fNofPixelsInRow;
Int_t sectorPhi=0;
map<pair<Int_t, Int_t>, HRich700PmtData>::iterator it = fPmtDataMapXY.find(make_pair(pmtIndX, pmtIndY));
if (it != fPmtDataMapXY.end()){
Float_t x = (Float_t)it->second.fX;
Float_t y = (Float_t)it->second.fY;
sectorPhi = getSector(x, y);
}
return sectorPhi;
}
Int_t HRich700DigiPar::getSectorPMTInd(Int_t xind,Int_t yind)
{
Int_t sectorPhi=0;
map<pair<Int_t, Int_t>, HRich700PmtData>::iterator it = fPmtDataMapXY.find(make_pair(xind,yind));
if (it != fPmtDataMapXY.end()){
sectorPhi = it->second.fSector;
}
return sectorPhi;
}
Int_t HRich700DigiPar::getSectorPMTId(Int_t pmtid)
{
Int_t sectorPhi=0;
map<Int_t,HRich700PmtData> ::iterator it = fPmtDataMapPmtId.find(pmtid);
if (it != fPmtDataMapPmtId.end()){
sectorPhi = it->second.fSector;
}
return sectorPhi;
}
Int_t HRich700DigiPar::getSectorPhiThetaDegPixels(Int_t col,Int_t row,Float_t& phiDeg,Float_t& thetaDeg)
{
Int_t pmtIndX = (Int_t) col / fNofPixelsInRow;
Int_t pmtIndY = (Int_t) row / fNofPixelsInRow;
Int_t sectorPhi = 0;
phiDeg = 0;
thetaDeg = 0;
map<pair<Int_t, Int_t>, HRich700PmtData>::iterator it = fPmtDataMapXY.find(make_pair(pmtIndX, pmtIndY));
if (it != fPmtDataMapXY.end()){
Float_t x = (Float_t)it->second.fX;
Float_t y = (Float_t)it->second.fY;
sectorPhi = getInterpolatedSectorThetaPhi(x,y,thetaDeg,phiDeg);
}
return sectorPhi;
}
Int_t HRich700DigiPar::getSectorPhiThetaDegPMTInd(Int_t xind,Int_t yind,Float_t& phiDeg,Float_t& thetaDeg)
{
Int_t sectorPhi = 0;
phiDeg = 0;
thetaDeg = 0;
map<pair<Int_t, Int_t>, HRich700PmtData>::iterator it = fPmtDataMapXY.find(make_pair(xind,yind));
if (it != fPmtDataMapXY.end()){
sectorPhi = it->second.fSector;
phiDeg = (Float_t)it->second.fPhi;
thetaDeg = (Float_t)it->second.fTheta;
}
return sectorPhi;
}
Int_t HRich700DigiPar::getSectorPhiThetaDegPMTId(Int_t pmtid,Float_t& phiDeg,Float_t& thetaDeg)
{
Int_t sectorPhi = 0;
phiDeg = 0;
thetaDeg = 0;
map<Int_t,HRich700PmtData>::iterator it = fPmtDataMapPmtId.find(pmtid);
if (it != fPmtDataMapPmtId.end()){
sectorPhi = it->second.fSector;
phiDeg = (Float_t)it->second.fPhi;
thetaDeg = (Float_t)it->second.fTheta;
}
return sectorPhi;
}