Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

JIdealFocalSurface - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: JIdealFocalSurface.cc,v 1.11 2005/10/02 14:17:55 thea Exp $
// 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();
}

About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005