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

// $Id: SphericalIdealFocalSurface.cc 2789 2008-02-11 07:26:05Z naumov $
// Author: A.Thea   2005/01/16

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

//_____________________________________________________________________________
//
// SphericalIdealFocalSurface
//
//   Describes a spherical focal surface.
//
//   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.
//
//   fSphRadius [mm] : radius of the sphere the focal surface is part of.
//
//   fPrec [mm] : precision required to claim that the photon has hit the FS.
//

#include "SphericalIdealFocalSurface.hh"
#include <TMath.h>

using namespace sou;

ClassImp(SphericalIdealFocalSurface)

const Int_t kMaxIter=50;

//_____________________________________________________________________________
SphericalIdealFocalSurface::SphericalIdealFocalSurface() {
    //
    // Constructor
    //

        fR         = Conf()->GetNum("SphericalIdealFocalSurface.fR")*mm;
        fPos       = EVector(0,0,Conf()->GetNum("SphericalIdealFocalSurface.fPos.Z"))*mm;
        fSphRadius = Conf()->GetNum("SphericalIdealFocalSurface.fSphRadius")*mm;
    fPrec      = Conf()->GetNum("KIdealFocalSurface.fPrec")*mm;
    fDZdown    = TMath::Abs(Zfs(fR)); //TOCHECK

}

//_____________________________________________________________________________
SphericalIdealFocalSurface::~SphericalIdealFocalSurface() {
    //
    // Destructor
    //
}

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

    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 min_dist=TMath::Abs((bottom_hit-ax_pos).Dot(dummy))/dummy.Mag();

    //interaction point on top of the cylinder
    //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]) < fPrec) {
            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 SphericalIdealFocalSurface::Zfs(Double_t r) {
    // 
    // Zfs returns the z value of the focal surface as a function of radius in 
    // local coordinates
    // 
    Double_t x = r/fSphRadius;
        return -fSphRadius*(1-TMath::Sqrt(1.-x*x)); 
}

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