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

// $Id: ParamOpticalSystem.cc 2922 2011-06-12 14:21:23Z mabl $
// Author: A.Thea   2004/07/19

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

//_____________________________________________________________________________
//  
//   Parameterized Optical System
//   ============================
//
//   Parameterized optical system. This Optical System is done to reproduce the
//   behaviour of an optics design through a set of parameters.
//   The main purpose is to simulate and OS with a given Triggering efficacy
//   without the need of the full MC simulation.
//
//   ParamOpticalSystem behaves as an ideal optical system with a given focal 
//   surface.
//
//   The position an incoming photon hits the focal surface is calculated 
//   accordingly the relation
//
//   d = (Dmax/ThetaMax)*theta
//
//   where:
//   d is the distance of the impact point on the FS from the center of the FS
//   calculated along the FS.
//
//   Dmax is the maximum distance.
//
//   theta is the angle between the incoming photon direction and the optical
//   axis.
//
//   ThetaMax is the maximum value of theta accepted by the optics.
//
//  
//   Config file parameters
//   ======================
//
//   fPos.Z [mm]: Z coordinate of the bottom base of the OpticalSystem in
//   the Detector Reference System.
//
//   fRadius [mm]: Radius of the Optical Sytem.
//
//   fRmax [mm]: Maximum radius of the Focal Surface.
//
//   fDZ [mm]: Heigth of the Optical System.
//
//   fThetaMax [deg]: Maximum valid incident angle of an incoming photon.
//
//   fFocalSurfaceZ [mm]: Z coordinate of the center (tip) of the focal
//   surface.
//
//   fDataDirectory: path where data files are stored. If it begins with '/'
//   path is assumed to be absolute, otherwise relative to this cfg file current
//   position
//   
//   fTotalEfficacy.fFilename [string]: two column file containing the
//   total efficacy as a function of the incident angle.
//   
//   fAnglesOfIncidence.fFilename [string]: two column file containing the chief
//   ray incident angle on the focal surface as a function of the incident
//   angle on the pupil.
//
//   fPsfType [string]: type of PSF to use.
//   * options
//      - point: No spread.
//      - gauss: Gaussian psf. RMS and angles of incidence listed in
//               fGaussSpread.filename
//
//   fGaussSpread.fFilename [string]: three column file containing RMS of the
//   gaussian psf and focal surface incident photons cone aperture.
//
//BEGIN_HTML <!--
/* -->
<!--*/
// -->END_HTML
//

#include "ParamOpticalSystem.hh"
#include "Config.hh"
#include "Etypes.hh"
#include "EConst.hh"
#include "EsafRandom.hh"
#include "IdealFocalSurface.hh"
#include "KIdealFocalSurface.hh"
#include "JIdealFocalSurface.hh"

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

const Double_t kThetaTolerance = 0.1*deg;

ClassImp(ParamOpticalSystem)

