Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

EAtmosphereHistoPainter - source file

// $Id: EAtmosphereHistoPainter.cc,v 1.41 2005/11/14 10:24:38 moreggia Exp $
// Author: Anne Stutz   2004/10/19

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  Id: EAtmosphereHistoPainter                                                           *
 *  Package: <packagename>                                                   *
 *  Coordinator: <coordinator>                                               *
 *                                                                           *
 *****************************************************************************/

//_____________________________________________________________________________
//
// EAtmosphereHistoPainters
//
// Several groups of plots. Every plot of a group are generated together.
// Possible to build only a group of plots
//
//   Config file parameters
//   ======================
//
//   <parameter name>: <parameter description>
//   -Valid options: <available options>
//

#include "EAtmosphereHistoPainter.hh"
#include "EAtmosphere.hh"
#include "EShower.hh"
#include "EShowerStep.hh"
#include "EBunchPhotons.hh"
#include "ESystemOfUnits.hh"
#include "ETruth.hh"

#include <TH1F.h>
#include <TH2F.h>
#include <TF1.h>
#include <TTimer.h>
#include <TPad.h>
#include <TPaveText.h>
#include <TGaxis.h>
#include <TVector3.h>
#include <TMath.h>
#include <TLegend.h>

ClassImp(EAtmosphereHistoPainter)

extern Double_t Zv(const TVector3&);
using namespace TMath;
using namespace sou;

//_____________________________________________________________________________
 EAtmosphereHistoPainter::EAtmosphereHistoPainter(EAtmosphere* atm, Double_t GTUlength, ETruth* truth, EShower* show) {
    //
    // Constructor
    // ETruth needed to get shower earth impact
    //

    fShower = 0;
    fAtmosphere = 0;

    if ( !atm ) {
        MakeZombie();
        return;
    }
    fAtmosphere = new EAtmosphere(*atm);
    fTruth = new ETruth(*truth);
    if(show) fShower = new EShower(*show);
    fHistos = new TList();
    
    fGTUlength = GTUlength*ns;
    //fEUSO = TVector3(0*km,0*km,0); //GRNDetec
    fEUSO = TVector3(0,0,430*km);
    fRadius = 1250*mm;
    fTimeMin = 0.;
    fTimeMax = 0.;
    Info("EAtmosphereHistoPainter::Build()","USE pupil radius = %.2f mm, euso altitude = %.2f km, GTUlength = %.2f microsec",fRadius/mm,fEUSO.Z()/km,fGTUlength);

    fTimeBins = 0;
    fAltBins = 0;
    fLongBins = 0; 
    fGramBins = 0;
    // initialize fIndex_list
    Int_t ns = fAtmosphere->GetNumSingles();
    fMaxScatOrder = fAtmosphere->GetMaxScatOrder();
    Printf("Simu has been run with a max scattering order = %d",fMaxScatOrder);
    fIndex_list.Set(ns);
    for(Int_t i=0; i<ns; i++) fIndex_list[i] = i;
    
    BuildAll();
}

//_____________________________________________________________________________
 EAtmosphereHistoPainter::~EAtmosphereHistoPainter() {
    //
    // Destructor
    //

    if (IsZombie()) return;

    if ( fTimeBins ) delete[] fTimeBins;fTimeBins=0;
    if ( fAltBins ) delete[] fAltBins;fAltBins=0;
    if ( fLongBins ) delete[] fLongBins;fLongBins=0;
    if ( fGramBins ) delete[] fGramBins;fGramBins=0;
    if ( fHistos ) fHistos->Delete();
}

//_____________________________________________________________________________
 void EAtmosphereHistoPainter::BuildAll() {
    //
    // Build
    // 

    // histos of photons vs date of creation
    BuildDateHistos();
    BuildAltHistos();
    if(fShower) BuildLongitudinalHistos();
    BuildTofHistos();
    BuildWlHistos();
    BuildFovHistos();
    BuildAngularPupilHistos();
    BuildMScattHistos();
} 

//_____________________________________________________________________________
 void EAtmosphereHistoPainter::MScattCut(Double_t BDFmin, Double_t BDFmax, Double_t SNR) {
    //
    // defines a cut in space-time coordinates to get a more realistic view of Multiple Scattering contribution  
    // then replot all the histos in taking the cut into account
    //     - BDF min and max values can be configured (then linear interpolation w.r.t. distance to Nadir)
    //     - Signal To Noise ratio can be configured (SNR = S / BDF ---> NOT logarithmic here)
    //       otherwise, it is fixed to 3 sigmas of bdf fluctuations (poissonian assumed)
    //
    
    // if no data available
    if ( !fAtmosphere ) {
        Info("BuildDateHistos()","EAtmosphere object is NULL. Painter made zombie.");
        MakeZombie();
        return;
    }
    
    // determine the nb of GTU which allow to consider all the photons
    Int_t nbofGTU = 0;
    if(fTimeMin*fTimeMax <= 0 ) BuildAll();
    else nbofGTU = Int_t((fTimeMax-fTimeMin) / fGTUlength) + 1;    

    // delete all existing histograms     
    if(fHistos) fHistos->Delete();
    
    // list of histos (one for each GTU) + 1 histo integrated over all GTU (before the cut)
    TList tempHistos;
    TString prefix("temp");
    TString name(prefix);
    TH2F* temp = 0;
    for(Int_t i=0; i<nbofGTU; i++) {
        name = prefix;
	name += i;
	temp = new TH2F(name.Data(),name.Data(),625,-250,250,625,-250,250);
	tempHistos.Add(temp);
    }
    TH2F* integ = GetFovHisto("2D-FoV_nocut","FoV view, integrated over all GTU (NO CUT)");
    
    // fill histos
    TString histonb(prefix);
    Int_t ns = fAtmosphere->GetNumSingles();
    for (Int_t i=0; i<ns; i++) { 
        ESinglePhoton* s = fAtmosphere->GetSingle(i);
	if(s->IsAbsorbed()) continue;
        Double_t tof = (s->GetTof() + s->GetDate() - fTimeMin) / fGTUlength;
	Int_t gtunb  = Int_t(tof);
	Double_t x   = s->GetPos().X()/km;
	Double_t y   = s->GetPos().Y()/km;
	histonb = prefix;
	histonb += gtunb;
	temp = (TH2F*)tempHistos.FindObject(histonb.Data());
	if(temp) temp->Fill(x,y);
	else Printf("<MScattCut> 1. Wrong histo identification  %s",histonb.Data());
	integ->Fill(x,y);
    }
    
    // photon per photon analysis : does it make part of a "fired" pixel ?
    // if yes -> it is kept in memory for future plots
    Double_t bdf = 0.;
    Double_t sig_convol = 0.;
    Double_t theta_inc = 0.;
    Int_t bin = 0;
    Int_t listpos = 0;
    for (Int_t i=0; i<ns; i++) { 
        ESinglePhoton* s = fAtmosphere->GetSingle(i);
	if(s->IsAbsorbed()) continue;
        Double_t tof = (s->GetTof() + s->GetDate() - fTimeMin) / fGTUlength;
	Int_t gtunb  = Int_t(tof);
	Double_t x   = s->GetPos().X()/km;
	Double_t y   = s->GetPos().Y()/km;
	
	// bdf and SNR definition
	// bdf interpolation vs. position in the FoV
	bdf = BDFmax - (sqrt(x*x + y*y)) / 250. * (BDFmax - BDFmin);
	if(SNR <= 0) SNR = 3 / sqrt(bdf);
        
	// signal convolution factor (by an "ideal detector") ---> translated in pe-
	if(fEUSO.Z()>0) theta_inc = atan(sqrt(x*x + y*y) * km / fEUSO.Z());
	else Printf("<MScattCut> CUT not available for ground detector");
	// parametrized euso optics
	//sig_convol = 3.86e-2 * (1 - 0.258/1.581*theta_inc*theta_inc - 2.131/1.581*theta_inc*theta_inc*theta_inc*theta_inc); 
	// only cos(theta) effect
	sig_convol = 3.86e-2 * cos(theta_inc); 
	
	// then fill relevant histos
	histonb = prefix;
	histonb += gtunb;
	temp = (TH2F*)tempHistos.FindObject(histonb.Data());
	if(temp) {
	    bin = temp->FindBin(x,y);
	    if( (temp->GetBinContent(bin) * sig_convol / bdf) >= SNR ) fIndex_list[listpos++] = i;
	}
	else Printf("<MScattCut> 2. Wrong histo identification  %s",histonb.Data());
	integ->Fill(x,y);
    }
    
    // fill the rest of index list with flag values (-1)
    for(Int_t i=listpos; i<fIndex_list.GetSize(); i++) fIndex_list[i] = -1;
    
    tempHistos.Delete();
    
    // rebuild all the histograms with the defined cut
    BuildAll();
}

//_____________________________________________________________________________
 void EAtmosphereHistoPainter::BuildDateHistos() {
    //
    // generate histos of photons at creation as function of creation time
    //
    
    // if no data available
    if ( !fAtmosphere ) {
        Info("BuildDateHistos()","EAtmosphere object is NULL. Painter made zombie.");
        MakeZombie();
        return;
    }
        
    Double_t NphMax = 0; 
    Int_t n = fAtmosphere->GetNumBunches();

    // if photons data are empty
    if ( n <= 0 ) {
        Info("Build()","No Bunches in Atmosphere ( NumBunch = 0 ). Painter made zombie.");
        MakeZombie();
        return;
    }

    // determine bins for photons histos as function of creation time
    Double_t temp[n+1];
    UInt_t Nbins = 0;
    for (Int_t i=0; i<n; i++) {
        EBunchPhotons* b = fAtmosphere->GetBunch(i);
        if (b->GetType() == 0 ) {
            temp[Nbins] = b->GetDatei()/ms;  
            temp[Nbins+1] = b->GetDatef()/ms;  
            if( NphMax < b->GetWeight() ) NphMax = b->GetWeight();
            Nbins++;
        }
    }
    // if only cerenkov simulation
    if(!Nbins) {
	for (Int_t i=0; i<n; i++) {
            EBunchPhotons* b = fAtmosphere->GetBunch(i);
	    if (b->GetType() == 1 ) {
                temp[Nbins] = b->GetDatei()/ms;  
                temp[Nbins+1] = b->GetDatef()/ms;  
                if( NphMax < b->GetWeight() ) NphMax = b->GetWeight();
                Nbins++;
	    }
	}
    }

    Printf("<BuildDateHistos> %d bins in time to draw",Nbins);
    
    // array for X-axis
    fTimeBins = new Double_t[Nbins+1];
    for (UInt_t i=0; i<Nbins+1; i++) {
        fTimeBins[i] = temp[i];
    }
    
    // build histograms from bunches vs date of creation 
    TH1F *th0 = GetTimeHisto("Bunch_Nph", "Nph vs time",Nbins);
    TH1F *th1 = GetTimeHisto("Bunch_Nph_F", "Nph_Fluo vs time",Nbins);
    TH1F *th2 = GetTimeHisto("Bunch_Nph_C", "Nph_Cerenkov vs time",Nbins);
    TH1F *th3 = GetTimeHisto("Bunch_Yield_F","Fluo Yield vs time",Nbins);
    TH1F *th4 = GetTimeHisto("Bunch_Yield_C","Cerenkov Yield vs time",Nbins);
    TH1F *th5 = GetTimeHisto("Bunch_Nph_F_SA", "Nph_Fluo_in_Solide_Angle vs time",Nbins);
    th0->SetYTitle("Nph");
    th1->SetYTitle("Nph Fluo");
    th2->SetYTitle("Nph Cerenkov");
    th3->SetYTitle("Nph/m Fluo");
    th4->SetYTitle("Nph/m Cerenkov");
    th5->SetYTitle("Nph Fluo");

    
    // build histograms from singles in Detector solid angle _BEFORE LAST TRANSMISSION_ vs date of creation 
    TH1F *th33 = GetTimeHisto("Single_Nph_F", "Nph_Fluo_in_Solide_Angle vs time",Nbins);    
    TH1F *th35 = GetTimeHisto("Single_Nph_CB", "Nph_Cerenkov_backscttered_in_Solide_Angle vs time",Nbins);
    TH1F *th35_1 = GetTimeHisto("Single_Nph_ckov_AeroScat", "Nph_Cerenkov_AeroScat_in_Solide_Angle vs time",Nbins);
    TH1F *th35_2 = GetTimeHisto("Single_Nph_ckov_CloudBackscat", "Nph_Cerenkov_CloudBackscat_in_Solide_Angle vs time",Nbins);
    TH1F *th35_tot = GetTimeHisto("Single_Nph_CB_tot", "Nph_Cerenkov_Total_backscattered_in_Solide_Angle vs time",Nbins);  // not used so far
    
    th33->SetYTitle("Nph Fluo");
    th35->SetYTitle("Nph Cerenkov backscattered");
    th35_1->SetYTitle("Nph Cerenkov backscattered by aerosols");
    th35_2->SetYTitle("Nph Cerenkov backscattered by clouds");
    th35_tot->SetYTitle("Nph Cerenkov total backscattered");
    
    // fill histograms from bunches
    for (Int_t i=0; i<n; i++) { 
        EBunchPhotons* b = fAtmosphere->GetBunch(i);
        Double_t time =  (b->GetDatei() + b->GetDatef()) / 2 / ms;
        Double_t Nph =    b->GetWeight(); 
        Double_t Yield =  b->GetYield()*m;

	th0->Fill(time,Nph); 

        if (b->GetType()==0) {
	   // fluorescence bunches
           th1->Fill(time,Nph);  
           th3->Fill(time,Yield);

           TVector3 dist = EUSO() - (b->GetShowerPosi()+ b->GetShowerPosf())*0.5;
           Double_t omega = TMath::Pi() *DetectorRadius()*DetectorRadius() *cos( dist.Theta() ) / dist.Mag2();
           th5->Fill(time,Nph*omega/(4*TMath::Pi()));
        } 

        if (b->GetType()==1) { 
	   // Cerenkov bunches
           th2->Fill(time,Nph); 
           th4->Fill(time,Yield);
        }  
    } 
    
    
    // fill histograms from singles
    for (Int_t i=0; i<fIndex_list.GetSize(); i++) {
        if(fIndex_list[i] < 0) break;
	ESinglePhoton* s = fAtmosphere->GetSingle(fIndex_list[i]);
        Double_t time =  s->GetDate()/ms;

        if (s->GetType() == 0 ) { 
           // for fluorescence photons
	   th33->Fill(time,1); 
        }
        if (s->GetType() == 1 ) { 
	    // for Backscattered Cerenkov photons 
            if (s->GetHistory()==2) {
		th35->Fill(time,1);
		th35_tot->Fill(time,1);
            }
            if ( s->GetHistory()==4) {
		th35_1->Fill(time,1);
		th35_tot->Fill(time,1);
            }
            if ( s->GetHistory()==3) {
		th35_2->Fill(time,1);
		th35_tot->Fill(time,1);
            }
        }
    }
}

