source: JEM-EUSO/esaf_lal/tags/v1_r0/esaf/packages/simulation/detector/optics/src/JIdealFocalSurface.cc @ 117

// ESAF : Euso Simulation and Analysis Framework
// $Id: JIdealFocalSurface.cc 2136 2005-10-02 14:18:08Z thea $
// A.Thea created Apr,  3 2003

#include "JIdealFocalSurface.hh"
#include "EVector.hh"
#include "Photon.hh"
#include "TMath.h"

using namespace sou;

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();
}