Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

ETreePainter - source file

// $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;
}
About Us | EUSO Official Website | Web pages created by Roberto Pesce and Alessandro Thea - Last Update Wed Nov 16 16:57:39 2005 Wed Nov 16 16:29:22 2005