//_____________________________________________________________________________
 void EAtmosphereHistoPainter::BuildAltHistos() {
    //
    // generate histos of YIELD as function of altitude
    //
    
    // if no data available
    if ( !fAtmosphere ) {
        Info("BuildAltHistos()","EAtmosphere object is NULL. Painter made zombie.");
        MakeZombie();
        return;
    }
        
    Int_t n = fAtmosphere->GetNumBunches();

    // if photons data are empty
    if ( n <= 0 ) {
        Info("Build()","No Bunches in Atmosphere ( NumBunch = 0 ). Painter made zombie.");
        MakeZombie();
        return;
    }

    // determine bins for development w.r.t altitude
    Double_t temp[n+1];
    Int_t Nbins = 0;
    for (Int_t i=0; i<n; i++) {
        EBunchPhotons* b = fAtmosphere->GetBunch(i);
        if (b->GetType() == 0 ) {
            temp[Nbins] = Zv(b->GetShowerPosi())/km;
	    temp[Nbins+1] = Zv(b->GetShowerPosf())/km; 
	    Nbins++;
        }
    }
    // if only cerenkov simulation
    if(!Nbins) {
        for (Int_t i=0; i<n; i++) {
            EBunchPhotons* b = fAtmosphere->GetBunch(i);
            if (b->GetType() == 1 ) {
                temp[Nbins] = Zv(b->GetShowerPosi())/km;
	        temp[Nbins+1] = Zv(b->GetShowerPosf())/km; 
	        Nbins++;
            }
        }
    }
    Printf("<BuildAltHistos> %d bins in time to draw",Nbins);
    
    // fill X-axis array
    fAltBins = new Double_t[Nbins+1];
    Int_t incrm = 0;
    for (Int_t i=0; i<Nbins+1; i++) {
        fAltBins[Nbins - incrm] =  temp[i]; 
	incrm++; 
    }
    
    
    // build histo vs altitude
    TH1F *alt0 = GetAltHisto("alt_Bunch_Nph", "Nph vs altitude",Nbins);
    TH1F *alt1 = GetAltHisto("alt_Bunch_Nph_F", "Nph_Fluo vs altitude",Nbins);
    TH1F *alt2 = GetAltHisto("alt_Bunch_Nph_C", "Nph_Cerenkov vs altitude",Nbins);
    TH1F *alt3 = GetAltHisto("alt_Bunch_Yield_F","Fluo Yield vs altitude",Nbins);
    TH1F *alt4 = GetAltHisto("alt_Bunch_Yield_C","Cerenkov Yield vs altitude",Nbins);
    
    TH1F* hs1 = GetSingleAltHisto("alt_directfluo","Direct fluo and cerenkov v.s. altitude (not transmitted)"); 
    TH1F* hs2 = GetSingleAltHisto("alt_directckov","Direct fluo and cerenkov v.s. altitude (not transmitted)"); 
    TH1F* hs3 = GetSingleAltHisto("alt_scatfluo","Scattered fluo and cerenkov v.s. altitude (not transmitted)"); 
    TH1F* hs4 = GetSingleAltHisto("alt_scatckov","Scattered fluo and cerenkov v.s. altitude (not transmitted)"); 
    TH1F* hs5 = GetSingleAltHisto("alt_airscat","Air Scattered photons v.s. altitude (not transmitted)"); 
    TH1F* hs6 = GetSingleAltHisto("alt_cloudscat","Clouds Scattered photons v.s. altitude (not transmitted)"); 
    TH1F* hs7 = GetSingleAltHisto("alt_aeroscat","Aerosols Scattered photons v.s. altitude (not transmitted)"); 
    TH1F* hs11   = GetSingleAltHisto("alt_fluo_Ntrans","Fluo transmission v.s. altitude"); 
    TH1F* hs11_t = GetSingleAltHisto("alt_fluo_trans","Fluo transmission v.s. altitude"); 
    TH1F* hs12   = GetSingleAltHisto("alt_ckov_Ntrans","Cerenkov transmission v.s. altitude"); 
    TH1F* hs12_t = GetSingleAltHisto("alt_ckov_trans","Cerenkov transmission v.s. altitude"); 
    TH1F* trans1 = GetSingleAltHisto("alt_trans_fluo","Fluo transmission v.s. altitude");      // transmission value within [0,1]
    TH1F* trans2 = GetSingleAltHisto("alt_trans_ckov","Cerenkov transmission v.s. altitude");  // transmission value within [0,1]
    
    TH2F* tottrans_c   = GetAlt2DHisto("tottrans_c","Total transmission v.s. altitude");
    TH2F* tottrans_f   = GetAlt2DHisto("tottrans_f","Total transmission v.s. altitude");
    TH2F* rayltrans_c  = GetAlt2DHisto("rayltrans_c","Rayleigh transmission v.s. altitude");
    TH2F* rayltrans_f  = GetAlt2DHisto("rayltrans_f","Rayleigh transmission v.s. altitude");
    TH2F* ozonetrans_c = GetAlt2DHisto("ozonetrans_c","Ozone transmission v.s. altitude");
    TH2F* ozonetrans_f = GetAlt2DHisto("ozonetrans_f","Ozone transmission v.s. altitude");
    TH2F* aerotrans_c  = GetAlt2DHisto("aerotrans_c","Aerosols transmission v.s. altitude");
    TH2F* aerotrans_f  = GetAlt2DHisto("aerotrans_f","Aerosols transmission v.s. altitude");
    TH2F* cloudtrans_c = GetAlt2DHisto("cloudtrans_c","Clouds transmission v.s. altitude");
    TH2F* cloudtrans_f = GetAlt2DHisto("cloudtrans_f","Clouds transmission v.s. altitude");
    
    alt0->SetYTitle("Nph");
    alt1->SetYTitle("Nph Fluo");
    alt2->SetYTitle("Nph Cerenkov");
    alt3->SetYTitle("Nph/m Fluo");
    alt4->SetYTitle("Nph/m Cerenkov");
    
    
    // fill histograms from bunches
    for (Int_t i=0; i<n; i++) { 
	EBunchPhotons* b = fAtmosphere->GetBunch(i);
	Double_t alt =   ( Zv(b->GetShowerPosi())  
                         + Zv(b->GetShowerPosf()) )/2/km; 
        Double_t Yield = b->GetYield()*m;
        Double_t Nph =   b->GetWeight(); 

	alt0->Fill(alt,Nph);
        if (b->GetType()==0) {
	    alt3->Fill(alt,Yield);
	    alt1->Fill(alt,Nph);
	}
        if (b->GetType()==1) {
	    alt2->Fill(alt,Nph);
	    alt4->Fill(alt,Yield); 
	}
    }
    
    // fill histo with singlephotons attributes
    for (Int_t i=0; i<fIndex_list.GetSize(); i++) { 
        if(fIndex_list[i] < 0) break;
        ESinglePhoton* s = fAtmosphere->GetSingle(fIndex_list[i]);
	Double_t alt = Zv(s->GetPos())/km;
	Double_t totaltrans = s->GetLastTrans("tot");
	Double_t rayltrans = s->GetLastTrans("rayl");
	Double_t ozonetrans = s->GetLastTrans("ozone");
	Double_t aerotrans = s->GetLastTrans("aero");
	Double_t cloudtrans = s->GetLastTrans("cloud");
	// fluo
	if(s->GetType() == 0) {
	    tottrans_f->Fill(alt,totaltrans);
	    rayltrans_f->Fill(alt,rayltrans);
	    ozonetrans_f->Fill(alt,ozonetrans);
	    aerotrans_f->Fill(alt,aerotrans);
	    cloudtrans_f->Fill(alt,cloudtrans);
	    hs11->Fill(alt);
	    if(!s->IsAbsorbed())           hs11_t->Fill(alt);
	    if(s->GetHistory() == 0)       hs1->Fill(alt);
	    else if(s->GetHistory() == 2)  {
	        hs3->Fill(alt);
		hs5->Fill(alt);
	    }
	    else if(s->GetHistory() == 3)  hs6->Fill(alt);
	    else if(s->GetHistory() == 4)  hs7->Fill(alt);
	}
	// ckov
	else if(s->GetType() == 1) {
	    tottrans_c->Fill(alt,totaltrans);
	    rayltrans_c->Fill(alt,rayltrans);
	    ozonetrans_c->Fill(alt,ozonetrans);
	    aerotrans_c->Fill(alt,aerotrans);
	    cloudtrans_c->Fill(alt,cloudtrans);
	    hs12->Fill(alt);
	    if(!s->IsAbsorbed())           hs12_t->Fill(alt);
	    if(s->GetHistory() == 0)       hs2->Fill(alt);
	    else if(s->GetHistory() == 2)  {
	        hs4->Fill(alt);
		hs5->Fill(alt);
	    }
	    else if(s->GetHistory() == 3)  hs6->Fill(alt);
	    else if(s->GetHistory() == 4)  hs7->Fill(alt);
	}
    }
    
    // fill transmission histos
    for(Int_t i=0; i<trans1->GetNbinsX();i++) {
        if(hs11->GetBinContent(i+1) > 10) trans1->SetBinContent(i+1,hs11_t->GetBinContent(i+1)/hs11->GetBinContent(i+1));
	else trans1->SetBinContent(i+1,0);
        if(hs12->GetBinContent(i+1) > 10) trans2->SetBinContent(i+1,hs12_t->GetBinContent(i+1)/hs12->GetBinContent(i+1));
	else trans2->SetBinContent(i+1,0);
    }
}

//_____________________________________________________________________________
 void EAtmosphereHistoPainter::BuildLongitudinalHistos() {
    //
    // generate histos of photons at creation along the track (Longitudinal Profiles)
    // 
    // generate also histos for SinglePhotons : longitudinal and lateral distribution along track (in km only)
    //
    
    // if no atmosphere data available
    if ( !fAtmosphere ) {
        Info("BuildLongitudinalHistos()","EAtmosphere object is NULL. Painter made zombie.");
        MakeZombie();
        return;
    }
    // if no shower data available
    if ( !fShower ) {
        Info("BuildLongitudinalHistos()","EShower object is NULL. It is needed for this method ");
        //MakeZombie();
        return;
    }
        
    Int_t n = fAtmosphere->GetNumBunches();
    Int_t nsh = fShower->GetNumSteps();
    if(n != 2*nsh) Warning("BuildLongitudinalHistos()","Nb of bunches and Nb of showersteps SHOULD BE THE SAME");

    // if photons data are empty
    if ( n <= 0 ) {
        Info("Build()","No Bunches in Atmosphere ( NumBunch = 0 ). Painter made zombie.");
        MakeZombie();
        return;
    }

    // determine bins for longitudinal development (both in km and g/cm2)
    UInt_t Nbins = UInt_t(nsh);
    UInt_t Nbins_dist = 0;
    Double_t temp[Nbins+1];
    fGramBins = new Double_t[Nbins+1];
    TVector3 initpos(1,0,0);
    Bool_t flag = false;  // turn true when starts light track
    EShowerStep* shstep = 0;
    for (Int_t i=0; i<nsh; i++) {
        shstep = fShower->GetStep(i);
	if(fAtmosphere->GetBunch(2*i)->GetFate()!=2 || fAtmosphere->GetBunch(2*i+1)->GetFate()!=2) {
	    if(!flag) {
	        initpos = shstep->GetPosi();
		flag = true;
            }
	    temp[Nbins_dist] = (shstep->GetPosi() - initpos).Mag()/km;  
	    temp[Nbins_dist+1] = (shstep->GetPosf() - initpos).Mag()/km;
	    Nbins_dist++;
	}
	fGramBins[i] = shstep->GetXi()*cm2/g;
	fGramBins[i+1] = shstep->GetXf()*cm2/g;
    }
    fLongBins = new Double_t[Nbins_dist+1];
    for(UInt_t m=0; m<Nbins_dist+1; m++) fLongBins[m] = temp[m];

    // get shower geometry for SinglePhotons histos
    TVector3 first_interac = fShower->GetInitPos();
    TVector3 earth_impact  = fTruth->GetTrueEarthImpact();
    TVector3 track = earth_impact - first_interac;

    // build histo along track
    TH1F *h10 = GetLongitudinalHisto("Fluo_Longit_prof_km","Fluo_Longitudinal_profile_km",false,Nbins_dist);
    TH1F *h11 = GetLongitudinalHisto("Ckov_Longit_prof_km","Ckov_Longit_prof_km",false,Nbins_dist);
    TH1F *h20 = GetLongitudinalHisto("Fluo_Longit_prof_gram","Fluo_Longit_prof_gram",true,Nbins);
    TH1F *h21 = GetLongitudinalHisto("Ckov_Longit_prof_gram","Ckov_Longit_prof_gram",true,Nbins);
    TH1F *hs1 = GetSingleLongitudinalHisto("fluo_direct_longi","Direct fluo and cerenkov photons along track (not transmitted)",track.Mag()/km);
    TH1F *hs2 = GetSingleLongitudinalHisto("ckov_direct_longi","Direct fluo and cerenkov photons along track (not transmitted)",track.Mag()/km);
    TH1F *hs3 = GetSingleLongitudinalHisto("fluo_scat_longi","Scattered fluo and cerenkov photons along track (not transmitted)",track.Mag()/km);
    TH1F *hs4 = GetSingleLongitudinalHisto("ckov_scat_longi","Scattered fluo and cerenkov photons along track (not transmitted)",track.Mag()/km);
    TH1F *hs11 = GetSingleLateralHisto("fluo_direct_lat","Direct fluo and cerenkov photons lateral distribution (not transmitted)");
    TH1F *hs12 = GetSingleLateralHisto("ckov_direct_lat","Direct fluo and cerenkov photons lateral distribution (not transmitted)");
    TH1F *hs13 = GetSingleLateralHisto("fluo_scat_lat","Scattered fluo and cerenkov photons lateral distribution (not transmitted)");
    TH1F *hs14 = GetSingleLateralHisto("ckov_scat_lat","Scattered fluo and cerenkov photons lateral distribution (not transmitted)");
    TH2F *hs21 = GetSingle2DHisto("fluo_direct_2D","Direct fluo and cerenkov photons 2D distribution along track (not transmitted)",track.Mag()/km);
    TH2F *hs22 = GetSingle2DHisto("ckov_direct_2D","Direct fluo and cerenkov photons 2D distribution along track (not transmitted)",track.Mag()/km);
    TH2F *hs23 = GetSingle2DHisto("fluo_scat_2D","Scattered fluo and cerenkov photons 2D distribution along track (not transmitted)",track.Mag()/km);
    TH2F *hs24 = GetSingle2DHisto("ckov_scat_2D","Scattered fluo and cerenkov photons 2D distribution along track (not transmitted)",track.Mag()/km);
   
    h10->SetYTitle("Fluo Longitudinal profile in km");
    h20->SetYTitle("Fluo Longitudinal profile in g/cm2");
    h11->SetYTitle("Cerenkov Longitudinal profile in km");
    h21->SetYTitle("Cerenkov Longitudinal profile in g/cm2");
    
    // build one histo per interaction order
    TString prefix("lateral_nbinter");
    TString name(prefix);
    TH1F* temphisto = 0;
    for(Int_t i=0; i<=fMaxScatOrder; i++) {
        name = prefix;
	name += i;
	temphisto = GetSingleLateralHisto(name.Data(),"Lateral distribution - scattering order contribution");
    }
    
    
    // fill histograms from fluo bunches
    for (Int_t i=0; i<n; i++) {
        EBunchPhotons* b = fAtmosphere->GetBunch(i);
	if(!(i%2)) shstep =     fShower->GetStep(i/2);
	Double_t dist =      ((0.5*(shstep->GetPosf() + shstep->GetPosi())) - initpos).Mag()/km;
	Double_t grammage =  (0.5*(shstep->GetXf() + shstep->GetXi()))*cm2/g;
        Double_t nph =       b->GetWeight();
	UChar_t fate = b->GetFate();  // if bunch weight not too small

        if (b->GetType()==0) {
	    if(fate !=2) h10->Fill(dist,nph);
	    h20->Fill(grammage,nph);
        }
        if (b->GetType()==1) {
	    if(fate !=2) h11->Fill(dist,nph);
	    h21->Fill(grammage,nph);
        }
    }
    
    // fill SinglePhoton histograms
    TString histonb(prefix);
    for (Int_t i=0; i<fIndex_list.GetSize(); i++) { 
        if(fIndex_list[i] < 0) break;
        ESinglePhoton* s    = fAtmosphere->GetSingle(fIndex_list[i]);
	TVector3 temp       = s->GetPos() - first_interac;
	Double_t longimag   = temp.Dot(track.Unit())/km;                          // projection along track
	Double_t lateralmag = (temp - longimag*km*track.Unit()).Mag()/km;  // orthogonal to the track
	
	// for MScatt details
	Int_t nbinter  = s->GetNbinter();
	histonb = prefix;
	histonb += nbinter;
	temphisto = (TH1F*)fHistos->FindObject(histonb.Data());
	if(temphisto) temphisto->Fill(lateralmag);
	else Printf("<BuildTofHistos> Wrong histo identification  %s",histonb.Data());
	
	// fluo
	if(s->GetType() == 0) {
	    if(s->GetHistory() == 0) {
	        hs1->Fill(longimag);
	        hs11->Fill(lateralmag);
	        hs21->Fill(longimag,lateralmag);
	    }
	    else if(s->GetHistory() == 2) {
	        hs3->Fill(longimag);
	        hs13->Fill(lateralmag);
	        hs23->Fill(longimag,lateralmag);
	    }
	}
	// ckov
	else if(s->GetType() == 1) {
	    if(s->GetHistory() == 0) {
	        hs2->Fill(longimag);
	        hs12->Fill(lateralmag);
	        hs22->Fill(longimag,lateralmag);
	    }
	    else if(s->GetHistory() == 2) {
	        hs4->Fill(longimag);
	        hs14->Fill(lateralmag);
	        hs24->Fill(longimag,lateralmag);
	    }
	}
    }
}

