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

// $Id: LOpticalSystem.cc 2849 2010-11-11 11:59:41Z biktem $
// Author: Dmitry V.Naumov   2007/11/13

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

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

#include "LOpticalSystem.hh"
#include "Ltrace_optF1v4.hh"
#include "Config.hh"
#include "EsafRandom.hh"
#include "EEvent.hh"
#include <TMath.h>

#define CHECK_TIMING
#ifdef CHECK_TIMING
#include <TBenchmark.h>
#include <TTree.h>
#include <TFile.h>
#include <TH1D.h>
#endif

using namespace sou;
ClassImp(LOpticalSystem)
tel_param  fLOpticalSystemParam;

//#define DEBUG
// extern  TBenchmark* theBenchMark;
// extern  TTree* theBenchTree;
// extern  double theTime;
 extern  TH1D* theh;
//_____________________________________________________________________________
LOpticalSystem::LOpticalSystem(): fLogFile(0) {
    //
    // Constructor
    //

    #ifdef CHECK_TIMING
    CreateHisto();
    #endif

    if(Lread_tel_param("config/Optics/LOpticalSystem/telparm_F1v4.dat",&fLOpticalSystemParam,fLogFile)!=0) exit(-1);

    int code;
    cout<<"init"<<endl;
    if( (code=LReadLensData(&fLOpticalSystemParam,fLogFile))<0 ) {
       fprintf(stderr,"reading lens surface data files failed. Code %i\n", code);
       exit(code);
    }
    if( (code=LReadDiffractData(fLOpticalSystemParam.Llens_dir,fLogFile))<0 ) {
        fprintf(stderr,"reading diffractive opt. data files failed (%s). Code %i\n", fLOpticalSystemParam.Llens_dir,code);
        exit(code);
    }

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

    fDZdown=1*mm;//1
    fDZup=2072.*mm;//2073
    //fFirstLensDZ=15*mm;
    //fFirstLensDZ=412.942*mm;
    fFirstLensDZ=412.942*mm;
    //fSecondLensDZ=10*mm;
    //fSecondLensDZ=210.469*mm;
    fSecondLensDZ=210.469*mm;

}

//_____________________________________________________________________________
LOpticalSystem::~LOpticalSystem() {
    //
    // Destructor
    //
    cout<<"Destructor ************************************"<<endl;
    #define CHECK_TIMING
    #ifdef CHECK_TIMING

    DeleteHisto();
    #endif

}
//______________________________________________________________________________
Photon* LOpticalSystem::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<<"LOpticalSystem"<<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 = Ltrace_main(&fLOpticalSystemParam,ph_in, ph_out, fLogFile);

    // FIXME: add to ltrace 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<<"LOpticalSystem 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 ltrace returns photons out of the boundaries
    if (p->pos.Perp() > fLOpticalSystemParam.Lr_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 LOpticalSystem::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 LOpticalSystem::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;
}