Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

PmtGeometry - source file

// $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() );
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005