Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

EShowerPainter - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: EShowerPainter.cc,v 1.20 2005/06/01 00:16:26 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;
using namespace sou;

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;
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005