// ESAF : Euso Simulation and Analysis Framework
// class KOpticalSystem
// $Id: KOpticalSystem.cc,v 1.35 2005/04/16 15:24:08 thea Exp $
// J.Watts created , Sep 21 2003
//
#include "SystemOfUnits.hh"
#include "PhysicalConstants.hh"
#include "Ktrace_optF1v1.hh"
#include "KEUSO_optF1v1.hh"
#include "KOpticalSystem.hh"
#include "Config.hh"
#include "EsafRandom.hh"
#include "EEvent.hh"
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 */
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/c_light;
}
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/c_light;
}
#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;
}