// $Id: KIdealFocalSurface.cc,v 1.14 2005/04/15 18:54:41 thea Exp $
// Author: J.Watts 2003/09/21
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: KIdealFocalSurface *
* Package: Optics *
* Coordinator: Alessandro.Thea *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// KIdealFocalSurface
//
// KIdealFocalSurface is the optimized focal surface for the KOpticalSystem
// Its shape is described by the radial profile:
//
// r^2
// z(r) = --------------------------- + A*r^4 + B*r^6 + C*r^8 + D*r^10
// (1+(sqrt( (1+k)*(r/rho)^2) ))
//
// 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.
//
// fRho, fK, fA, fB, fC, fD : focal surface parameters.
//
// fPrec [mm] : precision required to claim that the photon has hit the FS.
//
#include "KIdealFocalSurface.hh"
#include <TVector3.h>
#include "Photon.hh"
#include "TMath.h"
ClassImp(KIdealFocalSurface)
const Int_t kMaxIter=50;
using namespace TMath;
//______________________________________________________________________________
KIdealFocalSurface::KIdealFocalSurface() {
//
// Constructor
//
fR = Conf()->GetNum("KIdealFocalSurface.fR")*mm;
fPos = TVector3(0,0,Conf()->GetNum("KIdealFocalSurface.fPos.Z"))*mm;
fRho = Conf()->GetNum("KIdealFocalSurface.rho");
fK = Conf()->GetNum("KIdealFocalSurface.k");
fA = Conf()->GetNum("KIdealFocalSurface.a");
fB = Conf()->GetNum("KIdealFocalSurface.b");
fC = Conf()->GetNum("KIdealFocalSurface.c");
fD = Conf()->GetNum("KIdealFocalSurface.d");
fPrec = Conf()->GetNum("KIdealFocalSurface.precision")*mm;
fDZdown = Abs(Zfs(fR)); //TOCHECK
}
//______________________________________________________________________________
KIdealFocalSurface::~KIdealFocalSurface() {
//
// Destructor
//
}
//______________________________________________________________________________
Bool_t KIdealFocalSurface::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;
TVector3 ax_pos(0,0,-fDZdown), ax_dir(0,0,1);
TVector3 bottomHit = in_pos+p->dir*((-fDZdown-in_pos[Z])/p->dir[Z]);
if (bottomHit.Perp() > fR)
return kFALSE;
//the minimum distance between the photon and the Zfs axis
TVector3 dummy=ax_dir.Cross(p->dir);
//double min_dist=Abs((bottomHit-ax_pos).Dot(dummy))/dummy.Mag();
//interaction point on top of the cylinder
//limit inclination of the photon
Double_t dt = NextInt( bottomHit+fPos, p->dir);
TVector3 topHit = bottomHit+p->dir*dt;
/*
#ifdef DEBUG
cout << "p->dir " << p->dir << endl;
cout << "topHit= " << topHit << " bottomHit=" << bottomHit << endl;
#endif
*/
TVector3 step = (topHit-bottomHit)*.5;
dummy=bottomHit+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 (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 )
return kFALSE;
p->posOnIfs=dummy+fPos;
p->hitIfs=kTRUE;
return kTRUE;
}
//______________________________________________________________________________
Double_t KIdealFocalSurface::Zfs(Double_t r) {
//
// Zfs returns the z value of the focal surface as a function of radius in
// local coordinates
//
Double_t rr=r*r;
return rr/((1+sqrt(1-((1+fK)*rr/(fRho*fRho))))*fRho)+
fA*rr*rr+fB*rr*rr*rr+fC*rr*rr*rr*rr+
fD*rr*rr*rr*rr*rr;
}
//______________________________________________________________________________
Double_t KIdealFocalSurface::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();
}