// $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 );
}