// $Id: PmtGeometry.cc,v 1.17 2004/12/23 00:43:33 thea Exp $
// Author: D.Demarco, M.Pallavicini
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: PmtGeometry *
* Package: Optics *
* Coordinator: Marco.Pallavicini *
* *
*****************************************************************************/
//______________________________________________________________________________
//
// Photomultiplier Geometry
// ========================
//
// Geometry descriptor of a PMT this object can be attached to any
// Photomultiplier object In the constructor it retrieves the relevant
// parameters from the Config object and store them in class variables Local
// Coordinates: x,y,z z is normal to the Photocathode (0,0,0) is the upper
// left corner of the pmt photocathode surface up is defined by the y axis
// right is defined by the x axis direction left is -x
//
#include <iostream>
#include "PmtGeometry.hh"
#include "Photomultiplier.hh"
#include "OpticalAdaptor.hh"
#include "EusoMapping.hh"
ClassImp(PmtGeometry)
#define MAX_PIXELS 500000
// class variables with geometry parameters
Double_t PmtGeometry::fgSide = -1.; // physical total width of the pmt
Double_t PmtGeometry::fgDeadBorder = 0.; // thickness of external dead space
Double_t PmtGeometry::fgDeadInner = 0.; // thickness of internal space
// ( between pads )
Double_t PmtGeometry::fgPadSide = 0.; // width of one pad
Int_t PmtGeometry::fgNumPads = 0; // number of active pads of this mapmt
Int_t PmtGeometry::fgNumRows = 0; // number of rows and columns
TVector3* PmtGeometry::fgCorners[10][10];
//______________________________________________________________________________
PmtGeometry::PmtGeometry(TVector3 x, TVector3 n, TVector3 m)
: EsafConfigurable(), fPosition(x), fZAxis(n.Unit()) {
//
// Construtor
//
fXAxis = ( m.Cross(fZAxis)).Unit();
fYAxis = fZAxis.Cross(fXAxis);
Double_t dx, dy, dz;
if ( (dz=(fXAxis.Dot(fYAxis))) > PMT_TOLERANCE ||
(dx=(fYAxis.Dot(fZAxis))) > PMT_TOLERANCE ||
(dy=(fZAxis.Dot(fXAxis))) > PMT_TOLERANCE )
throw runtime_error(Form("SOmething's wrong in pmt's axis dx = %f, dy = %f, dz = %f",dx, dy, dz));
fRight = NULL;
fLeft = NULL;
fFront = NULL;
fBack = NULL;
pPmt = NULL;
pOA = NULL;
if ( NumPads() <= 0 )
SetPmtGeometry();
}
//______________________________________________________________________________
PmtGeometry::~PmtGeometry() {
//
// Destructor
//
delete pOA;
}
//______________________________________________________________________________
void PmtGeometry::ResetClass() {
//
// Reset static data members
//
fgSide = -1.;
fgDeadBorder = 0.;
fgDeadInner = 0.;
fgPadSide = 0.;
fgNumPads = 0;
fgNumRows = 0;
}
//______________________________________________________________________________
void PmtGeometry::SetNearest(PmtGeometry* r, PmtGeometry* l,
PmtGeometry* f, PmtGeometry* b) {
//
// Set pointers to nearest neighbors
//
fRight = r;
fLeft = l;
fFront = f;
fBack = b;
#ifdef DEBUG
// cout << "Set nearest called for PMT = " << Pmt()->Id() << endl;
// if ( b )
// cout << " Back = " << b->Pmt()->Id() << endl;
// if ( f )
// cout << " Front = " << f->Pmt()->Id() << endl;
// if ( r )
// cout << " Right = " << r->Pmt()->Id() << endl;
// if ( l )
// cout << " Left = " << l->Pmt()->Id() << endl;
#endif
}
//______________________________________________________________________________
void PmtGeometry::InsertOA( OpticalAdaptor* oa ){
//
// Insert oa in front of the pmt and displace geometry of oa->GetHeight()
//
// push the pmt back
fPosition -= GetZ()*oa->GetThickness();
oa->SetGeometry(this);
SetOA(oa);
}
//______________________________________________________________________________
void PmtGeometry::SetPmtGeometry() {
//
// Set local geometry infos and global parameters only once
//
if ( fgNumPads == 0 ) {
ConfigFileParser* conf = Config::Get()->GetCF("Electronics","Photomultiplier");
fgNumRows = (Int_t)conf->GetNum("Photomultiplier.PmtSide");
fgNumPads = fgNumRows*fgNumRows;
fgSide = conf->GetNum("Photomultiplier.PmtSize");
fgDeadBorder = conf->GetNum("Photomultiplier.PmtDeadLateral");
fgDeadInner = conf->GetNum("Photomultiplier.PmtDeadInner");
fgPadSide = ( fgSide - 2*fgDeadBorder - (fgNumRows-1)*fgDeadInner ) / fgNumRows;
// compute the corner position of each pad in local coordinate
for(Int_t row=0; row<=fgNumRows; row++) {
Double_t row_shift = fgDeadBorder + row * fgPadSide + row * fgDeadInner;
for(Int_t col=0; col<=fgNumRows; col++) {
Double_t col_shift = fgDeadBorder + col * fgPadSide + col * fgDeadInner;
fgCorners[row][col] = new TVector3( col_shift, row_shift, 0. );
}
}
}
}
//______________________________________________________________________________
Bool_t PmtGeometry::IsInside( const Photon& ph ) const {
// check whether a Photon is inside the geometrical acceptance of this pmt
// it returns true when the Photon is within the external side and its
// distance from the surface is less than PMT_TOLERANCE
// check if the photon is on the cathode surface
TVector3 a = ( ph.pos - Position());
Double_t b = a.Dot(Normal());
if ( b < 0. ) b = -b;
if ( b > PMT_TOLERANCE ) {
cerr<<"PmtGeometry::IsInside: photon not on PMT surface, " << b << " mm away" <<endl;
return false;
}
// check if the photon direction is going to intersect the cathode
b = Normal().Dot( ph.dir );
if ( b > -PMT_TOLERANCE ) {
cerr<<"PmtGeometry::IsInside: photon not intersecting PMT surface"<<endl;
return false;
}
// check if it is inside
Double_t dx = a.Dot(GetX().Unit()); // along x is positive
if ( dx < PMT_TOLERANCE || dx > (fgSide-PMT_TOLERANCE) )
return false;
Double_t dy = -a.Dot(GetY().Unit()); // along y must be negative
if ( dy < PMT_TOLERANCE || dy > (fgSide-PMT_TOLERANCE) )
return false;
return true;
}
//______________________________________________________________________________
Double_t PmtGeometry::IsHit( const Photon& ph ) const {
//
// check whether a Photon
//
// ph is going away from cathode
if ( ph.dir.Dot(Normal()) > -PMT_TOLERANCE ) return -1;
// calculate the hit point between ph and cathode plane
Photon dummy;
dummy.pos=ph.pos-ph.dir*((ph.pos-Position()).Dot(Normal())/ph.dir.Dot(Normal()));
dummy.dir=ph.dir;
if ( IsInside(dummy) ){
cout << "PmtGeometry::IsHitn"
<< "ph.pos = " << ph.pos << " ph.dir = " << ph.dir << endl;
cout << "Position() = " << Position() << " Normal() = " << Normal() << endl;
cout << "DL = " << (dummy.pos-ph.pos).Mag() << endl;
return (dummy.pos-ph.pos).Mag();
}
return -1;
}
//______________________________________________________________________________
TVector3 PmtGeometry::Center() const {
//
// Return the position of the center of the pmt
//
TVector3 a = Position() + (fgSide/2.) * GetX().Unit();
a += (-fgSide/2.) * GetY().Unit();
return a;
}
// return local coordinate position of a point in global coord
// the transofrmation works for point on the pmt surface
TVector3 PmtGeometry::LocalCoord( const TVector3& pos) const{
TVector3 diff = pos - Position();
Double_t deltax = diff.Dot(GetX().Unit());
Double_t deltay = diff.Dot(GetY().Unit());
return TVector3( deltax, deltay, 0.0);
}
// return the number of a Pad hit by a Photon
Int_t PmtGeometry::Pad( const Photon& ph ) const{
return Pad( LocalCoord( ph.pos ) );
}
// return the pad number corresponding to local coordinate x
Int_t PmtGeometry::Pad( const TVector3& x) const{
Int_t row,col;
for(row=0; row < Rows(); row++ ) {
Double_t yy = -x.y();// y axis is negative!
Double_t yy_low = Corner(row,0).y();
Double_t yy_high = Corner(row+1,0).y();
if ( yy > yy_low && yy < yy_high)
break;
}
if ( row == Rows())
return -1;
for(col=0; col < Rows(); col++ ) {
Double_t xx = x.x();
Double_t xx_low = Corner(0,col).x();
Double_t xx_high = Corner(0,col+1).x();
if ( xx > xx_low && xx < xx_high)
break;
}
if ( col == Rows())
return -1;
return ( col + ( row*Rows() ) );
}
//______________________________________________________________________________
const PmtGeometry* PmtGeometry::Nearest( OAFace side ) const {
//
// Nearest neighbor
//
switch ( side ) {
case RIGHT:
return fRight; // x > 0
break;
case LEFT:
return fLeft; // x < 0
break;
case FRONT:
return fFront; // y < 0
break;
case BACK:
return fBack; // y > 0
break;
case TOP:
case BOTTOM:
default:
break;
}
cerr << "Invalid side type in PmtGeometry::Nearest" << endl;
return NULL;
}
//______________________________________________________________________________
TVector3 PmtGeometry::Position(Int_t r, Int_t c) const {
//
// Return the position in EUSO coordinate of the center of a pad
// r and c are not checked
//
TVector3 p = Position();
p += (fgDeadBorder + fgPadSide/2 + c*(fgPadSide+fgDeadInner)) * GetX();
p -= (fgDeadBorder + fgPadSide/2 + r*(fgPadSide+fgDeadInner)) * GetY();
return p;
}
//______________________________________________________________________________
TVector3 PmtGeometry::Position(Int_t ch) const {
//
// return the position in EUSO coordinate of the center of a pad
// ch is the channel number defined as ch = r*n+c
//
Int_t r = ch / Rows();
Int_t c = ch - r*Rows();
return Position(r,c);
}
//______________________________________________________________________________
Double_t PmtGeometry::ThetaFOV(Int_t ch) const {
//
// Theta in field of view (FOV) corresponding to the center of this pad
//
return EusoMapping::Get()->GetThetaFOV( GetUniqueId(ch) );
}
//______________________________________________________________________________
Double_t PmtGeometry::PhiFOV(Int_t ch) const {
//
// Theta in field of view (FOV) corresponding to the center of this pad
//
return EusoMapping::Get()->GetPhiFOV( GetUniqueId(ch) );
}
//______________________________________________________________________________
Int_t PmtGeometry::GetChannel( ChannelUniqueId chid ) const {
//
// Returns channel associated to unique id
//
if ( chid < GetStartUniqueId() ) {
return -1;
}
if ( chid > GetLastUniqueId() ) {
return -1;
}
return ( chid - GetStartUniqueId() );
}