//////////////////////////////////////////////////////////////////////////////
//
// $Id: $
//
//*-- Author  : S. Lebedev
//
//_HADES_CLASS_DESCRIPTION
//////////////////////////////////////////////////////////////////////////////
//
//  HRich700DigiPar
//
//  This class handles parameters for RICH digitizer
//
//////////////////////////////////////////////////////////////////////////////

#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()
{
  // destructor
}

void HRich700DigiPar::clear()
{
    fArrayPmt.Reset();
    fArrayZVertices.Reset();
    for (UInt_t i = 0; i < fArrayThetaMean.size(); i++) {
	fArrayThetaMean[i].Reset();
    }
    fArrayThetaMean.clear();
    fhxyThetaMean.clear();// =0;
    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)
{
    // intitializes the container from an input
  Bool_t rc = HParCond::init(inp,set);
  if (rc && changed){
      // make hists new
      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)
{
    // In PMT geometry there are no sectors any more ... lets take it for phi
    // phi calculated from x,y on PMT plane
    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;   // get sector from phi angle
    if (sectorPhi == -1) sectorPhi = 5;           // 1st sector is along y axis!
    return sectorPhi;
}

void HRich700DigiPar::getLocation(Int_t pmtId, Float_t x, Float_t y, Int_t *loc, Bool_t silent)
{
    // x, y inside pmt (NOT PMT PLANE x,y!) . Used in sim -> HGeantRichPhoton, Direct -> location

    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.;

     // not in the sensitive area
    if (x < -halfPmtSensSize || x > halfPmtSensSize || y < -halfPmtSensSize || y > halfPmtSensSize) {
        return;
    }

    // local pixels address
    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;
    }

    //cout << "x:" << x << " y:"<< y << " indX:" << indX << " indY:" << indY << endl;

    Int_t sectorPhi = getSector(x,y);

    //----------------------------------------------------------------------------------
    // location of the pixel
    loc[0] = sectorPhi;
    loc[1] = indX + pmtData.fIndX * ufNofPixelsInRow; // col
    loc[2] = indY + pmtData.fIndY * ufNofPixelsInRow; // row
}

void HRich700DigiPar::pmtIdPixelToColRowSec(Int_t pmtId,Int_t pixel,Int_t& sec,Int_t& col,Int_t& row, Bool_t silent)
{
    // pmtID
    // pixel (1-64)
    // sec   (0-5)
    // col   (0-192)
    // row   (0-192)

    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;

    //----------------------------------------------------------------------------------
    // pmtID, pixel -> col,row
    UInt_t pix_x = pixel % fNofPixelsInRow;
    UInt_t pix_y = pixel / fNofPixelsInRow;

    col = pix_x + pmtData.fIndX * ufNofPixelsInRow; // col
    row = pix_y + pmtData.fIndY * ufNofPixelsInRow; // row

    //----------------------------------------------------------------------------------
    // pmtID, pixel -> x,y
    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 padding = (fPmtSize - fPmtSensSize) / 2.;
    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;
    // cout << "locIndX:" << locIndX << " xLoc:" << xLoc << " locIndY:"<< locIndY <<  "  yLoc:" << yLoc << endl;
    // cout << "indX:" << indX << " x: " << xLoc + pmtData.fX << " indY:" << indY << " y:" << yLoc + pmtData.fY << endl;

    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)
{

    Int_t nofNoisePixelstmp = gRandom->Gaus(nofNoisePixels, nofNoisePixels/3);

    while (nofNoisePixelstmp < 0 ) {
	nofNoisePixelstmp = gRandom->Gaus(nofNoisePixels, nofNoisePixels/3);
    }

    nofNoisePixels = (UInt_t) nofNoisePixelstmp;
    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);
    //cout <<"x: " << x << " y:" << y << " pmtX:" << pmtIndX << " pmtY:" << pmtIndY << endl;
    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 (x,y) = (%f,%f).",x,y);
    	return -1;
    }
}

Int_t HRich700DigiPar::getPMTId(Int_t col, Int_t row, Bool_t silent)
{

    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 {
        if(!silent)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) ; // do not use Z vertex

    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 {
        // interpolation between two Z positions
	Int_t iZ0 = 0;
        Int_t iZ1 = 0;
	// we assume that Z values are sorted in the array
	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);

    //if Z-Vertex is too far out of range use target center as z-position
    Double_t z = (zv < -120 || zv > 20)?meanTargetWRTMirror + TargetCenter : meanTargetWRTMirror + zv;

    //Parameters for which the polynomials were derived are:
    // r = 872;target0 = 533.3;k_diff = 128;k0 = target0 - 122.2;k1 = k0 - k_diff;
    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;
}

// This is inverse transformation Theta, Phi -> XY Pmt
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;

    //calculate the intersection of planes
    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;

    //calculate the radial position as if the ring would be on the inner plane
    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;
    }

    //calculate the radial position as if the ring would be on the outer plane
    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;
    }

    //Ambiguous position on detector
    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;

    // indX 0 at min X, indX 23 at max X
    // indY 0 at min Y, indY 23 at max Y
    // get neighbour pmts for interpolation
    //
    //          iy+1            2     1
    //     ix+1 ix,iy ix-1         c
    //          iy-1            3     4
    //
    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){  // right upper (1)
        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){ // left upper (2)
        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){ // left lower (3)
        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){ // right lower (4)
        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;
}