//_____________________________________________________________________________
ParamOpticalSystem::ParamOpticalSystem() : fTotalEfficacy(0), fGaussSpread(0),
    fIntegral(0) {
    //
    // Constructor
    //

    ConfigFileParser *pC = Conf();

    string config = Conf()->GetStr("ParamOpticalSystem.fUseConfig");

    // path of the config dir
    string path;
    path = path+ClassType()+'/'+ClassName()+'/'+ config;

    ConfigFileParser *myParser = Config::Get()->GetCF(ClassType(),config,path+".cfg");

    pC = myParser;

    fPos           = EVector(0,0,pC->GetNum("ParamOpticalSystem.fPos.Z")*mm);
    fR             = pC->GetNum("ParamOpticalSystem.fR")*mm;

    fDZup          = pC->GetNum("ParamOpticalSystem.fDZup")*mm;

    fRmax          = pC->GetNum("ParamOpticalSystem.fRmax")*mm;
    fThetaMax      = pC->GetNum("ParamOpticalSystem.fThetaMax")*deg;
    fOutLensBorder = pC->GetNum("ParamOpticalSystem.fOutLensBorder")*mm;

    fDmax          = 0*mm;

    fNbins         = 1000;

    string dummy;

    // set the descriptor of the focal surface descriptor
    dummy = pC->GetStr("ParamOpticalSystem.fFSDescriptor") ;
    if ( dummy == "Ideal" ) {
        // in case of ideal focal surface descriptor build it
        fFSDescriptor = kIdeal;
        dummy = pC->GetStr("ParamOpticalSystem.fIdealFS.fType") ;
        if ( dummy == "KIdeal" )
            fIdealFocalSurface = new KIdealFocalSurface();
        else if ( dummy == "JIdeal" )
            fIdealFocalSurface = new JIdealFocalSurface();
        else
            FatalError("Ideal focal surface type "+dummy+" unknown");
    } else
        FatalError("fFSDescriptor: "+dummy+"\" is not a valid option.");

    string datadir   = pC->GetCfgDir()+'/'+path+'/';  

    // define point spread function
    string type     = pC->GetStr("ParamOpticalSystem.fPsfType");
    if ( type == "gauss" ) { 
        fPsfType    = kGauss;
        fGaussAmpl  = pC->GetNum("ParamOpticalSystem.fGaussAmpl");
        fGaussSigma = pC->GetNum("ParamOpticalSystem.fGaussSigma");
    } else if ( type == "gaussVsTheta" ) {
        fPsfType = kGaussVsTheta;
        // load gaussian shapes from datafile
        string gaussfile = datadir+pC->GetStr("ParamOpticalSystem.fGaussSpread.fFilename");
        fGaussSpread     = new Interpolate( gaussfile ,2);
    } else if ( type == "point" ) {
        fPsfType = kPoint;
    }
    else
        FatalError("fPsfType: "+type+" is not a valid option.");
    
    string toteffile = datadir+pC->GetStr("ParamOpticalSystem.fTotalEfficacy.fFilename");
    fTotalEfficacy   = new Interpolate( toteffile );

    fEPD             = pC->GetNum("ParamOpticalSystem.fEPD");

    fSpreadFun       = new TF2("spreadgauss", "xygaus", -10, 10, -10, 10);

    fSpreadFun->SetNpx(200);
    fSpreadFun->SetNpy(200);

    fSpreadFun->SetParameter(0,1.); // normalization
    fSpreadFun->SetParameter(1,0.); // psf center - X 
    fSpreadFun->SetParameter(3,0.); // psf center - Y
    fSpreadFun->SetParameter(2,1.); // psf rms - X
    fSpreadFun->SetParameter(4,1.); // psf rms - Y

    InitIntegral();
}

//_____________________________________________________________________________
ParamOpticalSystem::~ParamOpticalSystem() {
    //
    // Destructor
    //

    if ( fIdealFocalSurface ) delete fIdealFocalSurface;
    if ( fIntegral )          delete [] fIntegral;
    if ( fGaussSpread)        delete fGaussSpread;
    if ( fSpreadFun )         delete fSpreadFun;

}

//______________________________________________________________________________
Bool_t ParamOpticalSystem::InitIntegral() {
    //
    // Initialize integral tables
    //

    if ( fIntegral ) delete [] fIntegral;

    fIntegral = new Double_t[fNbins+1];

    fIntegral[0] = 0;
    Double_t dz, step, r;
    Double_t dr = (fRmax/fNbins);

    for( Int_t i(1); i < fNbins+1; i++ ) {
        r = dr*i;
        dz = FocalSurfProfile(r)-FocalSurfProfile(r-dr); 
        step = Sqrt( dr*dr+dz*dz );
        fIntegral[i] = fIntegral[i-1]+step;
    }

    // fDmax is the length from center to the border of the focal surface

    fDmax = fIntegral[fNbins];
    
    return kTRUE;
}
    

