// ESAF : Euso Simulation and Analysis Framework
// $Id: JIdealFocalSurface.cc,v 1.10 2005/04/14 21:49:58 thea Exp $
// A.Thea created Apr, 3 2003
#include "JIdealFocalSurface.hh"
#include "EVector.hh"
#include "Photon.hh"
#include "TMath.h"
ClassImp(JIdealFocalSurface)
//______________________________________________________________________________
JIdealFocalSurface::JIdealFocalSurface() {
//
// Constructor
//
fR = Conf()->GetNum("JIdealFocalSurface.fR")*mm;
fPos = EVector(0,0,Conf()->GetNum("JIdealFocalSurface.fPos.Z"))*mm;
_rho = Conf()->GetNum("JIdealFocalSurface.rho");
_k = Conf()->GetNum("JIdealFocalSurface.k");
_a = Conf()->GetNum("JIdealFocalSurface.a");
_b = Conf()->GetNum("JIdealFocalSurface.b");
_c = Conf()->GetNum("JIdealFocalSurface.c");
_d = Conf()->GetNum("JIdealFocalSurface.d");
_prec = Conf()->GetNum("JIdealFocalSurface.precision")*mm;
fDZdown = TMath::Abs(Zfs(fR)); //TOCHECK
#ifdef DEBUG
cout << "JIdealFocalSurface created" << endl;
cout << "rho=" << _rho << "\n" <<
"k=" << _k << "\n" <<
"a=" << _a << "\n" <<
"b=" << _b << "\n" <<
"c=" << _c << "\n" <<
"d=" << _d << "\n" <<
"prec=" << _prec << "\n" <<
"fDZdown" << fDZdown << "\n" << endl;
cout << "Zfs(0)=" << Zfs(0) << "\n" <<
"Zfs(100)=" << Zfs(100) << "\n" <<
"Zfs(200)=" << Zfs(200) << "\n" <<
"Zfs(300)=" << Zfs(300) << "\n" <<
"Zfs(400)=" << Zfs(400) << "\n" <<
"Zfs(500)=" << Zfs(500) << "\n" <<
"Zfs(600)=" << Zfs(600) << "\n" <<
"Zfs(700)=" << Zfs(700) << "\n" <<
"Zfs(800)=" << Zfs(800) << "\n" <<
"Zfs(900)=" << Zfs(900) << "\n" <<
"Zfs(1000)=" << Zfs(1000) << "\n" <<
"Zfs(1100)=" << Zfs(1100) << "\n" <<
"Zfs(1200)=" << Zfs(1200) << "\n" << endl;
#endif
}
//______________________________________________________________________________
JIdealFocalSurface::~JIdealFocalSurface() {
//
// Destructor
//
}
const Int_t kMaxIter=50;
//______________________________________________________________________________
Bool_t JIdealFocalSurface::HitPosition(Photon *p)
{
TVector3 in_pos = p->pos-fPos;
EVector ax_pos(0,0,-fDZdown), ax_dir(0,0,1);
// check if the ph is in the IFS boundaries
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_t min_dist=TMath::Abs((bottom_hit-ax_pos).Dot(dummy))/dummy.Mag();
//interaction point on top of the cilinder
//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]) < _prec) {
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 JIdealFocalSurface::Zfs(Double_t r) {
//
// Zfs returns the z value of the focal surface as a function of radius
//
Double_t rr=r*r;
return rr/((1+sqrt(1-((1+_k)*rr/(_rho*_rho))))*_rho)+
_a*rr*rr+_b*rr*rr*rr+_c*rr*rr*rr*rr+
_d*rr*rr*rr*rr*rr;
}
//______________________________________________________________________________
Double_t JIdealFocalSurface::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();
}