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

// $Id: LIdealFocalSurface.cc 2824 2009-03-16 14:36:44Z naumov $
// Author: Dmitry V.Naumov  2007/11/16

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: NIdealFocalSurface                                                   *
 *  Package: Optics                                                          *
 *  Coordinator: Alessandro.Thea                                             *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// NIdealFocalSurface
//
//   NIdealFocalSurface is the optimized focal surface for the MOpticalSystem
//   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 coordinate 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 "NIdealFocalSurface.hh"
#include "Ntrace_optics.hh"
#include <TVector3.h>
#include "Photon.hh"
#include "TMath.h"
#include <cmath>

ClassImp(NIdealFocalSurface)

const Int_t kMaxIter=50;

using namespace sou;
using namespace TMath;
using namespace NTraceLens;



//______________________________________________________________________________
NIdealFocalSurface::NIdealFocalSurface() {
    //
    // Constructor
    //
    const char* lens_dir = Conf()->GetStr("NOpticalSystem.lens_dir").data();
    const char* tel_par = Conf()->GetStr("NOpticalSystem.tel_par").data();
    ReadTelParm(lens_dir,tel_par);
    fR   = fOpticalSystemParam->r_fs;
    Double_t FS_Offset = fOpticalSystemParam->surface[0].Sz+fOpticalSystemParam->FS_OFFSET;
    fPos = TVector3(0,0,FS_Offset*mm);
    fDZdown = Zfs(0)-Zfs(fR*mm)-1.*mm;
}

//______________________________________________________________________________
NIdealFocalSurface::~NIdealFocalSurface() {
    //
    // Destructor
    //

}

//______________________________________________________________________________
Bool_t NIdealFocalSurface::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;
    TVector3 ax_pos(0,0,-fDZdown), ax_dir(0,0,1);

    TVector3 bottomHit = in_pos+p->dir*((fPos[Z]-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);
    //cout <<" dummy[Z]=" << dummy[Z] << endl;
    //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, p->dir);
    //cout << "dt " << dt << endl;
    TVector3 topHit = bottomHit+p->dir*dt;
   // cout << "p->dir " << p->dir << endl;
   // cout << "topHit= " << topHit << " bottomHit=" << bottomHit << endl;
/*
#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;
    }
     double r_cut = 950.;
    if ( dummy.Perp() > fR || fabs(dummy[Y])>r_cut)
        return kFALSE;



    p->posOnIfs=dummy;
    //cout <<" dummy[Z]=" << dummy[Z] <<" dummy[Y]=" << dummy[Y]  <<" dummy[X]=" << dummy[X] << endl;
    p->hitIfs=kTRUE;

    return kTRUE;
}


//______________________________________________________________________________
Double_t NIdealFocalSurface::Zfs(Double_t r) {
    //
    // Zfs returns the z value of the focal surface as a function of radius in
    // local coordinates
    //

    return FocalSurface(r,fOpticalSystemParam);
}

//______________________________________________________________________________
Double_t NIdealFocalSurface::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();
}