Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

ETreePainter - source file

// $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(); 
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update 14-May-2005 21:31