// $Id: DetectorTransportManager.cc,v 1.72 2005/05/08 21:02:44 thea Exp $
// Author: D.DeMarco
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: DetectorTransportManager *
* Package: Optics *
* Coordinator: Alessandro.Thea *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// DetectorTransportManager
// ========================
//
// The raytracer of the inner detector.
//
#include <math.h>
#include <algorithm>
#include "DetectorTransportManager.hh"
#include "utils.hh"
#include "EEvent.hh"
#include "Etypes.hh"
#include "PhotonsOnPupil.hh"
#include "ParentPhoton.hh"
#include "EDetectorPhotonAdder.hh"
ClassImp(DetectorTransportManager)
using namespace TMath;
//______________________________________________________________________________
DetectorTransportManager::DetectorTransportManager() : DetectorPhotonTransporter(),
fDebug(kFALSE), fEnableFocalPlane(kTRUE) {
//
// Constructor
//
// size of the cyl that includes detector and the baffle
fDZup = 10000*mm;
fDZdown = 7000*mm;
fR = 2000*mm;
fPos = EVector(0,0,0);
fInnerRadius=Conf()->GetNum("DetectorTransportManager.fInnerRadius")*mm;
of = OpticsFactory::Get();
of->Build();
fOptics = of->GetOpticalSystem();
fFocalPlane = of->GetFocalPlane();
fIdealFocalSurface=of->GetIdealFocalSurface();
fWalls = of->GetWalls();
fBaffle = of->GetBaffle();
if ( fFocalPlane->Radius() > fInnerRadius ||
fOptics->Radius() > fInnerRadius )
FatalError("FocalPlane radius or Optics radius is bigger than fInnerRadius. Aborting");
else if ( fBaffle->Radius() < fInnerRadius )
FatalError("Baffle cannot be smaller than fInnerRadius. Aborting");
else
fOuterRadius = fBaffle->Radius();
fMaxIterations=(int)Conf()->GetNum("DetectorTransportManager.fMaxIterations");
}
//______________________________________________________________________________
void DetectorTransportManager::Go( PhotonsOnPupil* pupil ) {
//
// Get and trace all the photons in PhotonsOnPupil
//
if ( !pupil ) {
MsgForm(EsafMsg::Warning,"PhotonsOnPupil is NULL. Photon tracing disabled.");
return;
}
EEvent *ev = EEvent::GetCurrent();
Photon *p=0;
Photon *p_old=0;
int nTracked=0, nTotal(pupil->GetNphotons());
fTimeFirstPhoton = kHuge;
fTimeLastPhoton = -kHuge;
if (!fEnableFocalPlane) Msg(EsafMsg::Info) << "FocalPlane disabled" << MsgDispatch;
Msg(EsafMsg::Info) << "Tracking Photons: " << MsgFlush;
while(1) {
// get a new photon
// if p==0 and p_old not, photon is arrived to destination
// use the copy p_old to draw graphics
if(p_old && ev ) {
EDetectorPhotonAdder a(p_old,0,kFALSE);
ev->Fill(a);
}
// get a new photon
p = pupil->Get();
// build an empty EDetPhoton in the event
if ( !p ) break;
if ( ev ) {
EDetectorPhotonAdder a(p,0,kTRUE);
ev->Fill(a);
}
p_old=p; // keep a copy for later use
// keep track of first and last photon times
if ( p->time > fTimeLastPhoton )
fTimeLastPhoton = p->time;
if ( p->time < fTimeFirstPhoton )
fTimeFirstPhoton = p->time;
// update display
if ( nTracked*50 % nTotal == 0 || (nTracked-1)*50 % nTotal > nTracked*50 % nTotal )
if ( nTracked*10 % nTotal == 0 || (nTracked-1)*10 % nTotal > nTracked*10 % nTotal ){
Msg(EsafMsg::Info) << Form("%d",nTracked*100/nTotal) << '%' << MsgFlush;
}
else
Msg(EsafMsg::Info) << "." << MsgFlush;
nTracked++;
Transport(p);
}
Msg(EsafMsg::Info) << " -done" << MsgDispatch;
Msg(EsafMsg::Info) << "Tracking complete. " << nTracked <<" traced." << MsgDispatch;
}
//______________________________________________________________________________
Photon* DetectorTransportManager::Transport( Photon* p) const {
//
// Trace one photon in the detector.
//
EVector intPoint;
if ( fDebug ) Msg(EsafMsg::Debug) << "@ new photon" << MsgDispatch;
// check if the photon is in front of the detector
// should be improved
if ( p->pos.Perp()-fBaffle->Radius() > -kTolerance ) {
p->history=kNotInside;
return 0;
}
while( p ){
if((p->iterations)++ > fMaxIterations) {
MsgForm(EsafMsg::Warning,"Photon reached maximum number of iterations. killed.");
Msg(EsafMsg::Warning) << " p->id = " << p->parent->Id() << MsgDispatch;
Msg(EsafMsg::Warning) << " p->pos = " << p->pos << MsgDispatch;
Msg(EsafMsg::Warning) << " p->pos.Perp() = " << p->pos.Perp() << MsgDispatch;
Msg(EsafMsg::Warning) << " p->dir = " << p->dir << MsgDispatch;
Msg(EsafMsg::Warning) << " p->history = " << p->history << MsgDispatch;
Msg(EsafMsg::Warning) << " p->iterations = " << p->iterations << MsgDispatch;
p = 0;
continue;
}
if(p->fate != 0) {
FatalError("Shouldn't happen p->fate != 0 while photon is still alive.");
}
// reality check on the photon position
// FIXME these are debug conditions
Bool_t isInEuso;
isInEuso = p->pos[Z]-fBaffle->Bottom() > -kTolerance &&
p->pos[Z]-fFocalPlane->Top() < kTolerance &&
p->pos.Perp()-fInnerRadius < kTolerance;
if ( !(isInEuso || fBaffle->IsInside( p )) ) {
Msg(EsafMsg::Warning) << "This photon seems to be where it shouldn't:" << MsgDispatch;
Msg(EsafMsg::Warning) << " p->pos = " << p->pos << MsgDispatch;
Msg(EsafMsg::Warning) << " p->pos.Perp() = " << p->pos.Perp() << MsgDispatch;
Msg(EsafMsg::Warning) << " p->dir = " << p->dir << MsgDispatch;
Msg(EsafMsg::Warning) << " p->history = " << p->history << MsgDispatch;
Msg(EsafMsg::Warning) << " p->iterations = " << p->iterations << MsgDispatch;
Msg(EsafMsg::Warning) << " isInEuso = " << isInEuso << " isInsideBaffle = " << fBaffle->IsInside( p ) << MsgDispatch;
p->history=kNotInside;
p = 0;
continue;
}
// photon before the optics
if ( p->pos[Z] <= fOptics->FirstLensTop() ){
if ( fDebug ) Msg(EsafMsg::Debug) << "--first zone" << MsgDispatch;
if ( Abs(p->pos[Z]-fBaffle->Bottom()) < kTolerance && p->dir[Z] < 0) {
if ( fDebug ) Msg(EsafMsg::Debug) << " ~ back to the earth" << MsgDispatch;
p->history=kEarth;
p = 0;
continue;
} else
if ( p->pos.Perp() > fInnerRadius ) {
// photon outside the inner cylinder and inside the baffle
if ( fDebug ) Msg(EsafMsg::Debug) << " - OUTSIDE" << MsgDispatch;
// check if gets throught inner cylinder
Double_t DL = CylinderIntersection( p , fInnerRadius, fBaffle->Top(), fBaffle->Bottom(), kTRUE );
if ( DL > 0 ) {
// bring the photon to the borders of the inner cyl
if ( fDebug ) Msg(EsafMsg::Debug) << " ~ going in" << MsgDispatch;
intPoint = p->pos+p->dir.Unit()*DL;
p->time += DL/c_light;
p->pos = intPoint;
} else {
if ( fDebug ) Msg(EsafMsg::Debug) << " ~ baffle" << MsgDispatch;
p->history = kBaffle;
p = fBaffle->Transport( p );
}
} else {//if ( (p->pos.Perp()-fInnerRadius) < kTolerance )
if ( Abs((p->pos-fPos).Perp()-fInnerRadius) < kTolerance){
if ( fDebug ) Msg(EsafMsg::Debug) << " - BORDERS ";
if ((p->pos-fPos).XYvector()*p->dir.XYvector() <0) {
// inward, nothing to do
if ( fDebug ) Msg(EsafMsg::Debug) << "inward" << MsgDispatch;
} else {
// going outward
if ( fDebug ) Msg(EsafMsg::Debug) << "outward" << MsgDispatch;
if (p->pos[Z] < fBaffle->Top() ){
p->history = kBaffle;
if ( fDebug ) Msg(EsafMsg::Debug) << " ~ baffle" << MsgDispatch;
p = fBaffle->Transport( p );
} else {
p->history = kWalls;
if ( fDebug ) Msg(EsafMsg::Debug) << " ~ walls" << MsgDispatch;
p = fWalls->Transport( p );
}
continue;
}
}
// photon inside the inner cylinder
if ( fDebug ) Msg(EsafMsg::Debug) << " - INSIDE" << MsgDispatch;
intPoint=InteractionPoint(p);
if ( Abs((p->pos-intPoint).Mag()) < kTolerance ){
// if (fLastCylInt.second[1] < 0)
Msg(EsafMsg::Warning) << "TraceOnePhoton: got confused in bringing ph in the inner cyl" << MsgDispatch;
// try with the second option
// intPoint = p->pos+p->dir.Unit()*fLastCylInt.second[1];
}
if ( intPoint[Z] >= fOptics->Bottom() ){
if ( fDebug ) Msg(EsafMsg::Debug) << " ~ optics" << MsgDispatch;
Double_t DZ = fOptics->Bottom()-p->pos[Z];
intPoint=NewInteractionPoint(p,DZ);
p->time+=(intPoint-p->pos).Mag()/c_light;
p->pos=intPoint;
p->history=kOptics;
p=fOptics->Transport( p );
} else if ( intPoint[Z] < fBaffle->Bottom()) {
if ( fDebug ) Msg(EsafMsg::Debug) << " ~ back to the earth" << MsgDispatch;
Double_t DL = fBaffle->CylinderIntersection( p );
p->time+=DL/c_light;
p->pos+=p->dir.Unit()*DL;
p->history=kEarth;
p=0;
} else {
if ( fDebug ) Msg(EsafMsg::Debug) << " ~out of the inner cyl" << MsgDispatch;
p->time+=(intPoint-p->pos).Mag()/c_light;
p->pos=intPoint;
}
}
} else if (p->pos[Z] > Max(fOptics->FirstLensTop(), fBaffle->Top())
&& p->pos[Z] < fOptics->SecondLensBottom()) {
Msg(EsafMsg::Debug) << "--optics forbidden zone" << MsgDispatch;
Msg(EsafMsg::Debug) << "- p->pos="<< p->pos << " - p->dir=" << p->dir << MsgDispatch;
Msg(EsafMsg::Debug) << "- p->pos.Perp()=" << p->pos.Perp() << MsgDispatch;
p->history=kNotInside;
p=0;
} else {
// calculate possible impact point of the photon on the cylinder
if ( fDebug ) Msg(EsafMsg::Debug) << "--second zone" << MsgDispatch;
intPoint=InteractionPoint(p);
// photon is in the second zone:
// between the optics and the focal plane
if(intPoint[Z] > fFocalPlane->Bottom()) {
// photon on focal plane
if ( fDebug ) Msg(EsafMsg::Debug) << " ~focalplane" << MsgDispatch;
fIdealFocalSurface->HitPosition(p);
p->history=kFocalPlane;
if (!fEnableFocalPlane) {
p=0;
continue;
}
// DZ distance along Z from photon to focal plane
double DZ = fFocalPlane->Bottom() - p->pos[Z];
// calculate impact point on the focal plane
intPoint=NewInteractionPoint(p, DZ);
// calculate interaction of photon with focal plane
p->time+=(p->pos - intPoint).Mag()/c_light;
p->pos=intPoint;
p=fFocalPlane->Transport(p);
} else if(intPoint[Z] < fOptics->Top() && p->dir[Z] < 0) {
// photon on optics
if ( fDebug ) Msg(EsafMsg::Debug) << " ~optics from above" << MsgDispatch;
if ( p->pos[Z] > fOptics->Top() ) {
// DZ distance along Z from photon to optics
double DZ = fOptics->Top() - p->pos[Z];
// calculate impact point on the optics surface
intPoint=NewInteractionPoint(p, DZ);
// calculate interaction of photon with optics
p->time+=(p->pos - intPoint).Mag()/c_light;
p->pos=intPoint;
p->history=kOptics;
p=fOptics->Transport(p);
} else if (p->history == kOptics ){
p->history=kWalls;
p=fWalls->Transport(p);
} else {
p->history=kOptics;
p=fOptics->Transport(p);
}
} else {
// photon on wall: reflection or absorption
if ( fDebug ) Msg(EsafMsg::Debug) << " ~walls" << MsgDispatch;
p->time+=(p->pos - intPoint).Mag()/c_light;
p->pos=intPoint;
p->history=kWalls;
p=fWalls->Transport(p);
}
}
}
return 0;
}
//______________________________________________________________________________
Double_t DetectorTransportManager::EusoCylIntersection( Photon *p ) const {
//
// Returns the distance (if any) of p from the next intersection with the
// EUSO shell.
//
return CylinderIntersection( p , fInnerRadius, fPos[Z]+fDZup, fPos[Z]-fDZdown);
}
//______________________________________________________________________________
EVector DetectorTransportManager::InteractionPoint(Photon *p, Double_t radius) const {
//
// Returns the next interaction point of p with a cylinder with a given
// radius. The cylinder axis is the optical axis.
//
// intesection point
EVector intPoint(p->pos);
if ( radius < 0 ) radius = fInnerRadius;
Double_t DL;
DL = CylinderIntersection( p , radius, fPos[Z]+fDZup, fPos[Z]-fDZdown, kTRUE);
if ( DL < 0 ){
Msg(EsafMsg::Warning) << "InteractionPoint" << MsgDispatch;
// MsgForm(EsafMsg::Warning,"DL = %f - LastCyInt[0] = %f - LastCylInt[1] = %f",
// DL, fLastCylInt.second[0], fLastCylInt.second[1]);
Msg(EsafMsg::Warning) << "p->pos: " << p->pos << " - p->pos.Perp(): " << p->pos.Perp() << MsgDispatch;
MsgForm(EsafMsg::Warning,"R: - %f pos*dir: %f - pos*dir:",
radius, p->pos.XYvector()*p->dir.XYvector());
p->Dump();
FatalError("should't happen DL<0. Photon killed");
}
intPoint=p->pos+p->dir.Unit()*DL;
return intPoint;
}