// $Id: DetectorPhotonTransporter.cc,v 1.15 2005/04/17 15:50:08 thea Exp $
// Author: D.Demarco, M.Pallavicini
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: DetectorPhotontransporter *
* Package: Optics *
* Coordinator: Alessandro.Thea *
* *
*****************************************************************************/
//______________________________________________________________________________
//
// DetectorPhotonTransporter
// =========================
//
// Abstarct base class providing the interface of a generic photon
// transporter. Every element in the detector that is able to interact
// with photons is a DetectorPhotonTransporter.
// The DetectorPhotonTransporters are represented as cylinders with position
// and height (fDZdown+fDZup).
//
#include "DetectorPhotonTransporter.hh"
#include "Etypes.hh"
ClassImp(DetectorPhotonTransporter)
using namespace TMath;
//______________________________________________________________________________
DetectorPhotonTransporter::DetectorPhotonTransporter(): EsafConfigurable(),
fPos(0,0,0), fDZdown(0), fDZup(0), fR(0) {
//
// Constructor
//
// fLastCylInt.first=0;
// fLastCylInt.second=new Double_t[2];
}
//______________________________________________________________________________
DetectorPhotonTransporter::~DetectorPhotonTransporter(){
//
// Destructor
//
// delete [] fLastCylInt.second;
}
//______________________________________________________________________________
Bool_t DetectorPhotonTransporter::IsInside( Photon *p ) const {
//
// True if the photon is inside the boundaries of the transporter
//
return ((p->pos[Z]-Bottom()) > -kTolerance && (p->pos[Z]-Top()) < kTolerance &&
(p->pos.Perp()-Radius()) < kTolerance);
}
//______________________________________________________________________________
Double_t DetectorPhotonTransporter::CylinderIntersection( const TVector3 &pos,
const TVector3 &dir, Bool_t zero) const {
//
// Calculate the next intersection of p with this transporter boundaries
//
return CylinderIntersection( pos, dir, fR, fPos[Z]+fDZup, fPos[Z]-fDZdown, zero);
}
//______________________________________________________________________________
Double_t DetectorPhotonTransporter::CylinderIntersection( const Photon *p,
Bool_t zero) const {
//
// Calculate the next intersection of p with this transporter boundaries
//
return CylinderIntersection( p, fR, fPos[Z]+fDZup, fPos[Z]-fDZdown, zero );
}
//______________________________________________________________________________
Double_t DetectorPhotonTransporter::CylinderIntersection( const Photon *p,
Double_t radius, Double_t zup, Double_t zdown, Bool_t zero ) const {
//
// Finds, if exists, the next interaction of p over cyl. In any case saves
// the distances between p->pos and the two int points in fLastCylInt
//
return CylinderIntersection(p->pos, p->dir, radius, zup, zdown, zero );
}
//______________________________________________________________________________
Double_t DetectorPhotonTransporter::CylinderIntersection( const TVector3 &pos,
const TVector3 &dir, Double_t radius, Double_t zup, Double_t zdown,
Bool_t zero ) const {
//
// Finds, if exists, the next interaction of p over cyl. In any case saves
// the distances between p->pos and the two int points in fLastCylInt
//
Double_t dist[4] = {-kHuge,-kHuge,-kHuge,-kHuge};
// array of intersections, 0, 1 basis, 3, 4 lateral surface
// int with the basis
if ( Abs(dir[Z]) > kTolerance ) {
dist[0] = (zup-pos[Z])/dir[Z];
dist[1] = (-zdown-pos[Z])/dir[Z];
}
Double_t a, b, c;
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]-radius*radius;
Double_t delta = (b*b)-4*a*c;
if ( a > 0 )
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 threshold = zero ? -kTolerance : kTolerance;
Double_t dt = MaxElement(4,dist);
for(Int_t i(0); i<4; i++)
if ( dist[i] > threshold && dist[i] < dt ) dt = dist[i];
// inbetween +/- ktolerance dt is basically 0
if ( zero && dt > -kTolerance && dt < kTolerance) dt = 0;
//MsgForm(EsafMsg::Info,"dt %3ft,d(%8f,%8f,%8f,%8f)",dt,dist[0],dist[1],dist[2],dist[3]);
return dt;
/***************************************************************************
*
* Old code, kept for the time being as reference
*
* *************************************************************************
// intesection point
TVector3 intPoint(0,0,0), uDir;
uDir = dir.Unit();
Double_t a, b, c, dup, ddown, dummy[2];
fLastCylInt.first=0;
fLastCylInt.second[0]=fLastCylInt.second[1]=0;
// distances to the upper and lower planes
if ( TMath::Abs(uDir[Z]) >= kTolerance ){
dup = (zUp-pos[Z])/uDir[Z];
ddown = (zDown-pos[Z])/uDir[Z];
}
if((TMath::Abs(uDir.Theta()) < kTolerance) ||
(TMath::Abs(uDir.Theta() - TMath::Pi()) < kTolerance)) {
// photon directed along Z: special case
// doesn't hit the cyl.
if(pos.Perp()-radius > kTolerance) {
return -1*kHuge;
}
// intersection with the bases
if(uDir[Z] > 0 && pos[Z] < zDown ||
uDir[Z] < 0 && pos[Z] > zDown && pos[Z] < zUp ){
// lower base
fLastCylInt.second[0]= ddown;
fLastCylInt.second[1]= dup;
}
else if(uDir[Z] < 0 && pos[Z] > zUp ||
uDir[Z] > 0 && pos[Z] > zDown && pos[Z] < zUp ){
// upper base
fLastCylInt.second[0]= dup;
fLastCylInt.second[1]= ddown;
} else {
MsgForm(EsafMsg::Warning,"CylinderIntersection(): photon lost");
return -2*kHuge;
}
return fLastCylInt.second[0];
}
// fill the coefficients of the equation
a=uDir[X]*uDir[X]+uDir[Y]*uDir[Y];
b=2*(uDir[X]*pos[X]+uDir[Y]*pos[Y]);
c=pos[X]*pos[X]+pos[Y]*pos[Y]-radius*radius;
// find the first intersection with the cylinder if it exists
pair<int, double* > res;
res=findRoots(a,b,c);
if (res.first == 0) {
// photon doesn't cross the cylinder
// return something;
return -3*kHuge;
} else if (res.first == 1) {
// one solution, tangent photon
intPoint=pos+uDir*res.second[0];
// check the sol to be inside cylinder
if ( intPoint[Z] > zUp || intPoint[Z] < zDown )
return -4*kHuge;
fLastCylInt.first=1;
fLastCylInt.second[0]=fLastCylInt.second[1]=res.second[0];
return fLastCylInt.second[0];
} else if ( TMath::Abs(uDir[Z]) < kTolerance ){
// two crossings but no intersection with the bases
// horizontal photon
intPoint=pos+uDir*res.second[0];
// check the sol to be inside cylinder
if ( intPoint[Z] > zUp || intPoint[Z] < zDown )
return -5*kHuge;
}
// 4 candidates left, just 2 of them on the borders of the cyl; in fact 3 if
// the photons crosses the intersection between the side surface and one of
// the bases
fLastCylInt.first=2;
// distances of the 4 intpoints from pos
// 0: upper base plane
// 1: lower base plane
// 2: 1st lateral surface
// 3: 2nd lateral surface
Double_t intDist[4];
intDist[0]=dup;
intDist[1]=ddown;
intDist[2]=res.second[0];
intDist[3]=res.second[1];
Int_t nTrueInts(0);
Bool_t onSide, onBase;
for(Int_t i(0); i<4; i++ ){
intPoint = pos+uDir*intDist[i];
onSide = TMath::Abs(intPoint.Perp()-radius) < kTolerance &&
(intPoint[Z] < zUp && intPoint[Z] > zDown);
onBase = (TMath::Abs(intPoint[Z]-zUp) < kTolerance ||
TMath::Abs(intPoint[Z]-zDown) < kTolerance) &&
intPoint.Perp() <= radius;
if ( onSide || onBase )
// if intDist[i] < kTolerance, save just 0
dummy[nTrueInts++]=( TMath::Abs(intDist[i]) > kTolerance ? intDist[i] : 0 );
}
if ( nTrueInts == 3 ) {
// check if two of the ints are the same point (one with the bases and
// one with the side)
if ( TMath::Abs(intDist[0]-intDist[2]) < kTolerance ||
TMath::Abs(intDist[0]-intDist[3]) < kTolerance ||
TMath::Abs(intDist[1]-intDist[2]) < kTolerance ||
TMath::Abs(intDist[1]-intDist[3]) < kTolerance) {
nTrueInts--;
}
}
if ( nTrueInts == 2) {
// last thing to do, sort the sols
if ( dummy[0]*dummy[1] < 0 || dummy[0]+dummy[1] < 0 ){
// opposite sign sols or sols < 0
fLastCylInt.second[0]=TMath::Max(dummy[0],dummy[1]);
fLastCylInt.second[1]=TMath::Min(dummy[0],dummy[1]);
} else {
// sols >= 0
fLastCylInt.second[0]=TMath::Min(dummy[0],dummy[1]);
fLastCylInt.second[1]=TMath::Max(dummy[0],dummy[1]);
}
if ( fLastCylInt.second[0] < 0 )
return -6*kHuge;
else
return fLastCylInt.second[0];
} else if ( nTrueInts==0 ) {
// non of the sols was on the cyl
return -7*kHuge;
} else {
// dump to screen all the intersections
Msg(EsafMsg::Warning) << "pos" << pos << " dir" << uDir << endl;
for(Int_t i(0); i<4; i++ ){
MsgForm(EsafMsg::Warning, "Intersection %d:",i);
intPoint = pos+uDir*intDist[i];
MsgForm(EsafMsg::Warning,"intPoint.PerP() = %f intPoint[Z] = %f",
intPoint.Perp(),intPoint[Z]);
MsgForm(EsafMsg::Warning,"Radius = %f Zup = %f Zdown = %f",
radius, zUp, zDown);
onSide = TMath::Abs(intPoint.Perp()-radius) < kTolerance &&
(intPoint[Z] < zUp && intPoint[Z] > zDown);
onBase = (TMath::Abs(intPoint[Z]-zUp) < kTolerance ||
TMath::Abs(intPoint[Z]-zDown) < kTolerance) &&
intPoint.Perp()-radius <= kTolerance;
MsgForm(EsafMsg::Warning," side = %f base = %f", onSide, onBase);
}
Msg(EsafMsg::Warning) << " nTrueInts = " << nTrueInts << MsgDispatch;
MsgForm(EsafMsg::Panic,"CylinderIntersection(): something gone wrong!");
return -8*kHuge;
}
*******************************************************************************/
}