//_____________________________________________________________________________
 void EAtmosphereHistoPainter::BuildTofHistos() {
    //
    // generate histos of photons on pupil as function of time since PRIMARY cosmic ray first interaction
    //
    
    // if no data available
    if ( !fAtmosphere ) {
        Info("BuildTofHistos()","EAtmosphere object is NULL. Painter made zombie.");
        MakeZombie();
        return;
    }
        
    // determine bins for time on pupil
    if(!(fTimeMin*fTimeMax)) { fTimeMin = 3000.*microsecond;
        fTimeMax = 0;
        for (Int_t i=0; i<fIndex_list.GetSize(); i++) { 
            if(fIndex_list[i] < 0) break;
            ESinglePhoton* s = fAtmosphere->GetSingle(fIndex_list[i]);
            if ( !(s->IsAbsorbed()) ) {
                if ( fTimeMin > (s->GetTof()+s->GetDate()) )
                    fTimeMin = s->GetTof()+s->GetDate();
                if ( fTimeMax < (s->GetTof()+s->GetDate()) )
                    fTimeMax = s->GetTof()+s->GetDate();
            }
        }
    }
    
    // build histogram vs time on pupil
    TH1F *th20 = GetTofHisto("Single_Nph_P_t","Nph_Total_on_Pupil vs time");
    TH1F *th24 = GetTofHisto("Single_Nph_F_P_t","Direct fluo on pupil vs time");
    TH1F *th25 = GetTofHisto("Single_Nph_direct_ckov_P_t","Direct cerenkov  on pupil vs time");
    TH1F *th26 = GetTofHisto("Single_Nph_CB_P_t","Cerenkov air scattered on pupil vs time");
    TH1F *th36 = GetTofHisto("Single_Nph_fluo_airscat_P_t","Fluo air scattered on pupil vs time");
    TH1F *th26_1 = GetTofHisto("Single_Nph_ckov_AeroScat_P_t","Cerenkov aerosol scattered on pupil vs time");
    TH1F *th36_1 = GetTofHisto("Single_Nph_fluo_AeroScat_P_t","Fluo aerosol scattered on pupil vs time");
    TH1F *th26_2 = GetTofHisto("Single_Nph_ckov_CloudBackscat_P_t","Cerenkov clouds scattered on pupil vs time");
    TH1F *th36_2 = GetTofHisto("Single_Nph_fluo_CloudBackscat_P_t","Fluo clouds scattered on pupil vs time");
    TH1F *th26_tot = GetTofHisto("Single_Nph_ckov_Total_Backscat_P_t","Cerenkov total scattered on pupil vs time");
    TH1F *th36_tot = GetTofHisto("Single_Nph_fluo_Total_Backscat_P_t","Fluo total scattered on pupil vs time");
    TH1F *th28 = GetTofHisto("Single_Nph_CR_P_t","Nph_Cerenkov refl on pupil vs time");
    TH1F *th38 = GetTofHisto("Single_Nph_fluo_refl_P_t","Nph_fluo refl on pupil vs time");
    
    th20->SetYTitle("Nph ");
    th24->SetYTitle("Nph Fluo");
    th25->SetYTitle("Nph direct ckov");
    th26->SetYTitle("Nph Airscattered Cerenkov");
    th36->SetYTitle("Nph Airscattered Fluo");
    th26_1->SetYTitle("Nph AeroScattered Cerenkov");
    th26_2->SetYTitle("Nph CloudBackscattered Cerenkov");
    th26_tot->SetYTitle("Nph total Backscattered Cerenkov");
    th28->SetYTitle("Nph Refelected Cerenkov");
    th38->SetYTitle("Nph Refelected Fluo");
    
       
    // build one histo per interaction order
    TString prefix("pupil_nbinter");
    TString name(prefix);
    TH1F* temp = 0;
    for(Int_t i=0; i<=fMaxScatOrder; i++) {
        name = prefix;
	name += i;
	temp = GetTofHisto(name.Data(),"Photons on pupil - scattering order contribution");
    }
    
    
    
    // fill histograms from singles
    TString histonb(prefix);
    for (Int_t i=0; i<fIndex_list.GetSize(); i++) { 
        if(fIndex_list[i] < 0) break;
        ESinglePhoton* s = fAtmosphere->GetSingle(fIndex_list[i]);
	
	// for MScatt details
	Int_t nbinter  = s->GetNbinter();
	histonb = prefix;
	histonb += nbinter;
	temp = (TH1F*)fHistos->FindObject(histonb.Data());
        Double_t tof = (s->GetTof() + s->GetDate() - fTimeMin) / fGTUlength;
        
	// 
	if ( !(s->IsAbsorbed()) )  {
	    if(temp) temp->Fill(tof);
	    else Printf("<BuildTofHistos> Wrong histo identification  %s",histonb.Data());
	    th20->Fill(tof,1);
            if ( s->GetType() == 0 ) {
        	if(s->GetHistory()==0) th24->Fill(tof,1);
		if(s->GetHistory()==2) th36->Fill(tof,1);
		if(s->GetHistory()==1) th38->Fill(tof,1);
		if(s->GetHistory()==4) th36_1->Fill(tof,1);
		if(s->GetHistory()==3) th36_2->Fill(tof,1);
        	if(s->GetHistory() > 0) th36_tot->Fill(tof,1);
	    }
            if ( s->GetType() == 1 ) { 
                if ( s->GetHistory()==0) {
		    th25->Fill(tof,1);
		}
                if ( s->GetHistory()==2) {
        	   th26->Fill(tof,1);
        	   th26_tot->Fill(tof,1);
                }
                if ( s->GetHistory()==4) {
                   th26_1->Fill(tof,1);
                   th26_tot->Fill(tof,1);
                }
                if ( s->GetHistory()==3) {
                   th26_2->Fill(tof,1);
                   th26_tot->Fill(tof,1);
                }
                if ( s->GetHistory()==1 )
                   th28->Fill(tof,1);
            }
	}
    }
}

//_____________________________________________________________________________
 void EAtmosphereHistoPainter::BuildWlHistos() {
    //
    // generate wavelength spectra of photons on pupil BEFORE and AFTER last transmission
    // generate also wavelength spectra of photons in atmosphere at creation
    //

    // if no data available
    if ( !fAtmosphere ) {
        Info("BuildWlHistos()","EAtmosphere object is NULL. Painter made zombie.");
        MakeZombie();
        return;
    }
        
    Int_t n = fAtmosphere->GetNumBunches();

    // if photons data are empty
    if ( n <= 0 ) {
        Info("Build()","No Bunches in Atmosphere ( NumBunch = 0 ). Painter made zombie.");
        MakeZombie();
        return;
    }
    
    //build histogram vs wavelenght
    TH1F *th11 = GetWlHisto("Bunch_Wl_F", "Nph_F vs Lambda");
    TH1F *th12 = GetWlHisto("Bunch_Wl_C", "Nph_C vs Lambda");

    TH1F *th15 = GetWlHisto("Single_Wl_CB", "Nph_Cerenkov_back_in_Solid_Angle vs Lambda");
    TH1F *th15_tot = GetWlHisto("Single_Wl_CB_tot", "Nph_Cerenkov_Total_back_in_Solid_Angle vs Lambda");
    TH1F *th16 = GetWlHisto("Single_Wl_CB_P", "Nph_Cerenkov_back_on_pupil vs Lambda");
    TH1F *th16_1 = GetWlHisto("Single_Wl_ckov_AeroScat_P", "Nph_Cerenkov_AeroScat_on_pupil vs Lambda");
    TH1F *th16_2 = GetWlHisto("Single_Wl_ckov_CloudBackscat_P", "Nph_Cerenkov_Cloudback_on_pupil vs Lambda");
    TH1F *th17 = GetWlHisto("Single_Wl_CR","Nph_Cerenkov refl in Solid Angle vs lambda");
    TH1F *th18 = GetWlHisto("Single_Wl_CR_P","Nph_Cerenkov refl on pupil vs lambda");
    
    TH1F *th13 = GetWlHisto("Single_Wl_F", "Fluo transmission vs Lambda");
    TH1F *th14 = GetWlHisto("Single_Wl_F_P", "Fluo transmission vs Lambda");
    TH1F *hs2 = GetWlHisto("wl_ckov_Ntrans","Cerenkov transmission vs Lambda"); // all the cerenkov
    TH1F *hs2_t = GetWlHisto("wl_ckov_trans","Cerenkov transmission vs Lambda");
    TH1F *trans1 = GetWlHisto("wl_trans_fluo","Cerenkov transmission");
    TH1F *trans2 = GetWlHisto("wl_trans_ckov","Cerenkov transmission");
    
    TH2F* tottrans_c   = GetWl2DHisto("WLtottrans_c","Total transmission v.s. Lambda");
    TH2F* tottrans_f   = GetWl2DHisto("WLtottrans_f","Total transmission v.s. Lambda");
    TH2F* rayltrans_c  = GetWl2DHisto("WLrayltrans_c","Rayleigh transmission v.s. Lambda");
    TH2F* rayltrans_f  = GetWl2DHisto("WLrayltrans_f","Rayleigh transmission v.s. Lambda");
    TH2F* ozonetrans_c = GetWl2DHisto("WLozonetrans_c","Ozone transmission v.s. Lambda");
    TH2F* ozonetrans_f = GetWl2DHisto("WLozonetrans_f","Ozone transmission v.s. Lambda");
    TH2F* aerotrans_c  = GetWl2DHisto("WLaerotrans_c","Aerosols transmission v.s. Lambda");
    TH2F* aerotrans_f  = GetWl2DHisto("WLaerotrans_f","Aerosols transmission v.s. Lambda");
    TH2F* cloudtrans_c = GetWl2DHisto("WLcloudtrans_c","Clouds transmission v.s. Lambda");
    TH2F* cloudtrans_f = GetWl2DHisto("WLcloudtrans_f","Clouds transmission v.s. Lambda");
    
    
    th11->SetYTitle("Nph Fluo");
    th12->SetYTitle("Nph Cerenkov");
    th13->SetYTitle("Nph Fluo");
    th14->SetYTitle("Nph Fluo");
    th15->SetYTitle("Nph Cerenkov");
    th15_tot->SetYTitle("Nph total Cerenkov");
    th16->SetYTitle("Nph Cerenkov");
    th16_1->SetYTitle("Nph aerosol Cerenkov");
    th16_2->SetYTitle("Nph clouds Cerenkov");
    th17->SetYTitle("Nph Cerenkov");
    th18->SetYTitle("Nph Cerenkov");
    
    // fill histograms from bunches
    for (Int_t i=0; i<n; i++) { 
        EBunchPhotons* b = fAtmosphere->GetBunch(i);
        Double_t Nph =      b->GetWeight(); 
        Double_t NumWl =    b->GetNumWavelengths();
        const Float_t* lambda = 0;
        const Float_t* wlweight = 0;
        wlweight = b->GetTable("weight");
        lambda = b->GetTable("lambda");

        if (b->GetType()==0) {
	   // fluorescence bunches
           for ( Int_t j=0; j<NumWl; j++ ) {
	       th11->Fill( lambda[j]/nm , wlweight[j]*Nph );  
	   }
        } 

        if (b->GetType()==1) { 
	   // Cerenkov bunches
           for ( Int_t j=0; j<NumWl; j++ ) {
	       th12->Fill( lambda[j]/nm , wlweight[j]*Nph );  
	   }
        }  
    } 

    // fill histograms from singles
    for (Int_t i=0; i<fIndex_list.GetSize(); i++) { 
        if(fIndex_list[i] < 0) break;
        ESinglePhoton* s = fAtmosphere->GetSingle(fIndex_list[i]);
        Double_t wl = s->GetWl()/nm;
	Double_t totaltrans = s->GetLastTrans("tot");
	Double_t rayltrans = s->GetLastTrans("rayl");
	Double_t ozonetrans = s->GetLastTrans("ozone");
	Double_t aerotrans = s->GetLastTrans("aero");
	Double_t cloudtrans = s->GetLastTrans("cloud");
        // for fluorescence photons
        if (s->GetType() == 0 ) { 
	    tottrans_f->Fill(wl,totaltrans);
	    rayltrans_f->Fill(wl,rayltrans);
	    ozonetrans_f->Fill(wl,ozonetrans);
	    aerotrans_f->Fill(wl,aerotrans);
	    cloudtrans_f->Fill(wl,cloudtrans);
            th13->Fill(wl,1);
            if ( !(s->IsAbsorbed()) ) th14->Fill(wl,1);
        }
        // for cerenkov photons
        if (s->GetType() == 1 ) {
	    tottrans_c->Fill(wl,totaltrans);
	    rayltrans_c->Fill(wl,rayltrans);
	    ozonetrans_c->Fill(wl,ozonetrans);
	    aerotrans_c->Fill(wl,aerotrans);
	    cloudtrans_c->Fill(wl,cloudtrans);
	    hs2->Fill(wl);
	    if(!(s->IsAbsorbed())) hs2_t->Fill(wl);  
	    // for rayleigh Cerenkov photons 
            if (s->GetHistory()==2) {
		th15->Fill(wl,1);
		th15_tot->Fill(wl,1);
        	if ( !(s->IsAbsorbed()) ) th16->Fill(wl,1);
            }
            if (s->GetHistory()==4) {
		th15_tot->Fill(wl,1);
        	if (!(s->IsAbsorbed())) th16_1->Fill(wl,1);
            }
            if (s->GetHistory()==3) {
		th15_tot->Fill(wl,1);
        	if (!(s->IsAbsorbed())) th16_2->Fill(wl,1);
            }
            // for reflected Cerenkov photons
            if (s->GetHistory()==1 ) {
		th17->Fill(wl,1);
        	if (!(s->IsAbsorbed())) th18->Fill(wl,1);
            }
        }
    }
    
    // fill transmission histos
    for(Int_t i=0; i<trans1->GetNbinsX();i++) {
        if(th13->GetBinContent(i+1) > 10) trans1->SetBinContent(i+1,th14->GetBinContent(i+1)/th13->GetBinContent(i+1));
	else trans1->SetBinContent(i+1,0);
        if(hs2->GetBinContent(i+1) > 10) trans2->SetBinContent(i+1,hs2_t->GetBinContent(i+1)/hs2->GetBinContent(i+1));
	else trans2->SetBinContent(i+1,0);
    }

}

//_____________________________________________________________________________
 void EAtmosphereHistoPainter::BuildFovHistos() {
    //
    // generate histos of photons on pupil as function of time since PRIMARY cosmic ray first interaction
    //
    
    // if no data available
    if ( !fAtmosphere ) {
        Info("BuildTofHistos()","EAtmosphere object is NULL. Painter made zombie.");
        MakeZombie();
        return;
    }
        
    // build histogram vs time on pupil
    TH2F *th = GetFovHisto("2D-FoV","FoV view, integrated over all GTU (WITH CUT)");
    
    // fill histograms from singles
    for (Int_t i=0; i<fIndex_list.GetSize(); i++) { 
        if(fIndex_list[i] < 0) break;
        ESinglePhoton* s = fAtmosphere->GetSingle(fIndex_list[i]);
	if(s->IsAbsorbed()) continue;
	Double_t x     = s->GetPos().X()/km;
	Double_t y     = s->GetPos().Y()/km;
	th->Fill(x,y);
    }
}

