// $Id: RecoShowerTrack.cc,v 1.9 2005/11/01 11:08:59 naumov Exp $
// Author: Dmitry.Naumov 2005/10/11
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: RecoShowerTrack *
* Package: Energy *
* Coordinator: <Dmitry.Naumov> *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// RecoShowerTrack
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// <parameter name>: <parameter description>
// -Valid options: <available options>
//
#include "RecoShowerTrack.hh"
#include "EConst.hh"
#include "SlastShowerGenerator.hh"
#include "ShowerTrack.hh"
#include "FluoCalculator.hh"
#include "Config.hh"
#include "EsafSpectrum.hh"
#include "LightSourceFactory.hh"
#include "Ground.hh"
#include "O1_ClearSkyPropagator.hh"
#include "DetectorGeometry.hh"
#include "RadiativeFactory.hh"
#include "EnergyModule.hh"
#include "RecoEvent.hh"
#include <TMath.h>
#include <TCanvas.h>
#include <TProfile.h>
using namespace sou;
using namespace TMath;
using namespace EConst;
ClassImp(RecoShowerTrack)
// constructors, destructor
//_____________________________________________________________________________
RecoShowerTrack::RecoShowerTrack(): EsafMsgSource() {
//
// Constructor
//
Constructor();
}
//_____________________________________________________________________________
RecoShowerTrack::RecoShowerTrack(Float_t t, Float_t p, Float_t h, Float_t tfov, Float_t pfov): EsafMsgSource() {
//
// Constructor
//
Constructor();
SetThetaPhi(t,p);
SetThetaPhiFoV(tfov,pfov);
SetHmax(h);
}
//_____________________________________________________________________________
void RecoShowerTrack::Constructor() {
//
// This is default ctor
fHmax = -100;
fShowerDir.SetXYZ(0,0,0);
fEusoDir.SetXYZ(0,0,0);
string fluoname = Config::Get()->GetCF("Reco","EnergyModule")->GetStr("EnergyModule.fFluoCalculator");
fFluoCalculator = LightSourceFactory::Get()->GetFluoCalculator( fluoname );
fEusoAlt = Config::Get()->GetCF("General","Euso")->GetNum("Euso.fAltitude")*km;
fEnergyDistributionName = Config::Get()->GetCF("Reco","EnergyModule")->GetStr("EnergyModule.fEnergyDistributionName");
fEsafSpectrum = new EsafSpectrum();
fShowerGenerator = new SlastShowerGenerator(kTRUE);
fShowerGenerator->SetQuiet();
fDetectorGeometry = new DetectorGeometry();
Double_t radius = Config::Get()->GetCF("Optics","ConicBaffle")->GetNum("ConicBaffle.fTopRadius")*mm;
fDetectorGeometry->SetRadius(radius);
fDetectorGeometry->SetPos(TVector3(0,0,fEusoAlt));
fGround = RadiativeFactory::Get()->GetGround();
fClearSkyPropagator = new O1_ClearSkyPropagator(fGround);
fClearSkyPropagator->CopyDetectorGeometry(fDetectorGeometry,false,false);
fEnergyModule = NULL;
}
//_____________________________________________________________________________
RecoShowerTrack::~RecoShowerTrack() {
//
// Destructor
//
SafeDelete(fShowerGenerator);
SafeDelete(fDetectorGeometry);
SafeDelete(fFluoCalculator);
SafeDelete(fClearSkyPropagator);
}
//_____________________________________________________________________________
void RecoShowerTrack::Clear() {
//
//
//
for (UInt_t i(0); i < fSteps.size(); i++) {
RecoShowerStep& rs = fSteps[i];
rs.fSpectrum->Delete();
}
fSteps.clear();
}
//______________________________________________________________________________
const RecoShowerStep& RecoShowerTrack::GetStepThetaPhi(Float_t Theta, Float_t Phi) const {
//
// Return the step of the shower corresponding to the given time
Int_t index = 0;
Double_t mmin = -1.;
TVector3 nFoV;
nFoV.SetMagThetaPhi(1.,Theta,Phi);
for (UInt_t i(0); i < fSteps.size(); i++) {
const RecoShowerStep& rs = fSteps[i];
EarthVector dir = (rs.fDir).Unit();
Double_t mag = dir*nFoV;
if (mag > mmin) {mmin = mag; index = i;}
}
return (fSteps[index]);
}
//_____________________________________________________________________________
void RecoShowerTrack::MakeTrack() {
//
//
//
Int_t progress = 0;
Clear();
if (fHmax == -100 || !fShowerDir.Mag() || !fEusoDir.Mag()) {
if (fHmax == -100) Msg(EsafMsg::Warning) << "Hmax is not setupn" << MsgDispatch;
if (!fShowerDir.Mag()) Msg(EsafMsg::Warning) << "Shower Direction is not setupn" << MsgDispatch;
if (!fEusoDir.Mag()) Msg(EsafMsg::Warning) << "Shower Direction is not setupn" << MsgDispatch;
Msg(EsafMsg::Panic) << "Can not do a showern" << MsgDispatch;
}
// find 3D point of shower maximum
Float_t AB = EarthRadius() + fEusoAlt;
Float_t AC = EarthRadius() + fHmax;
Float_t cs = Cos(fEusoDir.Theta());
Float_t sn = Sin(fEusoDir.Theta());
Float_t z = fEusoAlt - AB*cs*cs + cs*Sqrt(AC*AC - AB*AB*sn*sn);
Float_t rho = (fEusoAlt-z)*Tan(fEusoDir.Theta());
Float_t x = -rho*Cos(fEusoDir.Phi());
Float_t y = -rho*Sin(fEusoDir.Phi());
Float_t tMin = 0, tMax = 0;
Float_t age, TotalYield, Attenuation;
EarthVector MaximumPosition(x,y,z);
EarthVector ISS(0,0,fEusoAlt);
Float_t Jump = 60.*km;
EarthVector InitPos = MaximumPosition + Jump*fShowerDir;
fShowerGenerator->SetShower(1.e21,fShowerDir.Theta(),fShowerDir.Phi(),
InitPos);
fShowerGenerator->SetQuiet(kTRUE);
fShowerGenerator->SetDepthStep(0.2);
ShowerTrack *fShowerTrack = (ShowerTrack*)fShowerGenerator->Get();
// Define tMin and tMax of shower evolution
{
const ShowerStep& s = fShowerTrack->GetStep(0);
Float_t time = 0.5*(s.GetTimei() + s.GetTimef());
EarthVector dir = (ISS - 0.5*(s.GetXYZi() + s.GetXYZf()));
tMin = (time + dir.Mag()/Clight())*1.e-3; // in microsecond
}
{
const ShowerStep& s = fShowerTrack->GetStep(fShowerTrack->Size()-1);
Float_t time = 0.5*(s.GetTimei() + s.GetTimef());
EarthVector dir = (ISS - 0.5*(s.GetXYZi() + s.GetXYZf()));
tMax = (time + dir.Mag()/Clight())*1.e-3; // in microsecond
}
Float_t GtuStep = 2.5;
if (fEnergyModule)
GtuStep=fEnergyModule->GetRecoEvent()->GetHeader().GetGtuLength()*1.e-3;
Int_t Nbins = (Int_t)Floor((tMax - tMin)/GtuStep);
TH1F* h[10];
for (Int_t i(0); i < 10; i++) {
TString title = Form("EnergyModule_Step%d",i);
TH1F* hh = NULL;
if ((hh = (TH1F*)gDirectory->FindObject(title))) hh->Delete();
h[i] = new TH1F(title,"",Nbins,0,tMax - tMin);
}
for (UInt_t i(0); i < fShowerTrack->Size(); i++) {
const ShowerStep& s = fShowerTrack->GetStep(i);
// get position and compute for it: fluo yield + attenuation + time at EUSO
EarthVector pos = 0.5*(s.GetXYZi() + s.GetXYZf());
Float_t time = 0.5*(s.GetTimei() + s.GetTimef());
Float_t X = 0.5*(s.GetXi() + s.GetXf());
EarthVector dir = ISS - pos;
Float_t EusoTime = (time + dir.Mag()/Clight())*1.e-3 - tMin; // in microsecond
if (!(i % 50)) {
age = 0.5*(s.GetAgei() + s.GetAgef());
TF2* EnergyDistribution = (TF2*)ShowerTrack::GetEnergyDistribution(fEnergyDistributionName);
TF12 EnergyDistributionProj("EDS_name",EnergyDistribution,age,"X");
TotalYield = fFluoCalculator->GetFluoYield(pos.Zv(),&EnergyDistributionProj,fEsafSpectrum);
Attenuation=fClearSkyPropagator->GoToDetector(*fEsafSpectrum,pos.Zv());
}
Float_t Length = (s.GetXYZi() - s.GetXYZf()).Mag();
Double_t SolidAngle = 1./(4.*Pi()*dir.Mag2());
Double_t flux = SolidAngle*s.GetNelectrons()*Attenuation*TotalYield*Length;
h[0]->Fill(EusoTime,s.GetNelectrons());
h[1]->Fill(EusoTime,flux);
h[2]->Fill(EusoTime,pos.X()*flux);
h[3]->Fill(EusoTime,pos.Y()*flux);
h[4]->Fill(EusoTime,pos.Z()*flux);
h[5]->Fill(EusoTime,EusoTime*flux);
h[6]->Fill(EusoTime,TotalYield*flux);
h[7]->Fill(EusoTime,Attenuation*flux);
h[8]->Fill(EusoTime,age*flux);
h[9]->Fill(EusoTime,X*flux);
if (100.*(i+1)/fShowerTrack->Size() >= progress) {
Msg(EsafMsg::Info).SetProgress(progress);
Msg(EsafMsg::Info) << "MakeTrack():" << MsgCount;
progress += 10;
}
}
if (progress <= 100) Msg(EsafMsg::Info) << MsgDispatch;
for (Int_t i = 1; i <= Nbins; i++) {
Double_t flux = h[1]->GetBinContent(i);
if (!flux) continue;
Double_t ne = h[0]->GetBinContent(i);
Double_t x = h[2]->GetBinContent(i)/flux;
Double_t y = h[3]->GetBinContent(i)/flux;
Double_t z = h[4]->GetBinContent(i)/flux;
Double_t time = h[5]->GetBinContent(i)/flux;
Double_t yield = h[6]->GetBinContent(i)/flux;
Double_t atten = h[7]->GetBinContent(i)/flux;
Double_t age = h[8]->GetBinContent(i)/flux;
Double_t X = h[9]->GetBinContent(i)/flux;
EarthVector pos(x,y,z);
EarthVector dir = ISS - pos;
TVector3 nFoV;
nFoV.SetMagThetaPhi(1.,dir.Theta(),dir.Phi());
Float_t length = 1.*Clight()*GtuStep*1.e+3/(1.+fShowerDir*nFoV);
RecoShowerStep rs;
rs.fEusoTime = time;
rs.fPos = pos;
rs.fDir = dir;
rs.fFluoYield = yield;
rs.fAttenuation = atten;
rs.fNe = ne;
rs.fFlux = flux;
rs.fLength = length;
rs.fX = X*cm2/g;
TF2 *EnergyDistribution = (TF2*)ShowerTrack::GetEnergyDistribution(fEnergyDistributionName);
TF12 EnergyDistributionProj("EDS_name",EnergyDistribution,age,"X");
fFluoCalculator->GetFluoYield(pos.Zv(),&EnergyDistributionProj,
fEsafSpectrum);
fClearSkyPropagator->GoToDetector(*fEsafSpectrum,pos.Zv());
rs.fSpectrum = new TH1F(Form("RecoShowerTrackSpectrum%d",i),"",
100,300*nm,400*nm);
for (Int_t nt = 0; nt < 1000; nt++)
rs.fSpectrum->Fill(fEsafSpectrum->GetLambda());
Add(rs);
}
for (Int_t i(0); i < 10; i++) h[i]->Delete();
fShowerTrack->Reset();
fShowerGenerator->Reset();
}
//_____________________________________________________________________________
void RecoShowerTrack::Debug(TCanvas *c, Int_t id) {
//
//
//
if (!c) {
Warning("RecoShowerTrack::Debug","Please specify the pad where I have to draw");
return;
}
if (id < 0 || id > 10) {
Warning("RecoShowerTrack::Debug","Please specify id between 0 and 9");
return;
}
Float_t tMin, tMax;
// Set tMin and tMax
{
RecoShowerStep& rs = fSteps[0];
tMin = rs.fEusoTime;
}
{
RecoShowerStep& rs = fSteps[fSteps.size()-1];
tMax = rs.fEusoTime;
}
TH1* h[11];
for (Int_t i(0); i < 11; i++) {
TString title = Form("debug%d",i);
TProfile* hh = NULL;
if (( hh = (TProfile*)gDirectory->FindObject(title))) hh->Delete();
if (i == 10) continue;
h[i] = new TProfile(title,"",fSteps.size(),tMin,tMax);
}
h[10] = new TH1F("debug10","",fSteps.size(),tMin,tMax);
for (UInt_t i(0); i < fSteps.size(); i++) {
RecoShowerStep& rs = fSteps[i];
h[0]->Fill(rs.fEusoTime,rs.fNe);
h[1]->Fill(rs.fEusoTime,rs.fPos.X()/km);
h[2]->Fill(rs.fEusoTime,rs.fPos.Y()/km);
h[3]->Fill(rs.fEusoTime,rs.fPos.Z()/km);
h[4]->Fill(rs.fEusoTime,rs.fFluoYield*m);
h[5]->Fill(rs.fEusoTime,rs.fAttenuation);
h[6]->Fill(rs.fEusoTime,rs.fFlux);
h[7]->Fill(rs.fEusoTime,rs.fLength/km);
h[8]->Fill(rs.fEusoTime,rs.fX);
h[10]->Fill(rs.fEusoTime,rs.fLength/km);
}
h[id]->Draw();
}