source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/electronics/src/PmtGeometry.cc @ 117

Last change on this file since 117 was 117, checked in by moretto, 11 years ago

ESAF version compilable on mac OS
File size: 11.3 KB
Line 
1// $Id: PmtGeometry.cc 2927 2011-06-16 18:02:41Z mabl $
2// Author: D.Demarco, M.Pallavicini
3
4/*****************************************************************************
5 * ESAF: Euso Simulation and Analysis Framework                              *
6 *                                                                           *
7 *  Id: PmtGeometry                                                          *
8 *  Package: Optics                                                          *
9 *  Coordinator: Marco.Pallavicini                                           *
10 *                                                                           *
11 *****************************************************************************/
12
13//______________________________________________________________________________
14//
15//    Photomultiplier Geometry
16//    ========================
17//
18//    Geometry descriptor of a PMT this object can be attached to any
19//    Photomultiplier object In the constructor it retrieves the relevant
20//    parameters from the Config object and store them in class variables Local
21//    Coordinates: x,y,z
22//    z is normal to the Photocathode
23//    (0,0,0) is the upper left corner of the pmt photocathode
24//    surface up is defined by the y axis
25//    right is defined by the x axis direction left is -x
26//
27
28#include <iostream>
29#include "PmtGeometry.hh"
30#include "Photomultiplier.hh"
31#include "OpticalAdaptor.hh"
32#include "EusoMapping.hh"
33
34ClassImp(PmtGeometry)
35
36#define MAX_PIXELS 500000
37
38// class variables with geometry parameters
39Double_t PmtGeometry::fgSide        = -1.; // physical total width of the pmt
40Double_t PmtGeometry::fgDeadBorder  = 0.;  // thickness of external dead space
41Double_t PmtGeometry::fgDeadInner   = 0.;  // thickness of internal space
42                                           // ( between pads )
43Double_t PmtGeometry::fgPadSide     = 0.;          // width of one pad
44
45Int_t PmtGeometry::fgNumPads        = 0;   // number of active pads of this mapmt
46Int_t PmtGeometry::fgNumRows        = 0;   // number of rows and columns
47TVector3* PmtGeometry::fgCorners[10][10];
48
49//______________________________________________________________________________
50PmtGeometry::PmtGeometry(TVector3 x, TVector3 n, TVector3 m)
51    : EsafConfigurable(), fPosition(x), fZAxis(n.Unit()) {
52    //
53    // Construtor
54    //
55    fXAxis = ( m.Cross(fZAxis)).Unit();
56    fYAxis = fZAxis.Cross(fXAxis);
57
58    Double_t dx = fYAxis.Dot(fZAxis);
59    Double_t dy = fZAxis.Dot(fXAxis);
60    Double_t dz = fXAxis.Dot(fYAxis);
61
62    if ( dz > PMT_TOLERANCE || dx > PMT_TOLERANCE || dy > PMT_TOLERANCE )
63        throw runtime_error(Form("Something's wrong in pmt's axis dx = %f, dy = %f, dz = %f",dx, dy, dz));
64
65    fRight = NULL;
66    fLeft = NULL;
67    fFront = NULL;
68    fBack = NULL;
69    pPmt = NULL;
70    pOA = NULL;
71    if ( NumPads() <= 0 )
72        SetPmtGeometry();
73}
74
75//______________________________________________________________________________
76PmtGeometry::~PmtGeometry() {
77    //
78    // Destructor
79    //
80    SafeDelete(pOA);
81}
82
83//______________________________________________________________________________
84void PmtGeometry::ResetClass() {
85    //
86    // Reset static data members
87    //
88
89    fgSide            = -1.;
90    fgDeadBorder      = 0.;
91    fgDeadInner       = 0.;
92
93    fgPadSide         = 0.;
94
95    fgNumPads         = 0;
96    fgNumRows         = 0;
97}
98
99//______________________________________________________________________________
100void PmtGeometry::SetNearest(PmtGeometry* r, PmtGeometry* l,
101                             PmtGeometry* f, PmtGeometry* b) {
102    //
103    // Set pointers to nearest neighbors
104    //
105
106    fRight = r;
107    fLeft = l;
108    fFront = f;
109    fBack = b;
110
111#ifdef DEBUG
112//  cout << "Set nearest called for PMT = " << Pmt()->Id() << endl;
113//  if ( b )
114//    cout << "  Back = " << b->Pmt()->Id() << endl;
115//  if ( f )
116//    cout << "  Front = " <<  f->Pmt()->Id() << endl;
117//  if ( r )
118//    cout << "  Right = " << r->Pmt()->Id() << endl;
119//  if ( l )
120//    cout << "  Left = " <<  l->Pmt()->Id() << endl;
121#endif
122}
123
124
125//______________________________________________________________________________
126void PmtGeometry::InsertOA( OpticalAdaptor* oa ){
127    //
128    // Insert oa in front of the pmt and displace geometry of oa->GetHeight()
129    //
130
131    // push the pmt back
132    fPosition -= GetZ()*oa->GetThickness();
133
134    oa->SetGeometry(this);
135
136    SetOA(oa);
137
138}
139//______________________________________________________________________________
140void PmtGeometry::SetPmtGeometry() {
141    //
142    // Set local geometry infos and global parameters only once
143    //
144
145    if ( fgNumPads == 0 ) {
146        ConfigFileParser* conf = Config::Get()->GetCF("Electronics","Photomultiplier");
147        fgNumRows = (Int_t)conf->GetNum("Photomultiplier.PmtSide");
148        fgNumPads = fgNumRows*fgNumRows;
149        fgSide = conf->GetNum("Photomultiplier.PmtSize");
150        fgDeadBorder = conf->GetNum("Photomultiplier.PmtDeadLateral");
151        fgDeadInner = conf->GetNum("Photomultiplier.PmtDeadInner");
152        fgPadSide = ( fgSide - 2*fgDeadBorder - (fgNumRows-1)*fgDeadInner ) / fgNumRows;
153
154        // compute the corner position of each pad in local coordinate
155        for(Int_t row=0; row<=fgNumRows; row++) {
156            Double_t row_shift = fgDeadBorder + row * fgPadSide + row * fgDeadInner;
157            for(Int_t col=0; col<=fgNumRows; col++) {
158                Double_t col_shift = fgDeadBorder + col * fgPadSide + col * fgDeadInner;
159                fgCorners[row][col] = new TVector3( col_shift, row_shift, 0. );
160            }
161        }
162    }
163}
164
165//______________________________________________________________________________
166Bool_t PmtGeometry::IsInside( const Photon& ph ) const {
167    // check whether a Photon is inside the geometrical acceptance of this pmt
168    // it returns true when the Photon is within the external side and its
169    // distance from the surface is less than PMT_TOLERANCE
170
171    // check if the photon is on the cathode surface
172    TVector3 a = ( ph.pos - Position());
173    Double_t b = a.Dot(Normal());
174    if ( b < 0. ) b = -b;
175    if ( b > PMT_TOLERANCE ) {
176        cerr<<"PmtGeometry::IsInside: photon not on PMT surface, " << b << " mm away" <<endl;
177        return false;
178    }
179
180    // check if the photon direction is going to intersect the cathode
181    b = Normal().Dot( ph.dir );
182    if ( b > -PMT_TOLERANCE ) {
183        cerr<<"PmtGeometry::IsInside: photon not intersecting PMT surface"<<endl;
184        return false;
185    }
186
187    // check if it is inside
188    Double_t dx = a.Dot(GetX().Unit());  // along x is positive
189    if ( dx < PMT_TOLERANCE || dx > (fgSide-PMT_TOLERANCE) )
190        return false;
191
192    Double_t dy = -a.Dot(GetY().Unit()); // along y must be negative
193    if ( dy < PMT_TOLERANCE || dy > (fgSide-PMT_TOLERANCE) )
194        return false;
195
196    return true;
197}
198
199//______________________________________________________________________________
200Double_t PmtGeometry::IsHit( const Photon& ph ) const {
201    //
202    // check whether a Photon
203    //
204
205    // ph is going away from cathode
206    if ( ph.dir.Dot(Normal()) > -PMT_TOLERANCE ) return -1;
207
208    // calculate the hit point between ph and cathode plane
209
210    Photon dummy;
211    dummy.pos=ph.pos-ph.dir*((ph.pos-Position()).Dot(Normal())/ph.dir.Dot(Normal()));
212    dummy.dir=ph.dir;
213    if ( IsInside(dummy) ){
214        cout << "PmtGeometry::IsHit\n"
215             << "ph.pos = " << ph.pos << "   ph.dir = " << ph.dir << endl;
216        cout << "Position() = " << EVector(Position()) << "   Normal() = " << EVector(Normal()) << endl;
217
218        cout << "DL = " << (dummy.pos-ph.pos).Mag() << endl;
219
220        return (dummy.pos-ph.pos).Mag();
221    }
222    return -1;
223}
224
225//______________________________________________________________________________
226TVector3 PmtGeometry::Center() const {
227    //
228    // Return the position of the center of the pmt
229    //
230    TVector3 a = Position() + (fgSide/2.) * GetX().Unit();
231    a += (-fgSide/2.) * GetY().Unit();
232    return a;
233}
234
235// return local coordinate position of a point in global coord
236// the transofrmation works for point on the pmt surface
237TVector3 PmtGeometry::LocalCoord( const TVector3& pos) const{
238    TVector3 diff = pos - Position();
239    Double_t deltax = diff.Dot(GetX().Unit());
240    Double_t deltay = diff.Dot(GetY().Unit());
241    return TVector3( deltax, deltay, 0.0);
242}
243
244// return the number of a Pad hit by a Photon
245Int_t PmtGeometry::Pad( const Photon& ph ) const{
246    return Pad( LocalCoord( ph.pos ) );
247}
248
249// return the pad number corresponding to local coordinate x
250Int_t PmtGeometry::Pad( const TVector3& x) const{
251    Int_t row,col;
252
253    for(row=0; row < Rows(); row++ ) {
254        Double_t yy = -x.y();// y axis is negative!
255        Double_t yy_low = Corner(row,0).y();
256        Double_t yy_high = Corner(row+1,0).y();
257        if ( yy > yy_low && yy < yy_high)
258            break;
259    }
260
261    if ( row == Rows())
262        return -1;
263
264    for(col=0; col < Rows(); col++ ) {
265        Double_t xx = x.x();
266        Double_t xx_low = Corner(0,col).x();
267        Double_t xx_high = Corner(0,col+1).x();
268        if ( xx > xx_low && xx < xx_high)
269            break;
270    }
271
272    if ( col == Rows())
273        return -1;
274
275    return ( col + ( row*Rows() ) );
276}
277
278//______________________________________________________________________________
279const PmtGeometry* PmtGeometry::Nearest( OAFace side ) const {
280    //
281    // Nearest neighbor
282    //
283    switch ( side ) {
284      case RIGHT:
285          return fRight;   // x > 0
286          break;
287      case LEFT:
288          return fLeft;    // x < 0
289          break;
290      case FRONT:
291          return fFront;   // y < 0
292          break;
293      case BACK:
294          return fBack;    // y > 0
295          break;
296      case TOP:
297      case BOTTOM:
298      default:
299          break;
300    }
301    cerr << "Invalid side type in PmtGeometry::Nearest" << endl;
302    return NULL;
303}
304
305
306//______________________________________________________________________________
307TVector3 PmtGeometry::Position(Int_t r, Int_t c) const {
308    //
309    // Return the position in EUSO coordinate of the center of a pad
310    // r and c are not checked
311    //
312
313    TVector3 p = Position();
314    p += (fgDeadBorder + fgPadSide/2 + c*(fgPadSide+fgDeadInner)) * GetX();
315    p -= (fgDeadBorder + fgPadSide/2 + r*(fgPadSide+fgDeadInner)) * GetY();
316    return p;
317}
318
319//______________________________________________________________________________
320TVector3 PmtGeometry::Position(Int_t ch) const {
321    //
322    // return the position in EUSO coordinate of the center of a pad
323    // ch is the channel number defined as ch = r*n+c
324    //
325    Int_t r = ch / Rows();
326    Int_t c = ch - r*Rows();
327    return Position(r,c);
328}
329
330//______________________________________________________________________________
331Double_t PmtGeometry::ThetaFOV(Int_t ch) const {
332    //
333    // Theta in field of view (FOV) corresponding to the center of this pad
334    //
335    return EusoMapping::Get()->GetThetaFOV( GetUniqueId(ch) );
336}
337
338//______________________________________________________________________________
339Double_t PmtGeometry::PhiFOV(Int_t ch) const {
340    //
341    // Theta in field of view (FOV) corresponding to the center of this pad
342    //
343    return EusoMapping::Get()->GetPhiFOV( GetUniqueId(ch) );
344}
345
346//______________________________________________________________________________
347Int_t PmtGeometry::GetChannel( ChannelUniqueId chid ) const {
348    //
349    // Returns channel associated to unique id
350    //
351    if ( chid < GetStartUniqueId() ) {
352        return -1;
353    }
354    if ( chid > GetLastUniqueId() ) {
355    return -1;
356    }
357    return ( chid - GetStartUniqueId() );
358}
Note: See TracBrowser for help on using the repository browser.