//_____________________________________________________________________________
 void EAtmosphereHistoPainter::BuildAngularPupilHistos() {
    //
    // generate histos of photons on pupil as function of their incident angle (w.r.t vertical axis)
    //
    
    // if no data available
    if ( !fAtmosphere ) {
        Info("BuildAngularPupilHistos()","EAtmosphere object is NULL. Painter made zombie.");
        MakeZombie();
        return;
    }
            
    // build histograms
    TH1F *th20 = GetAngularHisto("Single_Nph_P_t_angul","Photons on pupil vs entering angle");
    TH1F *th24 = GetAngularHisto("Single_Nph_F_P_t_angul","Direct fluo on pupil vs entering angle");
    TH1F *th25 = GetAngularHisto("Single_Nph_direct_ckov_P_t_angul","Direct cerenkov  on pupil vs entering angle");
    TH1F *th26 = GetAngularHisto("Single_Nph_CB_P_t_angul","Cerenkov air scattered on pupil vs entering angle");
    TH1F *th36 = GetAngularHisto("Single_Nph_fluo_airscat_P_t_angul","Fluo air scattered on pupil vs entering angle");
    TH1F *th26_1 = GetAngularHisto("Single_Nph_ckov_AeroScat_P_t_angul","Cerenkov aerosol scattered on pupil vs entering angle");
    TH1F *th36_1 = GetAngularHisto("Single_Nph_fluo_AeroScat_P_t_angul","Fluo aerosol scattered on pupil vs entering angle");
    TH1F *th26_2 = GetAngularHisto("Single_Nph_ckov_CloudBackscat_P_t_angul","Cerenkov clouds scattered on pupil vs entering angle");
    TH1F *th36_2 = GetAngularHisto("Single_Nph_fluo_CloudBackscat_P_t_angul","Fluo clouds scattered on pupil vs entering angle");
    TH1F *th26_tot = GetAngularHisto("Single_Nph_ckov_Total_Backscat_P_t_angul","Cerenkov total scattered on pupil vs entering angle");
    TH1F *th36_tot = GetAngularHisto("Single_Nph_fluo_Total_Backscat_P_t_angul","Fluo total scattered on pupil vs entering angle");
    TH1F *th28 = GetAngularHisto("Single_Nph_CR_P_t_angul","Nph_Cerenkov refl on pupil vs entering angle");
    TH1F *th38 = GetAngularHisto("Single_Nph_fluo_refl_P_t_angul","Nph_fluo refl on pupil vs entering angle");
    
    TH2F *th2D = Get2DAngularHisto("angul2D_pupil","Angular localisation at pupil");
           
    // build one histo per interaction order
    TString prefix("angular_nbinter");
    TString name(prefix);
    TH1F* temp = 0;
    for(Int_t i=0; i<=fMaxScatOrder; i++) {
        name = prefix;
	name += i;
	temp = GetAngularHisto(name.Data(),"Photons on pupil - scattering order contribution");
    }
    
    // fill histograms from singles
    TVector3 axis(1);
    if(fEUSO.Z() == 0) axis.SetXYZ(0,0,-1);
    else axis.SetXYZ(0,0,1);
    TVector3 dir(1);
    TString histonb(prefix);
    for (Int_t i=0; i<fIndex_list.GetSize(); i++) { 
        if(fIndex_list[i] < 0) break;
        ESinglePhoton* s = fAtmosphere->GetSingle(fIndex_list[i]);
	
	// for MScatt details
	Int_t nbinter  = s->GetNbinter();
	histonb = prefix;
	histonb += nbinter;
	temp = (TH1F*)fHistos->FindObject(histonb.Data());
	
	dir = (fEUSO - s->GetPos()).Unit();
	Double_t theta = dir.Angle(axis)*RadToDeg();
	Double_t phi   = ATan(s->GetPos().Y() / s->GetPos().X())*RadToDeg();
        
	// 
	if ( !(s->IsAbsorbed()) )  {
	    th2D->Fill(phi,Cos(theta*DegToRad()));
	    if(temp) temp->Fill(theta);
	    else Printf("<BuildAngularPupilHistos> Wrong histo identification  %s",histonb.Data());
	    th20->Fill(theta,1);
            if ( s->GetType() == 0 ) {
        	if(s->GetHistory()==0) th24->Fill(theta,1);
		if(s->GetHistory()==2) th36->Fill(theta,1);
		if(s->GetHistory()==1) th38->Fill(theta,1);
		if(s->GetHistory()==4) th36_1->Fill(theta,1);
		if(s->GetHistory()==3) th36_2->Fill(theta,1);
        	if(s->GetHistory() > 0) th36_tot->Fill(theta,1);
	    }
            if ( s->GetType() == 1 ) { 
                if ( s->GetHistory()==0) {
		    th25->Fill(theta,1);
		}
                if ( s->GetHistory()==2) {
        	   th26->Fill(theta,1);
        	   th26_tot->Fill(theta,1);
                }
                if ( s->GetHistory()==4) {
                   th26_1->Fill(theta,1);
                   th26_tot->Fill(theta,1);
                }
                if ( s->GetHistory()==3) {
                   th26_2->Fill(theta,1);
                   th26_tot->Fill(theta,1);
                }
                if ( s->GetHistory()==1 )
                   th28->Fill(theta,1);
            }
	}
    }
}

