// ESAF : Euso Simulation and Analysis Framework
// $Id: BunchOfPhotons.cc,v 1.37 2005/10/21 13:31:30 moreggia Exp $
// Anne Stutz created Dec, 1 2003
#include "BunchOfPhotons.hh"
#include <stdexcept>
#include "EsafRandom.hh"
#include "EConst.hh"
ClassImp(BunchOfPhotons)
using namespace sou;
using namespace TMath;
using namespace EConst;
UInt_t BunchOfPhotons::fCounter = 0;
//____________________________________________________________________________________________________
BunchOfPhotons::BunchOfPhotons(): fId(0), fParent(0) {
//
// ctor
//
Msg(EsafMsg::Panic) <<"Default ctor should not be used" <<MsgDispatch;
}
//____________________________________________________________________________________________________
BunchOfPhotons::BunchOfPhotons(const BunchOfPhotons& b) : fParent(0){
//
// copy ctor
//
fId = b.GetId();
fNbins = b.GetNbins();
fType = b.GetType();
fWeight = b.GetWeight();
fDate = b.GetDate();
fTof = b.GetTof();
fPos = b.GetPos();
fShowerPos = b.GetShowerPos();
fDir = b.GetDir();
fNextImpact = b.GetNextImpact();
fFinalImpact = b.GetFinalImpact();
b.GetWlSpectrum().Copy(fWlSpectrum);
fFate = b.GetFate();
fParent = new ParentBunch(*b.GetParent());
//b.GetTransSpectrum().Copy(fTransSpectrum);
//b.GetRaylSpectrum().Copy(fRaylSpectrum);
//b.GetOzoneSpectrum().Copy(fOzoneSpectrum);
//b.GetMieSpectrum().Copy(fMieSpectrum);
}
//____________________________________________________________________________________________________
BunchOfPhotons::BunchOfPhotons(const Double_t weight,const Double_t yield,
const EarthVector& showerposi, const EarthVector& showerposf,
const Double_t datei, const Double_t datef, const EsafSpectrum& spectrum,
const EarthVector& dir, PhotonType type) {
//
// ctor with creation of a Parent Bunch
//
fCounter++;
fId = fCounter;
fType = type;
fWeight = weight;
fDate = (datei+datef)/2;
fTof = 0;
fShowerPos = ( showerposi+showerposf )*0.5;
fPos = fShowerPos;
fDir = dir.Unit();
fNextImpact = fShowerPos;
fFinalImpact = fShowerPos;
spectrum.Copy(fWlSpectrum);
fNbins = spectrum.GetAxis().GetNbins();
fFate = 0;
fParent = new ParentBunch(weight,yield,showerposi,showerposf,datei,datef,spectrum,fDir);
//spectrum.Copy(fTransSpectrum);
//spectrum.Copy(fRaylSpectrum);
//spectrum.Copy(fOzoneSpectrum);
//spectrum.Copy(fMieSpectrum);
// following needed for convolution operations
//fTransSpectrum.Clear();
//fRaylSpectrum.Clear();
//fOzoneSpectrum.Clear();
//fMieSpectrum.Clear();
}
//____________________________________________________________________________________________________
BunchOfPhotons::~BunchOfPhotons() {
//
// dtor
//
SafeDelete(fParent);
}
//____________________________________________________________________________________________________
EarthVector BunchOfPhotons::RandomPosInShower() const {
//
// Returns a random position at creation, using lateral and longitudinal position distributions
// If no distributions, returns bunch showerpos
//
EarthVector rtn(1,0,0);
EarthVector longb = fParent->GetShowerPosf()- fParent->GetShowerPosi();
EarthVector axis = fParent->GetAxis();
TRandom* rndm = EsafRandom::Get();
rtn = fParent->GetShowerPosi() + rndm->Rndm()*longb;
// if working with histos
if (fParent->GetHistoLateral()) {
EarthVector v = axis.Orthogonal(); // vector orthogonal to direction
v.Rotate(rndm->Rndm() * 2 * TMath::Pi(),axis); // random rotation around direction
v.SetMag(fParent->GetHistoLateral()->GetRandom()*m); // lateral distribution in meter if NKG2 //TOFIX if NKG1
rtn += v ;
return rtn;
}
// if working with functions
else if(fParent->GetLateralDist()) {
EarthVector v = axis.Orthogonal(); // vector orthogonal to direction
v.Rotate(rndm->Rndm() * 2 * TMath::Pi(),axis); // random rotation around direction
v.SetMag(fParent->GetLateralDist()->GetRandom()*m); // lateral distribution in meter if NKG2 //TOFIX if NKG1
rtn += v ;
return rtn;
}
// if no distributions
return GetShowerPos();
}
//____________________________________________________________________________________________________
EarthVector BunchOfPhotons::PosRandomAngCorrec(const EarthVector& showerpos,const EarthVector& currentpos) const {
//
// Returns a corrected position at any propagation stage, using :
// - angular distributions at creation
// - the corrected SinglePhoton showerpos given in argument
// - the current SinglePhoton position given in argument
// NB : all corrections (in position and time) must be done before calling the present method
//
// correction due to corrected showerpos
EarthVector increment(currentpos - showerpos);
EarthVector rtn(currentpos);
// should not occur, but..
if(increment.Mag() == 0) return rtn;
// correction due to angular distributions
if (GetParent()->GetHistoAngular()) {
EarthVector v = increment.Orthogonal();
Double_t theta = fParent->GetHistoAngular()->GetRandom();
TRandom* rndm = EsafRandom::Get();
v.SetMag(increment.Mag() * TMath::Tan(theta));
v.Rotate(rndm->Rndm() * 2 * TMath::Pi(),increment);
rtn += v;
return rtn;
}
else if (GetParent()->GetAngularDist()) {
EarthVector v = increment.Orthogonal();
Double_t theta = fParent->GetAngularDist()->GetRandom();
TRandom* rndm = EsafRandom::Get();
v.SetMag(increment.Mag() * TMath::Tan(theta));
v.Rotate(rndm->Rndm() * 2 * TMath::Pi(),increment);
rtn += v;
return rtn;
}
// if no distibutions, return current position
return rtn;
}
//____________________________________________________________________________________________________
EarthVector BunchOfPhotons::RandomDirection() const {
//
// Returns a random direction at creation using the angular distributions
//
Double_t theta = 0;
EarthVector rtn(fDir);
TRandom* rndm = EsafRandom::Get();
// isotropic if fluo
if(fType == Fluo) {
Double_t x,y,z;
rndm->Sphere(x,y,z,1.);
rtn.SetXYZ(x,y,z);
return rtn;
}
// use angular distribution if cerenkov
if(fDir.Mag() == 0) Msg(EsafMsg::Warning) <<"<RandomDirection> not-Fluo bunch with a NULL direction -> SHOULD NOT" << MsgDispatch;
if(GetParent()->GetHistoAngular()) theta = fParent->GetHistoAngular()->GetRandom();
else if(GetParent()->GetAngularDist()) theta = fParent->GetAngularDist()->GetRandom();
// if no distribution, return bunch mean direction
else return fDir;
EarthVector v = fDir.Orthogonal();
v.SetMag(TMath::Tan(theta));
v.Rotate(rndm->Rndm() * 2 * TMath::Pi(),fDir);
rtn += v;
rtn = rtn.Unit();
return rtn;
}
//____________________________________________________________________________________________________
Double_t BunchOfPhotons::AngularDist_OverTwoPi(Double_t theta) const {
//
// returns 'angular distribution value' / 2pi, at theta
// theta taken w.r.t bunch axis, and zero means same direction as bunch axis
//
Double_t rtn = 0;
if (fParent->GetHistoAngular()) {
TH1F histo(*fParent->GetHistoAngular());
rtn = histo.GetBinContent(histo.FindBin(theta));
}
else if(fParent->GetAngularDist()) {
TF12 distrib(*fParent->GetAngularDist());
rtn = distrib.Eval(theta);
}
// if isotropic fluo bunch
else if(fType == Fluo) rtn = 0.5;
else rtn = 0;
return rtn/TwoPi();
}
//____________________________________________________________________________________________________
void BunchOfPhotons::AddToPosTof(const EarthVector& v) {
//
// Increment the position AND the tof
//
fTof += v.Mag() / Clight();
fPos += v;
}
//____________________________________________________________________________________________________
void BunchOfPhotons::Reset() {
//
// Reset method
//
fCounter = 0;
}
/*
//_____________________________________________________________________________
void BunchOfPhotons::ConvoluteTransSpectra(UInt_t nb, const vector<Double_t*>& coeff) {
//
// convolute the transmission spectra
//
if(nb != fNbins)
Msg(EsafMsg::Warning) <<"<ConvoluteTransSpectra> Number of wavelengths doesn't match what is foreseen" << MsgDispatch;
fTransSpectrum.Convolute(coeff[0]);
fRaylSpectrum.Convolute(coeff[1]);
fOzoneSpectrum.Convolute(coeff[2]);
fMieSpectrum.Convolute(coeff[3]);
}
*/