// $Id: SphericalIdealFocalSurface.cc,v 1.3 2005/04/14 21:49:58 thea Exp $
// Author: A.Thea 2005/01/16
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: SphericalIdealFocalSurface *
* Package: Optics *
* Coordinator: Alessandro.Thea *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// SphericalIdealFocalSurface
//
// Describes a spherical focal surface.
//
// Config file parameters
// ======================
//
// fR [mm] : Maximum extension of the ideal focal surface from the optical
// axis.
//
// fPos [mm] : z coodinate of the tip of the surface in detector coordinate
// system.
//
// fSphRadius [mm] : radius of the sphere the focal surface is part of.
//
// fPrec [mm] : precision required to claim that the photon has hit the FS.
//
#include "SphericalIdealFocalSurface.hh"
ClassImp(SphericalIdealFocalSurface)
const Int_t kMaxIter=50;
//_____________________________________________________________________________
SphericalIdealFocalSurface::SphericalIdealFocalSurface() {
//
// Constructor
//
fR = Conf()->GetNum("SphericalIdealFocalSurface.fR")*mm;
fPos = EVector(0,0,Conf()->GetNum("SphericalIdealFocalSurface.fPos.Z"))*mm;
fSphRadius = Conf()->GetNum("SphericalIdealFocalSurface.fSphRadius")*mm;
fPrec = Conf()->GetNum("KIdealFocalSurface.fPrec")*mm;
fDZdown = TMath::Abs(Zfs(fR)); //TOCHECK
}
//_____________________________________________________________________________
SphericalIdealFocalSurface::~SphericalIdealFocalSurface() {
//
// Destructor
//
}
//______________________________________________________________________________
Bool_t SphericalIdealFocalSurface::HitPosition(Photon *p)
{
//
// Find whether the photons hits the ideal surface, and in such a case
// calculated the final position
//
// reject photons going downward
if ( p->dir.Z() < 0 ) return kFALSE;
TVector3 in_pos = p->pos-fPos;
EVector ax_pos(0,0,-fDZdown), ax_dir(0,0,1);
EVector bottom_hit = in_pos+p->dir*((-fDZdown-in_pos[Z])/p->dir[Z]);
if (bottom_hit.Perp() > fR) return kFALSE;
//the minimum distance between the photon and the Zfs axis
EVector dummy=ax_dir.Cross(p->dir);
double min_dist=TMath::Abs((bottom_hit-ax_pos).Dot(dummy))/dummy.Mag();
//interaction point on top of the cylinder
//limit inclination of the photon
EVector top_hit;
if (p->dir[Z]/p->dir.Perp() < (fDZup+fDZdown)/(2*fR)) {
// photon hit the lateral surface of the cylinder
top_hit=bottom_hit+p->dir*
(fR*sqrt(1-min_dist*min_dist/(fR*fR))/p->dir.Perp());
} else {
// check if the photon hit the top of the cylinder
top_hit=bottom_hit+p->dir*((fDZup+fDZdown)/p->dir[Z]);
if (top_hit.Perp() > fR) {
top_hit=bottom_hit+p->dir*
(fR*sqrt(1-min_dist*min_dist/(fR*fR))/p->dir.Perp());
}
}
#ifdef DEBUG
cout << "p->dir " << p->dir << endl;
cout << "top_hit= " << top_hit << " bottom_hit=" << bottom_hit << endl;
#endif
EVector step=(top_hit-bottom_hit)*.5;
dummy=bottom_hit+step;
for(int i=0; i < kMaxIter ; i++) {
#ifdef DEBUG
cout << "i=" << i << " dummy[Z]=" << dummy[Z] << " Zfs(dummy.Perp())=" << Zfs(dummy.Perp()) << endl;
#endif
if (TMath::Abs(Zfs(dummy.Perp())-dummy[Z]) < fPrec) {
break;
}
step*=.5;
if (Zfs(dummy.Perp()) > dummy[Z]) dummy+=step;
else dummy-=step;
}
if ( dummy.Perp() < fR ) {
p->posOnIfs=dummy+fPos;
p->hitIfs=kTRUE;
return kTRUE;
}
return kFALSE;
}
//______________________________________________________________________________
Double_t SphericalIdealFocalSurface::Zfs(Double_t r) {
//
// Zfs returns the z value of the focal surface as a function of radius in
// local coordinates
//
Double_t x = r/fSphRadius;
return -fSphRadius*(1-TMath::Sqrt(1.-x*x));
}
//______________________________________________________________________________
Double_t SphericalIdealFocalSurface::Profile(Double_t r) {
//
// Profile returns the z value of the focal surface as a function of radius
// in detector coordinates
//
return Zfs(r)+fPos.Z();
}