//_____________________________________________________________________________
 void EAtmosphereHistoPainter::BuildMScattHistos() {
    //
    // Multiple Scattering details
    //
    
    // if no data available
    if ( !fAtmosphere ) {
        Info("BuildMScattHistos()","EAtmosphere object is NULL. Painter made zombie.");
        MakeZombie();
        return;
    }
    
    Bool_t addir = TH1::AddDirectoryStatus();
    TH1::AddDirectory(kFALSE);
    // build histos
    // nbinter
    TH1F *th = new TH1F("nbinter","Photon nb of interactions ",fMaxScatOrder+1,0,fMaxScatOrder+1);
    TH1F *th2 = new TH1F("nbinter_trans","Photon nb of interactions ",fMaxScatOrder+1,0,fMaxScatOrder+1);
    TH1F *ratio = new TH1F("nbinter_trans_ratio","Photon nb of interactions ",fMaxScatOrder+1,0,fMaxScatOrder+1);
    fHistos->Add(th);
    fHistos->Add(th2);
    fHistos->Add(ratio);
    
    // omegadiff - build one histo per interaction order
    TString prefix("omegadiff");
    TString name(prefix);
    TH1F* temp = 0;
    for(Int_t i=0; i<=fMaxScatOrder; i++) {
        name = prefix;
	name += i;
	//temp = new TH1F(name.Data(),"#cbar #Omega_{scatpos} - #Omega_{showerpos} #cbar / #Omega_{showerpos}",1000,0,1);
	temp = new TH1F(name.Data(),"Log10(Omega) distribution",100,-20,1);  //DELETE
	//temp = new TH1F(name.Data(),"Distance to detector distribution (km)",100000,0,100);  //DELETE
	fHistos->Add(temp);
    }
    
    
    
    // fill histos
    // nbinter, before and after transmission
    for (Int_t i=0; i<fIndex_list.GetSize(); i++) { 
        if(fIndex_list[i] < 0) break;
        ESinglePhoton* s = fAtmosphere->GetSingle(fIndex_list[i]);
	th->Fill(s->GetNbinter()+0.1);
	if(!s->IsAbsorbed()) th2->Fill(s->GetNbinter());
    }
    for(Int_t i=0; i<th->GetNbinsX();i++) {
        if(th->GetBinContent(i+1) > 1) ratio->SetBinContent(i+1,th2->GetBinContent(i+1)/th->GetBinContent(i+1));
	else ratio->SetBinContent(i+1,0);
    }
    
    // omega difference for scattered photons only
    // between showerpos and scatpos              //TOFIX : must be adapted if more precise (but CPU expensive) MScatt algorithm
    TString histonb(prefix);
    for (Int_t i=0; i<fIndex_list.GetSize(); i++) { 
        if(fIndex_list[i] < 0) break;
        ESinglePhoton* s = fAtmosphere->GetSingle(fIndex_list[i]);
	Int_t nbinter  = s->GetNbinter();
	if(nbinter == 0) continue; // only scattered photons
	histonb = prefix;
	histonb += nbinter;
	temp = (TH1F*)fHistos->FindObject(histonb.Data());
	// OMEGA      abs(scat - shower) / shower    becomes
	// DISTANCE2  abs(scat - shower) / scat
	//if(temp) temp->Fill( fabs((s->GetShowerPos() - fEUSO).Mag2() - (s->GetPos() - fEUSO).Mag2()) / (s->GetPos() - fEUSO).Mag2() );
	//DELETE
	//if(temp) temp->Fill(Log10(TwoPi() * (1 - (s->GetPos() - fEUSO).Mag() / Sqrt((s->GetPos() - fEUSO).Mag2() + fRadius*fRadius) ) * Cos((s->GetPos() - fEUSO).Theta())));
	if(temp) temp->Fill(Log10(TwoPi() * (1 - (s->GetPos() - fEUSO).Mag() / Sqrt((s->GetPos() - fEUSO).Mag2() + fRadius*fRadius) )));
	//if(temp) temp->Fill( (s->GetPos() - fEUSO).Mag() / km);
	else Printf("<BuildMScattHistos> Wrong histo identification  %s",histonb.Data());
    }
        
    
    // build a summary of photons interactions in the atmosphere
    TH1F* Reflected = GetMScattHistory("reflected","Summary of photons hisotry in atmosphere");
    TH1F* Rayleigh = GetMScattHistory("rayleigh","Summary of photons hisotry in atmosphere");
    TH1F* Clouds = GetMScattHistory("clouds","Summary of photons hisotry in atmosphere");
    TH1F* Aerosols = GetMScattHistory("aerosols","Summary of photons hisotry in atmosphere");
    TH1F* norm = GetMScattHistory("norm","norm");
    
    Char_t val = 0;
    for(Int_t j=0;j<fMaxScatOrder;j++) {
        for (Int_t i=0; i<fIndex_list.GetSize(); i++) { 
            if(fIndex_list[i] < 0) break;
            ESinglePhoton* s = fAtmosphere->GetSingle(i);
	    val = s->GetAllHistory(j);
	    switch(val) {
		case 0 : break;
		case 1 : Reflected->Fill(j+1.5); break;
		case 2 : Rayleigh->Fill(j+1.5); break;
		case 3 : Clouds->Fill(j+1.5); break;
		case 4 : Aerosols->Fill(j+1.5); break;
		default: Printf("nTH1::AddDirectory(addir);
}

//_____________________________________________________________________________
 TH1F *EAtmosphereHistoPainter::GetTimeHisto(const char *name, const char *title, UInt_t Nbins) {
    //
    // Get histograms versus time
    //

    if ( title == 0 ) title = name;
    TH1F* th = (TH1F*)fHistos->FindObject(name);

    if (!th && Nbins) {
        th = new TH1F( name, title, Nbins, fTimeBins);
        fHistos->Add(th);
        th->SetStats(0);
        th->SetXTitle("Development Time, #mus");
        th->SetFillColor(9);
        th->SetDirectory(0);
    }
    if(!th && !Nbins) Warning("GetTimeHisto()","Nb of bins must be specified");

    return th;
}

//_____________________________________________________________________________
 TH1F *EAtmosphereHistoPainter::GetAltHisto(const char *name, const char *title, UInt_t Nbins) {
    //
    // Get histograms versus altitude
    //

    if ( title == 0 ) title = name;
    TH1F* th = (TH1F*)fHistos->FindObject(name);

    if (!th && Nbins) {
        th = new TH1F( name, title, Nbins, fAltBins);
        fHistos->Add(th);
        th->SetStats(0);
        th->SetXTitle("Altitude (km)");
        th->SetFillColor(9);
        th->SetDirectory(0);
    }
    else if(!th && !Nbins) Warning("GetAltHisto()","Nb of bins must be specified");

    return th;
}

//_____________________________________________________________________________
 TH1F *EAtmosphereHistoPainter::GetSingleAltHisto(const char *name, const char *title) {
    //
    // Get SinglePhoton histograms versus altitude
    //

    if ( title == 0 ) title = name;
    TH1F* th = (TH1F*)fHistos->FindObject(name);
    if (!th) {
        th = new TH1F(name,title,150,0,30);  // 200m binwidth
        fHistos->Add(th);
        th->SetStats(0);
        th->SetXTitle("Altitude (km)");
        th->SetDirectory(0);
    }

    return th;
}

//_____________________________________________________________________________
 TH2F *EAtmosphereHistoPainter::GetAlt2DHisto(const char *name, const char *title) {
    //
    // Get SinglePhoton histograms versus altitude
    //

    if ( title == 0 ) title = name;
    TH2F* th = (TH2F*)fHistos->FindObject(name);
    if (!th) {
        th = new TH2F(name,title,300,0,30,1000,0,1.1);  // 100m binwidth
        fHistos->Add(th);
        th->SetXTitle("Altitude (km)");
        th->SetDirectory(0);
    }

    return th;
}

//_____________________________________________________________________________
 TH1F *EAtmosphereHistoPainter::GetLongitudinalHisto(const char *name, const char *title, Bool_t flag, UInt_t Nbins) {
    //
    // Get longitudinal profiles
    // if flag == true  -> in g/cm2
    // if flag == false -> in km
    //

    if ( title == 0 ) title = name;
    TH1F* th = (TH1F*)fHistos->FindObject(name);

    if (!th && Nbins) {
        if(flag) th = new TH1F( name, title, Nbins, fGramBins);
	if(!flag) th = new TH1F( name, title, Nbins, fLongBins);
        fHistos->Add(th);
        th->SetStats(0);
        if(flag) th->SetXTitle("Depth along track (g/cm2)");
        if(!flag) th->SetXTitle("Length along VISIBLE track (km)");
        th->SetFillColor(9);
        th->SetDirectory(0);
    }
    else if(!th && !Nbins) Warning("GetLongitudinalHisto()","Nb of bins must be specified");

    return th;
}

//_____________________________________________________________________________
 TH1F *EAtmosphereHistoPainter::GetSingleLongitudinalHisto(const char *name, const char *title, Double_t max) {
    //
    // Get SinglePhoton histograms along the track (in km)
    //

    if ( title == 0 ) title = name;
    TH1F* th = (TH1F*)fHistos->FindObject(name);
    
    if (!th) {
        th = new TH1F( name, title, Int_t(max / 0.2) + 3, 0,max);  // binwidth = 200m
        fHistos->Add(th);
        th->SetStats(0);
        th->SetXTitle("Length along track from first interaction (km)");
        th->SetFillColor(9);
        th->SetDirectory(0);
    }

    return th;
}

//_____________________________________________________________________________
 TH1F *EAtmosphereHistoPainter::GetSingleLateralHisto(const char *name, const char *title) {
    //
    // Get SinglePhoton distance from the track (in km)
    //

    if ( title == 0 ) title = name;
    TH1F* th = (TH1F*)fHistos->FindObject(name);
    
    if (!th) {
        th = new TH1F( name, title, 200, 0,20);  // binwidth = 100, max at 20km
        fHistos->Add(th);
        th->SetStats(0);
        th->SetXTitle("Distance from shower axis (km)");
        th->SetFillColor(9);
        th->SetDirectory(0);
    }

    return th;
}

//_____________________________________________________________________________
 TH1F *EAtmosphereHistoPainter::GetMScattHistory(const char *name, const char *title) {
    //
    // Summary of photons interactions in the atmosphere as function of scattering order
    //

    if ( title == 0 ) title = name;
    TH1F* th = (TH1F*)fHistos->FindObject(name);
        
    if (!th) {
        th = new TH1F(name,title,fMaxScatOrder,1,fMaxScatOrder+1);
        fHistos->Add(th);
        th->SetStats(0);
        th->SetXTitle("Scattering order");
        th->SetYTitle("% of each component at given order");
        th->SetDirectory(0);
    }

    return th;
}

//_____________________________________________________________________________
 TH2F *EAtmosphereHistoPainter::GetSingle2DHisto(const char *name, const char *title, Double_t max) {
    //
    // Get SinglePhoton 2D-histograms (longitudinal, lateral) coordinates
    //

    if ( title == 0 ) title = name;
    TH2F* th = (TH2F*)fHistos->FindObject(name);
    
    if (!th) {
        th = new TH2F( name, title, Int_t(max / 0.1) + 3, 0,max, 200, 0,20);
        fHistos->Add(th);
        th->SetStats(0);
        th->SetXTitle("Longitudinal distance (km)");
        th->SetYTitle("Lateral distance (km)");
        th->SetFillColor(9);
        th->SetDirectory(0);
    }

    return th;
}

//_____________________________________________________________________________
 TH1F *EAtmosphereHistoPainter::GetWlHisto(const char *name, const char *title) {
    //
    // Get wavelength histograms 
    //

    if ( title == 0 ) title = name;
    TH1F* th = (TH1F*)fHistos->FindObject(name);

    if (!th) {
      th = new TH1F( name, title, 150, 300., 450.);
        fHistos->Add(th);
        th->SetStats(0);
        th->SetXTitle("Wavelength, nm");
        th->SetFillColor(9);
        th->SetDirectory(0);
    }

    return th;
}

//_____________________________________________________________________________
 TH2F *EAtmosphereHistoPainter::GetWl2DHisto(const char *name, const char *title) {
    //
    // Get wavelength histograms 
    //

    if ( title == 0 ) title = name;
    TH2F* th = (TH2F*)fHistos->FindObject(name);

    if (!th) {
      th = new TH2F( name, title, 150, 300., 450., 1000,0,1.1);
        fHistos->Add(th);
        th->SetXTitle("Wavelength, nm");
        th->SetFillColor(9);
        th->SetDirectory(0);
    }

    return th;
}

//_____________________________________________________________________________
 TH1F *EAtmosphereHistoPainter::GetTofHisto(const char *name, const char *title) {
    //
    // Get wavelength histograms 
    //

    if ( title == 0 ) title = name;
    TH1F* th = (TH1F*)fHistos->FindObject(name);

    Int_t nbins = (Int_t)((fTimeMax-fTimeMin)/(fGTUlength)) + 20;    

    if (!th && fTimeMin) {
      th = new TH1F( name, title, nbins, 0, Float_t(nbins)*(fGTUlength/microsecond));
        fHistos->Add(th);
        th->SetStats(0);
        th->SetXTitle("Time  (GTU)");
        th->SetFillColor(9);
        th->SetDirectory(0);
    }
    else if(!th && !fTimeMin) Warning("GetTofHisto()","fTimeMin must be specified");

    return th;
}

//_____________________________________________________________________________
 TH2F *EAtmosphereHistoPainter::GetFovHisto(const char *name, const char *title) {
    //
    // Get 2D-histograms in (x,y) MES-coordinates
    //

    if ( title == 0 ) title = name;
    TH2F* th = (TH2F*)fHistos->FindObject(name);
    
    if (!th) {
        th = new TH2F( name, title,625,-250,250,625,-250,250);  // 0.8x0.8 cellwidth
        fHistos->Add(th);
        th->SetStats(0);
        th->SetXTitle("X (km)");
        th->SetYTitle("Y (km)");
        th->SetFillColor(9);
        th->SetDirectory(0);
    }

    return th;
}

//_____________________________________________________________________________
 TH1F *EAtmosphereHistoPainter::GetAngularHisto(const char *name, const char *title) {
    //
    // Get angular distribution of photon direction entering pupil
    //

    if ( title == 0 ) title = name;
    TH1F* th = (TH1F*)fHistos->FindObject(name);

    if (!th) {
      th = new TH1F( name, title,180,0,90);
        fHistos->Add(th);
        th->SetStats(0);
        th->SetXTitle("#theta  (#circ)");
        th->SetDirectory(0);
    }

    return th;
}

//_____________________________________________________________________________
 TH2F *EAtmosphereHistoPainter::Get2DAngularHisto(const char *name, const char *title) {
    //
    // Get 2D angular distribution of photon direction entering pupil
    // [ phi,cos(theta) ] coord.
    //

    if ( title == 0 ) title = name;
    TH2F* th = (TH2F*)fHistos->FindObject(name);

    if (!th) {
      th = new TH2F( name, title,270,-180,180,270,0,1);
        fHistos->Add(th);
        th->SetStats(0);
        th->SetYTitle("cos(#theta)");
        th->SetXTitle("#phi");
        th->SetDirectory(0);
    }

    return th;
}

//______________________________________________________________________________
 void EAtmosphereHistoPainter::Draw(Option_t *opt) {
    //
    //
    //
    
    if ( IsZombie() ) return;
    TString option(opt);
    TH1F* th(0);
    TH1F* th1(0);
    TH1F* th2(0);
    TH1F* th3(0);
    TH1F* th4(0);
    TH1F* th5(0);
    TH1F* th6(0);
    TH1F* th7(0);
    TH1F* th8(0);
    TH1F* th9(0);
    TH1F* th10(0);
    TH1* hcopy(0);
    TH2F* hcopy2D(0);
    TH2F* th2D(0);
    TH2F* prof1(0);
    TH2F* prof2(0);
    Double_t TimeMin = 0;
    if ( option == "" ) option = "Nph";


// ******** CREATION TIME ************* //
    // Bunches are plotted VS creation time
    // FLUO and CKOV on the same plot
    if ( option == "Nph" ) {
        th2 = GetTimeHisto("Bunch_Nph_C");
        if ( th2 ) {
            hcopy = (TH1F*)th2->DrawCopy();
            hcopy->SetFillColor(39);
        }       
        th1 = GetTimeHisto("Bunch_Nph_F");        
        if ( th1 ) {
           hcopy = (TH1F*)th1->DrawCopy("SAME"); 
           hcopy->SetFillColor(29);
        }  
    } 

    // FLUO Bunches are plotted VS creation time
    else if ( option == "Nph_F" ) {
        th = GetTimeHisto("Bunch_Nph_F");        
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(29);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           sprintf(txt,"Total:%e",th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
        }  
    }
    
    // CKOV Bunches are plotted VS creation time
    else if ( option == "Nph_C" ) {
        th = GetTimeHisto("Bunch_Nph_C");        
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(39);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           sprintf(txt,"Total:%e",th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
        } 
    }
  
    // FLUO in Solid angle (Singles are plotted)
    else if ( option == "Nph_F_SA" ) {
        th = GetTimeHisto("Bunch_Nph_F_SA");        
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(29);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           sprintf(txt,"Total:%e",th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
        }  
    }
    
    // FLUO Yield
    else if ( option == "Yield_F" ) {
        th = GetTimeHisto("Bunch_Yield_F");
        if ( th ) {
	   th->SetMinimum(0.9*th->GetMinimum());
           hcopy = (TH1F*)th->DrawCopy();
           hcopy->SetFillColor(29);
	}
    }            
    
    // CKOV Yield
    else if ( option == "Yield_C" ) {
        th = GetTimeHisto("Bunch_Yield_C");  
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
           hcopy->SetFillColor(39);
	}
    } 
    
    // FLUO BEFORE AND AFTER last transmission VS creation time
    else if ( option == "Nphs_F" ) {
        th = GetTimeHisto("Single_Nph_F"); 
        th1 = GetTimeHisto("Single_Nph_F_P");       
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(29);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(29);
           sprintf(txt,"Total:%d",(int)th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
        
           if ( th1 ) {
              hcopy = (TH1F*)th1->DrawCopy("SAME"); 
              hcopy->SetFillColor(2);
              t->SetTextColor(2);
              sprintf(txt,"Total:%d",(int)th1->Integral());
              t->DrawText(0.005,th1->GetMaximum()*0.94,txt);
          }  
        }  
    }

    // reflected CKOV photons BEFORE and AFTER propagation VS creation time
    else if ( option == "Nphs_CR" ) {
        th = GetTimeHisto("Single_Nph_CR"); 
        th1 = GetTimeHisto("Single_Nph_CR_P");       
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(39);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(39);
           sprintf(txt,"Total:%d",(int)th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
        
           if ( th1 ) {
              hcopy = (TH1F*)th1->DrawCopy("SAME"); 
              hcopy->SetFillColor(50);
              t->SetTextColor(50);
              sprintf(txt,"Total:%d",(int)th1->Integral());
              t->DrawText(0.005,th1->GetMaximum()*0.94,txt);
          }  
        }  
    }
    
    // scattered CKOV photons BEFORE and AFTER last transmission VS creation time
    else if ( option == "Nphs_CB" ) {
	th = GetTimeHisto("Single_Nph_CB_tot");                 // total ckov backscat before last transmission
        th1 = GetTimeHisto("Single_Nph_CB_P");                  // air backscat on pupil
        th3 = GetTimeHisto("Single_Nph_ckov_AeroScat_P");       // aerosols backscat on pupil
        th2 = GetTimeHisto("Single_Nph_ckov_CloudBackscat_P");  // clouds backscat on pupil
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(39);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(39);
           sprintf(txt,"Total:%d",(int)th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
           if ( th1 ) {
              hcopy = (TH1F*)th1->DrawCopy("SAME"); 
              hcopy->SetFillColor(50);
              t->SetTextColor(50);
              sprintf(txt,"Total:%d",(int)th1->Integral());
              t->DrawText(0.005,th1->GetMaximum()*0.94,txt);
           }
           if ( th3 ) {
              hcopy = (TH1F*)th3->DrawCopy("SAME"); 
              hcopy->SetFillColor(kBlue);
              t->SetTextColor(kBlue);
              sprintf(txt,"Total:%d",(int)th3->Integral());
              t->DrawText(0.005,th3->GetMaximum()*0.94,txt);
           }
           if ( th2 ) {
              hcopy = (TH1F*)th2->DrawCopy("SAME"); 
              hcopy->SetFillColor(5);
              t->SetTextColor(5);
              sprintf(txt,"Total:%d",(int)th2->Integral());
              t->DrawText(0.005,th2->GetMaximum()*0.94,txt);
          }  
  
        } 
    }
    
    
    
  
// ************************************ TIME ON PUPIL ***************************************//
    
    
    
    // Fluorescence photons vs time on pupil
    else if ( option == "Nphs_F_t" ) { 
        th = GetTofHisto("Single_Nph_F_P_t");
	TimeMin = (th->GetBinLowEdge(0)+1);
        Double_t tx = TimeMin+10;       
        if ( th ) {
	    Double_t TimeMax = (th->GetBinLowEdge(th->GetXaxis()->GetLast())+th->GetBinWidth(th->GetXaxis()->GetLast()));
            TF1* mygaus = new TF1("mygaus","gaus",TimeMin,TimeMax);
	    Double_t par[3];
	    th->Fit("mygaus","R");
	    mygaus->GetParameters(par);
// 	    Double_t binwidth = th->GetBinWidth(th->FindBin(par[1]));
            hcopy = (TH1F*)th->DrawCopy(""); 
            hcopy->SetFillColor(2);
            TText *t = new TText();
            char txt[200];
            t->SetTextFont(62);  
            t->SetTextSize(0.06);
            t->SetTextColor(2);
            sprintf(txt,"Total:%f",th->Integral());
            t->DrawText(tx,th->GetMaximum()*0.94,txt);
            t->SetTextColor(1);
	    sprintf(txt,"DTmax : %f microsec",par[1]-TimeMin);
            t->DrawText(tx,th->GetMaximum()*0.84,txt);
	    sprintf(txt,"Nph max (per GTU unit) : %f",par[0]); // integrated within 1GTU
            t->DrawText(tx,th->GetMaximum()*0.74,txt);
        }    
    }
    
    // direct and reflected CKOV vs time on pupil
    else if ( option == "Nphs_CR_t" ) {
        th1 = GetTofHisto("Single_Nph_CR_P_t"); 
        th2 = GetTofHisto("Single_Nph_direct_ckov_P_t"); 

	TimeMin = (th1->GetBinLowEdge(0)+1)*microsecond;
	if(TimeMin > (th2->GetBinLowEdge(0)+1)*microsecond) TimeMin = (th2->GetBinLowEdge(0)+1)*microsecond;
        Double_t tx = TimeMin/microsecond+10;
        if ( th1 ) {
           hcopy = (TH1F*)th1->DrawCopy(""); 
           hcopy->SetFillColor(50);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(50);
           sprintf(txt,"Reflected ckov:%d",(int)th1->Integral());
           if(th1->Integral()) t->DrawText(tx,th1->GetMaximum()*0.94,txt);       
        }  
        if ( th2 ) {
           hcopy = (TH1F*)th2->DrawCopy(""); 
           hcopy->SetFillColor(40);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(40);
           sprintf(txt,"Direct Cerenkov:%d",(int)th2->Integral());
           if(th2->Integral()) t->DrawText(tx,th2->GetMaximum()*0.94,txt);       
        }  
    }
    
    // direct fluo and ckov vs time on pupil
    else if ( option == "direct_on_pupil" ) {
        th1 = GetTofHisto("Single_Nph_F_P_t"); 
        th2 = GetTofHisto("Single_Nph_direct_ckov_P_t"); 

	TimeMin = (th1->GetBinLowEdge(0)+1)*microsecond;
	if(TimeMin > (th2->GetBinLowEdge(0)+1)*microsecond) TimeMin = (th2->GetBinLowEdge(0)+1)*microsecond;
        Double_t tx = TimeMin/microsecond+10;
        if ( th1 ) {
           hcopy = (TH1F*)th1->DrawCopy(""); 
           hcopy->SetFillColor(0);
           hcopy->SetLineColor(2);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(2);
           sprintf(txt,"Direct fluo:%d",(int)th1->Integral());
           if(th1->Integral()) t->DrawText(tx,th1->GetMaximum()*0.94,txt);       
        }  
        if ( th2 ) {
	   Double_t Ymax = hcopy->GetMaximum();
	   if(th2->GetMaximum() > Ymax) hcopy->SetMaximum(1.1*th2->GetMaximum());
           hcopy = (TH1F*)th2->DrawCopy("same"); 
           hcopy->SetLineColor(4);
           hcopy->SetFillColor(0);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(4);
           sprintf(txt,"Direct Cerenkov:%d",(int)th2->Integral());
           if(th2->Integral()) t->DrawText(tx,th2->GetMaximum()*0.94,txt);       
        }  
    }
    
    // scattered CKOV vs time on pupil
    else if ( option == "Nphs_CB_t" ) {
        th = GetTofHisto("Single_Nph_ckov_Total_Backscat_P_t"); // total
	th2 = GetTofHisto("Single_Nph_CB_P_t");                 // air scattered
	th3 = GetTofHisto("Single_Nph_ckov_AeroScat_P_t");      // aerosols scattered
        th1 = GetTofHisto("Single_Nph_ckov_CloudBackscat_P_t"); // cloud scattered
	
	TimeMin = (th->GetBinLowEdge(0)+1)*microsecond;
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(39);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(39);
           sprintf(txt,"Total:%d",(int)th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
           if ( th1 ) {
              hcopy = (TH1F*)th1->DrawCopy("SAME"); 
              hcopy->SetFillColor(5);
              t->SetTextColor(5);
              sprintf(txt,"Cloud scattered:%d",(int)th2->Integral());
              t->DrawText(0.005,th2->GetMaximum()*0.94,txt);
          }  
           if ( th2 ) {
              hcopy = (TH1F*)th2->DrawCopy("SAME"); 
              hcopy->SetFillColor(50);
              t->SetTextColor(50);
              sprintf(txt,"Air scattered:%d",(int)th2->Integral());
              t->DrawText(0.005,th2->GetMaximum()*0.94,txt);
          }  
           if ( th3 ) {
              hcopy = (TH1F*)th3->DrawCopy("SAME"); 
              hcopy->SetFillColor(kBlue);
              t->SetTextColor(kBlue);
              sprintf(txt,"Aerosol scattered:%d",(int)th3->Integral());
              t->DrawText(0.005,th3->GetMaximum()*0.94,txt);
          }  
        }  
    }
    
    // fluo and ckov air scattered
    else if ( option == "air_scat_pupil" ) {
	th = GetTofHisto("Single_Nph_fluo_airscat_P_t"); // airscat fluo
        th2 = GetTofHisto("Single_Nph_CB_P_t");          // airscat ckov
	TimeMin = (th->GetBinLowEdge(0)+1)*microsecond;
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy("");
	   hcopy->SetFillColor(0); 
	   hcopy->SetLineColor(2); 
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(2);
           sprintf(txt,"fluo air scat:%d",(int)th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
           if ( th2 ) {
	      Double_t Ymax = hcopy->GetMaximum();
	      if(th2->GetMaximum() > Ymax) hcopy->SetMaximum(1.1*th2->GetMaximum());
              hcopy = (TH1F*)th2->DrawCopy("SAME");
	      hcopy->SetFillColor(0); 
	      hcopy->SetLineColor(4); 
              t->SetTextColor(4);
              sprintf(txt,"Cerenkov air scat:%d",(int)th2->Integral());
              t->DrawText(0.005,th2->GetMaximum()*0.94,txt);
          }  
        }  
    }
    
    // fluo and ckov aerosol scattered
    else if ( option == "aero_scat_pupil" ) {
	th = GetTofHisto("Single_Nph_fluo_AeroScat_P_t");           // aeroscat fluo
        th2 = GetTofHisto("Single_Nph_ckov_AeroScat_P_t");          // aeroscat ckov
	TimeMin = (th->GetBinLowEdge(0)+1)*microsecond;
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
	   hcopy->SetFillColor(0); 
	   hcopy->SetLineColor(2); 
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(2);
           sprintf(txt,"fluo aerosol scat:%d",(int)th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
           if ( th2 ) {
	      Double_t Ymax = hcopy->GetMaximum();
	      if(th2->GetMaximum() > Ymax) hcopy->SetMaximum(1.1*th2->GetMaximum());
              hcopy = (TH1F*)th2->DrawCopy("SAME"); 
	      hcopy->SetFillColor(0); 
	      hcopy->SetLineColor(4); 
              t->SetTextColor(4);
              sprintf(txt,"Cerenkov aerosol scat:%d",(int)th2->Integral());
              t->DrawText(0.005,th2->GetMaximum()*0.94,txt);
          }  
        }  
    }
    
    // fluo and ckov cloud scattered
    else if ( option == "cloud_scat_pupil" ) {
	th = GetTofHisto("Single_Nph_fluo_CloudBackscat_P_t");           // cloud scat fluo
        th2 = GetTofHisto("Single_Nph_ckov_CloudBackscat_P_t");          // cloud scat ckov
	TimeMin = (th->GetBinLowEdge(0)+1)*microsecond;
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
	   hcopy->SetFillColor(0); 
	   hcopy->SetLineColor(2); 
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(2);
           sprintf(txt,"fluo cloud scat:%d",(int)th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
          if ( th2 ) {
	      Double_t Ymax = hcopy->GetMaximum();
	      if(th2->GetMaximum() > Ymax) hcopy->SetMaximum(1.1*th2->GetMaximum());
              hcopy = (TH1F*)th2->DrawCopy("SAME"); 
	      hcopy->SetFillColor(0); 
	      hcopy->SetLineColor(4); 
              t->SetTextColor(4);
              sprintf(txt,"Cerenkov cloud scat:%d",(int)th2->Integral());
              t->DrawText(0.005,th2->GetMaximum()*0.94,txt);
          }  
        }  
    }
    
    // fluo and ckov total scattered
    else if ( option == "total_scat_pupil" ) {
	th = GetTofHisto("Single_Nph_fluo_Total_Backscat_P_t");           // total scat fluo
        th2 = GetTofHisto("Single_Nph_ckov_Total_Backscat_P_t");          // total scat ckov
	TimeMin = (th->GetBinLowEdge(0)+1)*microsecond;
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
	   hcopy->SetFillColor(0); 
	   hcopy->SetLineColor(2); 
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(2);
           sprintf(txt,"fluo total scat:%d",(int)th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
           if ( th2 ) {
	      Double_t Ymax = hcopy->GetMaximum();
	      if(th2->GetMaximum() > Ymax) hcopy->SetMaximum(1.1*th2->GetMaximum());
              hcopy = (TH1F*)th2->DrawCopy("SAME"); 
	      hcopy->SetFillColor(0); 
	      hcopy->SetLineColor(4); 
              t->SetTextColor(4);
              sprintf(txt,"Cerenkov total scat:%d",(int)th2->Integral());
              t->DrawText(0.005,th2->GetMaximum()*0.94,txt);
          }  
        }  
    }
    
    // fluo and ckov reflected
    else if ( option == "reflected_pupil" ) {
	th = GetTofHisto("Single_Nph_fluo_refl_P_t");   // reflected fluo
        th1 = GetTofHisto("Single_Nph_CR_P_t");         // reflected ckov
	TimeMin = (th->GetBinLowEdge(0)+1)*microsecond;
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
	   hcopy->SetFillColor(0); 
	   hcopy->SetLineColor(2); 
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           t->SetTextColor(2);
           sprintf(txt,"fluo reflected:%d",(int)th->Integral());
           t->DrawText(0.005,th->GetMaximum()*0.94,txt);
           if ( th1 ) {
	      Double_t Ymax = hcopy->GetMaximum();
	      if(th1->GetMaximum() > Ymax) hcopy->SetMaximum(1.1*th1->GetMaximum());
              hcopy = (TH1F*)th1->DrawCopy("SAME"); 
	      hcopy->SetFillColor(0); 
	      hcopy->SetLineColor(4); 
              t->SetTextColor(4);
              sprintf(txt,"Cerenkov reflected:%d",(int)th1->Integral());
              t->DrawText(0.005,th1->GetMaximum()*0.94,txt);
          }  
        }  
    }
    
    // All the Photons on Pupil
    else if ( option == "Nphs_t" ) {
        th = GetTofHisto("Single_Nph_P_t");
        th1 = GetTofHisto("Single_Nph_F_P_t"); 
        th2 = GetTofHisto("Single_Nph_CR_P_t"); 
        th3 = GetTofHisto("Single_Nph_CB_P_t"); 
        th4 = GetTofHisto("Single_Nph_ckov_CloudBackscat_P_t"); 
        th5 = GetTofHisto("Single_Nph_ckov_AeroScat_P_t"); 
        th6 = GetTofHisto("Single_Nph_direct_ckov_P_t"); 
        th7 = GetTofHisto("Single_Nph_fluo_airscat_P_t"); 
        th8 = GetTofHisto("Single_Nph_fluo_refl_P_t"); 
        th9 = GetTofHisto("Single_Nph_fluo_CloudBackscat_P_t"); 
        th10 = GetTofHisto("Single_Nph_fluo_AeroScat_P_t"); 

	TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
	TString name("name");
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
           //hcopy->SetFillColor(kBlack);
           hcopy->SetLineColor(kBlack);
           hcopy->SetFillColor(0);
           char txt[200];
           sprintf(txt,"TOTAL: %d",(int)th->Integral());
	   name = TString(txt);
	   name += TString("  ");
           leg->AddEntry(hcopy,name.Data(),"lf");
           if ( th1 ) {
              hcopy = (TH1F*)th1->DrawCopy("SAME"); 
              //hcopy->SetFillColor(2);
              hcopy->SetLineColor(2);
              hcopy->SetFillColor(0);
              sprintf(txt,"Direct Fluo: %d",(int)th1->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th1->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           } 
           if ( th4 ) {
              hcopy = (TH1F*)th4->DrawCopy("SAME"); 
              //hcopy->SetFillColor(25);
              hcopy->SetFillColor(21);
              hcopy->SetLineColor(21);
              sprintf(txt,"Cloud scat Ckov: %d",(int)th4->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th4->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }       
           if ( th9 ) {
              hcopy = (TH1F*)th9->DrawCopy("SAME"); 
              //hcopy->SetFillColor(29);
              hcopy->SetFillColor(28);
              hcopy->SetLineColor(28);
              sprintf(txt,"Cloud scat Fluo: %d",(int)th9->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th9->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }       
           if ( th2 ) {
              hcopy = (TH1F*)th2->DrawCopy("SAME"); 
              //hcopy->SetFillColor(6);
              hcopy->SetFillColor(6);
              hcopy->SetLineColor(6);
              sprintf(txt,"Reflected Ckov: %d",(int)th2->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th2->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }    
           if ( th8 ) {
              hcopy = (TH1F*)th8->DrawCopy("SAME"); 
              //hcopy->SetFillColor(48);
              hcopy->SetFillColor(5);
              hcopy->SetLineColor(5);
              sprintf(txt,"Reflected Fluo: %d",(int)th8->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th8->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }    
           if ( th7 ) {
              hcopy = (TH1F*)th7->DrawCopy("SAME"); 
              hcopy->SetFillColor(0);
              //hcopy->SetFillColor(8);
              hcopy->SetLineColor(3);
              hcopy->SetLineWidth(2);
              sprintf(txt,"Air scat Fluo: %d",(int)th7->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th7->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           } 
           if ( th3 ) {
              hcopy = (TH1F*)th3->DrawCopy("SAME"); 
              hcopy->SetFillColor(7);
              hcopy->SetFillStyle(3010);
              hcopy->SetLineColor(7);
              sprintf(txt,"Air scat Ckov: %d",(int)th3->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th3->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }       
           if ( th6 ) {
              hcopy = (TH1F*)th6->DrawCopy("SAME"); 
              hcopy->SetFillColor(0);
              //hcopy->SetFillColor(9);
              hcopy->SetLineColor(9);
              sprintf(txt,"Direct Ckov: %d",(int)th6->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th6->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }
           if ( th5 ) {
              hcopy = (TH1F*)th5->DrawCopy("SAME"); 
              //hcopy->SetFillColor(35);
              hcopy->SetFillColor(0);
              hcopy->SetLineColor(35);
              sprintf(txt,"Aerosol scat Ckov: %d",(int)th5->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th5->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }
           if ( th10 ) {
              hcopy = (TH1F*)th10->DrawCopy("SAME"); 
              //hcopy->SetFillColor(39);
              hcopy->SetFillColor(0);
              hcopy->SetLineColor(39);
              sprintf(txt,"Aerosol scat Fluo: %d",(int)th10->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th10->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }
        }
	leg->Draw();  
    }
    
    // all photons as function of their nb of interaction in atmosphere
    else if ( option == "photon_pupil_nbinter" ) {
        TString prefix("pupil_nbinter");
        TString prefix2("order ");
        TString histonb(prefix);
        th2 = GetTofHisto("Single_Nph_P_t");
	TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
        if( th2 ) {
           hcopy = (TH1F*)th2->DrawCopy(""); 
           hcopy->SetFillColor(1);
           //hcopy->SetLineColor(1);
	   hcopy->SetTitle("Time on pupil : scattering order contributions");
	   histonb = TString("total -> ");
	   histonb += hcopy->Integral();
	   leg->AddEntry(hcopy,histonb.Data(),"f");
        } 
        for (Int_t i=0; i<=fMaxScatOrder; i++) { 
	    histonb = prefix;
	    histonb += i;
	    th = (TH1F*)fHistos->FindObject(histonb.Data());
            if( th ) {
	        hcopy = (TH1F*)th->DrawCopy("same"); 
		histonb = prefix2;
		histonb += i;
		histonb += TString(" -> ");
		histonb += hcopy->Integral();
		leg->AddEntry(hcopy,histonb.Data(),"f");
		hcopy->SetFillColor(i+2); 
	        //hcopy->SetLineColor(i+2); 
	    }
	    else Printf("<BuildTofHistos> Wrong histo identification  %s",histonb.Data());
	}
	leg->Draw();
    }




// ****************** ENTERING ANGLE IN DETECTOR ********************* //


    // All the Photons
    else if ( option == "angul_distrib_all" ) {
        th = GetAngularHisto("Single_Nph_P_t_angul");
        th1 = GetAngularHisto("Single_Nph_F_P_t_angul"); 
        th2 = GetAngularHisto("Single_Nph_CR_P_t_angul"); 
        th3 = GetAngularHisto("Single_Nph_CB_P_t_angul"); 
        th4 = GetAngularHisto("Single_Nph_ckov_CloudBackscat_P_t_angul"); 
        th5 = GetAngularHisto("Single_Nph_ckov_AeroScat_P_t_angul"); 
        th6 = GetAngularHisto("Single_Nph_direct_ckov_P_t_angul"); 
        th7 = GetAngularHisto("Single_Nph_fluo_airscat_P_t_angul"); 
        th8 = GetAngularHisto("Single_Nph_fluo_refl_P_t_angul"); 
        th9 = GetAngularHisto("Single_Nph_fluo_CloudBackscat_P_t_angul"); 
        th10 = GetAngularHisto("Single_Nph_fluo_AeroScat_P_t_angul"); 

	TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
	TString name("name");
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
           //hcopy->SetFillColor(kBlack);
           hcopy->SetLineColor(kBlack);
           hcopy->SetFillColor(0);
           char txt[200];
           sprintf(txt,"TOTAL: %d",(int)th->Integral());
	   name = TString(txt);
	   name += TString("  ");
           leg->AddEntry(hcopy,name.Data(),"lf");
           if ( th1 ) {
              hcopy = (TH1F*)th1->DrawCopy("SAME"); 
              //hcopy->SetFillColor(2);
              hcopy->SetLineColor(2);
              hcopy->SetFillColor(0);
              sprintf(txt,"Direct Fluo: %d",(int)th1->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th1->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           } 
           if ( th4 ) {
              hcopy = (TH1F*)th4->DrawCopy("SAME"); 
              //hcopy->SetFillColor(25);
              hcopy->SetFillColor(21);
              hcopy->SetLineColor(21);
              sprintf(txt,"Cloud scat Ckov: %d",(int)th4->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th4->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }       
           if ( th9 ) {
              hcopy = (TH1F*)th9->DrawCopy("SAME"); 
              //hcopy->SetFillColor(29);
              hcopy->SetFillColor(28);
              hcopy->SetLineColor(28);
              sprintf(txt,"Cloud scat Fluo: %d",(int)th9->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th9->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }       
           if ( th2 ) {
              hcopy = (TH1F*)th2->DrawCopy("SAME"); 
              //hcopy->SetFillColor(6);
              hcopy->SetFillColor(6);
              hcopy->SetLineColor(6);
              sprintf(txt,"Reflected Ckov: %d",(int)th2->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th2->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }    
           if ( th8 ) {
              hcopy = (TH1F*)th8->DrawCopy("SAME"); 
              //hcopy->SetFillColor(48);
              hcopy->SetFillColor(5);
              hcopy->SetLineColor(5);
              sprintf(txt,"Reflected Fluo: %d",(int)th8->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th8->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }    
           if ( th7 ) {
              hcopy = (TH1F*)th7->DrawCopy("SAME"); 
              hcopy->SetFillColor(0);
              //hcopy->SetFillColor(8);
              hcopy->SetLineColor(3);
              hcopy->SetLineWidth(2);
              sprintf(txt,"Air scat Fluo: %d",(int)th7->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th7->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           } 
           if ( th3 ) {
              hcopy = (TH1F*)th3->DrawCopy("SAME"); 
              hcopy->SetFillColor(7);
              hcopy->SetFillStyle(3010);
              hcopy->SetLineColor(7);
              sprintf(txt,"Air scat Ckov: %d",(int)th3->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th3->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }       
           if ( th6 ) {
              hcopy = (TH1F*)th6->DrawCopy("SAME"); 
              hcopy->SetFillColor(0);
              //hcopy->SetFillColor(9);
              hcopy->SetLineColor(9);
              sprintf(txt,"Direct Ckov: %d",(int)th6->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th6->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }
           if ( th5 ) {
              hcopy = (TH1F*)th5->DrawCopy("SAME"); 
              //hcopy->SetFillColor(35);
              hcopy->SetFillColor(0);
              hcopy->SetLineColor(35);
              sprintf(txt,"Aerosol scat Ckov: %d",(int)th5->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th5->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }
           if ( th10 ) {
              hcopy = (TH1F*)th10->DrawCopy("SAME"); 
              //hcopy->SetFillColor(39);
              hcopy->SetFillColor(0);
              hcopy->SetLineColor(39);
              sprintf(txt,"Aerosol scat Fluo: %d",(int)th10->Integral());
	      name = TString(txt);
	      name += TString("  ");
              if(th10->Integral()) leg->AddEntry(hcopy,name.Data(),"lf");
           }
        }
	leg->Draw();  
    }
    
    // all photons as function of their nb of interaction in atmosphere
    else if ( option == "photon_angul_nbinter" ) {
        TString prefix("angular_nbinter");
        TString prefix2("order ");
        TString histonb(prefix);
        th2 = GetAngularHisto("Single_Nph_P_t_angul");
	TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
        if( th2 ) {
           hcopy = (TH1F*)th2->DrawCopy(""); 
           hcopy->SetFillColor(1);
           //hcopy->SetLineColor(1);
	   hcopy->SetTitle("Angle on pupil : scattering order contributions");
	   histonb = TString("total -> ");
	   histonb += hcopy->Integral();
	   leg->AddEntry(hcopy,histonb.Data(),"f");
        } 
        for (Int_t i=0; i<=fMaxScatOrder; i++) { 
	    histonb = prefix;
	    histonb += i;
	    th = (TH1F*)fHistos->FindObject(histonb.Data());
            if( th ) {
	        hcopy = (TH1F*)th->DrawCopy("same"); 
		histonb = prefix2;
		histonb += i;
		histonb += TString(" -> ");
		histonb += hcopy->Integral();
		leg->AddEntry(hcopy,histonb.Data(),"f");
		hcopy->SetFillColor(i+2); 
	        //hcopy->SetLineColor(i+2); 
	    }
	    else Printf("<BuildAngularPupilHistos.Draw> Wrong histo identification  %s",histonb.Data());
	}
	leg->Draw();
    }
    
    else if ( option == "2D_direct" ) {
        th2D = Get2DAngularHisto("angul2D_pupil");
        if ( th2D ) hcopy2D = (TH2F*)th2D->DrawCopy("colz");   
    }
    
    
    
    
// ******** ALTITUDE ************** //
// **** for bunches and singles *** //

    // FLUO Bunches are plotted VS altitude
    else if ( option == "alt_Nph_F" ) {
        th = GetAltHisto("alt_Bunch_Nph_F");        
        if ( th ) {
           //TF1* mygaus = new TF1("mygaus","gaus",0,10);
           //th->Fit("mygaus","R");
           Double_t par[3];
           //mygaus->GetParameters(par);
	   //Double_t binwidth = th->GetBinWidth(th->FindBin(par[1]));
           Double_t binwidth =th->GetBinWidth(th->GetMaximumBin());
           par[0] = th->GetBinContent(th->GetMaximumBin());
	   par[1] = th->GetBinCenter(th->GetMaximumBin());
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(29);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           sprintf(txt,"Total:%e",th->Integral());
           t->DrawText(5,th->GetMaximum()*0.94,txt);
           sprintf(txt,"Alt max : %e",par[1]);
           t->DrawText(5,th->GetMaximum()*0.84,txt);
           sprintf(txt,"Normalized Nph max : %e km-1",par[0]/binwidth);
           t->DrawText(5,th->GetMaximum()*0.74,txt);
        }  
    }
    
    // CKOV Bunches are plotted VS altitude
    else if ( option == "alt_Nph_C" ) {
        th = GetTimeHisto("alt_Bunch_Nph_C");        
        if ( th ) {
           //TF1* mygaus = new TF1("mygaus","gaus",0,10);
	   //th->Fit("mygaus","R");
           Double_t par[3];
           //mygaus->GetParameters(par);
	   //Double_t binwidth = th->GetBinWidth(th->FindBin(par[1]));
           Double_t binwidth =th->GetBinWidth(th->GetMaximumBin());
           par[0] = th->GetBinContent(th->GetMaximumBin());
	   par[1] = th->GetBinCenter(th->GetMaximumBin());
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(39);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           sprintf(txt,"Total:%e",th->Integral());
           t->DrawText(5,th->GetMaximum()*0.94,txt);
           sprintf(txt,"Alt max : %e",par[1]);
           t->DrawText(5,th->GetMaximum()*0.84,txt);
           sprintf(txt,"Normalized Nph max : %e km-1",par[0]/binwidth);
           t->DrawText(5,th->GetMaximum()*0.74,txt);
        } 
    }
    
    // FLUO Yield
    else if ( option == "alt_Yield_F" ) {
        th = GetAltHisto("alt_Bunch_Yield_F");
        if ( th ) {
	   th->SetMinimum(0.9*th->GetMinimum());
           hcopy = (TH1F*)th->DrawCopy();
           hcopy->SetFillColor(29);
	}
    }            
    
    // CKOV Yield
    else if ( option == "alt_Yield_C" ) {
        th = GetAltHisto("alt_Bunch_Yield_C");  
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
           hcopy->SetFillColor(39);
	}
    } 
    
    // SinglePhoton vs altitude
    else if ( option == "alt_direct" ) {
        th = GetSingleAltHisto("alt_directfluo");
        th2 = GetSingleAltHisto("alt_directckov");
	Double_t Ymax = th->GetMaximum();
	if(th2->GetMaximum() > Ymax) th->SetMaximum(th2->GetMaximum());
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
	   hcopy->SetLineColor(2); 
	}
        if ( th2 ) {
           hcopy = (TH1F*)th2->DrawCopy("same");
	   hcopy->SetLineColor(4); 
	}
    }     
    else if ( option == "alt_scat" ) {
        th = GetSingleAltHisto("alt_scatfluo");
        th2 = GetSingleAltHisto("alt_scatckov");
	Double_t Ymax = th->GetMaximum();
	if(th2->GetMaximum() > Ymax) th->SetMaximum(th2->GetMaximum());
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
	   hcopy->SetLineColor(2); 
	}
        if ( th2) {
           hcopy = (TH1F*)th2->DrawCopy("same");
	   hcopy->SetLineColor(4); 
	}
    }
    else if ( option == "alt_scat_type" ) {
	TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
        th = GetSingleAltHisto("alt_airscat");
        th2 = GetSingleAltHisto("alt_cloudscat");
        th3 = GetSingleAltHisto("alt_aeroscat");
	Double_t Ymax = th->GetMaximum();
	if(th2->GetMaximum() > Ymax) th->SetMaximum(th2->GetMaximum());
	if(th3->GetMaximum() > Ymax) th->SetMaximum(th3->GetMaximum());
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
	   hcopy->SetLineColor(8); 
	   leg->AddEntry(hcopy,"air scattered","l");
	}
        if ( th2) {
           hcopy = (TH1F*)th2->DrawCopy("same");
	   hcopy->SetLineColor(2); 
	   leg->AddEntry(hcopy,"clouds scattered","l");
	}
        if ( th3) {
           hcopy = (TH1F*)th3->DrawCopy("same");
	   hcopy->SetLineColor(4); 
	   leg->AddEntry(hcopy,"aerosols scattered","l");
	}
	leg->Draw();
    }
    
    // SinglePhoton transmission OVERVIEW vs altitude    
    // fluo
    else if ( option == "alt_fluo_trans" ) {
        th = GetSingleAltHisto("alt_fluo_Ntrans");
        th2 = GetSingleAltHisto("alt_fluo_trans");
	th3 = GetSingleAltHisto("alt_trans_fluo");
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
           hcopy->SetFillColor(29);
	}
        if ( th2) {
           hcopy = (TH1F*)th2->DrawCopy("same");
           hcopy->SetFillColor(2);
	   gPad->Update();
	}
	if(th3) {
	    Double_t rightmax = 1.1*1;
	    Double_t scale = gPad->GetUymax()/rightmax;
	    th3->Rebin(2);
	    th3->Scale(0.5*scale);  // 0.5 because of rebinning
	    hcopy = (TH1F*)th3->DrawCopy("same");
	    hcopy->SetLineColor(4); 
	    TGaxis* axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),gPad->GetUxmax(),gPad->GetUymax(),0,rightmax,510,"+L");
	    axis->SetLineColor(4);
	    axis->SetLabelColor(4);
	    axis->SetTextColor(4);
	    axis->Draw();
	}
    }
    // ckov
    else if ( option == "alt_ckov_trans" ) {
        th = GetSingleAltHisto("alt_ckov_Ntrans");
        th2 = GetSingleAltHisto("alt_ckov_trans");
	th3 = GetSingleAltHisto("alt_trans_ckov");
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
           hcopy->SetFillColor(29);
	}
        if ( th2) {
           hcopy = (TH1F*)th2->DrawCopy("same");
           hcopy->SetFillColor(2);
	   gPad->Update();
	}
	if(th3) {
	    Double_t rightmax = 1.1*1;
	    Double_t scale = gPad->GetUymax()/rightmax;
	    th3->Rebin(2);
	    th3->Scale(0.5*scale);  // 0.5 because of rebinning
	    hcopy = (TH1F*)th3->DrawCopy("same");
	    hcopy->SetLineColor(4); 
	    TGaxis* axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),gPad->GetUxmax(),gPad->GetUymax(),0,rightmax,510,"+L");
	    axis->SetLineColor(4);
	    axis->SetLabelColor(4);
	    axis->SetTextColor(4);
	    axis->Draw();
	}
    }
    
    // transmission DETAILS v.s. altitude
    else if ( option == "alt_tottrans_details" ) {
        prof1 = GetAlt2DHisto("tottrans_c");
        prof2 = GetAlt2DHisto("tottrans_f");
        if ( prof1 ) {
           hcopy = (TH2F*)prof1->DrawCopy();
           hcopy->SetMarkerColor(4);
	}
        if ( prof2) {
           hcopy = (TH2F*)prof2->DrawCopy("same");
           hcopy->SetMarkerColor(2);
	   gPad->Update();
	}
    }
    else if ( option == "alt_rayltrans_details" ) {
        prof1 = GetAlt2DHisto("rayltrans_c");
        prof2 = GetAlt2DHisto("rayltrans_f");
        if ( prof1 ) {
           hcopy = (TH2F*)prof1->DrawCopy();
           hcopy->SetMarkerColor(4);
	}
        if ( prof2) {
           hcopy = (TH2F*)prof2->DrawCopy("same");
           hcopy->SetMarkerColor(2);
	   gPad->Update();
	}
    }
    else if ( option == "alt_ozonetrans_details" ) {
        prof1 = GetAlt2DHisto("ozonetrans_c");
        prof2 = GetAlt2DHisto("ozonetrans_f");
        if ( prof1 ) {
           hcopy = (TH2F*)prof1->DrawCopy();
           hcopy->SetMarkerColor(4);
	}
        if ( prof2) {
           hcopy = (TH2F*)prof2->DrawCopy("same");
           hcopy->SetMarkerColor(2);
	   gPad->Update();
	}
    }
    else if ( option == "alt_aerotrans_details" ) {
        prof1 = GetAlt2DHisto("aerotrans_c");
        prof2 = GetAlt2DHisto("aerotrans_f");
        if ( prof1 ) {
           hcopy = (TH2F*)prof1->DrawCopy();
           hcopy->SetMarkerColor(4);
	}
        if ( prof2) {
           hcopy = (TH2F*)prof2->DrawCopy("same");
           hcopy->SetMarkerColor(2);
	   gPad->Update();
	}
    }
    else if ( option == "alt_cloudtrans_details" ) {
        prof1 = GetAlt2DHisto("cloudtrans_c");
        prof2 = GetAlt2DHisto("cloudtrans_f");
        if ( prof1 ) {
           hcopy = (TH2F*)prof1->DrawCopy();
           hcopy->SetMarkerColor(4);
	}
        if ( prof2) {
           hcopy = (TH2F*)prof2->DrawCopy("same");
           hcopy->SetMarkerColor(2);
	   gPad->Update();
	}
    }
    

// ******** LONGITUDINAL PROFILES ************* //
    // FLUO Bunches longitudinal profile in km
    else if ( option == "Fluo_longi_prof_km" ) {
        th = GetLongitudinalHisto("Fluo_Longit_prof_km");        
        if ( th ) {
           Double_t par[3];
	   Double_t binwidth = th->GetBinWidth(th->GetMaximumBin());
           par[0] = th->GetBinContent(th->GetMaximumBin());
	   par[1] = th->GetBinCenter(th->GetMaximumBin());
	   //Double_t maxf = th->GetBinLowEdge(th->GetXaxis()->GetLast()) + binwidth;
           //TF1* mygaus = new TF1("mygaus","gaus",0,maxf);
           //th->Fit("mygaus","R");
           //mygaus->GetParameters(par);
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(29);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           sprintf(txt,"Total:%e",th->Integral());
           t->DrawText(0.2,th->GetMaximum()*0.94,txt);
           sprintf(txt,"Tracklength max : %e",par[1]);
           t->DrawText(0.2,th->GetMaximum()*0.84,txt);
           sprintf(txt,"Normalized Nph max : %e km-1",par[0]/binwidth);
           t->DrawText(0.2,th->GetMaximum()*0.74,txt);
        }  
    }
    
    // FLUO Bunches longitudinal profile in g/cm2
    else if ( option == "Fluo_longi_prof_gram" ) {
        th = GetLongitudinalHisto("Fluo_Longit_prof_gram");        
        if ( th ) {
	   Double_t binwidth = th->GetBinWidth(th->GetMaximumBin());
           /*
	   Double_t maxf = th->GetBinLowEdge(th->GetXaxis()->GetLast()) + binwidth;
	   TF1* mygaus = new TF1("mygaus","gaus",0,maxf);
           th->Fit("mygaus","R");
           Double_t par[3];
           mygaus->GetParameters(par);
           */
	   Double_t xmax = th->GetBinCenter(th->GetMaximumBin());
	   Double_t nphmax = th->GetBinContent(th->GetMaximumBin());
	   hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(29);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           sprintf(txt,"Total:%e",th->Integral());
           t->DrawText(0.2,th->GetMaximum()*0.94,txt);
           sprintf(txt,"Xmax : %e",xmax);
           t->DrawText(0.2,th->GetMaximum()*0.84,txt);
           sprintf(txt,"Normalized Nph max : %e cm2/g",nphmax/binwidth);
           t->DrawText(0.2,th->GetMaximum()*0.74,txt);
        }  
    }
    
    // CKOV Bunches longitudinal profile in km
    else if ( option == "Ckov_longi_prof_km" ) {
        th = GetLongitudinalHisto("Ckov_Longit_prof_km");        
        if ( th ) {
           Double_t par[3];
	   Double_t binwidth = th->GetBinWidth(th->GetMaximumBin());
           par[0] = th->GetBinContent(th->GetMaximumBin());
	   par[1] = th->GetBinCenter(th->GetMaximumBin());
	   //Double_t maxf = th->GetBinLowEdge(th->GetXaxis()->GetLast()) + binwidth;
           //TF1* mygaus = new TF1("mygaus","gaus",0,maxf);
           //th->Fit("mygaus","R");
           //mygaus->GetParameters(par);
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(29);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           sprintf(txt,"Total:%e",th->Integral());
           t->DrawText(0.2,th->GetMaximum()*0.94,txt);
           sprintf(txt,"Tracklength max : %e",par[1]);
           t->DrawText(0.2,th->GetMaximum()*0.84,txt);
           sprintf(txt,"Normalized Nph max : %e km-1",par[0]/binwidth);
           t->DrawText(0.2,th->GetMaximum()*0.74,txt);
        }  
    }
    
    // CKOV Bunches longitudinal profile in g/cm2
    else if ( option == "Ckov_longi_prof_gram" ) {
        th = GetLongitudinalHisto("Ckov_Longit_prof_gram");        
        if ( th ) {
	   Double_t binwidth = th->GetBinWidth(th->GetMaximumBin());
           /*
	   Double_t maxf = th->GetBinLowEdge(th->GetXaxis()->GetLast()) + binwidth;
	   TF1* mygaus = new TF1("mygaus","gaus",0,maxf);
           th->Fit("mygaus","R");
           Double_t par[3];
           mygaus->GetParameters(par);
	   */
	   Double_t xmax = th->GetBinCenter(th->GetMaximumBin());
	   Double_t nphmax = th->GetBinContent(th->GetMaximumBin());
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(29);
           TText *t = new TText();
           char txt[200];
           t->SetTextFont(62);  
           t->SetTextSize(0.06);
           sprintf(txt,"Total:%e",th->Integral());
           t->DrawText(0.2,th->GetMaximum()*0.94,txt);
           sprintf(txt,"Xmax : %e",xmax);
           t->DrawText(0.2,th->GetMaximum()*0.84,txt);
           sprintf(txt,"Normalized Nph max : %e cm2/g",nphmax/binwidth);
           t->DrawText(0.2,th->GetMaximum()*0.74,txt);
        }  
    }

    // SinglePhotons longitudinal distribution
    else if ( option == "longi_direct" ) {
        th = GetSingleLongitudinalHisto("fluo_direct_longi");
        th2 = GetSingleLongitudinalHisto("ckov_direct_longi");
	Double_t Ymax = th->GetMaximum();
	if(th2->GetMaximum() > Ymax) th->SetMaximum(th2->GetMaximum());
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
	   hcopy->SetLineColor(2); 
           hcopy->SetFillColor(0);
	}
        if ( th2 ) {
           hcopy = (TH1F*)th2->DrawCopy("same");
	   hcopy->SetLineColor(4); 
           hcopy->SetFillColor(0);
	}
    }     
    else if ( option == "longi_scat" ) {
        th = GetSingleLongitudinalHisto("fluo_scat_longi");
        th2 = GetSingleLongitudinalHisto("ckov_scat_longi");
	Double_t Ymax = th->GetMaximum();
	if(th2->GetMaximum() > Ymax) th->SetMaximum(th2->GetMaximum());
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
	   hcopy->SetLineColor(2); 
           hcopy->SetFillColor(0);
	}
        if ( th2 ) {
           hcopy = (TH1F*)th2->DrawCopy("same");
	   hcopy->SetLineColor(4); 
           hcopy->SetFillColor(0);
	}
    }     

    // SinglePhotons lateral distribution
    else if ( option == "lateral_direct" ) {
        th = GetSingleLongitudinalHisto("fluo_direct_lat");
        th2 = GetSingleLongitudinalHisto("ckov_direct_lat");
	Double_t Ymax = th->GetMaximum();
	if(th2->GetMaximum() > Ymax) th->SetMaximum(th2->GetMaximum());
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
	   hcopy->SetLineColor(2); 
           hcopy->SetFillColor(0);
	}
        if ( th2 ) {
           hcopy = (TH1F*)th2->DrawCopy("same");
	   hcopy->SetLineColor(4); 
	}
    }     
    else if ( option == "lateral_scat" ) {
        th = GetSingleLongitudinalHisto("fluo_scat_lat");
        th2 = GetSingleLongitudinalHisto("ckov_scat_lat");
	Double_t Ymax = th->GetMaximum();
	if(th2->GetMaximum() > Ymax) th->SetMaximum(th2->GetMaximum());
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
	   hcopy->SetLineColor(2); 
           hcopy->SetFillColor(0);
	}
        if ( th2 ) {
           hcopy = (TH1F*)th2->DrawCopy("same");
	   hcopy->SetLineColor(4); 
           hcopy->SetFillColor(0);
	}
    }     

    // SinglePhotons 2D-distribution
    else if ( option == "2D_direct" ) {
        th = GetSingleLongitudinalHisto("fluo_direct_2D");
        th2 = GetSingleLongitudinalHisto("ckov_direct_2D");
	Double_t Ymax = th->GetMaximum();
	if(th2->GetMaximum() > Ymax) th->SetMaximum(th2->GetMaximum());
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();
           hcopy->SetMarkerColor(2);
           hcopy->SetMarkerSize(0.1);
	}
        if ( th2 ) {
           hcopy = (TH1F*)th2->DrawCopy("same");
           hcopy->SetMarkerColor(4);
           hcopy->SetMarkerSize(0.1);
	}
    }     
    else if ( option == "2D_scat" ) {
        th = GetSingleLongitudinalHisto("fluo_scat_2D");
        th2 = GetSingleLongitudinalHisto("ckov_scat_2D");
	Double_t Ymax = th->GetMaximum();
	if(th2->GetMaximum() > Ymax) th->SetMaximum(th2->GetMaximum());
        if ( th ) {
           th->SetMarkerColor(2);
           th->SetMarkerSize(0.1);
           th->DrawCopy();
	}
        if ( th2 ) {
           th2->SetMarkerColor(4);
           th2->SetMarkerSize(0.1);
           th2->DrawCopy("same");
	}
    } 

// ******** WAVELENGTH ************* //
    
    // FLUO spectrum at creation
    else if ( option == "Sp_F" ) {
        th = GetWlHisto("Bunch_Wl_F");        
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(29);
        }  
    }

    // CKOV spectrum at creation
    else if ( option == "Sp_C" ) {
        th = GetWlHisto("Bunch_Wl_C");        
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy(""); 
           hcopy->SetFillColor(39);
        }  
    }

    // ALL fluorescence transmission spectrum
    else if ( option == "Sps_F" ) {
        th = GetWlHisto("Single_Wl_F"); 
        th1 = GetWlHisto("Single_Wl_F_P");
	th2 = GetWlHisto("wl_trans_fluo");
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();   
           hcopy->SetFillColor(29);
           if( th1 ) { 
              hcopy = (TH1F*)th1->DrawCopy("SAME"); 
              hcopy->SetFillColor(2); 
	      gPad->Update();
           }
	   if(th2) {
	       Double_t rightmax = 1.1*1;
	       Double_t scale = gPad->GetUymax()/rightmax;
	       th2->Rebin(2);
	       th2->Scale(0.5*scale);
	       hcopy = (TH1F*)th2->DrawCopy("same");
	       hcopy->SetLineColor(4);
	       hcopy->SetFillColor(0);
	       TGaxis* axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),gPad->GetUxmax(),gPad->GetUymax(),0,rightmax,510,"+L");
	       axis->SetLineColor(4);
	       axis->SetLabelColor(4);
	       axis->SetTextColor(4);
	       axis->Draw();
	   }
        }  
    }
    // ALL cerenkov transmission spectrum
    else if ( option == "Sps_C" ) {
        th = GetWlHisto("wl_ckov_Ntrans"); 
        th1 = GetWlHisto("wl_ckov_trans");       
	th2 = GetWlHisto("wl_trans_ckov");
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();   
           hcopy->SetFillColor(29);
           if( th1 ) { 
              hcopy = (TH1F*)th1->DrawCopy("SAME"); 
              hcopy->SetFillColor(2); 
	      gPad->Update();
           }
	   if(th2) {
	       Double_t rightmax = 1.1*1;
	       Double_t scale = gPad->GetUymax()/rightmax;
	       th2->Rebin(2);
	       th2->Scale(0.5*scale);
	       hcopy = (TH1F*)th2->DrawCopy("same");
	       hcopy->SetLineColor(4);
	       hcopy->SetFillColor(0);
	       TGaxis* axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),gPad->GetUxmax(),gPad->GetUymax(),0,rightmax,510,"+L");
	       axis->SetLineColor(4);
	       axis->SetLabelColor(4);
	       axis->SetTextColor(4);
	       axis->Draw();
	   }
        }  
    }

    // reflected Cerenkov spectrum
    else if ( option == "Sps_CR" ) {
        th = GetWlHisto("Single_Wl_CR"); 
        th1 = GetWlHisto("Single_Wl_CR_P"); 
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();   
           hcopy->SetFillColor(39);
           if( th1 ) { 
              hcopy = (TH1F*)th1->DrawCopy("SAME"); 
              hcopy->SetFillColor(50); 
           }
        }  
    }
    
    // scattered Cerenkov spectrum
    else if ( option == "Sps_CB" ) {
        //th = GetWlHisto("Single_Wl_CB"); 
        th = GetWlHisto("Single_Wl_CB_tot"); 
        th1 = GetWlHisto("Single_Wl_CB_P"); 
        th3 = GetWlHisto("Single_Wl_ckov_AeroScat_P"); 
        th2 = GetWlHisto("Single_Wl_ckov_CloudBackscat_P"); 
        if ( th ) {
           hcopy = (TH1F*)th->DrawCopy();   
           hcopy->SetFillColor(39);
           if( th1 ) { 
              hcopy = (TH1F*)th1->DrawCopy("SAME"); 
              hcopy->SetFillColor(50); 
           }
           if( th3 ) { 
              hcopy = (TH1F*)th3->DrawCopy("SAME"); 
              hcopy->SetFillColor(kBlue); 
           }
           if( th2 ) { 
              hcopy = (TH1F*)th2->DrawCopy("SAME"); 
              hcopy->SetFillColor(5); 
           }
        }  
    }    
    
    // transmission DETAILS v.s. lamda
    else if ( option == "wl_tottrans_details" ) {
        prof1 = GetWl2DHisto("WLtottrans_c");
        prof2 = GetWl2DHisto("WLtottrans_f");
        if ( prof1 ) {
           hcopy = (TH2F*)prof1->DrawCopy();
           hcopy->SetMarkerColor(4);
	}
        if ( prof2) {
           hcopy = (TH2F*)prof2->DrawCopy("same");
           hcopy->SetMarkerColor(2);
	   gPad->Update();
	}
    }
    else if ( option == "wl_rayltrans_details" ) {
        prof1 = GetWl2DHisto("WLrayltrans_c");
        prof2 = GetWl2DHisto("WLrayltrans_f");
        if ( prof1 ) {
           hcopy = (TH2F*)prof1->DrawCopy();
           hcopy->SetMarkerColor(4);
	}
        if ( prof2) {
           hcopy = (TH2F*)prof2->DrawCopy("same");
           hcopy->SetMarkerColor(2);
	   gPad->Update();
	}
    }
    else if ( option == "wl_ozonetrans_details" ) {
        prof1 = GetWl2DHisto("WLozonetrans_c");
        prof2 = GetWl2DHisto("WLozonetrans_f");
        if ( prof1 ) {
           hcopy = (TH2F*)prof1->DrawCopy();
           hcopy->SetMarkerColor(4);
	}
        if ( prof2) {
           hcopy = (TH2F*)prof2->DrawCopy("same");
           hcopy->SetMarkerColor(2);
	   gPad->Update();
	}
    }
    else if ( option == "wl_aerotrans_details" ) {
        prof1 = GetWl2DHisto("WLaerotrans_c");
        prof2 = GetWl2DHisto("WLaerotrans_f");
        if ( prof1 ) {
           hcopy = (TH2F*)prof1->DrawCopy();
           hcopy->SetMarkerColor(4);
	}
        if ( prof2) {
           hcopy = (TH2F*)prof2->DrawCopy("same");
           hcopy->SetMarkerColor(2);
	   gPad->Update();
	}
    }
    else if ( option == "wl_cloudtrans_details" ) {
        prof1 = GetWl2DHisto("WLcloudtrans_c");
        prof2 = GetWl2DHisto("WLcloudtrans_f");
        if ( prof1 ) {
           hcopy = (TH2F*)prof1->DrawCopy();
           hcopy->SetMarkerColor(4);
	}
        if ( prof2) {
           hcopy = (TH2F*)prof2->DrawCopy("same");
           hcopy->SetMarkerColor(2);
	   gPad->Update();
	}
    }


