Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

OpticalAdaptor - source file

// $Id: OpticalAdaptor.cc,v 1.25 2005/10/02 14:17:55 thea Exp $
// Author: D.Demarco, M.Pallavicini

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: OpticalAdaptor                                                       *
 *  Package: Optics                                                          *
 *  Coordinator: Alessandro.Thea                                             *
 *                                                                           *
 *****************************************************************************/

//______________________________________________________________________________
//
//   Optical Adapter abstract interface
//   ==================================
//   
//   Abstarct base class providing the interface of an Optical Adaptor.
//   It provides some general methods that can be used in the actual
//   implementation.
//

#include <math.h>
#include "TVector.h"
#include "OpticalAdaptor.hh"
#include "euso.hh"

using namespace TMath;
using namespace sou;

ClassImp(OpticalAdaptor)

Interpolate* OpticalAdaptor::fgBG3 = 0;
Interpolate* OpticalAdaptor::fgMultiLayer = 0;

//______________________________________________________________________________
OpticalAdaptor::OpticalAdaptor() : pmt(0), fSide(0), fHeight(0), fReflectivity(0), fRadiomEff(1.) {
    //
    // Constructor
    //

    fFilterType = kSquare; //default

    if (!fgBG3) {
        fgBG3 = new Interpolate("config/Optics/OpticalAdaptor/BG3.dat",2);
        fgBG3->SetXUnit(nm);
    }
    
    if (!fgMultiLayer) {
        fgMultiLayer = new Interpolate("config/Optics/OpticalAdaptor/Multilayer.dat",2);
        fgMultiLayer->SetXUnit(nm);
    }
}

//______________________________________________________________________________
 void OpticalAdaptor::SetGeometry( const PmtGeometry *g ) {
    //
    // Associate this with PmtGeometry g
    //

    pmt = g;
    fXaxis=pmt->GetX();
    fYaxis=pmt->GetY();
    fZaxis=pmt->GetZ();
    fPos=pmt->Position();

    // fRot is the rotation matrix to go from global to local
    // fInvRot is the rotation matrix to go from global to local
    fRot=fRot.RotateAxes(fXaxis, fYaxis, fZaxis);
    fInvRot=fRot.Inverse();

}


//______________________________________________________________________________
 Bool_t OpticalAdaptor::IsAbsorbed(Photon *p) const {
    //
    // Applies the filter to the photon 
    //

    Double_t value = 0;
    switch ( fFilterType ) {
        case kSquare:
            // really stupid filter
            return ((p->wl < 300*nm || p->wl > 400*nm) &&
                   EsafRandom::Get()->Rndm() > fRadiomEff);
        case kBG3:
            if ( p->wl < fgBG3->GetXmin() || p->wl > fgBG3->GetXmax() )
                return kTRUE;

//            value = Exp(fFilterThickness*Log(fgBG3->GetValue(p->wl))/(1.*mm));      
            value = Power(fgBG3->GetValue(p->wl),
                    fFilterThickness/(1.*mm));      

            return ( EsafRandom::Get()->Rndm() > value );
        case kMultilayer:
            if ( p->wl < fgMultiLayer->GetXmin() || p->wl > fgMultiLayer->GetXmax() )
                return kTRUE;

//            value = Exp(fFilterThickness*Log(fgMultiLayer->GetValue(p->wl))/(1.*mm));      
            value = Power(fgMultiLayer->GetValue(p->wl), 
                    fFilterThickness/(1.*mm));      

            return ( EsafRandom::Get()->Rndm() > value );
        default:
            FatalError("Unknown filter type.");

    }

    return kTRUE;
}

//______________________________________________________________________________
 Double_t OpticalAdaptor::IsHit( const Photon& ph) const {
    //
    // checks if ph is going to hit the top face of the OA and, in that case,
    // it returns the distance it has to travel. Otherwise returns -1
    //
    if (fZaxis.Dot(ph.dir) > -kTolerance) {
        cerr << "OpticalAdaptor::IsHit: photon is not going toward this oa" << endl;
        return -1;
    }
    EVector topface_pos = fPos+fZaxis.Unit()*GetThickness();
    EVector intPoint = ph.pos-ph.dir*((ph.pos-topface_pos).Dot(fZaxis)/ph.dir.Dot(fZaxis));
    if ( isInside(goLocal(intPoint-fPos), TOP) )
        return (intPoint-ph.pos).Mag();
    return -1;
}

//______________________________________________________________________________
 int OpticalAdaptor::whichFace(Photon *p) const {
    // lpos is the position of photon in the local reference
    EVector lpos=p->pos - fPos;
    lpos=goLocal(lpos);

    if(fabs(lpos[Z]-fHeight)<kTolerance) {
        // photon on the top surface
#ifdef DEBUG
        cout<<"TOP"<<endl;
#endif /* DEBUG */
        return TOP;
    } else if(fabs(lpos[X]-fSide)<kTolerance) {
        // photon on the right surface
#ifdef DEBUG
        cout<<"RIGHT"<<endl;
#endif /* DEBUG */
        return RIGHT;
    } else if(fabs(lpos[X])<kTolerance) {
        // photon on the left surface
#ifdef DEBUG
        cout<<"LEFT"<<endl;
#endif /* DEBUG */
        return LEFT;
    } else if(fabs(lpos[Y]+fSide)<kTolerance) {
        // lpos[Y] < 0
        // photon on the front surface
#ifdef DEBUG
        cout<<"FRONT"<<endl;
#endif /* DEBUG */
        return FRONT;
    } else if(fabs(lpos[Y])<kTolerance) {
        // photon on the back surface
#ifdef DEBUG
        cout<<"BACK"<<endl;
#endif /* DEBUG */
        return BACK;
    } else if(fabs(lpos[Z])<kTolerance) {
        // photon on the bottom surface
#ifdef DEBUG
        cout<<"BOTTOM"<<endl;
#endif /* DEBUG */
        return BOTTOM;
    } else throw runtime_error("whichFace: not of surface");
}

