// $Id: DetectorTransportManager.cc,v 1.83 2005/10/31 14:43:56 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 <iomanip>
#include "DetectorTransportManager.hh"
#include "utils.hh"
#include "EEvent.hh"
#include "Etypes.hh"
#include "PhotonsOnPupil.hh"
#include "ParentPhoton.hh"
#include "EDetectorPhotonAdder.hh"
#include "EConst.hh"
ClassImp(DetectorTransportManager)
using namespace TMath;
using namespace sou;
using namespace EConst;
//______________________________________________________________________________
DetectorTransportManager::DetectorTransportManager() : DetectorPhotonTransporter(),
fDebug(kFALSE), fScaleFactor(1.), 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");
Reset();
}
//______________________________________________________________________________
void DetectorTransportManager::SetScaleFactor( Double_t scale ) {
//
// Set the detector size scaling factor
//
if ( scale > 0 )
fScaleFactor = scale;
}
//______________________________________________________________________________
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());
Int_t progress = 0;
fTimeFirstPhoton = kHuge;
fTimeLastPhoton = -kHuge;
if (!fEnableFocalPlane) Msg(EsafMsg::Info) << "FocalPlane disabled" << MsgDispatch;
while(1) {
if ( p_old ) {
//
// Update the statistics
//
switch ( p_old->history ) {
case kNotInside:
fNout++;
break;
case kBaffle:
fNbaffle++;
break;
case kOptics:
fNoptics++;
break;
case kWalls:
fNwalls++;
break;
case kFocalPlane:
fNfocalplane++;
break;
case kEarth:
fNreflected++;
break;
default:
fNlost++;
}
//
if( ev ) {
EDetectorPhotonAdder a(p_old,0,kFALSE);
ev->Fill(a);
}
}
// get a new photon
p = pupil->GetPhoton();
// build an empty EDetPhoton in the event
if ( !p ) break;
// scale ph position in respect of the pupil center
// simulate greater or smaller experiments
if ( Abs(fScaleFactor-1) > kTolerance )
p->pos *= 1./fScaleFactor;
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
nTracked++;
if (100.*nTracked/nTotal >= progress) {
Msg(EsafMsg::Info).SetProgress(progress);
Msg(EsafMsg::Info) << "Tracking Photons: " << MsgCount;
progress+=10;
}
Transport(p);
}
// fill a summary of photons history
if ( ev ) {
TArrayI photonHistory(7);
photonHistory.AddAt(fNlost, kUndefined);
photonHistory.AddAt(fNout, kNotInside);
photonHistory.AddAt(fNbaffle, kBaffle);
photonHistory.AddAt(fNoptics, kOptics);
photonHistory.AddAt(fNfocalplane, kFocalPlane);
photonHistory.AddAt(fNwalls, kWalls);
photonHistory.AddAt(fNreflected, kEarth);
EDetectorPhotonAdder a(photonHistory);
ev->Fill(a);
}
// debug staistics
Msg(EsafMsg::Debug) << "TransportManager Photon's transport statistics" << MsgDispatch;
Msg(EsafMsg::Debug) << "NotInside : " << fNout << MsgDispatch;
Msg(EsafMsg::Debug) << "Baffle : " << fNbaffle << MsgDispatch;
Msg(EsafMsg::Debug) << "Optics : " << fNoptics << MsgDispatch;
Msg(EsafMsg::Debug) << "Walls : " << fNwalls << MsgDispatch;
Msg(EsafMsg::Debug) << "FocalPlane : " << fNfocalplane << MsgDispatch;
Msg(EsafMsg::Debug) << "Reflected : " << fNreflected << MsgDispatch;
Msg(EsafMsg::Debug) << "Lost : " << fNlost << MsgDispatch;
Msg(EsafMsg::Info) << "Tracking complete. " << nTracked <<" traced." << MsgDispatch;
}
//______________________________________________________________________________
void DetectorTransportManager::Reset() {
//
// Reset the manager.
//
fNout = 0;
fNbaffle = 0;
fNoptics = 0;
fNwalls = 0;
fNfocalplane = 0;
fNreflected = 0;
fNlost = 0;
}
//______________________________________________________________________________
Photon* DetectorTransportManager::Transport( Photon* p) const {
//
// Trace one photon in the detector.
//
// check if the photon has been already rejected by the pupil
if ( p->history == kNotInside ) return 0;
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/Clight();
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()/Clight();
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/Clight();
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()/Clight();
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()/Clight();
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()/Clight();
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()/Clight();
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;
}