// $Id: ParamOpticalSystem.cc,v 1.27 2005/05/06 17:34:27 thea Exp $
// 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 pourpose 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 distancs.
//
// 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]: Maximmum 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.filename [string]: two column file containing the
// total efficacy as a function of the incident angle.
//
// fAnglesOfIncidence.filename [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.filename [string]: three column file containing RMS of the
// gaussian psf and focal surface incident photons cone aperture.
//
//
//
#include "ParamOpticalSystem.hh"
#include "Config.hh"
#include "Etypes.hh"
#include "EsafRandom.hh"
#include "IdealFocalSurface.hh"
#include "KIdealFocalSurface.hh"
#include "JIdealFocalSurface.hh"
using namespace TMath;
const Double_t kThetaTolerance = 0.1*deg;
ClassImp(ParamOpticalSystem)
//_____________________________________________________________________________
ParamOpticalSystem::ParamOpticalSystem() : 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.type") ;
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.");
// define point spread function
string type = pC->GetStr("ParamOpticalSystem.fPsfType");
if ( type == "gauss" )
fPsfType = kGaussian;
else if ( type == "point" )
fPsfType = kPointLike;
else
FatalError("fPsfType: "+type+" is not a valid option.");
string datadir = pC->GetCfgDir()+'/'+path+'/';
string gaussfile = datadir+pC->GetStr("ParamOpticalSystem.fGaussSpread.filename");
string toteffile = datadir+pC->GetStr("ParamOpticalSystem.fTotalEfficacy.filename");
fGaussSpread = new Interpolate( gaussfile ,2);
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->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 >> 1;
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)
//
switch ( fPsfType ) {
case kPointLike:
x = 0;
y = 0;
a = 0;
b = 0;
break;
case kGaussian:
Double_t theta = th*RadToDeg();
TRandom *r = EsafRandom::Get();
// 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;
}