//______________________________________________________________________________
vector<EVector> OpticalAdaptor::intPoints(Photon *p) const {
    EVector lpos, ldir;
    lpos=p->pos - fPos;
    ldir=p->dir;
    lpos=goLocal(lpos);
    ldir=goLocal(ldir);

    int face=whichFace(p);

#ifdef DEBUG
    cout<<"lpos: "<<lpos<<endl;
    cout<<"ldir: "<<ldir<<endl;
#endif /* DEBUG */

    vector<EVector> ips(6);
    double dist;
    EVector out(1e6*mm, 1e6*mm, 1e6*mm);

    // TOP
    if(face==TOP || ldir[Z] < 0) ips[TOP]=out;
    else {
        dist=fHeight-lpos[Z];
        ips[TOP]=lpos + ldir*(dist/ldir[Z]);
#ifdef DEBUG
        cout<<"dist, ips[TOP]: "<<dist<<", "<<ips[TOP]<<endl;
#endif /* DEBUG */
    }



    // BOTTOM
    if(face==BOTTOM || ldir[Z] > 0) ips[BOTTOM]=out;
    else {
        dist=-lpos[Z];
        ips[BOTTOM]=lpos + ldir*(dist/ldir[Z]);
#ifdef DEBUG
        cout<<"dist, ips[BOTTOM]: "<<dist<<", "<<ips[BOTTOM]<<endl;
#endif /* DEBUG */
    }



    // FRONT
    if(face==FRONT || ldir[Y] > 0) ips[FRONT]=out;
    else {
        dist=-lpos[Y]-fSide;
        ips[FRONT]=lpos + ldir*(dist/ldir[Y]);
#ifdef DEBUG
        cout<<"dist, ips[FRONT]: "<<dist<<", "<<ips[FRONT]<<endl;
#endif /* DEBUG */
    }



    // BACK
    if(face==BACK || ldir[Y] < 0) ips[BACK]=out;
    else {
        dist=-lpos[Y];
        ips[BACK]=lpos + ldir*(dist/ldir[Y]);
#ifdef DEBUG
        cout<<"dist, ips[BACK]: "<<dist<<", "<<ips[BACK]<<endl;
#endif /* DEBUG */
    }



    // RIGHT
    if(face==RIGHT || ldir[X] < 0) ips[RIGHT]=out;
    else {
        dist=fSide-lpos[X];
        ips[RIGHT]=lpos + ldir*(dist/ldir[X]);
#ifdef DEBUG
        cout<<"dist, ips[RIGHT]: "<<dist<<", "<<ips[RIGHT]<<endl;
#endif /* DEBUG */
    }



    // LEFT
    if(face==LEFT || ldir[X] > 0) ips[LEFT]=out;
    else {
        dist=0-lpos[X];
        ips[LEFT]=lpos + ldir*(dist/ldir[X]);
#ifdef DEBUG
        cout<<"dist, ips[LEFT]: "<<dist<<", "<<ips[LEFT]<<endl;
#endif /* DEBUG */
    }

    return ips;
}

//______________________________________________________________________________
 bool OpticalAdaptor::isInside(const EVector &v, int face) const {
    switch(face) {
        case TOP:
            if( v[X] > kTolerance && v[X] < fSide-kTolerance &&
                    -v[Y] > kTolerance && -v[Y] < fSide-kTolerance &&
                    fabs(v[Z]-fHeight) < kTolerance ) return true;
            break;

        case BOTTOM:
            if( v[X] > kTolerance && v[X] < fSide-kTolerance &&
                    -v[Y] > kTolerance && -v[Y] < fSide-kTolerance &&
                    fabs(v[Z]) < kTolerance ) return true;
            break;

        case RIGHT:
            if( v[Z] > kTolerance && v[Z] < fHeight-kTolerance &&
                    -v[Y] > kTolerance && -v[Y] < fSide-kTolerance &&
                    fabs(v[X]-fSide) < kTolerance ) return true;
            break;

        case LEFT:
            if( v[Z] > kTolerance && v[Z] < fHeight-kTolerance &&
                    -v[Y] > kTolerance && -v[Y] < fSide-kTolerance &&
                    fabs(v[X]) < kTolerance ) return true;
            break;

        case FRONT:
            if( v[Z] > kTolerance && v[Z] < fHeight-kTolerance &&
                    v[X] > kTolerance && v[X] < fSide-kTolerance &&
                    fabs(v[Y]+fSide) < kTolerance ) return true;
            break;

        case BACK:
            if( v[Z] > kTolerance && v[Z] < fHeight-kTolerance &&
                    v[X] > kTolerance && v[X] < fSide-kTolerance &&
                    fabs(-v[Y]) < kTolerance ) return true;
            break;
    }
        return false;
}
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