// ESAF : Euso Simulation and Analysis Framework
// $Id: EShowerHistoPainter.cc,v 1.8 2005/01/18 16:52:10 moreggia Exp $
// Author: Naumov Dmitry Jul, 26 2004
#include "EShowerHistoPainter.hh"
#include "EShower.hh"
#include "EShowerStep.hh"
#include "ESystemOfUnits.hh"
#include <TTimer.h>
#include <TPad.h>
#include <TH1F.h>
#include <TPaveText.h>
ClassImp(EShowerHistoPainter)
// TOOL function : give local altitude
Double_t Zv(const TVector3& v) {
return sqrt( pow(6.370949e6*m + v.Z(),2) + v.Perp2() ) - 6.370949e6*m;
}
//_____________________________________________________________________________
EShowerHistoPainter::EShowerHistoPainter( EShower* sh ) {
//
// Constructor
//
fShower = sh;
fHistos = new TList();
fTimeBins = 0;
fAltBins = 0;
fNumFrames = 100;
fFrame = -1;
Build();
}
//_____________________________________________________________________________
EShowerHistoPainter::~EShowerHistoPainter() {
//
// Destructor
//
if ( fTimeBins ) delete [] fTimeBins;
if ( fAltBins ) delete [] fAltBins;
if ( fHistos ) fHistos->Delete();
}
//_____________________________________________________________________________
void EShowerHistoPainter::Build() {
//
// Build
//
if ( !fShower ) {
Info("Build()","EShower object is NULL. Painter made zombie.");
MakeZombie();
return;
}
fNeMax = 0;
fAgeMax = 0;
fDepthMax = 0;
Int_t n = fShower->GetNumSteps();
if ( n <= 0 ) {
Info("Build()","Empty Shower ( NumSteps = 0 ). Painter made zombie.");
MakeZombie();
return;
}
fTimeBins = new Double_t[n+1];
fAltBins = new Double_t[n+1];
for (Int_t i=0; i<n; i++) {
fTimeBins[i] = fShower->GetStep(i)->GetTimei()/ms;
fTimeBins[i+1] = fShower->GetStep(i)->GetTimef()/ms;
// if ( i > 0 && fTimeBins[i-1] > fTimeBins[i] ) {
// Printf("Porcapupazza i = %d, fTimeBins[i-1] = %f, fTimeBins[i] = %f",fTimeBins[i-1], fTimeBins[i]);
// }
if( fNeMax < fShower->GetStep(i)->GetNumElectrons() )
fNeMax = fShower->GetStep(i)->GetNumElectrons();
if ( fAgeMax < 0.5*(fShower->GetStep(i)->GetAgei()
+ fShower->GetStep(i)->GetAgef()))
fAgeMax = 0.5*(fShower->GetStep(i)->GetAgei()
+ fShower->GetStep(i)->GetAgef());
if ( fDepthMax < (fShower->GetStep(i)->GetXi()
+ fShower->GetStep(i)->GetXf())/2/g*cm2)
fDepthMax = (fShower->GetStep(i)->GetXi()
+ fShower->GetStep(i)->GetXf())/2/g*cm2;
fAltBins[n - i] = Zv(fShower->GetStep(i)->GetPosi())/km;
fAltBins[n-1 - i] = Zv(fShower->GetStep(i)->GetPosf())/km;
}
// build histograms
TH1F *th0 = GetTimeHisto("Shower_Ne", "Ne vs time");
TH1F *th1 = GetTimeHisto("Shower_Age", "Age vs time");
TH1F *th2 = GetTimeHisto("Shower_Depth", "Depth vs time");
th0->SetYTitle("Ne");
th1->SetYTitle("Age");
th2->SetYTitle("Slant depth");
TH1F *alt0 = GetAltHisto("Shower_Ne_alt", "Ne vs altitude");
alt0->SetYTitle("Ne");
for (Int_t i=0; i<n; i++) {
Double_t time = (fShower->GetStep(i)->GetTimei()
+ fShower->GetStep(i)->GetTimef())/2/ms;
Double_t Ne = fShower->GetStep(i)->GetNumElectrons();
Double_t Age = (fShower->GetStep(i)->GetAgei()
+ fShower->GetStep(i)->GetAgef())/2;
Double_t Depth = (fShower->GetStep(i)->GetXi()
+ fShower->GetStep(i)->GetXf())/2/g*cm2;
Double_t alt = (Zv(fShower->GetStep(i)->GetPosi())
+ Zv(fShower->GetStep(i)->GetPosf()))/2/km;
th0->Fill(time,Ne);
th1->Fill(time,Age);
th2->Fill(time,Depth);
alt0->Fill(alt,Ne);
}
}
//_____________________________________________________________________________
TH1F *EShowerHistoPainter::GetTimeHisto(const char *name, const char *title) {
//
//
//
if ( title == 0 ) title = name;
TH1F* th = (TH1F*)fHistos->FindObject(name);
if (!th) {
th = new TH1F( name, title, fShower->GetNumSteps(), fTimeBins);
fHistos->Add(th);
th->SetStats(0);
th->SetXTitle("Development Time, #{m}s");
th->SetFillColor(9);
}
return th;
}
//_____________________________________________________________________________
TH1F *EShowerHistoPainter::GetAltHisto(const char *name, const char *title) {
//
//
//
if ( title == 0 ) title = name;
TH1F* th = (TH1F*)fHistos->FindObject(name);
if (!th) {
th = new TH1F( name, title, fShower->GetNumSteps(), fAltBins);
fHistos->Add(th);
th->SetStats(0);
th->SetXTitle("Altitude (km)");
th->SetFillColor(9);
}
return th;
}
//______________________________________________________________________________
void EShowerHistoPainter::Draw(Option_t *opt) {
//
//
//
if ( IsZombie() ) return;
TString option(opt);
TH1F* th(0);
if ( option == "" ) option = "Ne";
if(option == "Ne_alt") {
th = GetAltHisto("Shower_Ne_alt");
if(th) th->Draw();
}
if ( option == "Ne" ) {
th = GetTimeHisto("Shower_Ne");
if ( th ) th->Draw();
}
else if ( option == "Age" ) {
th = GetTimeHisto("Shower_Age");
if ( th ) th->Draw();
}
else if ( option == "Depth" ) {
th = GetTimeHisto("Shower_Depth");
if ( th ) th->Draw();
}
else if (option == "Truth") {
TPaveText* hdemo = new TPaveText(.05,.05,.95,.95);
hdemo->SetTextAlign(12);
hdemo->SetTextFont(52);
hdemo->SetTextColor(4);
hdemo->AddText("");
hdemo->AddText("This event was simulated with:");
hdemo->AddText("");
hdemo->AddText(Form("Energy %.4e MeV",fShower->GetEnergy()));
hdemo->AddText(Form("Theta %.2f deg",fShower->GetTheta()*TMath::RadToDeg()));
hdemo->AddText(Form("Phi %.2f deg",fShower->GetPhi()*TMath::RadToDeg()));
hdemo->AddText(Form("X1 %.2f g/cm^{2}",fShower->GetX1()*cm2/g));
hdemo->AddText(Form("Xinit %.2f km",fShower->GetInitPos().X()/km));
hdemo->AddText(Form("Yinit %.2f km",fShower->GetInitPos().Y()/km));
hdemo->AddText(Form("Zinit %.2f km",fShower->GetInitPos().Z()/km));
hdemo->AddText("");
hdemo->Draw();
}
}
//_____________________________________________________________________________
void EShowerHistoPainter::NextFrame() {
//
// Switch to the next frame and update histograms
//
fFrame++;
UpdateHistos();
}
//______________________________________________________________________________
void EShowerHistoPainter::UpdateHistos() {
//
// Update hitograms at the current frame.
//
if ( IsZombie() ) return;
Int_t n;
if ( fFrame == -1 )
n = fShower->GetNumSteps();
else
n = (Int_t)(((Float_t)fShower->GetNumSteps()/fNumFrames)*fFrame);
TH1F* thNe = GetTimeHisto("Shower_Ne");
TH1F* thAge = GetTimeHisto("Shower_Age");
TH1F* thDepth = GetTimeHisto("Shower_Depth");
for (Int_t i=0; i<fShower->GetNumSteps() ; i++) {
Double_t Ne = fShower->GetStep(i)->GetNumElectrons();
Double_t Age = (fShower->GetStep(i)->GetAgei()
+ fShower->GetStep(i)->GetAgef())/2;
Double_t Depth = (fShower->GetStep(i)->GetXi()
+ fShower->GetStep(i)->GetXf())/2/g*cm2;
thNe->SetBinContent(i+1, i < n ? Ne : 0 );
thAge->SetBinContent(i+1, i < n ? Age : 0 );
thDepth->SetBinContent(i+1, i < n ? Depth : 0 );
}
thNe->SetMaximum( fNeMax*1.1 );
thAge->SetMaximum( fAgeMax*1.1 );
thDepth->SetMaximum( fDepthMax*1.1 );
}