source: JEM-EUSO/esaf_cc_at_lal/packages/simulation/detector/optics/src/KIdealFocalSurface.cc @ 114

// $Id: KIdealFocalSurface.cc 2136 2005-10-02 14:18:08Z thea $
// 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 sou;
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();
}