// ESAF : Euso Simulation and Analysis Framework
// $Id: EShowerPainter.cc,v 1.19 2005/04/18 22:50:21 thea Exp $
// Naumov Dmitry created Jun, 28 2004
#include "EShowerPainter.hh"
#include "ESystemOfUnits.hh"
#include "EShower.hh"
#include "EShowerStep.hh"
#include "Etypes.hh"
#include "Riostream.h"
#include <TGeoBBox.h>
#include <TGeoVolume.h>
#include <TGeoMatrix.h>
#include <TVector3.h>
#include <TGeoManager.h>
#include <TGeoTube.h>
#include <TPad.h>
#include <TTimer.h>
#include <TSeqCollection.h>
#include <TH2F.h>
#include <TH3F.h>
using namespace TMath;
ClassImp(EShowerPainter)
//_____________________________________________________________________________
EShowerPainter::EShowerPainter( EShower *sh ) {
//
// Constructor
//
fShower = sh;
fNumFrames = 100;
fFrame = -1;
fWorldBuild = kFALSE;
fBoxDefined = kFALSE;
fHistoContainer = new TList();
fThreshold = 1.e-3;
if ( !fShower || fShower->GetNumSteps() <= 0 ) {
Info("Build()","Empty Shower ( NumSteps = 0 ). Painter made zombie.");
MakeZombie();
return;
}
}
//_____________________________________________________________________________
EShowerPainter::~EShowerPainter() {
//
// Destructor
//
fHistoContainer->Delete();
}
//______________________________________________________________________________
void EShowerPainter::DefineBox() {
if (IsBoxDefined()) return;
Int_t n = fShower->GetNumSteps();
Double_t Xmin(kHuge), Xmax(-kHuge), Ymin(kHuge), Ymax(-kHuge), Zmin(kHuge), Zmax(-kHuge), Nmax(1);
for (Int_t i=0; i<n; i++) {
if (Nmax < fShower->GetStep(i)->GetNumElectrons() )
Nmax = fShower->GetStep(i)->GetNumElectrons();
}
for (Int_t i=0; i<n; i++) {
Double_t xmean = (fShower->GetStep(i)->GetPosi().X() + fShower->GetStep(i)->GetPosf().X())/2/km;
Double_t ymean = (fShower->GetStep(i)->GetPosi().Y() + fShower->GetStep(i)->GetPosf().Y())/2/km;
Double_t zmean = (fShower->GetStep(i)->GetPosi().Z() + fShower->GetStep(i)->GetPosf().Z())/2/km;
Double_t ne = fShower->GetStep(i)->GetNumElectrons();
if (ne/Nmax > fThreshold) {
Xmin = (Xmin < xmean ? Xmin : xmean);
Xmax = (Xmax > xmean ? Xmax : xmean);
Ymin = (Ymin < ymean ? Ymin : ymean);
Ymax = (Ymax > ymean ? Ymax : ymean);
Zmin = (Zmin < zmean ? Zmin : zmean);
Zmax = (Zmax > zmean ? Zmax : zmean);
}
}
fXmin = Xmin;
fXmax = Xmax;
fYmin = Ymin;
fYmax = Ymax;
fZmin = Zmin;
fZmax = Zmax;
fNmax = Nmax;
fBoxDefined = kTRUE;
}
//______________________________________________________________________________
void EShowerPainter::Build() {
//
// Find the size of the box that contains the shower
//
if ( IsZombie() ) return;
DefineBox();
Int_t n = fShower->GetNumSteps();
fDx = 1.;
fDy = 1.;
fDz = (fShower->GetStep(0)->GetPosi() - fShower->GetStep(n-1)->GetPosf()).Mag()/2/km;
TVector3 initPos = (fShower->GetStep(0)->GetPosi() + fShower->GetStep(0)->GetPosf())*0.5;
TVector3 lastPos = (fShower->GetStep(n-1)->GetPosi() + fShower->GetStep(n-1)->GetPosf())*0.5;
fCenter = (initPos + lastPos)*0.5*(1/km);
// build the shower track
Double_t kScale = 1./fNmax;
fShowerVolume = gGeoManager->MakeBox("fShowerVolume",fVacuum,fDx,fDy,fDz);
Double_t dz_shift = 0;
for (Int_t i=0; i<n; i++) {
TVector3 stepCenter = (fShower->GetStep(i)->GetPosi() + fShower->GetStep(i)->GetPosf())*0.5;
double dz = (fShower->GetStep(i)->GetPosi() - fShower->GetStep(i)->GetPosf()).Mag()/2/km;
double rmax = fShower->GetStep(i)->GetNumElectrons()*kScale;
if (fShower->GetStep(i)->GetNumElectrons()/fNmax > fThreshold) {
rmax = rmax < 0.001 ? 0.001 : rmax;
TString name = Form("tube_%d",i);
dz_shift = (stepCenter-initPos).Mag()/km-fDz;
TGeoVolume *s = gGeoManager->MakeTube(name, fVacuum, 0.9*rmax, rmax, dz);
s->SetLineColor(50+(Int_t)(50.*log(fShower->GetStep(i)->GetNumElectrons()/fNmax)));
fShowerVolume->AddNode(s,1, new TGeoTranslation(0,0,-dz_shift));
}
}
}
//______________________________________________________________________________
void EShowerPainter::BuildWorld() {
//
// Build World based on the shower geometry
//
if ( IsZombie() ) return;
if( IsWorldBuilt() ) return;
new TGeoManager("ShowerViewer","Shower Viewer");
TGeoMaterial *matVacuum = new TGeoMaterial("Al", 26.98,13,2.7);
fVacuum = new TGeoMedium("Al", 1, matVacuum);
fTOP = gGeoManager->MakeBox("fTOP",fVacuum,230,230,20);
gGeoManager->SetTopVolume(fTOP);
Build();
Double_t OffSet = 1.1;
Double_t x1(-5),x2(5),y1(-5),y2(5); // make 10x10x30 kms box for default
if (fXmin < x1) x1 = fXmin*OffSet;
if (fXmax > x2) x2 = fXmax*OffSet;
if (fYmin < y1) y1 = fYmin*OffSet;
if (fYmax > y2) y2 = fYmax*OffSet;
fFoVVolume = gGeoManager->MakeBox("fFoVVolume",fVacuum,(x2-x1)/2,(y2-y1)/2,15);
Double_t fTheta = fShower->GetTheta()*RadToDeg();
Double_t fPhi = fShower->GetPhi()*RadToDeg();
TGeoTranslation tr;
TGeoRotation rot;
tr.SetTranslation(fCenter.X(),fCenter.Y(),fCenter.Z());
rot.SetAngles(fTheta,Min(fPhi,360-fPhi),0);
TGeoHMatrix h = tr*rot;
fShowerVolume->SetLineColor(kGreen);
fFoVVolume->SetLineColor(kBlue);
fTOP->AddNode(fShowerVolume,1,new TGeoHMatrix(h));
fTOP->AddNode(fFoVVolume,1, new TGeoTranslation((x1+x2)/2,(y1+y2)/2,15));
gGeoManager->CloseGeometry();
}
//______________________________________________________________________________
void EShowerPainter::UpdateVisibility( Option_t *opt ) {
//
// Fast update ov fisibility status of the objects
//
if ( !gGeoManager ) return;
TGeoVolume *top;
if ( !(top = gGeoManager->GetTopVolume()) ) return;
TGeoVolume *shower = top->FindNode("fShowerVolume_1")->GetVolume();
if ( !shower ) return;
Int_t n;
if ( fFrame == -1)
n = shower->GetNdaughters();
else
n = (Int_t)(((Float_t)shower->GetNdaughters()/fNumFrames)*fFrame);
shower->SetVisibility(kFALSE);
shower->VisibleDaughters(kFALSE);
for (Int_t i=0; i<shower->GetNdaughters(); i++) {
Bool_t vis = ( i <= n );
if ( shower->GetNode(i)->GetVolume()->IsVisible() != vis ) {
shower->GetNode(i)->GetVolume()->SetVisibility(vis);
}
}
shower->SetVisibility(kTRUE);
shower->VisibleDaughters(kTRUE);
return;
}
//______________________________________________________________________________
void EShowerPainter::Draw( Option_t* ) {
//
// Draw the shower in its world volume
//
BuildWorld();
UpdateVisibility();
if ( !IsZombie() ) {
fTOP->Draw();
}
}
//______________________________________________________________________________
void EShowerPainter::DrawXY() {
// Draw the shower projection on XY Earth surface
// build histograms
TH2F *th0 = GetShowerHisto2D("ShowerXY", "XY projection");
for (Int_t i=0; i<fShower->GetNumSteps(); i++) {
double xmean = (fShower->GetStep(i)->GetPosi().X() + fShower->GetStep(i)->GetPosf().X())/2/km;
double ymean = (fShower->GetStep(i)->GetPosi().Y() + fShower->GetStep(i)->GetPosf().Y())/2/km;
double Ne = fShower->GetStep(i)->GetNumElectrons();
th0->Fill(xmean,ymean,Ne);
}
th0->Draw("colz");
}
//______________________________________________________________________________
void EShowerPainter::DrawXYZ() {
// Draw the shower track in XYZ
// build histograms
TH3F *th0 = GetShowerHisto3D("ShowerXYZ", "XYZ development");
for (Int_t i=0; i<fShower->GetNumSteps(); i++) {
double xmean = (fShower->GetStep(i)->GetPosi().X() + fShower->GetStep(i)->GetPosf().X())/2/km;
double ymean = (fShower->GetStep(i)->GetPosi().Y() + fShower->GetStep(i)->GetPosf().Y())/2/km;
double zmean = (fShower->GetStep(i)->GetPosi().Z() + fShower->GetStep(i)->GetPosf().Z())/2/km;
double Ne = fShower->GetStep(i)->GetNumElectrons();
th0->Fill(xmean,ymean,zmean,Ne/fNmax);
}
th0->Draw("box");
}
//_____________________________________________________________________________
void EShowerPainter::Animate() {
//
// Move to next animation step and increase the counter
//
if ( IsEnded() ) {
// leave the last frame to be painted
fFrame = fNumFrames-1;
Stop();
return;
}
NextFrame();
}
//_____________________________________________________________________________
void EShowerPainter::NextFrame() {
//
// Move to the next animation frame
//
if ( !gGeoManager ) return;
TGeoVolume *top = gGeoManager->GetTopVolume();
TGeoVolume *shower = top->FindNode("fShowerVolume_1")->GetVolume();
Int_t n = (Int_t)(((Float_t)shower->GetNdaughters()/fNumFrames)*fFrame);
for (Int_t i=0; i<shower->GetNdaughters(); i++) {
Bool_t vis = ( i <= n );
if ( shower->GetNode(i)->GetVolume()->IsVisible() != vis ) {
shower->GetNode(i)->GetVolume()->SetVisibility(vis);
}
}
fFrame++;
}
//_____________________________________________________________________________
TH2F *EShowerPainter::GetShowerHisto2D(const char *name, const char *title) {
// Draw XY shower track projection. X and Y of the histogram have always center
// at zero for a better display
//
DefineBox();
if ( title == 0 ) title = name;
TH2F* th = (TH2F*)fHistoContainer->FindObject(name);
if (!th) {
Double_t OffSet = 1.1;
Double_t x1(-5),x2(5),y1(-5),y2(5); // make 5 kms box for default
if (fXmin < x1) x1 = fXmin*OffSet;
if (fXmax > x2) x2 = fXmax*OffSet;
if (fYmin < y1) y1 = fYmin*OffSet;
if (fYmax > y2) y2 = fYmax*OffSet;
// define number of steps
Int_t nx = (Int_t)(x2-x1);
Int_t ny = (Int_t)(y2-y1);
th = new TH2F( name, title,nx,x1,x2,ny,y1,y2);
fHistoContainer->Add(th);
th->SetStats(0);
th->SetXTitle("X [km]");
th->SetYTitle("Y [km]");
}
return th;
}
//_____________________________________________________________________________
TH3F *EShowerPainter::GetShowerHisto3D(const char *name, const char *title) {
//
//
//
DefineBox();
if ( title == 0 ) title = name;
TH3F* th = (TH3F*)fHistoContainer->FindObject(name);
if (!th) {
Double_t OffSet = 1.1;
Double_t x1(-5),x2(5),y1(-5),y2(5),z1(0),z2(30); // make 10x10x30 kms box for default
if (fXmin < x1) x1 = fXmin*OffSet;
if (fXmax > x2) x2 = fXmax*OffSet;
if (fYmin < y1) y1 = fYmin*OffSet;
if (fYmax > y2) y2 = fYmax*OffSet;
th = new TH3F( name, title, 30,x1,x2,30,y1,y2,30,z1,z2);
fHistoContainer->Add(th);
th->SetStats(0);
th->SetXTitle("X [km]");
th->SetYTitle("Y [km]");
th->SetZTitle("Z [km]");
}
return th;
}