// ************************* FoV *****************************//
    
    // 2D-projection for FoV visualization - no MScattCUT
    else if ( option == "total_Fov_nocut" ) {
        th2D = GetFovHisto("2D-FoV_nocut"); 
        if ( th2D ) {
            hcopy2D = (TH2F*)th2D->DrawCopy("colz");   
	}
    }
    
    // idem but with the MScattCUT - if doesn't exist -> same plot as above
    else if ( option == "total_Fov" ) {
        th2D = GetFovHisto("2D-FoV"); 
        if ( th2D ) {
            hcopy2D = (TH2F*)th2D->DrawCopy("colz");
	}
    }


// ************************* Multiple Scattering details *****************************//
    else if ( option == "nbinteractions" ) {
        th = (TH1F*)fHistos->FindObject("nbinter");
        if(th) hcopy = (TH1F*)th->DrawCopy();
	hcopy->GetYaxis()->SetRangeUser(0,th->GetMaximum()*1.1);
        th2 = (TH1F*)fHistos->FindObject("nbinter_trans");
	th3 = (TH1F*)fHistos->FindObject("nbinter_trans_ratio");
        if(th2) {
	    hcopy = (TH1F*)th2->DrawCopy("same");
	    hcopy->SetLineColor(2);
	}
	if(th3) {
	    Double_t rightmax = 1.1*1;
	    Double_t scale = th->GetMaximum()*1.1/rightmax;
	    th3->Scale(scale);
	    th3->Draw("same");
	    th3->SetLineColor(4); 
	    th3->SetLineStyle(2);
	    TAxis* x_axis = th->GetXaxis();
	    TAxis* y_axis = th->GetYaxis();
	    TGaxis* axis = new TGaxis(x_axis->GetXmax(),y_axis->GetXmin(),x_axis->GetXmax(),th->GetMaximum()*1.1,0,rightmax,510,"+L");
	    axis->SetLineColor(4);
	    axis->SetLabelColor(4);
	    axis->SetTextColor(4);
	    axis->Draw();
	}
    }
    
    // omegadiff (per interaction order - direct not included)
    else if ( option == "omegadiff" ) {
        TString prefix("omegadiff");
        TString name(prefix);
	TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
	for (Int_t i=1; i<=fMaxScatOrder; i++) { 
            name = prefix;
	    name += i;
	    th = (TH1F*) fHistos->FindObject(name.Data());
	    if(!th) {Printf("<Draw> Wrong histo identification  %s",name.Data()); return;}
	    if(i == 1) {
	        th2 = (TH1F*)th->DrawCopy();  // first histo kept in memory
	        name = TString("order ");
	        name += i;
	        name += TString(": mean, rms = ");
	        name += Int_t(th2->GetMean()*100);
	        name += TString(", ");
	        name += Int_t(th2->GetRMS()*100);
	        name += TString(" %");
	        leg->AddEntry(th2,name.Data(),"l");
	    }
	    else {
	        Double_t Ymax = 1.1*th2->GetMaximum();
	        if(th->GetMaximum() > Ymax) th2->SetMaximum(1.1*th->GetMaximum());
	        hcopy = (TH1F*)th->DrawCopy("same");
	        hcopy->SetLineColor(i+2);
	        name = TString("order ");
	        name += i;
	        name += TString(": mean, rms = ");
	        name += Int_t(hcopy->GetMean()*100);
	        name += TString(", ");
	        name += Int_t(hcopy->GetRMS()*100);
	        name += TString(" %");
	        leg->AddEntry(hcopy,name.Data(),"l");
	    }
	}
	leg->Draw();
    }
    
    // Lateral distribution of all photons as function of their nb of interaction in atmosphere
    else if ( option == "lateral_distrib_nbinter" ) {
        TString prefix("lateral_nbinter");
        TString name(prefix);
	TLegend* leg = new TLegend(0.6,0.6,0.9,0.9);
	for (Int_t i=1; i<=fMaxScatOrder; i++) { 
            name = prefix;
	    name += i;
	    th = (TH1F*) fHistos->FindObject(name.Data());
	    if(!th) {Printf("<Draw> Wrong histo identification  %s",name.Data()); return;}
	    if(i == 1) {
	        th2 = (TH1F*)th->DrawCopy();  // first histo kept in memory
	        th2->SetFillColor(0);
		name = TString("order ");
	        name += i;
	        leg->AddEntry(th2,name.Data(),"l");
	    }
	    else {
	        Double_t Ymax = 1.1*th2->GetMaximum();
	        if(th->GetMaximum() > Ymax) th2->SetMaximum(1.1*th->GetMaximum());
	        hcopy = (TH1F*)th->DrawCopy("same");
	        hcopy->SetLineColor(i+2);
		hcopy->SetFillColor(0);
	        name = TString("order ");
	        name += i;
	        leg->AddEntry(hcopy,name.Data(),"l");
	    }
	}
	leg->Draw();
    }
    
     // Lateral distribution of all photons as function of their nb of interaction in atmosphere
    else if ( option == "mscatt_history" ) {
        th = GetMScattHistory("rayleigh"); 
        th1 = GetMScattHistory("reflected"); 
        th2 = GetMScattHistory("clouds"); 
        th3 = GetMScattHistory("aerosols"); 
        th4 = GetMScattHistory("norm"); 
        if ( th && th4) {
	    th->Divide(th4);
	    th->SetBarWidth(0.225);
	    th->SetBarOffset(0.275);
	    th->SetFillColor(30);
	    th->Draw("bar2");
	    th->SetStats(0);
	    th->GetXaxis()->SetNdivisions(fMaxScatOrder);
	    th->GetYaxis()->SetRangeUser(0,1.1);
	    if(th1) {
	        th1->Divide(th4);
	        th1->SetBarWidth(0.225);
	        th1->SetBarOffset(0.05);
	        th1->SetFillColor(46);
	        th1->Draw("bar2same");
	    }
	    if(th2) {
		th2->Divide(th4);
		th2->SetBarWidth(0.225);
		th2->SetBarOffset(0.5);
		th2->SetFillColor(41);
		th2->Draw("bar2same");
            }
	    if(th3) {
		th3->Divide(th4);
		th3->SetBarWidth(0.225);
		th3->SetBarOffset(0.725);
		th3->SetFillColor(38);
		th3->Draw("bar2same");
	    }
	    TLegend *legend = new TLegend(0.55,0.65,0.76,0.82);
	    legend->AddEntry(th1,"Reflected","f");
	    legend->AddEntry(th,"Rayleigh","f");
	    legend->AddEntry(th2,"Clouds","f");
	    legend->AddEntry(th3,"Aerosols","f");
	    legend->Draw();
	} 
    }
    
}

