Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

BunchOfPhotons - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: BunchOfPhotons.cc,v 1.37 2005/10/21 13:31:30 moreggia Exp $
// Anne Stutz created Dec,  1 2003

#include "BunchOfPhotons.hh"
#include <stdexcept>
#include "EsafRandom.hh"
#include "EConst.hh"

ClassImp(BunchOfPhotons)

using namespace sou;
using namespace TMath;
using namespace EConst;

UInt_t BunchOfPhotons::fCounter = 0;

//____________________________________________________________________________________________________
 BunchOfPhotons::BunchOfPhotons(): fId(0), fParent(0) {
    //
    // ctor
    //
    Msg(EsafMsg::Panic) <<"Default ctor should not be used" <<MsgDispatch;
}

//____________________________________________________________________________________________________
 BunchOfPhotons::BunchOfPhotons(const BunchOfPhotons& b) : fParent(0){
    //
    // copy ctor
    //
    fId = b.GetId();
    fNbins = b.GetNbins();
    fType = b.GetType();
    fWeight = b.GetWeight();
    fDate = b.GetDate();
    fTof = b.GetTof();
    fPos = b.GetPos();
    fShowerPos = b.GetShowerPos();
    fDir = b.GetDir();
    fNextImpact = b.GetNextImpact();
    fFinalImpact = b.GetFinalImpact();
    b.GetWlSpectrum().Copy(fWlSpectrum);
    fFate = b.GetFate();
    fParent = new ParentBunch(*b.GetParent());
    
    //b.GetTransSpectrum().Copy(fTransSpectrum);
    //b.GetRaylSpectrum().Copy(fRaylSpectrum);
    //b.GetOzoneSpectrum().Copy(fOzoneSpectrum);
    //b.GetMieSpectrum().Copy(fMieSpectrum);
}

//____________________________________________________________________________________________________
 BunchOfPhotons::BunchOfPhotons(const Double_t weight,const Double_t yield,
    const EarthVector& showerposi, const EarthVector& showerposf,
    const Double_t datei, const Double_t datef, const EsafSpectrum& spectrum,
    const EarthVector& dir, PhotonType type) {
    //
    // ctor with creation of a Parent Bunch 
    //
    fCounter++;
    fId = fCounter;
    fType = type;
    fWeight = weight;
    fDate = (datei+datef)/2;
    fTof = 0;
    fShowerPos = ( showerposi+showerposf )*0.5;
    fPos = fShowerPos;
    fDir = dir.Unit();
    fNextImpact = fShowerPos;
    fFinalImpact = fShowerPos;
    spectrum.Copy(fWlSpectrum);
    fNbins = spectrum.GetAxis().GetNbins();
    fFate = 0;
    fParent = new ParentBunch(weight,yield,showerposi,showerposf,datei,datef,spectrum,fDir);
    
    //spectrum.Copy(fTransSpectrum);
    //spectrum.Copy(fRaylSpectrum);
    //spectrum.Copy(fOzoneSpectrum);
    //spectrum.Copy(fMieSpectrum);
    // following needed for convolution operations
    //fTransSpectrum.Clear();
    //fRaylSpectrum.Clear();
    //fOzoneSpectrum.Clear();
    //fMieSpectrum.Clear();
}

//____________________________________________________________________________________________________
 BunchOfPhotons::~BunchOfPhotons() {
    //
    // dtor
    //
    SafeDelete(fParent);
}

//____________________________________________________________________________________________________
 EarthVector BunchOfPhotons::RandomPosInShower() const {
    //
    // Returns a random position at creation, using lateral and longitudinal position distributions
    // If no distributions, returns bunch showerpos
    //
    
    EarthVector rtn(1,0,0);
    EarthVector longb =  fParent->GetShowerPosf()- fParent->GetShowerPosi();
    EarthVector axis = fParent->GetAxis();
    TRandom* rndm = EsafRandom::Get();
    rtn = fParent->GetShowerPosi() + rndm->Rndm()*longb;
    // if working with histos
    if (fParent->GetHistoLateral()) {
	EarthVector v = axis.Orthogonal();                   // vector orthogonal to direction
        v.Rotate(rndm->Rndm() * 2 * TMath::Pi(),axis);       // random rotation around direction
        v.SetMag(fParent->GetHistoLateral()->GetRandom()*m);  // lateral distribution in meter if NKG2 //TOFIX if NKG1
	rtn += v ;
	return rtn;
    }
    // if working with functions
    else if(fParent->GetLateralDist()) {
        EarthVector v = axis.Orthogonal();                   // vector orthogonal to direction
        v.Rotate(rndm->Rndm() * 2 * TMath::Pi(),axis);       // random rotation around direction
        v.SetMag(fParent->GetLateralDist()->GetRandom()*m);   // lateral distribution in meter if NKG2 //TOFIX if NKG1

        rtn += v ;
        return rtn;
    }

    // if no distributions
    return GetShowerPos();
}

