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

// ESAF : Euso Simulation and Analysis Framework
// class KOpticalSystem
// $Id: KOpticalSystem.cc 2836 2009-07-16 08:22:02Z biktem $
// J.Watts created , Sep 21 2003
//

#include "Ktrace_optF1v1.hh"
#include "KEUSO_optF1v1.hh"
#include "KOpticalSystem.hh"
#include "Config.hh"
#include "EsafRandom.hh"
#include "EEvent.hh"
#include <TMath.h>
using namespace sou;

ClassImp(KOpticalSystem)

Int_t Kread_para(char*);
Int_t KReadSD(char*);
int Ktrace(double*, double*);

//______________________________________________________________________________
KOpticalSystem::KOpticalSystem(): OpticalSystem() {
    //
    // Constructor
    //

    fDZdown=0*mm;
    fDZup=2002*mm;

    // take the size of the system from Ktrace
    fR=R_WALL;
    fPos=EVector(0,0,Conf()->GetNum("KOpticalSystem.fPos.Z")*mm);

    fFirstLensDZ=169*mm;
    fSecondLensDZ=226*mm;

#ifdef DEBUG
/*    cout << "FirstLensTop() "<< FirstLensTop() << endl;
    cout << "FirstLensBottom() "<< FirstLensBottom() << endl;
    cout << "SecondLensTop() "<< SecondLensTop() << endl;
    cout << "SecondLensBottom() "<< SecondLensBottom() << endl;*/
#endif /* DEBUG */

    if(Kread_para("config/Optics/KOpticalSystem")<0){
       cout<<"reading surface data files failed\n"<<endl;
       exit(-2);
        }
   if(KReadSD("config/Optics/KOpticalSystem")<0){
       cout<<"reading diffractive data files failed\n"<<endl;
       exit(-2);
        }
}

//______________________________________________________________________________
Photon *KOpticalSystem::Transport(Photon *p) const {
    double ph_in[8];
    double ph_out[8];

    EVector dummy=p->pos-fPos;
    ph_in[0]=dummy[X]/mm;
    ph_in[1]=dummy[Y]/mm;
    ph_in[2]=dummy[Z]/mm;

    ph_in[3]=p->dir[X];
    ph_in[4]=p->dir[Y];
    ph_in[5]=p->dir[Z];

    ph_in[6]=p->wl/nm;
    ph_in[7]=p->time/ns;


#ifdef DEBUG/*
    cout<<"KOpticalSystem"<<endl;
    cout<<"position = "<<p->pos<<endl;
    cout<<"direction = "<<p->dir<<endl;
    double theta=p->dir.Angle(EVector(0.,0.,1.));
    cout<<"theta = "<<theta/deg<<endl;
    cout<<"wl = "<<p->wl/nm<<endl;
    cout<<"time = "<<p->time/ns<<endl;*/
#endif /* DEBUG */
    int flag=Ktrace(ph_in, ph_out);
#ifdef DEBUG
    //cout << flag << endl;
#endif /* DEBUG */

    // FIXME: add to ktrace possibility to save latest photon information, for cases when it's absorbed
    if ( !flag ) {
        dummy[X]=ph_out[0]*mm;
        dummy[Y]=ph_out[1]*mm;
        dummy[Z]=ph_out[2]*mm;
        p->pos=dummy+fPos;

        p->dir[X]=ph_out[3];
        p->dir[Y]=ph_out[4];
        p->dir[Z]=ph_out[5];

        p->wl   = ph_out[6]*nm;
        p->time = ph_out[7]*ns;
    }
    p->fate = flag;

    if ( flag != 0 )
        return 0;

    Double_t DL;
    // sometime Ktrace returns photons out of the boundaries
    if (p->pos.Perp() > R_WALL || p->pos[Z] > Top() || p->pos[Z] < Bottom() ){
        // the photon is outside the cyl, trace it back to the borders
        EVector savedir = p->dir;
        p->dir = -p->dir;
        if ((DL = CylinderIntersection(p,kTRUE)) > 0){
            p->pos+=p->dir.Unit()*DL;
            p->time-=DL/TMath::C();
        }
        p->dir=savedir;
    } else if ((DL = CylinderIntersection(p,kTRUE)) > 0){
    // all the photon must be returned on the boundaries of the optics

        p->pos+=p->dir.Unit()*DL;
        p->time+=DL/TMath::C();
    }


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

    if (IsGoingOut(p))
        return p;
    else
        return 0;
}

//______________________________________________________________________________
Bool_t KOpticalSystem::IsGoingOut( Photon *p ) const {
    //
    // helper method to select outward going photons
    //

    Bool_t side, face;
    // check if photon is on the lat surface
    side = TMath::Abs((p->pos-fPos).Perp()-fR) < kTolerance  &&
        (p->pos-fPos).XYvector()*p->dir.XYvector() > 0 &&
        ((p->pos[Z] < FirstLensTop() && p->pos[Z] > FirstLensBottom()) ||
         (p->pos[Z] < SecondLensTop() && p->pos[Z] > SecondLensBottom()) );
    // or if it's on the top/bottom face
    face = (TMath::Abs(p->pos[Z]-FirstLensBottom()) < kTolerance && p->dir[Z]<0) ||
        (TMath::Abs(p->pos[Z]-SecondLensTop()) < kTolerance && p->dir[Z]>0);
    return side || face;
}

//______________________________________________________________________________
Bool_t KOpticalSystem::HitSide(Photon *p) const {
    //
    // Check if p hist the side of the Cylinder.
    //

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