//_____________________________________________________________________________
 void EAtmosphereHistoPainter::NextFrame() {
    //
    // Switch to the next frame and update histograms
    // 

    fFrame++;
    UpdateHistos();
}

//______________________________________________________________________________
 void EAtmosphereHistoPainter::UpdateHistos() {
    //
    // Update hitograms at the current frame.
    //

    if ( IsZombie() ) return;

    Int_t n;
    if ( fFrame == -1 )
        n = fAtmosphere->GetNumBunches();
    else
        n = (Int_t)(((Float_t)fAtmosphere->GetNumBunches()/fNumFrames)*fFrame);


    TH1F* thNph    = GetTimeHisto("Bunch_Nph");
    // TH1F* thNph_F   = GetTimeHisto("Bunch_Nph_F");
    // TH1F* thNph_C = GetTimeHisto("Bunch_Nph_C");

    for (Int_t i=0; i<fAtmosphere->GetNumBunches() ; i++) {
        Double_t Nph = fAtmosphere->GetBunch(i)->GetWeight();

        thNph->SetBinContent(i+1, i < n ? Nph : 0 );
	//	thNph_F->SetBinContent(i+1, i < n ? Nph_F : 0 );
        //thNph_C->SetBinContent(i+1, i < n ? Nph_C : 0 );
    }

    //thNph->SetMaximum( NphMax*1.1 );
    //thNph_F->SetMaximum( fNphFMax*1.1 );
    //thNph_C->SetMaximum( fNphCMax*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