// ESAF : Euso Simulation and Analysis Framework
// $Id: TestLightToEuso.cc,v 1.2 2005/04/14 12:43:25 moreggia Exp $
// M. Pallavicini created Jul, 2 2002
#include "TestLightToEuso.hh"
#include "EsafRandom.hh"
#include "ListPhotonsOnPupil.hh"
#include "utils.hh"
#include "SystemOfUnits.hh"
#include "MCTruth.hh"
#include "EEventTruthAdder.hh"
#include "EEvent.hh"
#include "EConst.hh"
#include "EsafSpectrum.hh"
#include "TMath.h"
ClassImp(TestLightToEuso)
using namespace TMath;
//______________________________________________________________________________
TestLightToEuso::TestLightToEuso() : LightToEuso("TEST"), fTruth(NULL) {
//
// Constructor
//
fPhOnPupil = new ListPhotonsOnPupil();//(&fPhotons);
}
//______________________________________________________________________________
TestLightToEuso::~TestLightToEuso() {
//
// Destructor
//
Reset();
delete fPhOnPupil;
if ( fTruth ) delete fTruth;
}
//______________________________________________________________________________
void TestLightToEuso::Reset() {
//
// Reset internal list of fPhotons
//
fPhOnPupil->Clear();
}
//______________________________________________________________________________
PhotonsOnPupil *TestLightToEuso::Get(const DetectorGeometry* dg) {
//
// Returns the list of fPhotons for the pupil.
// The type of photon distribution is determined
// by TestLightToEuso.Type
//
Reset();
string type = Conf()->GetStr("TestLightToEuso.Type");
Double_t th1 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Theta1");
Double_t th2 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Theta2");
Double_t ph1 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Phi1");
Double_t ph2 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Phi2");
Int_t nPhotons = (Int_t)Conf()->GetNum("TestLightToEuso.Photons");
Double_t tm = Conf()->GetNum("TestLightToEuso.Duration");
TRandom* rndm = EsafRandom::Get();
if ( type == "SPOT" ) { // spot in fixed position
return MakeSpot(th1,ph1,nPhotons,tm);
}
else if ( type == "RANDSPOT" ) { // spot in random position within Phi1,2 Theta1,2
Double_t cmax = Cos(th2);
Double_t cmin = Cos(th1);
if ( th1 > th2 ) {
Msg(EsafMsg::Info) << "In TestLightToEuso Theta1 and Theta2 will be switchedn";
cmax = cmin;
cmin = Cos(th2);
}
Double_t diff = cmax - cmin;
Double_t xc = cmin + rndm->Rndm()*diff;
Double_t th = acos(xc);
if ( ph2 > ph1 )
diff = ph2 - ph1;
else {
Msg(EsafMsg::Info) << "In TestLightToEuso Phi1 and Phi2 will be switchedn";
diff = ph1 - ph2;
Double_t temp = ph2;
ph2 = ph1;
ph1 = temp;
}
Double_t ph = rndm->Rndm()*diff + ph1;
return MakeSpot(th,ph,nPhotons,tm);
}
else if ( type == "EXTSPOT" ) {
return MakeExtSpot(th1,ph1,TMath::RadToDeg()*ph2,nPhotons,tm);
}
else if ( type == "CIRCLE" ) {
return MakeCircle(ph1,ph2,th1,nPhotons,tm);
}
else if ( type == "EXTCIRCLE" ) {
return MakeExtCircle(ph1,ph2,th1,th2,nPhotons,tm);
}
else if ( type == "DIFFUSE" ) {
return MakeDiffuse(th1,th2,TMath::RadToDeg()*ph1,nPhotons,tm);
}
else if ( type == "CFLUXDIFFUSE" ) {
return MakeCFluxDiffuse(th1,th2,TMath::RadToDeg()*ph1,nPhotons,tm);
}
else if ( type == "RADIUS" ) {
return MakeRadius(ph1,th1,th2,nPhotons,tm);
}
else if ( type == "TRACK" ) {
return MakeTrack(ph1,ph2,th1,th2,nPhotons,tm);
}
else if ( type == "SHOWERTRACK" ) {
return MakeShowerTrack();
}
throw runtime_error("Invalid parameter TestLightToEuso.Type = " + type );
return (PhotonsOnPupil*)0;
}
//______________________________________________________________________________
PhotonsOnPupil *TestLightToEuso::MakeSpot(const Double_t& th, const Double_t& ph, Int_t nPhotons,
const Double_t& tm) {
//
// Produce light in a single spot
//
Msg(EsafMsg::Info) << "SPOT CALLED th="<<th<<" ph="<<ph<<" Photons="<<nPhotons<< MsgDispatch;
Double_t wl = 357.*nm;
Double_t x[3] = {0.,0.,0.};
Double_t y[3] = {0.,0.,0.};
Double_t t;
TRandom* rndm = EsafRandom::Get();
for(Int_t i=0; i<nPhotons; i++) {
t = rndm->Rndm()*tm*microsecond;
ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
fPhOnPupil->Add( p );
}
return fPhOnPupil;
}
//______________________________________________________________________________
PhotonsOnPupil *TestLightToEuso::MakeExtSpot(const Double_t& th1, const Double_t& ph1, const Double_t rad,
Int_t nPhotons, const Double_t& tm) {
//
// Produce light in an extensive spot
//
Msg(EsafMsg::Info) << "EXTSPOT CALLED th1="<<th1<<" ph1="<<ph1<<" radius="<<rad<<" Photons="<<nPhotons<< MsgDispatch;
Double_t wl = 357.*nm;
Double_t x[3] = {0.,0.,0.};
Double_t y[3] = {0.,0.,0.};
Double_t t;
TRandom* rndm = EsafRandom::Get();
for(Int_t i=0; i<nPhotons; i++) {
t = rndm->Rndm()*tm*microsecond;
Double_t radius = Sqrt(rndm->Rndm())*rad;
Double_t th2 = rndm->Rndm()*TMath::TwoPi();
y[0] = radius*Cos(th2)*mm;
y[1] = radius*Sin(th2)*mm;
ParentPhoton *p = new ParentPhoton(i,x,y,th1,ph1,wl,t,Fluorescence);
fPhOnPupil->Add( p );
}
return fPhOnPupil;
}
//______________________________________________________________________________
PhotonsOnPupil *TestLightToEuso::MakeCircle(const Double_t& ph1, const Double_t& ph2, const Double_t& th,
Int_t nPhotons, const Double_t& tm) {
//
// Produce light in a circle
//
TRandom* rndm = EsafRandom::Get();
Msg(EsafMsg::Info) << "CIRCLE CALLED with ph1=" << ph1 << " ph2=" << ph2 << " th=" << th << MsgDispatch;
Double_t wl = 357.*nm;
Double_t x[3] = {0.,0.,0.};
Double_t y[3] = {0.,0.,0.};
Double_t t;
for(Int_t i=0; i<nPhotons; i++) {
t = rndm->Rndm()*tm*microsecond;
Double_t ph = ph1 + rndm->Rndm()*(ph2-ph1);
ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
fPhOnPupil->Add( p );
}
return fPhOnPupil;
}
//______________________________________________________________________________
PhotonsOnPupil *TestLightToEuso::MakeExtCircle(const Double_t& ph1, const Double_t& ph2, const Double_t& th1,
const Double_t th2, Int_t nPhotons, const Double_t& tm) {
//
// Produce light in a extensive circle (fPhotons with th1 < th < th2)
//
TRandom *rndm = EsafRandom::Get();
Msg(EsafMsg::Info) << "EXTCIRCLE CALLED with ph1=" << ph1 << " ph2=" << ph2 <<
" th1=" << th1 << " th2=" << th2 << MsgDispatch;
Double_t wl = 357.*nm;
Double_t x[3] = {0.,0.,0.};
Double_t y[3] = {0.,0.,0.};
Double_t t,theta1,theta2;
if ( th1 > th2 ) {
theta1=th2;
theta2=th1;
} else {
theta1=th1;
theta2=th2;
}
for(Int_t i=0; i<nPhotons; i++) {
t = rndm->Rndm()*tm*microsecond;
Double_t th = theta1+(theta2-theta1)*Sqrt((theta2-theta1)/(pi/6)*rndm->Rndm());
Double_t ph = ph1 + rndm->Rndm()*(ph2-ph1);
ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
fPhOnPupil->Add( p );
}
return fPhOnPupil;
}
//______________________________________________________________________________
PhotonsOnPupil *TestLightToEuso::MakeDiffuse(const Double_t& th1, const Double_t& th2,
const Double_t rad, Int_t nPhotons, const Double_t& tm) {
TRandom *rndm = EsafRandom::Get();
Msg(EsafMsg::Info) << "DIFFUSE CALLED with th1=" << th1 << " th2=" << th2 <<
" rad=" << rad << MsgDispatch;
Double_t wl = 357.*nm;
Double_t x[3] = {0.,0.,0.};
Double_t y[3] = {0.,0.,0.};
Double_t t(0),theta1,theta2;
if ( th1 > th2 ) {
theta1=th2;
theta2=th1;
} else {
theta1=th1;
theta2=th2;
}
Double_t c1 = Cos(theta1);
Double_t c2 = Cos(theta2);
for(Int_t i=0; i<nPhotons; i++) {
// Double_t th = theta1+(theta2-theta1)*Sqrt((theta2-theta1)/(pi/6)*rndm->Rndm());
Double_t th = ACos(rndm->Uniform(c1,c2));
Double_t ph = TMath::TwoPi()*rndm->Rndm();
Double_t rho = Sqrt(rndm->Rndm())*rad;
Double_t alpha = rndm->Rndm()*TMath::TwoPi();
y[0] = rho*Cos(alpha)*mm;
y[1] = rho*Sin(alpha)*mm;
ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
fPhOnPupil->Add( p );
}
return fPhOnPupil;
}
//______________________________________________________________________________
PhotonsOnPupil *TestLightToEuso::MakeCFluxDiffuse(const Double_t& th1, const Double_t& th2,
const Double_t rad, Int_t nPhotons, const Double_t& tm) {
TRandom *rndm = EsafRandom::Get();
Msg(EsafMsg::Info) << "C FLUX-DIFFUSE CALLED with th1=" << th1 << " th2=" << th2 <<
" rad=" << rad << MsgDispatch;
Double_t wl = 357.*nm;
Double_t x[3] = {0.,0.,0.};
Double_t y[3] = {0.,0.,0.};
Double_t t(0),theta1(0),theta2(0);
if ( th1 > th2 ) {
theta1=th2;
theta2=th1;
} else {
theta1=th1;
theta2=th2;
}
Double_t c1 = Cos(2*theta1);
Double_t c2 = Cos(2*theta2);
for(Int_t i(0); i<nPhotons; i++) {
Double_t th = 0.5*ACos(rndm->Uniform(c1,c2));
Double_t ph = TMath::TwoPi()*rndm->Rndm();
Double_t rho = Sqrt(rndm->Rndm())*rad;
Double_t alpha = rndm->Rndm()*TMath::TwoPi();
y[0] = rho*Cos(alpha)*mm;
y[1] = rho*Sin(alpha)*mm;
ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
fPhOnPupil->Add( p );
}
return fPhOnPupil;
}
//______________________________________________________________________________
PhotonsOnPupil *TestLightToEuso::MakeRadius(const Double_t& ph, const Double_t& th1, const Double_t& th2,
Int_t nPhotons, const Double_t& tm) {
//
// Produce light along a focal surface radius.
//
TRandom* rndm = EsafRandom::Get();
Msg(EsafMsg::Info) << "RADIUS CALLED" << MsgDispatch;
Double_t wl = 357.*nm;
Double_t x[3] = {0.,0.,0.};
Double_t y[3] = {0.,0.,0.};
Double_t t;
for(Int_t i=0; i<nPhotons; i++) {
t = rndm->Rndm()*tm*microsecond;
Double_t th = th1 + rndm->Rndm()*(th2-th1); // uniform in theta, not costheta
ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
fPhOnPupil->Add( p );
}
return fPhOnPupil;
}
//______________________________________________________________________________
PhotonsOnPupil *TestLightToEuso::MakeTrack(const Double_t& ph1, const Double_t& ph2, const Double_t& th1,
const Double_t& th2, Int_t nPhotons, const Double_t& tm) {
//
// Produce light along a track
//
//TRandom* rndm = EsafRandom::Get();
return 0;
}
//______________________________________________________________________________
PhotonsOnPupil *TestLightToEuso::MakeShowerTrack() {
// produce light along a track
Msg(EsafMsg::Info) << "SHOWERTRACK CALLED" << MsgDispatch;
TRandom* rndm = EsafRandom::Get();
Double_t iX,iY, iZ;
Double_t tLength = Conf()->GetNum("TestLightToEuso.ShowerTrack.fTrackLength")*km;
Double_t tTheta = Conf()->GetNum("TestLightToEuso.ShowerTrack.fTheta")*TMath::DegToRad();
Double_t tPhi = Conf()->GetNum("TestLightToEuso.ShowerTrack.fPhi")*TMath::DegToRad();
string sInitPos = Conf()->GetStr("TestLightToEuso.ShowerTrack.fPosInit");
if ( sInitPos == "fixed") {
// fixed int point
iX = Conf()->GetNum("TestLightToEuso.ShowerTrack.fXInit")*km;
iY = Conf()->GetNum("TestLightToEuso.ShowerTrack.fYInit")*km;
iZ = Conf()->GetNum("TestLightToEuso.ShowerTrack.fZInit")*km;
} else if (sInitPos == "random") {
// random interacion point
Double_t rInitMax = Conf()->GetNum("TestLightToEuso.ShowerTrack.fRInitMax")*km;
Double_t radInit = TMath::Sqrt(rndm->Rndm())*rInitMax;
Double_t phiInit = rndm->Rndm()*TMath::TwoPi();
iX = radInit*TMath::Cos(phiInit);
iY = radInit*TMath::Sin(phiInit);
// Z init between 0 and 50 km
iZ = rndm->Rndm()*Conf()->GetNum("TestLightToEuso.ShowerTrack.fZInitMax")*km;
} else if (sInitPos=="zfixed") {
// Z coordinate of the interacion point is kept fixed
Double_t rInitMax = Conf()->GetNum("TestLightToEuso.ShowerTrack.fRInitMax")*km;
Double_t radInit = TMath::Sqrt(rndm->Rndm())*rInitMax;
Double_t phiInit = rndm->Rndm()*TMath::TwoPi();
Msg(EsafMsg::Info) << rInitMax << " " << phiInit << " " << radInit << MsgDispatch;
iX = radInit*TMath::Cos(phiInit);
iY = radInit*TMath::Sin(phiInit);
iZ = Conf()->GetNum("TestLightToEuso.ShowerTrack.fZInit")*km;
} else {
Msg(EsafMsg::Info) << "TestLightToEuso.cfg:"
<< " wrong value in TestLightToEuso.ShowerTrack.fPosInit:"
<< sInitPos << MsgDispatch;
throw runtime_error("ESAF Error: Wrong value in TestLightToEuso.cfg");
}
Int_t nPhotons = (Int_t)Conf()->GetNum("TestLightToEuso.ShowerTrack.fPhotons");
TVector3 eusoPos(0,0,Config::Get()->GetCF("General", "Euso")->GetNum("Euso.fAltitude")*km);
Double_t eusoRad = Config::Get()->GetCF("General","Euso")->GetNum("Euso.fRadius")*mm;
TVector3 initPos(iX, iY, iZ);
TVector3 diffPos;
diffPos.SetMagThetaPhi(tLength,tTheta,tPhi);
TVector3 endPos = diffPos + initPos;
TVector3 dir = (endPos-initPos).Unit();
TVector3 dL = (endPos-initPos)*(1./nPhotons);
fTruth = new MCTruth;
fTruth->SetFirstInt(initPos-eusoPos, 0);
fTruth->SetThetaPhi(dir.Theta(), dir.Phi());
fTruth->SetShowerMax(0.5*(endPos-initPos)+initPos-eusoPos, 0);
if ( EEvent::GetCurrent() ) {
EEventTruthAdder a(fTruth);
EEvent::GetCurrent()->Fill( a );
}
EsafSpectrum *shSpec = new EsafSpectrum("nitro10km");
Double_t wl;
Double_t t;
Double_t x[3] = {0.,0.,0.};
Double_t y[3] = {0.,0.,0.};
Double_t l;
for(Int_t i(0); i<nPhotons; i++) {
// get a random number, gaussian distributed between 0 and 1 sigma
for(l = rndm->Gaus(0.5,1/4.); l < 0 || l > 1;l = rndm->Gaus(0.5,1/4.) );
TVector3 shPos = initPos+diffPos*l;
TVector3 inDir = -((shPos-eusoPos).Unit());
Double_t radius = Sqrt(rndm->Rndm())*eusoRad;
Double_t phi = rndm->Rndm()*TMath::TwoPi();
x[0] = (shPos-eusoPos)[0];
x[1] = (shPos-eusoPos)[1];
x[2] = (shPos-eusoPos)[2];
y[0] = radius*Cos(phi)*mm;
y[1] = radius*Sin(phi)*mm;
y[2] = 0*mm;
wl = shSpec->GetLambda();
t = (l*tLength+(shPos-eusoPos).Mag())/EConst::C();
ParentPhoton *p = new ParentPhoton(i,x,y,inDir.Theta(),inDir.Phi(),wl,t,Fluorescence);
fPhOnPupil->Add( p );
}
return fPhOnPupil;
}