//______________________________________________________________________________
Photon *ParamOpticalSystem::Transport(Photon *ph) const {
    //
    // Transport photon through the optics
    // 

    Double_t th  = ph->dir.Theta();
    Double_t phi = ph->dir.Phi();
    Double_t wl  = ph->wl;

    // check geometry
    if ( ph->pos.Perp() > fEPD/2. ) {
        ph->fate = kOutOfPupil;
        return 0;
    }

    // check incident angle
    if ( (fThetaMax-th) < kThetaTolerance ) {
        ph->fate = kOutOfFoV;
        return 0;
    }

    // check system absorption
    if ( IsAbsorbed( th, wl ) ) {
        ph->fate = kAbsorbed;
        return 0;
    }

    // distance from focal surface center
    Double_t d = (fDmax/fThetaMax)*th;

    // find the position in cyl coordinates
    Double_t r = DistanceToRadius( d );

    Double_t z = FocalSurfProfile( r );

    // r,z and phi are the coordinates on the focal surface
    TVector3 center(r, 0, z);

    Double_t step, tg;
    step = (fRmax/fNbins)*0.1;

    tg = (FocalSurfProfile(r+step)-FocalSurfProfile(r))/(step);
    TVector3 xAxis(1, 0, tg);
    TVector3 yAxis(0, 1, 0);
    TVector3 zAxis(-tg, 0, 1);
    TVector3 chiefDir;

    xAxis.SetMag(1);
    zAxis.SetMag(1);

    TVector3 pos;
    TVector3 fspos;
    TVector3 dir;

    // add spread
    Double_t dx, dy, a, b;
    if ( Spread( th, wl, dx, dy, a, b) ) {
        // the ph is inside the psf
        fspos = center+dx*xAxis+dy*yAxis; 
        /*
           dir.Rotate(a, yAxis);
           dir.Rotate(b, chiefDir);
         */
        dir.RotateZ( phi );
        fspos.RotateZ( phi );
    } else {
        // the photon is somewhere on the focal surface
        ph->pos[Z] = SecondLensTop();
        ph->fate = kOutOfPsf;
        return 0;
    }

    TRandom *rndm = EsafRandom::Get();
    Double_t rho = (fR-fOutLensBorder)*Sqrt(rndm->Rndm());
    Double_t alpha = rndm->Rndm()*TMath::TwoPi();

    pos.SetXYZ(rho*Cos(alpha),rho*Sin(alpha),SecondLensTop());
    dir = (fspos-pos).Unit();

    ph->time += (ph->pos-pos).Mag()/Clight();
    ph->pos = pos;
    ph->dir = dir;

    return ph;

}

//______________________________________________________________________________
Double_t ParamOpticalSystem::FocalSurfProfile( Double_t r ) const {
    // 
    // Profile of the focal surface Z(r)
    //

    if ( r < 0 || r > fRmax ) return 0;

    Double_t z;
    
    switch ( fFSDescriptor ) {
        case kIdeal:
        default:
            if (!fIdealFocalSurface) FatalError("No ideal focal surface defined");
            z = fIdealFocalSurface->Profile(r);        
    }
    return z;
}
    

//______________________________________________________________________________
Double_t ParamOpticalSystem::DistanceToRadius( Double_t d ) const {
    // 
    // Distance from the fs vertex to r using linear interpolation
    //

    if ( (d < 0) || (d > fDmax) ) return -kHuge;

    Int_t jl = 0;
    Int_t ju = fNbins;

    while ( ju-jl > 1 ) {
        Int_t jm = (jl+ju)/2;
        if ( d >= fIntegral[jm] )
            jl = jm;
        else
            ju = jm;
    }

    Double_t rl = (fRmax/fNbins)*jl;
    Double_t ru = (fRmax/fNbins)*ju;

    Double_t r = rl+(d-fIntegral[jl])/(fIntegral[ju]-fIntegral[jl])*(ru-rl);
    
    return r;
}

//______________________________________________________________________________
Bool_t ParamOpticalSystem::IsAbsorbed( Double_t th, Double_t wl) const {
    //
    // True if photon is absorbed
    //

    TRandom *r = EsafRandom::Get();
    
    Double_t area = Pi()*fEPD*fEPD/4.;
    return (r->Rndm() > fTotalEfficacy->GetValue(th*RadToDeg())/(area*Cos(th)));
    
}

//______________________________________________________________________________
Bool_t ParamOpticalSystem::Spread( Double_t th, Double_t wl, Double_t &x,
    Double_t &y, Double_t &a, Double_t &b ) const {
    //
    // Returns a random point (x,y) and angles (a,b) according psf(th,wl)
    //

    TRandom *r = EsafRandom::Get();
            
    switch ( fPsfType ) {
        case kPoint:
            x = 0;
            y = 0; 
            a = 0;
            b = 0;
            break;
        case kGauss:
            {
                // check if the ray ends up in the psf
                if ( r->Rndm() > fGaussAmpl ) return kFALSE;

                fSpreadFun->GetRandom2(x,y);
                x *= fGaussSigma;
                y *= fGaussSigma;

                break;
            }
        case kGaussVsTheta:
            {
                Double_t theta  = th*RadToDeg();

                // check if the ray ends up in the psf
                if ( r->Rndm() > fGaussSpread->GetValue(theta,kAmpl) ) return kFALSE;

                Double_t sig = fGaussSpread->GetValue(theta, kSigma)*mm;

                fSpreadFun->GetRandom2(x,y);
                x *= sig;
                y *= sig;

                break;
            }
    }

    return kTRUE;

}