//____________________________________________________________________________________________________
 EarthVector BunchOfPhotons::PosRandomAngCorrec(const EarthVector& showerpos,const EarthVector& currentpos) const {
    //
    // Returns a corrected position at any propagation stage, using : 
    //     - angular distributions at creation
    //     - the corrected SinglePhoton showerpos given in argument
    //     - the current SinglePhoton position given in argument
    // NB : all corrections (in position and time) must be done before calling the present method
    //
    
    // correction due to corrected showerpos
    EarthVector increment(currentpos - showerpos);
    EarthVector rtn(currentpos);
    
    // should not occur, but..
    if(increment.Mag() == 0) return rtn;
    
    // correction due to angular distributions
    if (GetParent()->GetHistoAngular()) {
        EarthVector v = increment.Orthogonal();
        Double_t theta = fParent->GetHistoAngular()->GetRandom();
       	TRandom* rndm = EsafRandom::Get();
        v.SetMag(increment.Mag() * TMath::Tan(theta));
        v.Rotate(rndm->Rndm() * 2 * TMath::Pi(),increment);
	rtn += v;
        return rtn;      
    }
    else if (GetParent()->GetAngularDist()) {
        EarthVector v = increment.Orthogonal();
        Double_t theta = fParent->GetAngularDist()->GetRandom();
       	TRandom* rndm = EsafRandom::Get();
        v.SetMag(increment.Mag() * TMath::Tan(theta));
        v.Rotate(rndm->Rndm() * 2 * TMath::Pi(),increment);
	rtn += v;
        return rtn;      
    } 
    
    // if no distibutions, return current position
    return rtn;
}

//____________________________________________________________________________________________________
 EarthVector BunchOfPhotons::RandomDirection() const {
    //
    // Returns a random direction at creation using the angular distributions
    //
    
    Double_t theta = 0;
    EarthVector rtn(fDir);
    TRandom* rndm = EsafRandom::Get();
    
    // isotropic if fluo
    if(fType == Fluo) {
        Double_t x,y,z;
	rndm->Sphere(x,y,z,1.);
	rtn.SetXYZ(x,y,z);
	return rtn;
    }
    
    // use angular distribution if cerenkov
    if(fDir.Mag() == 0) Msg(EsafMsg::Warning) <<"<RandomDirection> not-Fluo bunch with a NULL direction -> SHOULD NOT" << MsgDispatch;
    if(GetParent()->GetHistoAngular())     theta = fParent->GetHistoAngular()->GetRandom();
    else if(GetParent()->GetAngularDist()) theta = fParent->GetAngularDist()->GetRandom();
    // if no distribution, return bunch mean direction
    else return fDir;
    
    EarthVector v = fDir.Orthogonal();
    v.SetMag(TMath::Tan(theta));
    v.Rotate(rndm->Rndm() * 2 * TMath::Pi(),fDir);
    rtn += v;
    rtn = rtn.Unit();
    
    return rtn;
}

//____________________________________________________________________________________________________
 Double_t BunchOfPhotons::AngularDist_OverTwoPi(Double_t theta) const {
    //
    // returns 'angular distribution value' / 2pi, at theta
    // theta taken w.r.t bunch axis, and zero means same direction as bunch axis
    //
    
    Double_t rtn = 0;
    if (fParent->GetHistoAngular()) {
        TH1F histo(*fParent->GetHistoAngular());
	rtn = histo.GetBinContent(histo.FindBin(theta));
    }
    else if(fParent->GetAngularDist()) {
        TF12 distrib(*fParent->GetAngularDist());
	rtn = distrib.Eval(theta);
    }
    // if isotropic fluo bunch
    else if(fType == Fluo) rtn = 0.5;
    else rtn = 0;
    
    return rtn/TwoPi();
}

//____________________________________________________________________________________________________
 void BunchOfPhotons::AddToPosTof(const EarthVector& v) {
    //
    // Increment the position AND the tof
    //
    fTof += v.Mag() / Clight();
    fPos += v;
}

//____________________________________________________________________________________________________
 void BunchOfPhotons::Reset() {
    //
    // Reset method
    //
    fCounter = 0;
}

/*
//_____________________________________________________________________________
void BunchOfPhotons::ConvoluteTransSpectra(UInt_t nb, const vector<Double_t*>& coeff) {
    //
    // convolute the transmission spectra
    //
    if(nb != fNbins)
        Msg(EsafMsg::Warning) <<"<ConvoluteTransSpectra> Number of wavelengths doesn't match what is foreseen" << MsgDispatch;

    fTransSpectrum.Convolute(coeff[0]);
    fRaylSpectrum.Convolute(coeff[1]);
    fOzoneSpectrum.Convolute(coeff[2]);
    fMieSpectrum.Convolute(coeff[3]);
}
*/
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