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

// $Id: MOpticalSystem.cc 2828 2009-04-01 14:29:16Z biktem $
// Author: Dmitry V.Naumov   2007/11/13

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

//_____________________________________________________________________________
//
// MOpticalSystem
//
// <extensive class description>
//
//   Config file parameters
//   ======================
//
//   <parameter name>: <parameter description>
//   -Valid options: <available options>
//

#include <TMath.h>
#include "EEvent.hh"

#include "Config.hh"
#include "EsafRandom.hh"

#include "MOpticalSystem.hh"
#include "Mtrace_optF1v4.hh"
#include "Ntrace_optics.hh"
using namespace sou;

Mtel_param  fMOpticalSystemParam;

ClassImp(MOpticalSystem)
//_____________________________________________________________________________
MOpticalSystem::MOpticalSystem(): fLogFile(0) {
    //
    // Constructor
    //

    if(Mread_tel_param("config/Optics/MOpticalSystem/telparm_F1v6.dat",&fMOpticalSystemParam,fLogFile)!=0) exit(-1);

    if(ReadLensData(&fMOpticalSystemParam,fLogFile)<0) {
       fprintf(stderr,"reading lens surface data files failed\n");
       exit(-2);
    }
    if(ReadDiffractData(fMOpticalSystemParam.Mlens_dir,fLogFile)<0) {
       fprintf(stderr,"reading diffractive opt. data files failed\n");
       exit(-2);
    }
//Below fOpticalSystemParam for MOpics is filled
    const char* lens_dir = "config/Optics/NOpticalSystem/";
    const char* tel_par = "CPC_2007";
    ReadTelParm(lens_dir,tel_par);

    fDZdown=1*mm;
    fDZup=2197*mm;

    // take the size of the system from Ltrace
    fR=fMOpticalSystemParam.Mr_wall;
    fPos=EVector(0,0,Conf()->GetNum("MOpticalSystem.fPos.Z")*mm);

    //fFirstLensDZ=15.9*mm;
    //fFirstLensDZ=437.784*mm;
    fFirstLensDZ=437.784*mm;
    //fSecondLensDZ=10.6*mm;
    fSecondLensDZ=223.217*mm;

}

//_____________________________________________________________________________
MOpticalSystem::~MOpticalSystem() {
    //
    // Destructor
    //
}
//______________________________________________________________________________
Photon* MOpticalSystem::Transport(Photon *p) const {
    Double_t ph_in[8];
    Double_t 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<<"MOpticalSystem"<<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 = Mtrace_main(&fMOpticalSystemParam,ph_in, ph_out, fLogFile);
#ifdef DEBUG
#endif /* DEBUG */

    // FIXME: add to mtrace 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;
#ifdef DEBUG
    cout<<"MOpticalSystem OUT"<<endl;
    cout<<"position = "<<p->pos<<endl;
    cout<<"direction = "<<p->dir<<endl;
    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
    if ( flag != 0 ) return 0;

    Double_t DL;
    // sometime Ktrace returns photons out of the boundaries
    if (p->pos.Perp() > fMOpticalSystemParam.Mr_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 MOpticalSystem::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 MOpticalSystem::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;
}