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

// class MacroCellGeometry
// $Id: MacroCellGeometry.cc 2850 2010-11-11 12:01:39Z biktem $
// description of the geometry of a MacroCell
//

#include "MacroCellGeometry.hh"
#include "PmtGeometry.hh"
#include "OpticalAdaptor.hh"
#include "MacroCell.hh"
#include "ElementaryCell.hh"
#include <iostream>

using namespace sou;
using namespace std;

ClassImp(MacroCellGeometry)

//______________________________________________________________________________
MacroCellGeometry::MacroCellGeometry( MacroCell *parent) : pCell(parent) {
    //
    // ctor
    //
}

//______________________________________________________________________________
MacroCellGeometry::~MacroCellGeometry() {
    //
    // dtor
    //
}

//#define DEBUG
//______________________________________________________________________________
Double_t MacroCellGeometry::IsHit(const Photon &p) const {
    //
    //
    //

    // check if p is // to the MC plane
    if ( p.dir*GetNormal() > -TOLERANCE ) {
        return -1;
    }
    EVector dummy=p.pos-GetCenter();
    EVector intPoint=dummy-p.dir*((dummy*GetNormal())/(p.dir*GetNormal()));

    if ( intPoint*fXAxis < fXMax && intPoint*fXAxis > fXMin &&
         intPoint*fYAxis < fYMax && intPoint*fYAxis > fYMin ){

#ifdef DEBUG
        cout << "\nMacroCellGeometry::IsHit\n"
             << "p.pos = " << p.pos << "   p.dir = " << p.dir << endl;
        cout << "fCenter = " << GetCenter() << "   GetNormal() = " << GetNormal() << endl;
        cout << "fXAxis = " << fXAxis << "   fYAxis = " << fYAxis << endl;

        cout << "fXMax = " << fXMax << "   fXMin = " << fXMin
             << "   fYMax = " << fYMax << "   fYMin = " << fYMin << endl;
        cout << "intPoint = " << intPoint+GetCenter()<< "   intPoint*fXAxis = " << intPoint*fXAxis << "   intPoint*fYAxis = " << intPoint*fYAxis << endl;
        cout << "DL = " << (intPoint-dummy).Mag() << endl;
#endif /* DEBUG */

        //printf("%g %g %g dist=%g\n", intPoint.X(),intPoint.Y(),intPoint.Z(),(intPoint-dummy).Mag());
        return (intPoint-dummy).Mag();
    }
    return -1;

}

//______________________________________________________________________________
Photomultiplier* MacroCellGeometry::FindHitPmt(const Photon &p) const{
    //
    // finds which pmt has been hit, if exists
    //
    const vector<Photomultiplier*> &PMTs=*(pCell->Pmts());

    Photomultiplier *HitPMT=NULL;
    for(size_t j=0; j<PMTs.size(); j++) {
        // bring p on pmt cathode

        Double_t DL = PMTs[j]->Geometry()->GetOA()->IsHit(p);

        if( DL > 0.) {
                HitPMT=PMTs[j];
                break;
        }
    }
    return HitPMT;
}

//______________________________________________________________________________
void MacroCellGeometry::Compute() {
    //
    // Compute internal infos
    //

    if ( !Cell() ) {
        FatalError("Invalid MacroCell pointer in MacroCellGeometry");
    }

    // number of elementary cells
    Int_t nECs = Cell()->ECs()->size();
    EVector n(0.,0.,0.), y_axis(0., 0., 0.);

    // compute normal vector and the y axis
    for(Int_t i=0; i<nECs; i++) {
       const ElementaryCell *pEC = Cell()->GetEC(i);
       if ( pEC ) {
           n += pEC->GetNormal();
           y_axis += pEC->GetYAxis();
       } else {
           FatalError("Invalid Elementary Cell in MacroCellGeometry");
       }
    }

    fNormal = n.Unit();
    fXAxis = (y_axis.Cross( n )).Unit();
    fYAxis = fNormal.Cross( fXAxis );

    // compute center
    EVector center(0.,0.,0.);
    for(Int_t i=0; i < nECs; i++) {
        ElementaryCell *ec = Cell()->GetEC(i);
        center += ec->GetCenter();
    }
    center *= (1./(double)Cell()->GetNumECs());

    // compute corners and radius
    for(Int_t i=0; i < nECs; i++) {
        EVector d = Cell()->GetEC(i)->GetCenter() - center;
        // d projection on the MC plane
        if ( !i ){
            fXMin = d.Dot(fXAxis);
            fXMax = d.Dot(fXAxis);
            fYMin = d.Dot(fYAxis);
            fYMax = d.Dot(fYAxis);
        } else {
            if ( d.Dot(fXAxis) < fXMin ) fXMin = d.Dot(fXAxis);
            if ( d.Dot(fXAxis) > fXMax ) fXMax = d.Dot(fXAxis);
            if ( d.Dot(fYAxis) < fYMin ) fYMin = d.Dot(fYAxis);
            if ( d.Dot(fYAxis) > fYMax ) fYMax = d.Dot(fYAxis);
        }
    }


    // shift center of 50 mm along normal to keep the square in front of the macrocell
    fCenter = center + 50. * fNormal.Unit()*mm;

    // shift points out by two full PMT sides
    fXMin-=2. * Cell()->GetEC(0)->GetPmt(1)->Geometry()->Side();
    fXMax+=2. * Cell()->GetEC(0)->GetPmt(1)->Geometry()->Side();
    fYMin-=2. * Cell()->GetEC(0)->GetPmt(1)->Geometry()->Side();
    fYMax+=2. * Cell()->GetEC(0)->GetPmt(1)->Geometry()->Side();


}

void MacroCellGeometry::Draw() const{
//FIXME: to be done (call EC->Draw)
}