Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

EShowerHistoPainter - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: EShowerHistoPainter.cc,v 1.10 2005/10/26 16:32:36 thea 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)

using namespace sou;

// 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);
        th->SetDirectory(0);
    }

    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);
        th->SetDirectory(0);
    }

    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 );

}
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