// $Id: ETreePainter.cc,v 1.2 2005/04/07 15:22:12 stutz Exp $
// Author: Anne Stutz 2005/03/03
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: ETreePainter *
* Package: <packagename> *
* Coordinator: <coordinator> *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// ETreePainter
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// <parameter name>: <parameter description>
// -Valid options: <available options>
//
#include "ETreePainter.hh"
#include "ESystemOfUnits.hh"
#include "EAtmosphereHistoPainter.hh"
#include "EAtmosphere.hh"
#include "ETruth.hh"
#include "EShower.hh"
#include "EShowerStep.hh"
#include "EEvent.hh"
#include <TVirtualGeoTrack.h>
#include <TGeoVolume.h>
#include <TGeoManager.h>
#include <TGeoTube.h>
#include <TTree.h>
#include <TH1F.h>
#include <TF1.h>
#include <TNtuple.h>
#include <iostream>
ClassImp(ETreePainter)
//_____________________________________________________________________________
ETreePainter::ETreePainter( TTree* etree ) {
//
// Constructor
//
fTree = etree;
fEv = NULL;
fWorldBuild = kFALSE;
fCurrentTrack = 0;
fHistos = new TList;
Build();
}
//_____________________________________________________________________________
ETreePainter::~ETreePainter() {
//
// Destructor
//
}
//_____________________________________________________________________________
void ETreePainter::Build() {
//
// Build
//
if ( !fTree ) {
Info("Build()","ETree object is NULL. Painter made zombie.");
MakeZombie();
return;
}
fEv = new EEvent();
fEv->SetBranches(fTree);
fEvTot = (Int_t)fTree->GetEntries();
BuildWorld();
EAtmosphereHistoPainter* atmohistopainter = NULL;
for (Int_t i=0; i< fEvTot; i++) {
Printf( "Entry %d: %d bytes read", i, fTree->GetEntry(i) );
EAtmosphere *atmo = fEv->GetAtmosphere();
ETruth *truth = fEv->GetTruth();
EShower *shower = fEv->GetShower();
if ( atmohistopainter ) delete atmohistopainter;
atmohistopainter = new EAtmosphereHistoPainter( atmo );
FillAtmoHistos(atmohistopainter);
FillTruthHistos( truth );
AddShower(shower);
FillTree();
}
}
//_____________________________________________________________________________
void ETreePainter::BuildWorld() {
if(fWorldBuild) return;
new TGeoManager("TreeShowerViewer","Tree Shower Viewer");
TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
fVacuum = new TGeoMedium("Vacuum", 1, matVacuum);
fTop = gGeoManager->MakeBox("fTop", fVacuum, 120, 120, 40);
fFoVVolume = gGeoManager->MakeTube("fFoV",fVacuum, 0,230,20);
gGeoManager->SetTopVolume(fTop);
fTop->AddNode(fFoVVolume, 1, new TGeoTranslation(0,0,20));
gGeoManager->SetVisLevel(1);
fFoVVolume->SetLineColor(kBlue);
gGeoManager->CloseGeometry();
fWorldBuild = kTRUE;
}
//_____________________________________________________________________________
TTree* ETreePainter::BuildTree() {
//
// build a simple tree
//
TTree* stree = (TTree*)fHistos->FindObject("stree");
if ( !stree ) {
stree = new TTree("stree","a Tree with simple datas");
// add the branches
stree->Branch("Energy", &fTrack.fEnergy, "fEnergy/F");
stree->Branch("Theta", &fTrack.fTheta, "fTheta/F");
stree->Branch("Phi", &fTrack.fPhi, "fPhi/F");
stree->Branch("Z1", &fTrack.fZ1, "fZ1/F");
stree->Branch("Zmax", &fTrack.fZmax, "fZmax/F");
stree->Branch("X1", &fTrack.fX1, "fX1/F");
stree->Branch("Xmax", &fTrack.fXmax, "fXmax/F");
stree->Branch("Nph_fluo",&fTrack.fNph_f, "fNph_f/F");
stree->Branch("Nph_cer", &fTrack.fNph_c, "fNph_c/F");
stree->Branch("Nmax_fluo",&fTrack.fNmax_f, "fNmax_f/F");
stree->Branch("Nmax_cer", &fTrack.fNmax_c, "fNmax_c/F");
stree->Branch("Zmax_fluo",&fTrack.fZmax_f, "fZmax_f/F");
stree->Branch("Zmax_cer", &fTrack.fZmax_c, "fZmax_c/F");
stree->Branch("Nph_fluo_pupil",&fTrack.fNph_f_p, "fNph_f_p/F");
stree->Branch("Nph_cer_back_pupil", &fTrack.fNph_cb_p, "fNph_cb_p/F");
stree->Branch("Nph_cer_refl_pupil", &fTrack.fNph_cr_p, "fNph_cr_p/F");
stree->Branch("Tmax_fluo_pupil",&fTrack.fTmax_f_p, "fTmax_f_p/F");
stree->Branch("Nmax_fluo_pupil", &fTrack.fNmax_f_p, "fNmax_f_p/F");
stree->Branch("Sig_fluo_pupil", &fTrack.fSig_f_p, "fSig_f_p/F");
// add stree to the collection
fHistos->Add(stree);
}
return stree;
}
//_____________________________________________________________________________
TH1F* ETreePainter::BuildHistos(const char *name, const char *title) {
//
// build hitograms
//
TH1F* th = (TH1F*)fHistos->FindObject(name);
if ( !th ) {
if ( name == "Energy" ) {
th = new TH1F( name, title, 100, 0, 100 );
th->SetXTitle( "Energy / 1e19" );
}
if ( name == "Theta" ) {
th = new TH1F( name, title, 100, 0, 90 );
th->SetXTitle( " " );
}
if ( name == "Phi" ) {
th = new TH1F( name, title, 100, 0, 360 );
th->SetXTitle( " " );
}
if ( name == "Z1" ) {
th = new TH1F( name, title, 100, 0, 50 );
th->SetXTitle( "Alt in km" );
}
if ( name == "Zmax" ) {
th = new TH1F( name, title, 100, 0, 10 );
th->SetXTitle( "Alt in km" );
}
if ( name == "X1" ) {
th = new TH1F( name, title, 100, 0, 50 );
th->SetXTitle( "Depth in g/cm2" );
}
if ( name == "Xmax" ) {
th = new TH1F( name, title, 100, 0, 1500 );
th->SetXTitle( "Depth in g/cm2" );
}
if ( name == "Nph_Fluo" ) {
th = new TH1F( name, title, 100, 0, 100 );
th->SetXTitle( "Nb of Photons / 1e14" );
}
if ( name == "Nph_Cer" ) {
th = new TH1F( name, title, 100, 0, 100 );
th->SetXTitle( "Nb of Photons / 1e14" );
}
if ( name == "Nph_Fluo_P" ) {
th = new TH1F( name, title, 100, 0, 5000 );
th->SetXTitle( "Nb of Photons" );
}
if ( name == "Nph_CB_P" ) {
th = new TH1F( name, title, 100, 0, 2000 );
th->SetXTitle( "Nb of Photons" );
}
if ( name == "Nph_CR_P" ) {
th = new TH1F( name, title, 100, 0, 1000 );
th->SetXTitle( "Nb of Photons" );
}
if ( name == "Nph_tot_P" ) {
th = new TH1F( name, title, 100, 0, 5000 );
th->SetXTitle( "Nb of Photons" );
}
if ( name == "Zmax_Fluo" ) {
th = new TH1F( name, title, 100, 0, 10 );
th->SetXTitle( "Altitude(km)" );
}
if ( name == "Zmax_Cer" ) {
th = new TH1F( name, title, 100, 0, 10 );
th->SetXTitle( "Altitude(km)" );
}
if ( name == "Nmax_Fluo" ) {
th = new TH1F( name, title, 100, 0, 100 );
th->SetXTitle( "Nb of photons / 1e12" );
}
if ( name == "Nmax_Cer" ) {
th = new TH1F( name, title, 100, 0, 100 );
th->SetXTitle( "Nb of photons / 1e12" );
}
if ( name == "Tmax_Fluo_P" ) {
th = new TH1F( name, title, 100, 0, 100 );
th->SetXTitle( "Time " );
}
if ( name == "Nmax_Fluo_P" ) {
th = new TH1F( name, title, 100, 0, 200 );
th->SetXTitle( "Nb of photons" );
}
if ( name == "Sig_Fluo_P" ) {
th = new TH1F( name, title, 100, 0, 100 );
th->SetXTitle( "Time" );
}
fHistos->Add(th);
}
// th->SetStats(0);
th->SetFillColor(9);
return th;
}
//_____________________________________________________________________________
void ETreePainter::AddShower(EShower *fShower) {
Int_t track_index = gGeoManager->AddTrack(fCurrentTrack++,1);
TVirtualGeoTrack *current = gGeoManager->GetTrack(track_index);
Int_t n = fShower->GetNumSteps();
for (Int_t i=0; i<n; i++) {
double x = (fShower->GetStep(i)->GetPosi().X() + fShower->GetStep(i)->GetPosf().X())/2/km;
double y = (fShower->GetStep(i)->GetPosi().Y() + fShower->GetStep(i)->GetPosf().Y())/2/km;
double z = (fShower->GetStep(i)->GetPosi().Z() + fShower->GetStep(i)->GetPosf().Z())/2/km;
double t = (fShower->GetStep(i)->GetTimei() + fShower->GetStep(i)->GetTimef())/2;
current->AddPoint(x,y,z,t);
}
}
//_____________________________________________________________________________
void ETreePainter::FillAtmoHistos( EAtmosphereHistoPainter *atmohistopainter ) {
//
// fill hitograms and ntuple
//
TH1F *th1 = BuildHistos( "Nph_Fluo","Nb of Fluoresence photons");
TH1F *th2 = BuildHistos( "Nph_Cer","Nb of Cerenkov photons");
TH1F *th3 = BuildHistos( "Zmax_Fluo","Zmax for Fluoresence");
TH1F *th4 = BuildHistos( "Zmax_Cer","Zmax for Cerenkov");
TH1F *th5 = BuildHistos( "Nmax_Fluo","Nb of Fluoresence photons at max");
TH1F *th6 = BuildHistos( "Nmax_Cer","Nb of Cerenkov photons at max");
TH1F *th11 = BuildHistos( "Nph_Fluo_P","Nb of Fluoresence photons on pupil");
TH1F *th12 = BuildHistos( "Nph_CB_P","Nb of Backscattered Cerenkov photons on pupil");
TH1F *th13 = BuildHistos( "Nph_CR_P","Nb of Reflected Cerenkov photons on pupil");
TH1F *th14 = BuildHistos( "Nph_tot_P","Nb of total photons on pupil");
TH1F *th15 = BuildHistos( "Tmax_Fluo_P","Time of max Fluoresence photons on pupil");
TH1F *th16 = BuildHistos( "Nmax_Fluo_P","Nb of Fluoresence photons on pupil at max");
TH1F *th17 = BuildHistos( "Sig_Fluo_P","Sigma Nb of Fluoresence photons on pupil");
TH1F* th(0);
TF1* mygaus = new TF1("mygaus","gaus",0,10);
Double_t par[3];
th = atmohistopainter->GetAltHisto("alt_Bunch_Nph_F");
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
fTrack.fNph_f = th->Integral();
fTrack.fNmax_f = par[0];
fTrack.fZmax_f = par[1];
th = atmohistopainter->GetAltHisto("alt_Bunch_Nph_C");
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
fTrack.fNph_c = th->Integral();
fTrack.fNmax_c = par[0];
fTrack.fZmax_c = par[1];
mygaus = new TF1("mygaus","gaus",atmohistopainter->GetTMin(),atmohistopainter->GetTMax());
th = atmohistopainter->GetTofHisto("Single_Nph_F_P_t");
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
fTrack.fNph_f_p = th->Integral();
fTrack.fNmax_f_p = par[0];
fTrack.fTmax_f_p = atmohistopainter->GetTMax()-par[1];
fTrack.fSig_f_p = par[2];
th = atmohistopainter->GetTofHisto("Single_Nph_CB_P_t");
fTrack.fNph_cb_p = th->Integral();
th = atmohistopainter->GetTofHisto("Single_Nph_CR_P_t");
fTrack.fNph_cr_p = th->Integral();
th1->Fill( fTrack.fNph_f/1e14 );
th2->Fill( fTrack.fNph_c/1e14 );
th3->Fill( fTrack.fZmax_f );
th4->Fill( fTrack.fZmax_c );
th5->Fill( fTrack.fNmax_f/1e12 );
th6->Fill( fTrack.fNmax_c/1e12 );
th11->Fill( fTrack.fNph_f_p );
th12->Fill( fTrack.fNph_cb_p );
th13->Fill( fTrack.fNph_cr_p );
th14->Fill( fTrack.fNph_f_p + fTrack.fNph_cb_p + fTrack.fNph_cr_p );
th15->Fill( fTrack.fTmax_f_p );
th16->Fill( fTrack.fNmax_f_p );
th17->Fill( fTrack.fSig_f_p );
}
//_____________________________________________________________________________
void ETreePainter::FillTruthHistos( ETruth *truth ) {
//
// fill truth in ntuple
//
TH1F *th01 = BuildHistos( "Energy","True Energy");
TH1F *th02 = BuildHistos( "Theta", "True Theta");
TH1F *th03 = BuildHistos( "Phi", "True Phi");
TH1F *th04 = BuildHistos( "Z1", "True Z1");
TH1F *th05 = BuildHistos( "Zmax", "True Zmax");
TH1F *th06 = BuildHistos( "X1", "True X1");
TH1F *th07 = BuildHistos( "Xmax", "True Xmax");
fTrack.fEnergy = truth->GetTrueEnergy(); // en MeV
fTrack.fTheta = truth->GetTrueTheta()*TMath::RadToDeg();
fTrack.fPhi = truth->GetTruePhi()*TMath::RadToDeg();
fTrack.fZ1 = truth->GetTrueInitPos().Z()/ km;
fTrack.fZmax = truth->GetTrueShowerMaxPos().Z()/km;
fTrack.fX1 = truth->GetTrueX1()/g*cm2;
fTrack.fXmax = truth->GetTrueShowerXMax()/g*cm2;
th01->Fill( fTrack.fEnergy /1e13 );
th02->Fill( fTrack.fTheta );
th03->Fill( fTrack.fPhi );
th04->Fill( fTrack.fZ1 );
th05->Fill( fTrack.fZmax );
th06->Fill( fTrack.fX1 );
th07->Fill( fTrack.fXmax );
}
//_____________________________________________________________________________
void ETreePainter::FillTree() {
//
// fill stree
//
TTree *stree = BuildTree();
stree->Fill();
}
//_____________________________________________________________________________
void ETreePainter::DrawHistos( Option_t *opt ) {
//
// Build Histos
//
TString option(opt);
TH1F* th = 0;
if (option == "Energy")
th = BuildHistos("Energy","");
else if (option == "Theta")
th = BuildHistos("Theta","");
else if (option == "Phi")
th = BuildHistos("Phi","");
else if (option == "Z1")
th = BuildHistos("Z1","");
else if (option == "Zmax")
th = BuildHistos("Zmax","");
else if (option == "X1")
th = BuildHistos("X1","");
else if (option == "Xmax")
th = BuildHistos("Xmax","");
else if (option == "Nph_Fluo")
th = BuildHistos("Nph_Fluo","");
else if (option == "Nph_Cer")
th = BuildHistos("Nph_Cer","");
else if (option == "Nmax_Fluo")
th = BuildHistos("Nmax_Fluo","");
else if (option == "Nmax_Cer")
th = BuildHistos("Nmax_Cer","");
else if (option == "Zmax_Fluo")
th = BuildHistos("Zmax_Fluo","");
else if (option == "Zmax_Cer")
th = BuildHistos("Zmax_Cer","");
else if (option == "Nph_Fluo_P")
th = BuildHistos("Nph_Fluo_P","");
else if (option == "Nmax_Fluo_P")
th = BuildHistos("Nmax_Fluo_P","");
else if (option == "Tmax_Fluo_P")
th = BuildHistos("Tmax_Fluo_P","");
else if (option == "Sig_Fluo_P")
th = BuildHistos("Sig_Fluo_P","");
else if (option == "Nph_CB_P")
th = BuildHistos("Nph_CB_P","");
else if (option == "Nph_CR_P")
th = BuildHistos("Nph_CR_P","");
else if (option == "Nph_tot_P")
th = BuildHistos("Nph_tot_P","");
else
th = NULL;
th->Draw();
}
//_____________________________________________________________________________
void ETreePainter::Draw( Option_t *opt ) {
fTop->Draw();
gGeoManager->DrawTracks();
}