// $Id: ConicBaffle.cc,v 1.2 2005/04/15 18:54:40 thea Exp $
// Author: Alessandro Thea 2005/04/08
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: ConicBaffle *
* Package: <packagename> *
* Coordinator: <coordinator> *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// ConicBaffle
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// fTopRadius [mm]: radius of the optics the baffle is attached to.
// -Valid options: <available options>
//
// fAlpha [deg]: baffle's field of view.
// -Valid options: 0 to 60 degrees.
//
#include "ConicBaffle.hh"
#include "Etypes.hh"
#include "PhysicalConstants.hh"
using namespace TMath;
ClassImp(ConicBaffle)
//_____________________________________________________________________________
ConicBaffle::ConicBaffle() {
//
// Constructor
//
fRTop = Conf()->GetNum("ConicBaffle.fTopRadius");
fAlpha = Conf()->GetNum("ConicBaffle.fAlpha");
if ( fAlpha < 0 || fAlpha > 60 )
FatalError("fAlpha must be larger than 0 and smaller than 60 deg.");
fAlpha *= DegToRad();
fDZdown = Conf()->GetNum("ConicBaffle.DZdown");
fPos = EVector(0,0,Conf()->GetNum("ConicBaffle.fPos.Z"));
fR = fRTop + fDZdown*Tan(fAlpha);
if ( Abs(TopRadius()-BottomRadius()) < kTolerance ) {
fCotAlpha = 0;
fVPosZ = 0;
} else {
fVPosZ = -fDZdown+(fDZdown/(BottomRadius()-TopRadius()))*BottomRadius();
fCotAlpha = BottomRadius()/(-fDZdown-fVPosZ);
}
}
//_____________________________________________________________________________
ConicBaffle::~ConicBaffle() {
//
// Destructor
//
}
//______________________________________________________________________________
Double_t ConicBaffle::NextInteraction( Photon* ph ) const{
//
// finds the next interaction of the photon inside the cylinder
//
Bool_t fgDebug = kFALSE;
Double_t dist[4] = {-kHuge,-kHuge,-kHuge,-kHuge};
// array of intersections, 0, 1 basis, 3, 4 lateral surface
if (fgDebug) Printf("pos(%.8f,%.8f,%.8f) dir(%.8f,%.8f,%.8f)",
ph->pos[X],ph->pos[Y],ph->pos[Z],ph->dir[X],ph->dir[Y],ph->dir[Z]);
// improve read
TVector3& dir = ph->dir;
TVector3& pos = ph->pos;
// int with the basis
if ( Abs(dir[Z]) > kTolerance ) {
dist[0] = -pos[Z]/dir[Z]; // top basis
dist[1] = (-fDZdown-pos[Z])/dir[Z]; // bottom basis
}
Double_t a, b, c;
if ( fCotAlpha == 0 ) {
a = dir[X]*dir[X]+dir[Y]*dir[Y];
b = 2*(dir[X]*pos[X]+dir[Y]*pos[Y]);
c = pos[X]*pos[X]+pos[Y]*pos[Y]-TopRadius()*TopRadius();
} else {
// vertex
a = dir.Perp2()-fCotAlpha*fCotAlpha*dir[Z]*dir[Z];
b = 2*(dir[X]*pos[X]+dir[Y]*pos[Y]-fCotAlpha*fCotAlpha*dir[Z]*(pos[Z]-fVPosZ));
c = pos[X]*pos[X]+pos[Y]*pos[Y]-fCotAlpha*fCotAlpha*(pos[Z]-fVPosZ)*(pos[Z]-fVPosZ);
}
Double_t delta = (b*b)-4*a*c;
if (fgDebug) Printf("a=%.7e, b=%.7e, c=%.7e, delta=%.7e",a,b,c,delta);
if ( Abs(a) > kTolerance )
if ( Abs(delta) < kTolerance ) {
dist[2] = -b/(2*a);
} else if ( delta >= kTolerance ) {
dist[2] = (-b-Sqrt(delta))/(2*a);
dist[3] = (-b+Sqrt(delta))/(2*a);
}
Double_t dt = MaxElement(4,dist);
for(Int_t i(0); i<4; i++)
if ( dist[i] > kTolerance && dist[i] < dt ) dt = dist[i];
if (fgDebug) Printf("d(%e,%e,%e,%e) and dt=%.3f",dist[0],dist[1],dist[2],dist[3],dt);
return dt;
}
//______________________________________________________________________________
Photon* ConicBaffle::Transport( Photon* ph ) const {
//
// Transport
//
// to local
ph->pos -= fPos;
Double_t dt = NextInteraction( ph );
ph->pos += ph->dir*dt;
ph->time+= dt/c_light;
ph->pos += fPos;
if ( Abs(ph->pos[Z] - Bottom()) < kTolerance ||
Abs(ph->pos[Z] - Top()) < kTolerance )
return ph;
else
return 0;
}