// $Id: ETreePainter.cc,v 1.20 2005/10/19 12:10:26 moreggia Exp $
// Author: Anne Stutz 2005/03/03
// $Id: ETreePainter.cc,v 1.20 2005/10/19 12:10:26 moreggia Exp $
// Author: Anne Stutz 2005/03/03
/*****************************************************************************
* ESAF: Euso Simulation and Analysis Framework *
* *
* Id: ETreePainter *
* Package: <packagename> *
* Coordinator: <coordinator> *
* *
*****************************************************************************/
//_____________________________________________________________________________
//
// ETreePainter
//
// <extensive class description>
//
// Config file parameters
// ======================
//
// <parameter name>: <parameter description>
// -Valid options: <available options>
//
//
// WARNING : the relations between showerstep ID and corresponding bunchofphotons ID are :
// - array[i] for showerstep
// - array[2i] for fluo bunch -> bunchID = 2i+1
// - array[2i+1] for ckov bunch -> bunchID = 2i+2
//
//
//***************************** ONLY FOR ICRC_2005 PRODUCTIONS ***********************************
// *
// FLAG //X1pb -> to find line where X1 definition pb is involved *
// *
// "+fTrack.fX1" <-> slast _showerstep_ convention (X == X-X1) *
// (BUT true Xmax, true X1, same definitions for all config) *
// *
// "NOTHING" <-> unisim convention *
// *
//***************************** ONLY FOR ICRC_2005 PRODUCTIONS ***********************************
//
//
// //TOFIX : GTUmax treatment MUST be changed if fluo backscattering
//
//
#include <vector>
#include "ETreePainter.hh"
#include "ESystemOfUnits.hh"
#include "EAtmosphere.hh"
#include "ETruth.hh"
#include "EShower.hh"
#include "EShowerStep.hh"
#include "EBunchPhotons.hh"
#include "EEvent.hh"
#include <TVirtualGeoTrack.h>
#include <TGeoVolume.h>
#include <TGeoManager.h>
#include <TGeoTube.h>
#include <TTree.h>
#include <TChain.h>
#include <TH1F.h>
#include <TF1.h>
#include <TNtuple.h>
#include <TCut.h>
#include <iostream>
#include <TRefArray.h>
ClassImp(ETreePainter)
extern Double_t Zv(const TVector3&);
using namespace sou;
using namespace std;
// TOOL FUNCTION
//_____________________________________________________________________________
Double_t ETreePainter_my_gaussian_fit(Double_t *x, Double_t *par) {
//
// par[3] allows to fit a sub-range only
//
if(x[0] > par[3]) {
TF1::RejectPoint();
return 0;
}
return par[0]*exp(-0.5*pow((x[0] - par[1]) / par[2],2));
}
//_____________________________________________________________________________
ETreePainter::ETreePainter( TTree* etree ) {
//
// Constructor
//
// if stree already implemented, only Draw() method can be used
if(!etree) Fatal("ETreePainter ctor","TTree given in argument is NULL");
if(strncmp(etree->GetName(),"stree",5) == 0) {
fStree = etree;
Info("ETreePainter()","Better to CHECK VALUES used for [altEUSO, pupil radius, GTU length] at Stree building");
}
else {
Info("ETreePainter()","Building Stree");
fEUSO = TVector3(0,0,430*km);
fRadius = 1250*mm;
fGTU = 2.5*microsecond;
fTree = etree;
fChain = NULL;
fEv = NULL;
fStree = NULL;
Build();
Info("ETreePainter()","BUILT FROM A TREEnUSE pupil radius = %f mm, euso altitude = %f km, GTU length = %f ms",fRadius/mm,fEUSO.Z()/km,fGTU/microsecond);
}
}
//_____________________________________________________________________________
ETreePainter::ETreePainter( TChain* chain ) {
//
// Constructor
//
Info("ETreePainter()","Building Stree");
fEUSO = TVector3(0,0,430*km);
fRadius = 1250*mm;
fGTU = 2.5*microsecond;
fChain = chain;
fTree = chain;
fEv = NULL;
fStree = NULL;
Build();
Info("ETreePainter()","BUILT FROM A CHAINnUSE pupil radius = %f mm, euso altitude = %f km, GTU length = %f ms",fRadius/mm,fEUSO.Z()/km,fGTU/microsecond);
}
//_____________________________________________________________________________
ETreePainter::~ETreePainter() {
//
// Destructor
//
}
//_____________________________________________________________________________
void ETreePainter::Build() {
//
// Build
//
if ( !fTree ) {
Info("Build()","ETree object is NULL. Painter made zombie.");
MakeZombie();
return;
}
fEv = new EEvent();
if ( fChain ) {
SetBranchesStatus(fChain); // MUST be before setbranches()
fEv->SetBranches(fChain);
}
else {
SetBranchesStatus(fTree);
fEv->SetBranches(fTree);
}
fEvTot = (Int_t)fTree->GetEntries();
BuildTree();
for (Int_t i=0; i< fEvTot; i++) {
Printf( "Entry %d: %d bytes read", i, fTree->GetEntry(i) );
fAtmo = fEv->GetAtmosphere();
fTruth = fEv->GetTruth();
fShower = fEv->GetShower();
if(fAtmo && fTruth && fShower) FillTrack();
else Warning("Build()","EAtmosphere or ETruth or EShower object is NULL");
FillTree();
}
}
//_____________________________________________________________________________
void ETreePainter::SetBranchesStatus(TTree* t) {
//
// for fChain or fTree
//
t->SetBranchStatus("*",0);
t->SetBranchStatus("fTrue*",1);
// shower
t->SetBranchStatus("fDir*",1);
t->SetBranchStatus("fInit*",1);
t->SetBranchStatus("fHitGround",1);
t->SetBranchStatus("fNumSteps",1);
t->SetBranchStatus("fStep*",1);
// atmosphere
t->SetBranchStatus("fNumSingles",1);
t->SetBranchStatus("fNumBunches",1);
t->SetBranchStatus("fBunch*",1);
t->SetBranchStatus("fSingle*",1);
}
//_____________________________________________________________________________
void ETreePainter::BuildTree() {
//
// build a simple tree
//
if ( !fStree ) {
fStree = new TTree("stree","a Tree with simple datas");
// truth
fStree->Branch("Energy", &fTrack.fEnergy, "fEnergy/F");
fStree->Branch("Theta", &fTrack.fTheta, "fTheta/F");
fStree->Branch("Phi", &fTrack.fPhi, "fPhi/F");
fStree->Branch("H1", &fTrack.fH1, "fH1/F");
fStree->Branch("InitPosX", &fTrack.fInitPosX, "fInitPosX/F");
fStree->Branch("InitPosY", &fTrack.fInitPosY, "fInitPosY/F");
fStree->Branch("InitPosZ", &fTrack.fInitPosZ, "fInitPosZ/F");
fStree->Branch("X1", &fTrack.fX1, "fX1/F");
fStree->Branch("TrueImpactX", &fTrack.fTrueImpactX, "fTrueImpactX/F");
fStree->Branch("TrueImpactY", &fTrack.fTrueImpactY, "fTrueImpactY/F");
fStree->Branch("TrueImpactZ", &fTrack.fTrueImpactZ, "fTrueImpactZ/F");
fStree->Branch("TrueMaxPosX", &fTrack.fTrueMaxPosX, "fTrueMaxPosX/F");
fStree->Branch("TrueMaxPosY", &fTrack.fTrueMaxPosY, "fTrueMaxPosY/F");
fStree->Branch("TrueMaxPosZ", &fTrack.fTrueMaxPosZ, "fTrueMaxPosZ/F");
fStree->Branch("TrueHmax", &fTrack.fTrueHmax, "fTrueHmax/F");
fStree->Branch("TrueXmax", &fTrack.fTrueXmax, "fTrueXmax/F");
fStree->Branch("TrueDmax", &fTrack.fTrueDmax, "fTrueDmax/F");
fStree->Branch("Latitude", &fTrack.fLatitude, "fLatitude/F");
fStree->Branch("Longitude", &fTrack.fLongitude, "fLongitude/F");
fStree->Branch("Date", &fTrack.fDate, "fDate/F");
fStree->Branch("Hclouds", &fTrack.fHclouds, "fHclouds/F");
fStree->Branch("Clouds_thick", &fTrack.fCloudsThick, "fCloudsThick/F");
fStree->Branch("Clouds_OD", &fTrack.fCloudsOD, "fCloudsOD/F");
// shower
fStree->Branch("Ne", &fTrack.fNe, "fNe/F");
fStree->Branch("Ne2", &fTrack.fNe2, "fNe2/F");
fStree->Branch("Nemax", &fTrack.fNemax, "fNemax/F");
fStree->Branch("Xmax_show", &fTrack.fXmaxshow, "fXmaxshow/F");
fStree->Branch("Hmax_show", &fTrack.fHmaxshow, "fHmaxshow/F");
//fStree->Branch("Width_show", &fTrack.fWidthshow, "fWidthshow/F");
//fStree->Branch("Width2_show", &fTrack.fWidth2show, "fWidth2show/F");
// photons at creation
fStree->Branch("Nph_fluo", &fTrack.fNph_f, "fNph_f/F");
fStree->Branch("Nph2_fluo", &fTrack.fNph2_f, "fNph2_f/F");
fStree->Branch("Nph_cer", &fTrack.fNph_c, "fNph_c/F");
fStree->Branch("Nph2_cer", &fTrack.fNph2_c, "fNph2_c/F");
fStree->Branch("Nmax_fluo", &fTrack.fNmax_f, "fNmax_f/F");
fStree->Branch("Nmax_fluo_L", &fTrack.fNmax_f_L, "fNmax_f_L/F");
fStree->Branch("Nmax_cer", &fTrack.fNmax_c, "fNmax_c/F");
fStree->Branch("Hmax_fluo", &fTrack.fHmax_f, "fHmax_f/F");
fStree->Branch("Hmax_cer", &fTrack.fHmax_c, "fHmax_c/F");
fStree->Branch("Width_fluo", &fTrack.fWidth_f, "fWidth_f/F");
fStree->Branch("Width2_fluo", &fTrack.fWidth2_f, "fWidth2_f/F");
fStree->Branch("Xmaxph_fluo", &fTrack.fXmaxph_f, "fXmaxph_f/F");
fStree->Branch("Xmaxph_cer", &fTrack.fXmaxph_c, "fXmaxph_c/F");
fStree->Branch("Yield_f", &fTrack.fYieldmean, "fYieldmean/F");
fStree->Branch("Yield_fluomax", &fTrack.fYieldHmax, "fYieldHmax/F");
fStree->Branch("Wlmean_f", &fTrack.fWlmean_f, "fWlmean_f/F");
fStree->Branch("Wlmean_c", &fTrack.fWlmean_c, "fWlmean_c/F");
fStree->Branch("Wlmean_tot", &fTrack.fWlmean_tot, "fWlmean_tot/F");
// on pupil
fStree->Branch("FluoInOmega", &fTrack.fFluoInOmega, "fFluoInOmega/F");
fStree->Branch("CkovInOmega", &fTrack.fCkovInOmega, "fCkovInOmega/F");
fStree->Branch("Nph_fluo_pupil", &fTrack.fNph_f_p, "fNph_f_p/F");
fStree->Branch("Nph_cer_dir_pupil", &fTrack.fNph_cdirect_p, "fNph_cdirect_p/F");
fStree->Branch("Nph2_fluo_pupil", &fTrack.fNph2_f_p, "fNph2_f_p/F");
fStree->Branch("Nph_cer_air_pupil", &fTrack.fNph_airscat_p, "fNph_airscat_p/F");
fStree->Branch("Nph_cer_cloud_pupil",&fTrack.fNph_cloudscat_p,"fNph_cloudscat_p/F");
fStree->Branch("Nph_cer_refl_pupil", &fTrack.fNph_cr_p, "fNph_cr_p/F");
fStree->Branch("Nmax_fluo2_pupil", &fTrack.fNmax_f2_p, "fNmax_f2_p/F");
fStree->Branch("Nmax_fluo_pupil", &fTrack.fNmax_f_p, "fNmax_f_p/F");
fStree->Branch("Nmax_fluo_pupil_raw",&fTrack.fNmax_f_p_raw, "fNmax_f_p_raw/F");
fStree->Branch("Nmax_tot_pupil", &fTrack.fNmax_tot_p, "fNmax_tot_p/F");
fStree->Branch("Nmax_tot_pupil_raw", &fTrack.fNmax_tot_p_raw, "fNmax_tot_p_raw/F");
fStree->Branch("Nmax_cr_pupil", &fTrack.fNmax_cr_p, "fNmax_cr_p/F");
fStree->Branch("Nmax_crtot_pupil", &fTrack.fNmax_crtot_p, "fNmax_crtot_p/F");
fStree->Branch("Width_fluo2_pupil", &fTrack.fWidth_f2_p, "fWidth_f2_p/F");
fStree->Branch("Width_fluo_pupil", &fTrack.fWidth_f_p, "fWidth_f_p/F");
fStree->Branch("Width_tot_pupil", &fTrack.fWidth_tot_p, "fWidth_tot_p/F");
fStree->Branch("Yield_fluoGTUmax", &fTrack.fYieldGTUmax, "fYieldGTUmax/F");
fStree->Branch("NbBunchInGTUmax", &fTrack.fNbBunches, "fNbBunches/I");
fStree->Branch("X_GTUmax", &fTrack.fXGTUmax, "fXGTUmax/F");
fStree->Branch("X_GTUmax_tot", &fTrack.fXGTUmax_tot, "fXGTUmax_tot/F");
fStree->Branch("Trans_fluo_atMax", &fTrack.fTransMax, "fTransMax/F");
fStree->Branch("Trans_fluo_GTUmax", &fTrack.fTransGTUmax, "fTransGTUmax/F");
fStree->Branch("Wlmean_f_p", &fTrack.fWlmean_f_p, "fWlmean_f_p/F");
fStree->Branch("Wlmean_cb_p", &fTrack.fWlmean_cb_p, "fWlmean_cb_p/F");
fStree->Branch("Wlmean_cr_p", &fTrack.fWlmean_cr_p, "fWlmean_cr_p/F");
fStree->Branch("Wlmean_tot_p", &fTrack.fWlmean_tot_p, "fWlmean_tot_p/F");
}
}
//_____________________________________________________________________________
void ETreePainter::FillTrack() {
//
// fill fTrack
//
// fill truth data
if(fTruth) {
fTrack.fEnergy = fTruth->GetTrueEnergy(); // en MeV
fTrack.fTheta = fTruth->GetTrueTheta()*TMath::RadToDeg();
fTrack.fPhi = fTruth->GetTruePhi()*TMath::RadToDeg();
fTrack.fH1 = Zv(fTruth->GetTrueInitPos())/ km;
fTrack.fInitPosX = fTruth->GetTrueInitPos().X()/km;
fTrack.fInitPosY = fTruth->GetTrueInitPos().Y()/km;
fTrack.fInitPosZ = fTruth->GetTrueInitPos().Z()/km;
fTrack.fX1 = fTruth->GetTrueX1()/g*cm2;
fTrack.fTrueImpactX = fTruth->GetTrueEarthImpact().X()/km;
fTrack.fTrueImpactY = fTruth->GetTrueEarthImpact().Y()/km;
fTrack.fTrueImpactZ = fTruth->GetTrueEarthImpact().Z()/km;
fTrack.fTrueMaxPosX = fTruth->GetTrueShowerMaxPos().X()/km;
fTrack.fTrueMaxPosY = fTruth->GetTrueShowerMaxPos().Y()/km;
fTrack.fTrueMaxPosZ = fTruth->GetTrueShowerMaxPos().Z()/km;
fTrack.fTrueHmax = Zv(fTruth->GetTrueShowerMaxPos())/km;
fTrack.fTrueXmax = fTruth->GetTrueShowerXMax()/g*cm2;
fTrack.fTrueDmax = (fEUSO - fTruth->GetTrueShowerMaxPos()).Mag() / km;
fTrack.fLatitude = fTruth->GetTrueLatitude();
fTrack.fLongitude = fTruth->GetTrueLongitude();
fTrack.fDate = fTruth->GetTrueDate();
fTrack.fHclouds = fTruth->GetTrueHclouds() /km;
fTrack.fCloudsThick = fTruth->GetTrueCloudsthick() /km;
fTrack.fCloudsOD = fTruth->GetTrueCloudsOD();
}
else Info("FillTrack","ETruth object EMPTY");
if(fShower && fAtmo) {
TH1F* th(0);
TF1* mygaus = 0;
Double_t par[10];
Double_t binwidth = 0.;
Double_t TimeMin = 0.;
Double_t TimeMax = 0.;
Int_t maxbin = 0;
Double_t maxpos = 0;
Int_t leftbin = 0;
Int_t rightbin = 0;
Int_t NbinsAroundMax = 0; // to avoid step effects when searching Xmax, Hmax etc.. (fit around maximum)
////////////////////////////////////////////////////////////
////////////////////// SHOWER DATA /////////////////////////
////////////////////////////////////////////////////////////
// init
Int_t numsteps = fShower->GetNumSteps();
EShowerStep* shstep = 0;
Int_t indexofmax = 0; //TOFIX : also used for fTransMax
Float_t integrate = 0.;
Float_t integrate2 = 0.;
Float_t NORM_weight = 0.;
Bool_t flag = false; // to stop integrate2 calculation
fTrack.fYieldmean = 0.;
// gaussian fit of shower profile (longitudinal - in depth) around its maxumum -> Xmax_shower Nemax
th = BuildLongitudinalHisto("Shower_Longit_prof_gram");
leftbin = 1;
NbinsAroundMax = Int_t(30. / th->GetBinWidth(th->GetMaximumBin())) + 1; // 30 g/cm2 seems ok
if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
rightbin = th->GetNbinsX() + 1;
if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
mygaus = new TF1("mygaus","gaus",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
//fTrack.fXmaxshow = par[1] + fTrack.fX1; //X1pb
fTrack.fXmaxshow = par[1]; //X1pb
fTrack.fNemax = par[0];
delete th; th = 0;
delete mygaus; mygaus = 0;
// gaussian fit of shower profile (vs altitude - in km) around its maxumum -> Hmax_shower
th = BuildAltHisto("alt_Shower");
leftbin = 1;
if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
rightbin = th->GetNbinsX() + 1;
if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
mygaus = new TF1("mygaus","gaus",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
fTrack.fHmaxshow = par[1];
delete th; th = 0;
delete mygaus; mygaus = 0;
// integrate shower profile - find the showerstep of Ne maximum - calculate mean fluo yield (fYieldmean)
for(Int_t i=0; i<numsteps; i++) {
shstep = fShower->GetStep(i);
integrate += shstep->GetNumElectrons()*(shstep->GetXf() - shstep->GetXi())*cm2/g;
if(!flag) {
integrate2 += shstep->GetNumElectrons()*(shstep->GetXf() - shstep->GetXi())*cm2/g;
//if((shstep->GetXf()/g*cm2 + fTrack.fX1) > fTrack.fXmaxshow) { //X1pb
if((shstep->GetXf()/g*cm2) > fTrack.fXmaxshow) { //X1pb
indexofmax = i;
integrate2 += shstep->GetNumElectrons()*(fTrack.fXmaxshow - shstep->GetXi()*cm2/g);
flag = true;
}
}
fTrack.fYieldmean += fAtmo->GetBunch(2*i)->GetYield()*m * shstep->GetNumElectrons();
NORM_weight += shstep->GetNumElectrons();
}
if(NORM_weight) fTrack.fYieldmean /= NORM_weight;
shstep = fShower->GetStep(indexofmax);
fTrack.fYieldHmax = fAtmo->GetBunch(2*indexofmax)->GetYield()*m; // get fluo yield associated to shower max
fTrack.fNe = integrate;
fTrack.fNe2 = integrate2;
const Int_t indexofmax_cst = indexofmax;
////////////////////////////////////////////////////////////
////////////////////// ATMOSPHERE DATA /////////////////////
////////////////////////////////////////////////////////////
///////////////////// AT PHOTONS CREATION ///////////////////
// from data at creation VS atlitude
// gaussian fit profiles around their maxumum -> Hmax
th = BuildAltHisto("alt_Bunch_Nph_F");
leftbin = 1;
if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
rightbin = th->GetNbinsX() + 1;
if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
mygaus = new TF1("mygaus","gaus",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
fTrack.fHmax_f = par[1];
delete th; th = 0;
delete mygaus; mygaus = 0;
th = BuildAltHisto("alt_Bunch_Nph_C");
leftbin = 1;
if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
rightbin = th->GetNbinsX() + 1;
if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
mygaus = new TF1("mygaus","gaus",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
fTrack.fHmax_c = par[1];
delete th; th = 0;
delete mygaus; mygaus = 0;
// from longitudinal data in grammage
// - nb - HALF nb - nb at max - Xmax - of photons at creation
th = BuildLongitudinalHisto("Fluo_Longit_prof_gram");
leftbin = 1;
if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
rightbin = th->GetNbinsX() + 1;
if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
mygaus = new TF1("mygaus","gaus",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
//fTrack.fXmaxph_f = par[1] + fTrack.fX1; //X1pb
fTrack.fXmaxph_f = par[1]; //X1pb
maxbin = th->FindBin(fTrack.fXmaxph_f);
binwidth = th->GetBinWidth(maxbin);
fTrack.fNmax_f = par[0]/binwidth;
fTrack.fNph_f = th->Integral();
fTrack.fNph2_f = th->Integral(1,maxbin - 1);
//fTrack.fNph2_f += th->GetBinContent(maxbin) * (fTrack.fXmaxph_f - th->GetBinLowEdge(maxbin) - fTrack.fX1)/th->GetBinWidth(maxbin); //X1pb
fTrack.fNph2_f += th->GetBinContent(maxbin) * (fTrack.fXmaxph_f - th->GetBinLowEdge(maxbin))/th->GetBinWidth(maxbin); //X1pb
delete th; th = 0;
delete mygaus; mygaus = 0;
// Nmax_fluo calculated per km
th = BuildLongitudinalHisto("Fluo_Longit_prof_km");
leftbin = 1;
if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
rightbin = th->GetNbinsX() + 1;
if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
mygaus = new TF1("mygaus","gaus",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
Double_t maxx = par[1];
maxbin = th->FindBin(maxx);
binwidth = th->GetBinWidth(maxbin);
fTrack.fNmax_f_L = par[0]/binwidth;
delete th; th = 0;
delete mygaus; mygaus = 0;
th = BuildLongitudinalHisto("Ckov_Longit_prof_gram");
leftbin = 1;
if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
rightbin = th->GetNbinsX() + 1;
if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
mygaus = new TF1("mygaus","gaus",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
//fTrack.fXmaxph_c = par[1] + fTrack.fX1; //X1pb
fTrack.fXmaxph_c = par[1]; //X1pb
maxbin = th->FindBin(fTrack.fXmaxph_c);
binwidth = th->GetBinWidth(maxbin);
fTrack.fNmax_c = par[0]/binwidth;
fTrack.fNph_c = th->Integral();
fTrack.fNph2_c = th->Integral(1,maxbin - 1);
//fTrack.fNph2_c += th->GetBinContent(maxbin) * (fTrack.fXmaxph_c - th->GetBinLowEdge(maxbin) - fTrack.fX1)/th->GetBinWidth(maxbin); //X1pb
fTrack.fNph2_c += th->GetBinContent(maxbin) * (fTrack.fXmaxph_c - th->GetBinLowEdge(maxbin))/th->GetBinWidth(maxbin); //X1pb
delete th; th = 0;
delete mygaus; mygaus = 0;
// width - HALF width - of fluo profile (km)
th = BuildLongitudinalHisto("Fluo_Longit_prof_km");
leftbin = 1;
if(leftbin < (th->GetMaximumBin() - NbinsAroundMax)) leftbin = th->GetMaximumBin() - NbinsAroundMax;
rightbin = th->GetNbinsX() + 1;
if(rightbin > (th->GetMaximumBin() + NbinsAroundMax)) rightbin = th->GetMaximumBin() + NbinsAroundMax;
mygaus = new TF1("mygaus","gaus",th->GetBinLowEdge(leftbin),th->GetBinLowEdge(rightbin)+th->GetBinWidth(rightbin));
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
maxx = par[1];
maxbin = th->GetMaximumBin();
maxpos = th->GetBinContent(maxbin);
leftbin = 0;
rightbin = th->GetNbinsX();
Bool_t stopleft(false), stopright(false);
for(Int_t i=0; i<th->GetNbinsX(); i++) {
if((th->GetBinContent(i+1) < maxpos/2.) && (!stopleft)) leftbin = i+1;
else stopleft = true; // bin found
if((th->GetBinContent(th->GetNbinsX() - i) < maxpos/2.) && (!stopright)) rightbin = th->GetNbinsX() - i;
else stopright = true;
}
Double_t leftx = 0.; // for linear interpolation within a step
Double_t rightx = 0.;
if(rightbin > 1) rightx = th->GetBinLowEdge(rightbin) + (maxpos/2. - th->GetBinContent(rightbin-1))/(th->GetBinContent(rightbin-1) - th->GetBinContent(rightbin))*th->GetBinWidth(rightbin);
if(leftbin < th->GetNbinsX()) leftx = th->GetBinLowEdge(leftbin+1) + (maxpos/2. - th->GetBinContent(leftbin+1))/(th->GetBinContent(leftbin) - th->GetBinContent(leftbin+1))*th->GetBinWidth(leftbin);
fTrack.fWidth_f = rightx - leftx;
if(fTrack.fWidth_f < 0) {
Warning("FillTrack()","WIDTH PB rigthbin = %d, maxbin = %d, Nbins = %d",rightbin,maxbin,th->GetNbinsX());
fTrack.fWidth_f = th->GetBinCenter(maxbin) - th->GetBinCenter(leftbin);
}
fTrack.fWidth2_f = maxx - leftx;
delete th; th = 0;
////////////////////// ON PUPIL /////////////////////////////
// from data on pupil VS time
// fluo
Double_t tmax_pupil = 0.;
th = BuildTofHisto("Single_Nph_F_P_t");
fTrack.fNph_f_p = th->Integral();
fTrack.fNmax_f_p_raw = th->GetMaximum();
TimeMin = th->GetBinLowEdge(1);
TimeMax = th->GetBinLowEdge(th->GetXaxis()->GetLast())+th->GetBinWidth(th->GetXaxis()->GetLast());
mygaus = new TF1("mygaus","gaus",TimeMin,TimeMax);
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
fTrack.fNph2_f_p = th->Integral(1,th->FindBin(par[1])); // first HALF of profile
fTrack.fNmax_f_p = par[0]; // histobinwidth = GTU -> so Nmax is per GTU
tmax_pupil = par[1]*microsecond; // records absolute Tmax in ms for below search of fXGTUmax - fYieldGTUmax etc..
fTrack.fWidth_f_p = par[2]/(fGTU/microsecond);
delete th; th = 0;
delete mygaus; mygaus = 0;
// direct ckov
th = BuildTofHisto("Single_Nph_Cdir_P_t");
fTrack.fNph_cdirect_p = th->Integral();
delete th; th = 0;
// ground reflected ckov
th = BuildTofHisto("Single_Nph_CR_P_t");
Int_t localbin = th->GetMaximumBin();
Double_t time_cr_max = th->GetBinCenter(localbin); // used here below for fTrack.fNmax_crtot_p
fTrack.fNph_cr_p = th->Integral();
fTrack.fNmax_cr_p = th->GetBinContent(localbin)/(fGTU/microsecond);
delete th; th = 0;
// air scattered ckov
th = BuildTofHisto("Single_Nph_CB_P_t");
fTrack.fNph_airscat_p = th->Integral();
fTrack.fNmax_crtot_p = fTrack.fNmax_cr_p + th->GetBinContent(th->FindBin(time_cr_max))/(fGTU/microsecond);
delete th; th = 0;
// clouds scattered ckov
th = BuildTofHisto("Single_Nph_CB_cloud_P_t");
fTrack.fNph_cloudscat_p = th->Integral();
delete th; th = 0;
// direct fluo + scattered ckov ->
// 1. fit total
// 2. fit only the first profile half (to substract a part of ckov component)
Double_t tmax_tot_pupil = 0.;
th = BuildTofHisto("Single_Nph_F-CB_P_t");
fTrack.fNmax_tot_p_raw = th->GetMaximum();
TimeMin = th->GetBinLowEdge(1);
TimeMax = th->GetBinLowEdge(th->GetXaxis()->GetLast())+th->GetBinWidth(th->GetXaxis()->GetLast());
mygaus = new TF1("mygaus","gaus",TimeMin,TimeMax);
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
fTrack.fNmax_tot_p = par[0]; // histobinwidth = GTU -> so Nmax is per GTU
fTrack.fWidth_tot_p = par[2]/(fGTU/microsecond);
tmax_tot_pupil = par[1]*microsecond; // records absolute Tmax in ms for below search of fXGTUmax_tot
delete mygaus; mygaus = 0;
mygaus = new TF1("mygaus",ETreePainter_my_gaussian_fit,TimeMin,TimeMax,4);
mygaus->SetParameters(par[0],par[1],par[2],par[1]);
mygaus->SetParLimits(0,0,2*par[0]);
mygaus->SetParLimits(1,TimeMin,TimeMax);
mygaus->SetParLimits(2,0.,(TimeMax - TimeMin));
mygaus->SetParLimits(3,par[1],par[1]); // max_tot bound -> for fitting the FIRST half of total profile
th->Fit("mygaus","RQON");
mygaus->GetParameters(par);
fTrack.fNmax_f2_p = par[0]; // histobinwidth = GTU -> so Nmax is per GTU
fTrack.fWidth_f2_p = par[2]/(fGTU/microsecond);
delete th; th = 0;
delete mygaus; mygaus = 0;
///////////////// GTUmax RELATED + singles related (mean trans)/////////////////////////
// Calculate fYieldGTUmax - fXGTUmax - fXGTUmax_tot
// fFluoInOmega and fCkovInOmega : get nb of photons BEFORE last TRANSMISSION
vector<Int_t> list_bunchID; // ID list of bunches in GTUmax
vector<Int_t> list_bunchID_new; // ID list of bunches in GTUmax, one ID appear only once
list_bunchID.clear();
list_bunchID_new.clear();
vector<Int_t> list_bunchID_tot; // same object as here above, but for GTUmax_tot
vector<Int_t> list_bunchID_tot_new;
list_bunchID_tot.clear();
list_bunchID_tot_new.clear();
Int_t ns = fAtmo->GetNumSingles();
fTrack.fFluoInOmega = 0.;
fTrack.fCkovInOmega = 0.;
// only FLUO (scattered ckov used only to get shift in time on pupil (tmax_tot_pupil)
for (Int_t i=0; i<ns; i++) {
ESinglePhoton* s = fAtmo->GetSingle(i);
if(s->GetType() == 0) {
fTrack.fFluoInOmega++;
if ( !s->IsAbsorbed() ) {
// get Id of bunches which have generated pupil photons with arrival time within : tmax_pupil +/- (GTU/2)
if ( fabs((s->GetTof() + s->GetDate()) - tmax_pupil) <= 0.5*fGTU) {
list_bunchID.push_back(Int_t(s->GetBunchId()));
}
// idem but within : tmax_pupil_tot +/- (GTU/2)
if ( fabs((s->GetTof() + s->GetDate()) - tmax_tot_pupil) <= 0.5*fGTU) {
list_bunchID_tot.push_back(Int_t(s->GetBunchId()));
}
}
}
else if(s->GetType() == 1) fTrack.fCkovInOmega++;
}
// to get each ID only once in the ID-lists
Bool_t flagg = false; // true if value already exists in new list
for(size_t i=0; i<list_bunchID.size(); i++) {
flagg = false;
for(size_t j=0; j<list_bunchID_new.size(); j++) {
if(list_bunchID_new[j] == list_bunchID[i]) flagg = true;
}
if(!flagg) list_bunchID_new.push_back(list_bunchID[i]);
}
// idem for ID-lists_tot
for(size_t i=0; i<list_bunchID_tot.size(); i++) {
flagg = false;
for(size_t j=0; j<list_bunchID_tot_new.size(); j++) {
if(list_bunchID_tot_new[j] == list_bunchID_tot[i]) flagg = true;
}
if(!flagg) list_bunchID_tot_new.push_back(list_bunchID_tot[i]);
}
fTrack.fNbBunches = Int_t(list_bunchID_new.size());
// to get mean transmission (nb of bunches considered == fNbBunches, thus given by GTUmax <-> X translation)
// - at Xmax
// - at GTUmax
Float_t ph_atmax = 0.;
Float_t ph_trans_atmax = 0.;
Float_t ph_atGTUmax = 0.;
Float_t ph_trans_atGTUmax = 0.;
Int_t idty = 0;
/* // When TRefArray no implemented //DELETE
for (Int_t i=0; i<ns; i++) {
ESinglePhoton* s = fAtmo->GetSingle(i);
idty = Int_t(s->GetBunchId());
// fluo only
if(s->GetType() == 0) {
for(size_t j=0; j<list_bunchID_new.size(); j++) {
if( list_bunchID_new[j] == idty ) {
ph_atGTUmax++;
if( !s->IsAbsorbed() ) ph_trans_atGTUmax++;
break;
}
}
// if idty is around bunch(Xmax) ID <-> if idty is between maxID-2*(nb/2) and maxID+2*(nb/2 +1 if nb is odd)
// "2*" because one fluo bunch step further <-> two bunches ID step further
// "nb" is the nb of bunches within GTUmax"
// "2*indexofmax_cst" is the ID of showermax bunch
if( ((2*indexofmax_cst+1 - 2*(fTrack.fNbBunches/2)) <= idty) || (idty <= (2*indexofmax_cst+1 + 2*(fTrack.fNbBunches/2 + fTrack.fNbBunches%2)) ) ) {
ph_atmax++;
if ( !s->IsAbsorbed() ) ph_trans_atmax++;
}
}
}
*/
// GTUmax photons
EBunchPhotons* b = 0;
ESinglePhoton* s = 0;
// Bunch ID points here to FLUO bunches (look at above lines)
for(size_t j=0; j<list_bunchID_new.size(); j++) {
b = fAtmo->GetBunch(list_bunchID_new[j]-1); // j-1 because ID j <-> array[j-1]
if(b->GetType() != 0) Warning("FillTrack","1. SHOULD be FLUO bunch !!!");
ph_atGTUmax += b->GetNumDirect("fluo");
TRefArrayIter iter(b->GetSinglesArray());
while((s = (ESinglePhoton*)iter.Next()))
if(!s->IsAbsorbed()) ph_trans_atGTUmax++; //TOFIX : to be changed if fluo backscattering
}
// (shower-)Xmax photons (nb of bunches considered ('fNbBunches') comes from GTUmax studies, to get ~same statistics in both case)
// if idty is around bunch(Xmax) ID <-> if idty is between maxID-2*(nb/2) and maxID+2*(nb/2 +1 if nb is odd)
// "2*" because one fluo bunch step further <-> two bunches ID step further
// "nb" is the nb of bunches within GTUmax"
// "2*indexofmax_cst" is the ID of showermax bunch
idty = 2*indexofmax_cst +1 - fTrack.fNbBunches - fTrack.fNbBunches%2;
if(idty < 0) {
Warning("FillTrack()","TRONCATED shower ? -> indexofmax = %d, nbBunches = %d",indexofmax_cst,fTrack.fNbBunches);
while(idty <= 0) idty ++;
Warning("FillTrack()","so idty set to = %d",idty);
}
for(Int_t w=0; w<fTrack.fNbBunches; w++) {
idty += 2*w;
if((idty-1) > fAtmo->GetNumBunches()) {
Warning("FillTrack()","idty is too high -> idty = %d",idty);
fTrack.fNbBunches = w;
Warning("FillTrack()","so fTrack.fNbBunches set to = %d",fTrack.fNbBunches);
break;
}
b = fAtmo->GetBunch(idty-1); // idty-1 because ID idty <-> array[idty-1]
if(!b) { // occurs when end of list is reached
Warning("FillTrack","fAtmo->GetBunch(%d) returned NULL pointer",idty-1);
if(ph_atmax*ph_trans_atmax == 0) { ph_atmax=-1; ph_trans_atmax=-1; }
continue;
}
if(b->GetType() != 0) Warning("FillTrack","2. SHOULD be FLUO bunch !!!");
ph_atmax += b->GetNumDirect("fluo");
TRefArrayIter iter(b->GetSinglesArray());
while((s = (ESinglePhoton*)iter.Next()))
if(!s->IsAbsorbed()) ph_trans_atmax++; //TOFIX : to be changed if fluo backscattering
}
fTrack.fTransMax = ph_trans_atmax/ph_atmax;
fTrack.fTransGTUmax = ph_trans_atGTUmax/ph_atGTUmax;
fTrack.fYieldGTUmax = 0.;
fTrack.fXGTUmax = 0.;
for(size_t j=0; j<list_bunchID_new.size(); j++) {
// mean yield
fTrack.fYieldGTUmax += fAtmo->GetBunch(list_bunchID_new[j] - 1)->GetYield()*m;
// get grammage position from corresponding shower step (bunchID is necessarily of FLUO type)
// XGTUmax is taken at the middle between first and last bunch involved
//fTrack.fXGTUmax += 0.5*(fShower->GetStep((list_bunchID_new[j] - 1)/2)->GetXf() + fShower->GetStep((list_bunchID_new[j] - 1)/2)->GetXi())*cm2/g + fTrack.fX1; //X1pb
fTrack.fXGTUmax += 0.5*(fShower->GetStep((list_bunchID_new[j] - 1)/2)->GetXf() + fShower->GetStep((list_bunchID_new[j] - 1)/2)->GetXi())*cm2/g; //X1pb
}
if(list_bunchID_new.size()) {
fTrack.fYieldGTUmax /= Double_t(list_bunchID_new.size());
fTrack.fXGTUmax /= Double_t(list_bunchID_new.size());
}
// idem for fXGTUmax_tot
fTrack.fXGTUmax_tot = 0.;
for(size_t j=0; j<list_bunchID_tot_new.size(); j++) {
// get grammage position from corresponding shower step (bunchID is necessarily of FLUO type)
// XGTUmax is taken at the middle between first and last bunch involved
//fTrack.fXGTUmax_tot += 0.5*(fShower->GetStep((list_bunchID_tot_new[j] - 1)/2)->GetXf() + fShower->GetStep((list_bunchID_tot_new[j] - 1)/2)->GetXi())*cm2/g + fTrack.fX1; //X1pb
fTrack.fXGTUmax_tot += 0.5*(fShower->GetStep((list_bunchID_tot_new[j] - 1)/2)->GetXf() + fShower->GetStep((list_bunchID_tot_new[j] - 1)/2)->GetXi())*cm2/g; //X1pb
}
if(list_bunchID_tot_new.size()) fTrack.fXGTUmax_tot /= Double_t(list_bunchID_tot_new.size());
/////////////////// WAVELENGTH //////////////////////////
// from wavelength spectra BEFORE and AFTER transmission
// before
th = BuildWlHisto("Bunch_Wl_F");
fTrack.fWlmean_f = th->GetMean();
delete th; th = 0;
th = BuildWlHisto("Bunch_Wl_C");
fTrack.fWlmean_c = th->GetMean();
delete th; th = 0;
th = BuildWlHisto("Bunch_Wl_tot");
fTrack.fWlmean_tot = th->GetMean();
delete th; th = 0;
// after
th = BuildWlHisto("Single_Wl_F_P");
fTrack.fWlmean_f_p = th->GetMean();
delete th; th = 0;
th = BuildWlHisto("Single_Wl_CB_P");
fTrack.fWlmean_cb_p = th->GetMean();
delete th; th = 0;
th = BuildWlHisto("Single_Wl_CR_P");
fTrack.fWlmean_cr_p = th->GetMean();
delete th; th = 0;
th = BuildWlHisto("Single_Wl_tot_P");
fTrack.fWlmean_tot_p = th->GetMean();
delete th; th = 0;
delete mygaus; mygaus = 0;
}
else Info("FillTrack","EAtmosphere OR EShower object EMPTY");
}
//_____________________________________________________________________________
void ETreePainter::FillTree() {
//
// fill stree
//
BuildTree(); // should be useless
fStree->Fill();
}
//_____________________________________________________________________________
void ETreePainter::PrepareHistos(Option_t *opt , TCut& drawcut, Option_t *drawopt, Float_t min, Float_t max, UInt_t Nbins) {
//
// define histos features : range, binwidth..
//
if(!fStree) {
cout <<"NO STREE FOR PLOTS !!n";
return;
}
TString name("h_");
TString option(opt);
name += option.Data();
TString drawopt_str(drawopt);
drawcut += drawopt_str.Data();
if(Nbins > 0 && fabs(min*max) > 0) new TH1F(name.Data(),name.Data(),Nbins,min,max);
else cout <<"nHistogram '"<<option.Data()<<"' has not been configuredn";
}
//_____________________________________________________________________________
void ETreePainter::DrawHistos( Option_t *opt , TCut& drawcut, Option_t *drawopt, Float_t min, Float_t max, UInt_t Nbins) {
//
// local cpy of ETreePainter method
//
PrepareHistos(opt,drawcut,drawopt,min,max,Nbins);
TString option(opt);
TString drawopt_str(drawopt);
drawcut += drawopt_str.Data();
// **************************** RAW HISTOS ******************************
// truth histo
if(option == "Energy") fStree->Draw("TMath::Log10(Energy)>>h_Energy",drawcut,"");
else if(option == "Theta") fStree->Draw("Theta>>h_Theta",drawcut,"");
else if(option == "Phi") fStree->Draw("Phi>>h_Phi",drawcut,"");
else if(option == "H1") fStree->Draw("H1>>h_H1",drawcut,"");
else if(option == "Hmax") fStree->Draw("TrueHmax>>h_Hmax",drawcut,"");
else if(option == "X1") fStree->Draw("X1>>h_X1",drawcut,"");
else if(option == "Xmax") fStree->Draw("TrueXmax>>h_Xmax",drawcut,"");
else if(option == "Xmax_X1") fStree->Draw("TrueXmax - X1>>h_Xmax_X1",drawcut,"");
// shower histo
else if(option == "Nemax") fStree->Draw("Nemax>>h_Nemax",drawcut,"");
// atmosphere histos
else if(option == "Nph_Fluo") fStree->Draw("Nph_fluo>>h_Nph_Fluo",drawcut,"");
else if(option == "Nph_Cer") fStree->Draw("Nph_cer>>h_Nph_Cer",drawcut,"");
else if(option == "Nmax_Fluo") fStree->Draw("Nmax_fluo>>h_Nmax_Fluo",drawcut,"");
else if(option == "Nmax_Cer") fStree->Draw("Nmax_cer>>h_Nmax_Cer",drawcut,"");
else if(option == "Hmax_Fluo") fStree->Draw("Hmax_fluo>>h_Hmax_Fluo",drawcut,"");
else if(option == "Hmax_Cer") fStree->Draw("Hmax_cer>>h_Hmax_Cer",drawcut,"");
else if(option == "Xmaxph_Fluo") fStree->Draw("Xmaxph_fluo>>h_Xmaxph_Fluo",drawcut,"");
else if(option == "Xmaxph_Cer") fStree->Draw("Xmaxph_cer>>h_Xmaxph_Cer",drawcut,"");
else if(option == "Wlmean_f") fStree->Draw("Wlmean_f>>h_Wlmean_f",drawcut,"");
else if(option == "Wlmean_c") fStree->Draw("Wlmean_c>>h_Wlmean_c",drawcut,"");
else if(option == "Wlmean_tot") fStree->Draw("Wlmean_tot>>h_Wlmean_tot",drawcut,"");
// on pupil histos
else if(option == "Nph_Fluo_P") fStree->Draw("Nph_fluo_pupil>>h_Nph_Fluo_P",drawcut,"");
else if(option == "Nmax_Fluo_P") fStree->Draw("Nmax_fluo_pupil>>h_Nmax_Fluo_P",drawcut,"");
else if(option == "Sig_Fluo_P") fStree->Draw("Width_fluo_pupil>>h_Sig_Fluo_P",drawcut,"");
//else if(option == "Nph_CB_P") fStree->Draw("Nph_cer_back_pupil>>h_Nph_CB_P",drawcut,"");
else if(option == "Nph_CB_P") fStree->Draw("Nph_cer_air_pupil>>h_Nph_CB_P",drawcut,"");
else if(option == "Nph_CR_P") fStree->Draw("Nph_cer_refl_pupil>>h_Nph_CR_P",drawcut,"");
//else if(option == "Nph_tot_P") fStree->Draw("Nph_fluo_pupil+Nph_cer_back_pupil+Nph_cer_refl_pupil>>h_Nph_tot_P",drawcut,"");
else if(option == "Nph_tot_P") fStree->Draw("Nph_fluo_pupil+Nph_cer_air_pupil+Nph_cer_refl_pupil>>h_Nph_tot_P",drawcut,"");
else if(option == "Wlmean_f_p") fStree->Draw("Wlmean_f_p>>h_Wlmean_f_p",drawcut,"");
else if(option == "Wlmean_cb_p") fStree->Draw("Wlmean_cb_p>>h_Wlmean_cb_p",drawcut,"");
else if(option == "Wlmean_cr_p") fStree->Draw("Wlmean_cr_p>>h_Wlmean_cr_p",drawcut,"");
else if(option == "Wlmean_tot_p") fStree->Draw("Wlmean_tot_p>>h_Wlmean_tot_p",drawcut,"");
// **************************** CORRELATIONS ******************************
// show EAS crash on ground -> it is better to use half integrals
else if(option == "Nph2_Over_Nph_Hmax_fluo") fStree->Draw("Nph2_fluo/Nph_fluo:Hmax_show>>h_Nph2_Over_Nph_Hmax_fluo",drawcut,"");
else if(option == "Nph2_Over_Nph_Theta_fluo") fStree->Draw("Nph2_fluo/Nph_fluo:Theta>>h_Nph2_Over_Nph_Theta_fluo",drawcut,"");
else if(option == "Nph2_Over_Nph_Hmax_cer") fStree->Draw("Nph2_cer/Nph_cer:Hmax_show>>h_Nph2_Over_Nph_Hmax_cer",drawcut,"");
else if(option == "Nph2_Over_Nph_Theta_cer") fStree->Draw("Nph2_cer/Nph_cer:Theta>>h_Nph2_Over_Nph_Theta_cer",drawcut,"");
else if(option == "Ne2_Over_Ne_Hmax") fStree->Draw("Ne2/Ne:Hmax_show>>h_Ne2_Over_Ne_Hmax",drawcut,"");
else if(option == "Ne2_Over_Ne_Theta") fStree->Draw("Ne2/Ne:Theta>>h_Ne2_Over_Ne_Theta",drawcut,"");
else if(option == "Width2_Over_Width_Hmax_fluo") fStree->Draw("Width2_fluo/Width_fluo:Hmax_show>>h_Width2_Over_Width_Hmax_fluo",drawcut,"");
else if(option == "Width2_Over_Width_Theta_fluo") fStree->Draw("Width2_fluo/Width_fluo:Theta>>h_Width2_Over_Width_Theta_fluo",drawcut,"");
// shower correlations
// hmax -- Xmax
else if(option == "Hmax_Theta") fStree->Draw("Hmax_fluo:Theta>>h_Hmax_Theta",drawcut,"");
else if(option == "Xmax_Energy") fStree->Draw("TrueXmax:log10(Energy)>>h_Xmax_Energy",drawcut,"");
// energy
else if(option == "Nemax_Energy") fStree->Draw("log10(Nemax):log10(Energy)>>h_Nemax_Energy",drawcut,"profs");
else if(option == "Nemax_Energy_linear") fStree->Draw("Nemax:log10(Energy)>>h_Nemax_Energy_linear",drawcut,"goffprofs");
else if(option == "Ne2_Energy") fStree->Draw("log10(Ne2):log10(Energy)>>h_Ne2_Energy",drawcut,"profs");
else if(option == "Ne2_Energy_linear") fStree->Draw("Ne2:log10(Energy)>>h_Ne2_Energy_linear",drawcut,"goffprofs");
// theta
else if(option == "Ne_Theta") fStree->Draw("Ne/Energy:Theta>>h_Ne_Theta",drawcut,"");
else if(option == "Ne2_Theta") fStree->Draw("Ne2/Energy:Theta>>h_Ne2_Theta",drawcut,"");
else if(option == "Nemax_Theta") fStree->Draw("Nemax/Energy:Theta>>h_Nemax_Theta",drawcut,"");
// photons in atmosphere and shower correlations
else if(option == "Nphmaxfluo_Nemax") fStree->Draw("log10(Nmax_fluo):log10(Nemax)>>h_Nphmaxfluo_Nemax",drawcut,"profs");
else if(option == "Nphmaxfluo_Nemax_linear") fStree->Draw("Nmax_fluo*cos(Theta*TMath::Pi()/180):log10(Nemax)>>h_Nphmaxfluo_Nemax_linear",drawcut,"profs");
else if(option == "Nph2fluo_Ne2") fStree->Draw("log10(Nph2_fluo):log10(Ne2)>>h_Nph2fluo_Ne2",drawcut,"profs");
else if(option == "Nph2fluo_Ne2_linear") fStree->Draw("Nph2_fluo*cos(Theta*TMath::Pi()/180):log10(Ne2)>>h_Nph2fluo_Ne2_linear",drawcut,"profs");
else if(option == "Nph2ckov_Ne2") fStree->Draw("log10(Nph2_cer):log10(Ne2)>>h_Nph2ckov_Ne2",drawcut,"profs");
else if(option == "Nph2ckov_Ne2_linear") fStree->Draw("Nph2_cer:log10(Ne2)>>h_Nph2ckov_Ne2_linear",drawcut,"profs");
else if(option == "Nph2fluo_Nemax") fStree->Draw("log10(Nph2_fluo):log10(Nemax)>>h_Nph2fluo_Nemax",drawcut,"profs");
else if(option == "Nph2fluo_Nemax_linear") fStree->Draw("Nph2_fluo*cos(Theta*TMath::Pi()/180):log10(Nemax)>>h_Nph2fluo_Nemax_linear",drawcut,"profs");
else if(option == "Width2_fluo_Hmax_fluo") fStree->Draw("Width2_fluo:Hmax_fluo>>h_Width2_fluo_Hmax_fluo",drawcut,"");
else if(option == "Width2_fluo_Theta") fStree->Draw("Width2_fluo:Theta>>h_Width2_fluo_Theta",drawcut,"");
else if(option == "Xmaxfluo_Xmax_Theta") fStree->Draw("Xmaxph_fluo - Xmax_show:Theta>>h_Xmaxfluo_Xmax_Theta",drawcut,"");
else if(option == "Xmaxfluo_Xmax_Energy") fStree->Draw("Xmaxph_fluo - Xmax_show:log10(Energy)>>h_Xmaxfluo_Xmax_Energy",drawcut,"");
else if(option == "Hmaxfluo_Hmax_Theta") fStree->Draw("Hmax_fluo - Hmax_show:Theta>>h_Hmaxfluo_Hmax_Theta",drawcut,"");
else if(option == "Xmaxfluo_Xmax_X1") fStree->Draw("Xmaxph_fluo - Xmax_show:X1>>h_Xmaxfluo_Xmax_X1",drawcut,"");
// photons in atmosphere correlations VS Energy
else if(option == "Nmax_Energy_fluo") fStree->Draw("log10(Nmax_fluo*cos(Theta*TMath::Pi()/180)):log10(Energy)>>h_Nmax_Energy_fluo",drawcut,"profs");
else if(option == "Nmax_Energy_cer") fStree->Draw("log10(Nmax_cer):log10(Energy)>>h_Nmax_Energy_cer",drawcut,"profs");
else if(option == "Nph2_Energy_fluo") fStree->Draw("log10(Nph2_fluo*cos(Theta*TMath::Pi()/180)):log10(Energy)>>h_Nph2_Energy_fluo",drawcut,"profs");
else if(option == "Nph2_Energy_cer") fStree->Draw("log10(Nph2_cer):log10(Energy)>>h_Nph2_Energy_cer",drawcut,"profs");
else if(option == "Nmax_Energy_fluo_linear") fStree->Draw("Nmax_fluo*cos(Theta*TMath::Pi()/180):log10(Energy)>>h_Nmax_Energy_fluo_linear",drawcut,"goffprofs");
else if(option == "Nmax_Energy_cer_linear") fStree->Draw("Nmax_cer:log10(Energy)>>h_Nmax_Energy_cer_linear",drawcut,"goffprofs");
else if(option == "Nph2_Energy_fluo_linear") fStree->Draw("Nph2_fluo*cos(Theta*TMath::Pi()/180):log10(Energy)>>h_Nph2_Energy_fluo_linear",drawcut,"goffprofs");
else if(option == "Nph2_Energy_cer_linear") fStree->Draw("Nph2_cer:log10(Energy)>>h_Nph2_Energy_cer_linear",drawcut,"goffprofs");
// photons in atmosphere correlations VS theta
else if(option == "Nph2_Theta_fluo") fStree->Draw("Nph2_fluo/Energy:Theta>>h_Nph2_Theta_fluo",drawcut,"");
else if(option == "Nph2_Theta_cer") fStree->Draw("Nph2_cer/Energy:Theta>>h_Nph2_Theta_cer",drawcut,"");
else if(option == "Nmax_Theta_fluo") fStree->Draw("Nmax_fluo/Energy:Theta>>h_Nmax_Theta_fluo",drawcut,"");
else if(option == "Nmax_Nph2_fluo_Hmax") fStree->Draw("Nmax_fluo*1000/Nph2_fluo:Hmax_fluo>>h_Nmax_Nph2_fluo_Hmax",drawcut,"");
else if(option == "Nmax_Theta_cer") fStree->Draw("Nmax_cer/Energy:Theta>>h_Nmax_Theta_cer",drawcut,"");
else if(option == "Nph2_Hmax_cer") fStree->Draw("Nph2_cer/Energy:Hmax_show>>h_Nph2_Hmax_cer",drawcut,"");
else if(option == "Nph2_Hmax_fluo") fStree->Draw("Nph2_fluo/Energy:Hmax_show>>h_Nph2_Hmax_fluo",drawcut,"");
// photons on pupil correlations VS Energy
else if(option == "Nph_Energy_CB_P") fStree->Draw("log10(Nph_cer_air_pupil*cos(Theta*TMath::Pi()/180)):log10(Energy)>>h_Nph_Energy_CB_P",drawcut,"profs");
else if(option == "Nph_Energy_CR_P") fStree->Draw("log10(Nph_cer_refl_pupil):log10(Energy)>>h_Nph_Energy_CR_P",drawcut,"profs");
else if(option == "Nph2_Energy_Fluo_P") fStree->Draw("log10(Nph2_fluo_pupil*cos(Theta*TMath::Pi()/180)):log10(Energy)>>h_Nph2_Energy_Fluo_P",drawcut,"profs");
else if(option == "Nmax_Energy_Fluo_P") fStree->Draw("log10(Nmax_fluo_pupil*cos(Theta*TMath::Pi()/180)):log10(Energy)>>h_Nmax_Energy_Fluo_P",drawcut,"profs");
else if(option == "Nph_Energy_CB_P_linear") fStree->Draw("Nph_cer_air_pupil*cos(Theta*TMath::Pi()/180):log10(Energy)>>h_Nph_Energy_CB_P_linear",drawcut,"goffprofs");
else if(option == "Nph_Energy_CR_P_linear") fStree->Draw("Nph_cer_refl_pupil:log10(Energy)>>h_Nph_Energy_CR_P_linear",drawcut,"goffprofs");
else if(option == "Nph2_Energy_Fluo_P_linear") fStree->Draw("Nph2_fluo_pupil*cos(Theta*TMath::Pi()/180):log10(Energy)>>h_Nph2_Energy_Fluo_P_linear",drawcut,"goffprofs");
else if(option == "Nmax_Energy_Fluo_P_linear") fStree->Draw("Nmax_fluo_pupil*cos(Theta*TMath::Pi()/180):log10(Energy)>>h_Nmax_Energy_Fluo_P_linear",drawcut,"goffprofs");
// photons on pupil correlations VS theta
//else if(option == "Nph_Theta_CB_P") fStree->Draw("Nph_cer_back_pupil/Energy:Theta>>h_Nph_Theta_CB_P",drawcut,"");
else if(option == "Nph_Theta_CB_P") fStree->Draw("Nph_cer_air_pupil/Energy:Theta>>h_Nph_Theta_CB_P",drawcut,"");
else if(option == "Nph_Theta_CR_P") fStree->Draw("Nph_cer_refl_pupil/Energy:Theta>>h_Nph_Theta_CR_P",drawcut,"");
else if(option == "Nph_Theta_Fluo2_P") fStree->Draw("Nph2_fluo_pupil/Energy:Theta>>h_Nph_Theta_Fluo2_P",drawcut,"");
else if(option == "Nmax_fluo_Theta_P") fStree->Draw("Nmax_fluo_pupil/Energy:Theta>>h_Nmax_fluo_Theta_P",drawcut,"");
else if(option == "Nph_Theta_Fluo_P_crash") fStree->Draw("Nph2_fluo_pupil/Nph_fluo_pupil:Theta>>h_Nph_Theta_Fluo_P_crash",drawcut,"");
//else if(option == "Nph_airscat_effect") fStree->Draw("(Nph_fluo_pupil+Nph_cer_back_pupil+Nph_cer_refl_pupil)/(Nph_fluo_pupil+Nph_cer_refl_pupil):Theta>>h_Nph_airscat_effect",drawcut,"");
else if(option == "Nph_airscat_effect") fStree->Draw("(Nph_fluo_pupil+Nph_cer_air_pupil+Nph_cer_refl_pupil)/(Nph_fluo_pupil+Nph_cer_refl_pupil):Theta>>h_Nph_airscat_effect",drawcut,"");
else if(option == "TimeWidth2_Theta") fStree->Draw("Width_fluo_pupil:Theta>>h_TimeWidth2_Theta",drawcut,"");
// Transmissions fluo (+ckov) : ratios (inomega / onpupil)
else if(option == "Trans_fluo_Theta") fStree->Draw("Nph_fluo_pupil/FluoInOmega:Theta>>h_Trans_fluo_Theta",drawcut,"");
else if(option == "Trans_max_mean") fStree->Draw("(Trans_fluo_atMax - (FluoInOmega/Nph_fluo_pupil))/(FluoInOmega/Nph_fluo_pupil):Theta>>h_Trans_max_mean",drawcut,"");
else if(option == "Trans_GTUmax_max") fStree->Draw("(Trans_fluo_GTUmax - Trans_fluo_atMax)/Trans_fluo_atMax>>h_Trans_GTUmax_max",drawcut,"");
//else if(option == "Trans_ckov_Theta") fStree->Draw("(Nph_cer_back_pupil+Nph_cer_refl_pupil)/CkovInOmega:Theta>>h_Trans_ckov_Theta",drawcut,"");
else if(option == "Trans_ckov_Theta") fStree->Draw("(Nph_cer_air_pupil+Nph_cer_refl_pupil)/CkovInOmega:Theta>>h_Trans_ckov_Theta",drawcut,"");
// MISCELLANEOUS//////////////////////////
// fluo yield properties
else if(option == "Yield_max_mean") fStree->Draw("(Yield_fluomax - Yield_f)/Yield_f:Theta>>h_Yield_max_mean",drawcut,"");
else if(option == "Yield_GTUmax_max") fStree->Draw("(Yield_fluoGTUmax - Yield_fluomax)/Yield_fluomax:Theta>>h_Yield_GTUmax_max",drawcut,"");
// XGTUmax shift due to scattered ckov
else if(option == "Xmax_shift") fStree->Draw("(X_GTUmax_tot - X_GTUmax)/X_GTUmax:Theta>>h_Xmax_shift",drawcut,"");
// correlations between width of time profile on pupil VS hmax
else if(option == "TimeWidth2_Hmax") fStree->Draw("Width_fluo_pupil:Hmax_show>>h_TimeWidth2_Hmax",drawcut,"");
else Printf("<ETreePainter::DrawHistos()> : %s OPTION not allowed", option.Data());
}
//_____________________________________________________________________________
TH1F* ETreePainter::BuildAltHisto( Option_t *opt ) {
//
// generate histos as function of altitude (Xaxis in km)
//
// if no data available
if ( !fAtmo || !fShower) {
Info("BuildAltHistos()","EAtmosphere (or EShower) object is NULL. Painter made zombie.");
MakeZombie();
return 0;
}
Int_t n = fAtmo->GetNumBunches();
// if photons data are empty
if ( n <= 0 ) {
Info("BuildAltHisto()","No Bunches in Atmosphere ( NumBunch = 0 ). Painter made zombie.");
MakeZombie();
return 0;
}
// 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++) {
if (fAtmo->GetBunch(i)->GetType() == 0 ) {
temp[Nbins] = Zv(fAtmo->GetBunch(i)->GetShowerPosi())/km;
temp[Nbins+1] = Zv(fAtmo->GetBunch(i)->GetShowerPosf())/km;
Nbins++;
}
}
// fill X-axis array for bunches altitude
Double_t AltBins[Nbins+1];
Int_t incrm = 0;
for (Int_t i=0; i<Nbins+1; i++) {
AltBins[Nbins - incrm] = temp[i];
incrm++;
}
TH1F* h = 0;
if(opt == "alt_Bunch_Nph_F") h = new TH1F("alt_Bunch_Nph_F", "Nph_Fluo vs altitude", Nbins, AltBins);
else if(opt == "alt_Bunch_Nph_C" ) h = new TH1F("alt_Bunch_Nph_C", "Nph_Cerenkov vs altitude", Nbins, AltBins);
else if(opt == "alt_Shower" ) h = new TH1F("alt_Shower", "Shower vs altitude", Nbins, AltBins);
// fill histograms from bunches
EShowerStep* shstep = 0;
Int_t bin(-1);
for (Int_t i=0; i<n; i++) {
if(i%2 == 0) shstep = fShower->GetStep(i/2);
Double_t alt = ( Zv(fAtmo->GetBunch(i)->GetShowerPosi())
+ Zv(fAtmo->GetBunch(i)->GetShowerPosf()) )/2/km;
Double_t Nph = fAtmo->GetBunch(i)->GetWeight();
Double_t Nelec = shstep->GetNumElectrons();
if (fAtmo->GetBunch(i)->GetType()==0) {
if(opt == "alt_Bunch_Nph_F") {
bin = h->FindBin(alt);
h->SetBinContent(bin,Nph/h->GetBinWidth(bin));
}
if(opt == "alt_Shower") {
bin = h->FindBin(alt);
h->SetBinContent(bin,Nelec);
}
}
if (fAtmo->GetBunch(i)->GetType()==1) {
if(opt =="alt_Bunch_Nph_C" ) {
bin = h->FindBin(alt);
h->SetBinContent(bin,Nph/h->GetBinWidth(bin));
}
}
}
return h;
}
//_____________________________________________________________________________
TH1F* ETreePainter::BuildLongitudinalHisto( Option_t* opt ) {
//
// generate histos along the track (Longitudinal Profiles)
// (Xaxis in km or in g/cm2)
//
// if no data available
if ( !fAtmo || !fShower) {
Info("BuildLongitudinalHistos()","EAtmosphere or EShower (needed for this method) object is NULL. Painter made zombie.");
MakeZombie();
return 0;
}
Int_t n = fAtmo->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("BuildLongitudinalHisto()","No Bunches in Atmosphere ( NumBunch = 0 ). Painter made zombie.");
MakeZombie();
return 0;
}
// determine bins for longitudinal development (both in km and g/cm2)
UInt_t Nbins = UInt_t(nsh);
Double_t LongBins[Nbins+1];
Double_t GramBins[Nbins+1];
TVector3 initpos = fShower->GetStep(0)->GetPosi();
EShowerStep* shstep = 0;
for (Int_t i=0; i<nsh; i++) {
shstep = fShower->GetStep(i);
LongBins[i] = (shstep->GetPosi() - initpos).Mag()/km;
LongBins[i+1] = (shstep->GetPosf() - initpos).Mag()/km;
GramBins[i] = shstep->GetXi()*cm2/g;
GramBins[i+1] = shstep->GetXf()*cm2/g;
}
TH1F* h = 0;
if( opt == "Fluo_Longit_prof_km") h = new TH1F("Fluo_Longit_prof_km","Fluo_Longitudinal_profile_km",Nbins,LongBins);
else if(opt == "Ckov_Longit_prof_km") h = new TH1F("Ckov_Longit_prof_km","Ckov_Longit_prof_km",Nbins,LongBins);
else if(opt == "Fluo_Longit_prof_gram") h = new TH1F("Fluo_Longit_prof_gram","Fluo_Longit_prof_gram",Nbins,GramBins);
else if(opt == "Ckov_Longit_prof_gram" ) h = new TH1F("Ckov_Longit_prof_gram","Ckov_Longit_prof_gram",Nbins,GramBins);
else if(opt == "Shower_Longit_prof_gram" ) h = new TH1F("Shower_Longit_prof_gram","Shower_Longit_prof_gram",Nbins,GramBins);
// fill histograms from fluo bunches
for (Int_t i=0; i<n; i++) {
if(i%2 == 0) 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 = fAtmo->GetBunch(i)->GetWeight();
Double_t nelec = shstep->GetNumElectrons();
if (fAtmo->GetBunch(i)->GetType()==0) {
if(opt == "Fluo_Longit_prof_km") h->Fill(dist,nph);
if(opt == "Fluo_Longit_prof_gram") h->Fill(grammage,nph);
if(opt == "Shower_Longit_prof_gram") h->Fill(grammage,nelec);
}
if (fAtmo->GetBunch(i)->GetType()==1) {
if(opt == "Ckov_Longit_prof_km") h->Fill(dist,nph);
if(opt =="Ckov_Longit_prof_gram" ) h->Fill(grammage,nph);
}
}
return h;
}
//_____________________________________________________________________________
TH1F* ETreePainter::BuildTofHisto( Option_t* opt ) {
//
// generate histos of photons on pupil as function of time since PRIMARY cosmic ray first interaction
// (Xaxis in microsecond -> binwidth in GTU !!)
//
// if no data available
if ( !fAtmo ) {
Info("BuildTofHistos()","EAtmosphere object is NULL. Painter made zombie.");
MakeZombie();
return 0;
}
Int_t ns = fAtmo->GetNumSingles();
// if photons data are empty
if ( ns <= 0 ) {
Info("BuildTofHisto()","No SinglePhoton generated in Atmosphere");
}
// determine bins for time on pupil
Double_t TimeMin = 3000.*microsecond;
Double_t TimeMax = 0;
for (Int_t i=0; i<ns; i++) {
if ( !(fAtmo->GetSingle(i)->IsAbsorbed()) ) {
if ( TimeMin > (fAtmo->GetSingle(i)->GetTof()+fAtmo->GetSingle(i)->GetDate()) )
TimeMin = fAtmo->GetSingle(i)->GetTof()+fAtmo->GetSingle(i)->GetDate();
if ( TimeMax < (fAtmo->GetSingle(i)->GetTof()+fAtmo->GetSingle(i)->GetDate()) )
TimeMax = fAtmo->GetSingle(i)->GetTof()+fAtmo->GetSingle(i)->GetDate();
}
}
// build histogram vs time on pupil
Float_t min = (Int_t)(TimeMin/microsecond)-1;
Int_t nbins = (Int_t)((TimeMax-TimeMin)/(fGTU)) + 1;
TH1F* h = 0;
if(opt == "Single_Nph_F_P_t") h = new TH1F("Single_Nph_F_P_t","Nph_Fluo_on_Pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
else if(opt == "Single_Nph_Cdir_P_t") h = new TH1F("Single_Nph_CR_P_t","Nph_Cerenkov direct on pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
else if(opt == "Single_Nph_CB_P_t") h = new TH1F("Single_Nph_CB_P_t","Nph_Cerenkov back on pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
else if(opt == "Single_Nph_CR_P_t") h = new TH1F("Single_Nph_CR_P_t","Nph_Cerenkov refl on pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
else if(opt == "Single_Nph_CB_cloud_P_t") h = new TH1F("Single_Nph_CB_cloud_P_t","Nph_Cerenkov clouds scat on pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
else if(opt == "Single_Nph_F-CB_P_t") h = new TH1F("Single_Nph_F-CB_P_t","Nph_fluo_ckovBack on pupil vs time", nbins, min, min+nbins*(fGTU/microsecond));
// fill histograms from singles
for (Int_t i=0; i<ns; i++) {
Double_t tof = (fAtmo->GetSingle(i)->GetTof() + fAtmo->GetSingle(i)->GetDate()) / microsecond;
if ( !(fAtmo->GetSingle(i)->IsAbsorbed()) ) {
if ( fAtmo->GetSingle(i)->GetType() == 0 ) {
if(opt == "Single_Nph_F_P_t") h->Fill(tof,1);
if(opt == "Single_Nph_F-CB_P_t") h->Fill(tof,1);
}
if ( fAtmo->GetSingle(i)->GetType() == 1 ) {
if ( fAtmo->GetSingle(i)->GetHistory()==0) {
if(opt == "Single_Nph_Cdir_P_t") h->Fill(tof,1);
}
if ( fAtmo->GetSingle(i)->GetHistory()==2) {
if(opt == "Single_Nph_CB_P_t") h->Fill(tof,1);
if(opt == "Single_Nph_F-CB_P_t") h->Fill(tof,1);
}
if ( fAtmo->GetSingle(i)->GetHistory()==3) {
if(opt == "Single_Nph_CB_cloud_P_t") h->Fill(tof,1);
}
if ( fAtmo->GetSingle(i)->GetHistory()==1 )
if(opt == "Single_Nph_CR_P_t") h->Fill(tof,1);
}
}
}
return h;
}
//_____________________________________________________________________________
TH1F* ETreePainter::BuildWlHisto( Option_t* opt ) {
//
// generate wavelength spectra of photons on pupil BEFORE and AFTER last transmission
// generate also wavelength spectra of photons in atmosphere at creation
// (Xaxis in nm)
//
// if no data available
if ( !fAtmo ) {
Info("BuildWlHistos()","EAtmosphere object is NULL. Painter made zombie.");
MakeZombie();
return 0;
}
Int_t n = fAtmo->GetNumBunches();
Int_t ns = fAtmo->GetNumSingles();
// if photons data are empty
if ( n <= 0 ) {
Info("BuildWlHisto()","No Bunches in Atmosphere ( NumBunch = 0 ). Painter made zombie.");
MakeZombie();
return 0;
}
TH1F* h = 0;
// bunches at creation spectra
if(opt == "Bunch_Wl_F") h = new TH1F("Bunch_Wl_F", "Nph_F vs Lambda",150,300,450);
else if(opt == "Bunch_Wl_C") h = new TH1F("Bunch_Wl_C", "Nph_C vs Lambda",150,300,450);
else if(opt == "Bunch_Wl_tot") h = new TH1F("Bunch_Wl_tot", "Nph_tot vs Lambda",150,300,450);
// photons on pupil spectra
else if(opt == "Single_Wl_F_P") h = new TH1F("Single_Wl_F_P", "Nph_Fluo_on_pupil vs Lambda",150,300,450);
else if(opt == "Single_Wl_CB_P") h = new TH1F("Single_Wl_CB_P", "Nph_Cerenkov_back_on_pupil vs Lambda",150,300,450);
else if(opt == "Single_Wl_CR_P") h = new TH1F("Single_Wl_CR_P","Nph_Cerenkov refl on pupil vs lambda",150,300,450);
else if(opt == "Single_Wl_tot_P") h = new TH1F("Single_Wl_tot_P","Nph_tot on pupil vs lambda",150,300,450);
//TH1F *th16_2 = GetWlHisto("Single_Wl_ckov_CloudBackscat_P", "Nph_Cerenkov_Cloudback_on_pupil vs Lambda");
//TH1F *th16_tot = GetWlHisto("Single_Wl_CB_tot_P", "Nph_Cerenkov_Total_back_on_pupil vs Lambda"); // not used so far
// fill histograms from bunches
for (Int_t i=0; i<n; i++) {
Double_t Nph = fAtmo->GetBunch(i)->GetWeight();
Double_t NumWl = fAtmo->GetBunch(i)->GetNumWavelengths();
const Float_t* lambda = 0;
const Float_t* wlweight = 0;
wlweight = fAtmo->GetBunch(i)->GetTable("weight");
lambda = fAtmo->GetBunch(i)->GetTable("lambda");
if(opt == "Bunch_Wl_tot") {
for(Int_t j=0; j<NumWl; j++) h->Fill( lambda[j]/nm , wlweight[j]*Nph );
}
if (fAtmo->GetBunch(i)->GetType()==0 && opt == "Bunch_Wl_F") {
// fluorescence bunches
for(Int_t j=0; j<NumWl; j++) h->Fill( lambda[j]/nm , wlweight[j]*Nph );
}
if (fAtmo->GetBunch(i)->GetType()==1 && opt == "Bunch_Wl_C") {
// Cerenkov bunches
for(Int_t j=0; j<NumWl; j++) h->Fill( lambda[j]/nm , wlweight[j]*Nph );
}
}
// fill histograms from singles
for (Int_t i=0; i<ns; i++) {
Double_t wl = fAtmo->GetSingle(i)->GetWl()/nm;
if (!(fAtmo->GetSingle(i)->IsAbsorbed())) {
// all photons
if(opt == "Single_Wl_tot_P") h->Fill(wl,1);
// for fluorescence photons
if ( fAtmo->GetSingle(i)->GetType() == 0 && opt == "Single_Wl_F_P") {
h->Fill(wl,1);
}
// for Cerenkov photons
if ( fAtmo->GetSingle(i)->GetType() == 1 ) {
// for backscattered ones
if ( fAtmo->GetSingle(i)->GetHistory()==2 && opt == "Single_Wl_CB_P") {
h->Fill(wl,1);
}
// for reflected ones
if ( fAtmo->GetSingle(i)->GetHistory()==1 && opt == "Single_Wl_CR_P") {
h->Fill(wl,1);
}
}
}
}
return h;
}