source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/optics/src/OpticalSystem.cc @ 114

// ESAF : Euso Simulation and Analysis Framework
// class TestOpticalAdaptor
// $Id: OpticalSystem.cc 2835 2009-07-16 03:22:42Z biktem $
//
#include <TMath.h>
#include "TRotation.h"

#include "OpticalSystem.hh"
#include "EsafRandom.hh"
#include "EEvent.hh"
#include "EConst.hh"

using namespace sou;
using namespace TMath;
using namespace EConst;

#define DEBUG

ClassImp(OpticalSystem)
ClassImp(TestOpticalSystem)

//______________________________________________________________________________
TestOpticalSystem::TestOpticalSystem() : OpticalSystem() {
        fPos=EVector(0,0,Conf()->GetNum("TestOpticalSystem.pos.Z"))*mm;
        fDZup=Conf()->GetNum("TestOpticalSystem.DZup")*mm;
        fDZdown=Conf()->GetNum("TestOpticalSystem.DZdown")*mm;
        fFocalDistance=Conf()->GetNum("TestOpticalSystem.FocalDistance")*mm;
        fR=Conf()->GetNum("TestOpticalSystem.Radius")*mm;

        psf=new Interpolate("config/Optics/"+Conf()->GetStr("TestOpticalSystem.psf.filename"));
        tr=new Interpolate("config/Optics/"+Conf()->GetStr("TestOpticalSystem.tr.filename"));

}

//______________________________________________________________________________
Photon *TestOpticalSystem::Transport(Photon *p) const {

    // theta is the angle between the photon direction and the lens axis
    double theta=p->dir.Angle(EVector(0.,0.,1.));
    if(theta > TMath::PiOver2()) theta=TMath::Pi() - theta;
#ifdef DEBUG
    cout<<"position = "<<p->pos<<endl;
    cout<<"direction = "<<p->dir<<endl;
    cout<<"theta: "<<theta/TMath::Pi()*180.<<endl;
#endif /* DEBUG */

    // check if photon is in FOV (FOV=30 deg)
    if( theta > 30*deg) {
        // photon outside FOV
        // destroy it
#ifdef DEBUG
        cout<<"photon outside FOV: destroyed"<<endl;
#endif /* DEBUG */
        p=0;
        return p;
    }
    // photon inside FOV

    // check if it is absorbed or refracted
    TRandom* rndm = EsafRandom::Get();
    if(rndm->Rndm() > tr->GetValue(theta/deg)) {
#ifdef DEBUG
        cout<<"photon absorbed by opticalSystem: destroyed"<<endl;
#endif /* DEBUG */
        p=0;
        return p;
    }

    EVector dir=p->dir.Unit();

    // propagate photon to the center of the lens
    double DZ=(dir[Z]>0 ? fDZdown : -fDZup);
    EVector dummy=dir*(DZ/dir[Z]);
    p->time+=dummy.Mag()/Clight();
    p->pos+=dummy;

    if(p->pos.Perp() > fR) {
        // hit wall
#ifdef DEBUG
        cout<<"Hit Wall: "<<p->pos<<endl;
#endif /* DEBUG */
        p=0;
        return p;
    }

    // calculate new direction
    // This is for a plane focalplane
    // EVector newDir=dir*(fFocalDistance/fabs(dir[Z])) - (p->pos - fPos);

    // This is for a spherical focalplane with radius fFocalDistance
    // newDir is the point on the focal surface
    double thetain=p->dir.Theta();       // save original theta
                                            // used later for PSF
    EVector newDir=dir*fFocalDistance - (p->pos - fPos);
    p->dir=newDir.Unit();

    // Point Spread Function
    EVector PSF=getPSF(p, thetain);
    newDir+=PSF;
    p->dir=newDir.Unit();

    // propagate photon to the end of the lens
    dir=p->dir;
    DZ=(dir[Z]>0 ? fDZup : -fDZdown);
    dummy=dir*(DZ/dir[Z]);
    p->time+=dummy.Mag()/Clight();
    p->pos+=dummy;

    p->pos[Z]=(dir[Z]>0 ? Top() : Bottom());

    if(p->pos.Perp() > fR) {
        // hit wall
#ifdef DEBUG
        cout<<"Hit Wall: "<<p->pos<<endl;
#endif /* DEBUG */
        p=0;
        return p;
    }

#ifdef DEBUG
    cout<<"new position = "<<p->pos<<endl;
    cout<<"new direction = "<<p->dir<<endl;
#endif /* DEBUG */

    return p;
}

//______________________________________________________________________________
EVector TestOpticalSystem::getPSF(Photon *p, double thetain) const {
    // theta is the angle between the photon direction and the lens axis
    double theta=p->dir.Angle(EVector(0.,0.,1.));
    if(theta > TMath::PiOver2()) theta=TMath::Pi() - theta;

    double psf_t=psf->GetValue(thetain/deg);
#ifdef DEBUG
    cout<<"psf("<<thetain/deg<<"): "<<psf_t<<endl;
#endif /* DEBUG */

    TRandom* rndm=EsafRandom::Get();
    double r=rndm->Gaus(0., psf_t);
    double phi=rndm->Rndm()*360.*deg;
    EVector psf_v(1,0,0);
    psf_v.SetPhi(phi);
    psf_v.SetMag(r);
#ifdef DEBUG
    cout<<"psf_v: "<<psf_v<<endl;
    cout<<"p->dir: "<<p->dir<<endl;
#endif /* DEBUG */

    TRotation rot;
    EVector axis=p->dir.Cross(EVector(0.,0.,1.));
    rot=rot.Rotate(-asin(axis.Mag()), axis.Unit());
    psf_v.Transform(rot);

#ifdef DEBUG
    cout<<"psf_v = "<<psf_v<<endl;
    cout<<"dot = "<<psf_v.Dot(p->dir)<<endl;
#endif /* DEBUG */

    return psf_v;
}

//______________________________________________________________________________
Bool_t OpticalSystem::HitSide(Photon *p) const{
    // checks if p is going to hit the dead border of the optics when it's goint toward
    // the optics i.e. if it's going to hit the walls between lenses from outside


    Double_t dt = CylinderIntersection( p );
    if (dt < 0 || (p->pos+p->dir*dt)[Z]-fDZup < kTolerance
            || (p->pos+p->dir*dt)[Z]-fDZdown < kTolerance )
        return kFALSE;
    return kTRUE;
}