Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

TestLightToEuso - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: TestLightToEuso.cc,v 1.8 2005/11/07 23:24:05 thea Exp $
// M. Pallavicini created Jul,  2 2002

#include "TestLightToEuso.hh"
#include "EsafRandom.hh"
#include "ListPhotonsOnPupil.hh"
#include "utils.hh"
#include "MCTruth.hh"
#include "EEventTruthAdder.hh"
#include "EEvent.hh"
#include "EConst.hh"
#include "DetectorGeometry.hh"
#include "EsafSpectrum.hh"

#include "TMath.h"

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

ClassImp(TestLightToEuso)

//______________________________________________________________________________
//
// This class brings test fPhotons to the Euso Pupil
// The type of events supported are the following:
//     SPOT         light from a single point at the same time
//     EXTSPOT      light from a single point at the same time
//     RANDSPOT     single point randomly chosen
//     CIRCLE       light making a circle in focal surface (theta constant,
//                  phi 0-360 degrees)
//     EXTCIRCLE    light making a thick circle in focal surface (theta
//                  constant, phi 0-360 degrees)
//     RADIUS       light making a radius (phi constant, theta 0-30 degrees)
//     TRACK        light making a track in focal surface with peaked pattern
//     SHOWERTRACK  track in the atmosphere
//     
// The parameters needed are the following:
//     TestLightToEuso.Type       see above
//     TestLightToEuso.Duration   event duration in micro-seconds
//     TestLightToEuso.Photons    number of fPhotons generated per event
//     TestLightToEuso.Theta1     constant if CIRCLE, SPOT or EXTSPOT, 
//               				  starting point if EXTCIRCLE, RADIUS or TRACK
//     TestLightToEuso.Theta2     ignored if CIRCLE, SPOT or EXTSPOT,
//     							  end point if RADIU, EXTCIRCLE or TRACK
//     TestLightToEuso.Phi1       starting point if CIRCLE, EXTCIRCLE or TRACK,
//     							  constant if RADIUS, SPOT or EXTSPOT
//     TestLightToEuso.Phi2       end point if CIRCLE, EXTCIRCLE or TRACK,
//                                ignored if RADIUS
//     				  or SPOT, radius of the beam in EXTSPOT
//

//______________________________________________________________________________
 TestLightToEuso::TestLightToEuso() : LightToEuso("TEST"), fTruth(NULL) {
    // 
    // Constructor
    // 

    fPhOnPupil = new ListPhotonsOnPupil();//(&fPhotons);
}

//______________________________________________________________________________
 TestLightToEuso::~TestLightToEuso() {
    // 
    // Destructor
    //
    
    Reset();
    delete fPhOnPupil;
    if ( fTruth ) delete fTruth;
}

//______________________________________________________________________________
 void TestLightToEuso::Reset() {
    // 
    // Reset internal list of fPhotons
    // 

    fPhOnPupil->Clear();
}

//______________________________________________________________________________
 PhotonsOnPupil *TestLightToEuso::Get(const DetectorGeometry* dg) {
    // 
    // Returns the list of fPhotons for the pupil.
    // The type of photon distribution is determined 
    // by TestLightToEuso.Type
    //

    Reset();
    string type = Conf()->GetStr("TestLightToEuso.Type");
    Double_t th1 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Theta1");
    Double_t th2 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Theta2");
    Double_t ph1 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Phi1");
    Double_t ph2 = TMath::DegToRad()*Conf()->GetNum("TestLightToEuso.Phi2");
    Int_t nPhotons = (Int_t)Conf()->GetNum("TestLightToEuso.Photons");
    Double_t tm = Conf()->GetNum("TestLightToEuso.Duration");
    TRandom* rndm = EsafRandom::Get();
     
    if ( type == "SPOT" ) {              // spot in fixed position
        return MakeSpot(th1,ph1,nPhotons,tm);
    }
    else if ( type == "RANDSPOT" ) {     // spot in random position within Phi1,2 Theta1,2
	Double_t cmax = Cos(th2);
	Double_t cmin = Cos(th1);
	if ( th1 > th2 ) {
	    Msg(EsafMsg::Info) << "In TestLightToEuso Theta1 and Theta2 will be switchedn";
	    cmax = cmin;
	    cmin = Cos(th2);
	}
	Double_t diff = cmax - cmin;
	Double_t xc = cmin + rndm->Rndm()*diff;
	Double_t th = acos(xc);
	if ( ph2 > ph1 )
	    diff = ph2 - ph1;
	else {
	    Msg(EsafMsg::Info) << "In TestLightToEuso Phi1 and Phi2 will be switchedn";
            diff = ph1 - ph2;
	    Double_t temp = ph2;
	    ph2 = ph1;
	    ph1 = temp;
	}   
	Double_t ph = rndm->Rndm()*diff + ph1;
	return MakeSpot(th,ph,nPhotons,tm);
    }
    else if ( type == "EXTSPOT" ) {
        return MakeExtSpot(th1,ph1,TMath::RadToDeg()*ph2,nPhotons,tm);
    }
    else if ( type == "CIRCLE" ) {
        return MakeCircle(ph1,ph2,th1,nPhotons,tm);
    }
    else if ( type == "EXTCIRCLE" ) {
	    return MakeExtCircle(ph1,ph2,th1,th2,nPhotons,tm);
    }
    else if ( type == "DIFFUSE" ) {
	    return MakeDiffuse(th1,th2,TMath::RadToDeg()*ph1,nPhotons,tm);
    }
    else if ( type == "CFLUXDIFFUSE" ) {
	    return MakeCFluxDiffuse(th1,th2,TMath::RadToDeg()*ph1,nPhotons,tm);
    }
    else if ( type == "RADIUS" ) {
        return MakeRadius(ph1,th1,th2,nPhotons,tm);
    } 
    else if ( type == "EXTRADIUS" ) {
        return MakeExtRadius(ph1,th1,th2,nPhotons,tm);
    }
    else if ( type == "TRACK" ) {
        return MakeTrack(ph1,ph2,th1,th2,nPhotons,tm);
    }
    else if ( type == "SHOWERTRACK" ) {

        return MakeShowerTrack( dg );
    }
    throw runtime_error("Invalid parameter TestLightToEuso.Type = " + type );
    return (PhotonsOnPupil*)0;
}

//______________________________________________________________________________
 PhotonsOnPupil *TestLightToEuso::MakeSpot(const Double_t& th, const Double_t&  ph, Int_t nPhotons, 
		         const Double_t&  tm) {
    //
    // Produce light in a single spot
    //

    Msg(EsafMsg::Info) << "SPOT CALLED th="<<th<<" ph="<<ph<<" Photons="<<nPhotons<< MsgDispatch;
    Double_t wl = 357.*nm;
    Double_t x[3] = {0.,0.,0.};
    Double_t y[3] = {0.,0.,0.};
    Double_t t;
    TRandom* rndm = EsafRandom::Get();
    for(Int_t i=0; i<nPhotons; i++) {
       t = rndm->Rndm()*tm*microsecond;
       ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
       fPhOnPupil->Add( p );
    }    
    return fPhOnPupil;
}

//______________________________________________________________________________
 PhotonsOnPupil *TestLightToEuso::MakeExtSpot(const Double_t& theta, const Double_t&  phi, const Double_t radius,
	           	Int_t nPhotons, const Double_t&  tm) {
    //
    // Produce light in an extensive spot
    //

    Msg(EsafMsg::Info) << "EXTSPOT CALLED theta="<<theta<<" phi="<<phi<<" radius="<<radius<<" Photons="<<nPhotons<< MsgDispatch;
    Double_t wl = 357.*nm;
    Double_t x[3] = {0.,0.,0.};
    Double_t y[3] = {0.,0.,0.};
    Double_t t;
    TRandom* rndm = EsafRandom::Get();
    for(Int_t i=0; i<nPhotons; i++) {
       t = rndm->Rndm()*tm*microsecond;
       Double_t radius = Sqrt(rndm->Rndm())*rad;
       Double_t th2 = rndm->Rndm()*TMath::TwoPi();
       y[0] = radius*Cos(th2)*mm;
       y[1] = radius*Sin(th2)*mm;
       ParentPhoton *p = new ParentPhoton(i,x,y,theta,phi,wl,t,Fluorescence);
       fPhOnPupil->Add( p );
    }    
    return fPhOnPupil;
}

//______________________________________________________________________________
 PhotonsOnPupil *TestLightToEuso::MakeCircle(const Double_t& ph1, const Double_t&  ph2, const Double_t&  th, 
		           Int_t nPhotons, const Double_t&  tm) {
    // 
    // Produce light in a circle
    // 
    TRandom* rndm = EsafRandom::Get();
    Msg(EsafMsg::Info) << "CIRCLE CALLED with ph1=" << ph1 << " ph2=" << ph2 << " th=" << th << MsgDispatch;
    Double_t wl = 357.*nm;
    Double_t x[3] = {0.,0.,0.};
    Double_t y[3] = {0.,0.,0.};
    Double_t t;
    for(Int_t i=0; i<nPhotons; i++) {
       t = rndm->Rndm()*tm*microsecond;
       Double_t ph = ph1 + rndm->Rndm()*(ph2-ph1);
       ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
       fPhOnPupil->Add( p );
    }    
    return fPhOnPupil;
}

//______________________________________________________________________________
 PhotonsOnPupil *TestLightToEuso::MakeExtCircle(const Double_t& ph1, const Double_t&  ph2, const Double_t&  th1,
           const Double_t th2, Int_t nPhotons, const Double_t&  tm) {
    // 
    // Produce light in a extensive circle (fPhotons with th1 < th < th2)
    // 
    
    TRandom *rndm = EsafRandom::Get();
    Msg(EsafMsg::Info) << "EXTCIRCLE CALLED with ph1=" << ph1 << " ph2=" << ph2 << 
	    " th1=" << th1 << " th2=" << th2 << MsgDispatch;
    Double_t wl = 357.*nm;
    Double_t x[3] = {0.,0.,0.};
    Double_t y[3] = {0.,0.,0.};
    Double_t t,theta1,theta2;
    if ( th1 > th2 ) {
    theta1=th2;
    theta2=th1;
    } else {
    theta1=th1;
    theta2=th2;
    }
    for(Int_t i=0; i<nPhotons; i++) {
       t = rndm->Rndm()*tm*microsecond;
       Double_t th = theta1+(theta2-theta1)*Sqrt((theta2-theta1)/(Pi()/6)*rndm->Rndm());
       Double_t ph = ph1 + rndm->Rndm()*(ph2-ph1);
       ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
       fPhOnPupil->Add( p );
    }    
    return fPhOnPupil;

}

//______________________________________________________________________________
 PhotonsOnPupil *TestLightToEuso::MakeDiffuse(const Double_t& th1, const Double_t&  th2,
        const Double_t rad, Int_t nPhotons, const Double_t&  tm) {
    //
    //
    //
    
    TRandom *rndm = EsafRandom::Get();
    Msg(EsafMsg::Info) << "DIFFUSE CALLED with th1=" << th1 << " th2=" << th2 << 
        " rad=" << rad <<  MsgDispatch;
    Double_t wl = 357.*nm;
    Double_t x[3] = {0.,0.,0.};
    Double_t y[3] = {0.,0.,0.};
    Double_t t(0),theta1,theta2;
    if ( th1 > th2 ) {
        theta1=th2;
        theta2=th1;
    } else {
        theta1=th1;
        theta2=th2;
    }
    
    Double_t c1 = Cos(theta1);
    Double_t c2 = Cos(theta2);
    for(Int_t i=0; i<nPhotons; i++) {
//        Double_t th = theta1+(theta2-theta1)*Sqrt((theta2-theta1)/(Pi()/6)*rndm->Rndm());
        Double_t th = ACos(rndm->Uniform(c1,c2));
        Double_t ph = TwoPi()*rndm->Rndm();

        Double_t rho = Sqrt(rndm->Rndm())*rad;
        Double_t alpha = rndm->Rndm()*TwoPi();
        y[0] = rho*Cos(alpha)*mm;
        y[1] = rho*Sin(alpha)*mm;

        ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
        fPhOnPupil->Add( p );
    }    
    return fPhOnPupil;
}

//______________________________________________________________________________
 PhotonsOnPupil *TestLightToEuso::MakeCFluxDiffuse(const Double_t& th1, const Double_t&  th2,
        const Double_t rad, Int_t nPhotons, const Double_t&  tm) {
    TRandom *rndm = EsafRandom::Get();
    Msg(EsafMsg::Info) << "C FLUX-DIFFUSE CALLED with th1=" << th1 << " th2=" << th2 << 
        " rad=" << rad << MsgDispatch;
    Double_t wl = 357.*nm;
    Double_t x[3] = {0.,0.,0.};
    Double_t y[3] = {0.,0.,0.};
    Double_t t(0),theta1(0),theta2(0);
    if ( th1 > th2 ) {
        theta1=th2;
        theta2=th1;
    } else {
        theta1=th1;
        theta2=th2;
    }
    
    Double_t c1 = Cos(2*theta1);
    Double_t c2 = Cos(2*theta2);
    for(Int_t i(0); i<nPhotons; i++) {
        Double_t th = 0.5*ACos(rndm->Uniform(c1,c2));
        Double_t ph = TMath::TwoPi()*rndm->Rndm();

        Double_t rho = Sqrt(rndm->Rndm())*rad;
        Double_t alpha = rndm->Rndm()*TMath::TwoPi();
        y[0] = rho*Cos(alpha)*mm;
        y[1] = rho*Sin(alpha)*mm;

        ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
        fPhOnPupil->Add( p );
    }    
    return fPhOnPupil;
}

//______________________________________________________________________________
 PhotonsOnPupil *TestLightToEuso::MakeRadius(const Double_t& ph, const Double_t& th1, const Double_t&  th2, 
		           Int_t nPhotons, const Double_t&  tm) {
    //
    // Produce light along a focal surface radius.
    // 

    TRandom* rndm = EsafRandom::Get();
    Msg(EsafMsg::Info) << "RADIUS CALLED" << MsgDispatch;
    Double_t wl = 357.*nm;
    Double_t x[3] = {0.,0.,0.};
    Double_t y[3] = {0.,0.,0.};
    Double_t t;
    for(Int_t i=0; i<nPhotons; i++) {
       t = rndm->Rndm()*tm*microsecond;
       Double_t th = th1 + rndm->Rndm()*(th2-th1); // uniform in theta, not costheta
       ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
       fPhOnPupil->Add( p );
    }    
    return fPhOnPupil;
}

//______________________________________________________________________________
 PhotonsOnPupil *TestLightToEuso::MakeExtRadius(const Double_t& ph, const Double_t& th1, const Double_t&  th2, 
		           Int_t nPhotons, const Double_t&  tm) {
    //
    // Produce light along a focal surface radius.
    // 

    TRandom* rndm = EsafRandom::Get();
    Msg(EsafMsg::Info) << "EXTRADIUS CALLED" << MsgDispatch;
    Double_t wl = 357.*nm;
    Double_t x[3] = {0.,0.,0.};
    Double_t y[3] = {0.,0.,0.};
    Double_t t;
    for(Int_t i=0; i<nPhotons; i++) {
        t = rndm->Rndm()*tm*microsecond;
        Double_t th = th1 + rndm->Rndm()*(th2-th1); // uniform in theta, not costheta
        Double_t rho = Sqrt(rndm->Rndm())*rad;
        Double_t alpha = rndm->Rndm()*TMath::TwoPi();
        y[0] = rho*Cos(alpha)*mm;
        y[1] = rho*Sin(alpha)*mm;

        ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,wl,t,Fluorescence);
        fPhOnPupil->Add( p );
    }    
    return fPhOnPupil;
}


//______________________________________________________________________________
 PhotonsOnPupil *TestLightToEuso::MakeTrack(const Double_t& ph1, const Double_t&  ph2, const Double_t& th1, 
		          const Double_t&  th2, Int_t nPhotons, const Double_t&  tm) {
    // 
    // Produce light along a track
    //

    //TRandom* rndm = EsafRandom::Get();
    return 0;
}

//______________________________________________________________________________
 PhotonsOnPupil *TestLightToEuso::MakeShowerTrack(const DetectorGeometry* geo) {
    // produce light along a track
    
    Msg(EsafMsg::Info) << "SHOWERTRACK CALLED" << MsgDispatch;
    
    TRandom* rndm = EsafRandom::Get();
    
    Double_t iX,iY, iZ;
    Double_t tLength = Conf()->GetNum("TestLightToEuso.ShowerTrack.fTrackLength")*km;
    Double_t tTheta = Conf()->GetNum("TestLightToEuso.ShowerTrack.fTheta")*TMath::DegToRad();
    Double_t tPhi = Conf()->GetNum("TestLightToEuso.ShowerTrack.fPhi")*TMath::DegToRad();
    string sInitPos = Conf()->GetStr("TestLightToEuso.ShowerTrack.fPosInit");
    Int_t tProfileCode = Nint(Conf()->GetNum("TestLightToEuso.ShowerTrack.fProfileCode"));

    if ( sInitPos == "fixed") {
        // fixed int point
        iX = Conf()->GetNum("TestLightToEuso.ShowerTrack.fXInit")*km;
        iY = Conf()->GetNum("TestLightToEuso.ShowerTrack.fYInit")*km;
        iZ = Conf()->GetNum("TestLightToEuso.ShowerTrack.fZInit")*km;

    } else if (sInitPos == "random") {
        // random interacion point
        Double_t rInitMax = Conf()->GetNum("TestLightToEuso.ShowerTrack.fRInitMax")*km;
        Double_t radInit = TMath::Sqrt(rndm->Rndm())*rInitMax;
        Double_t phiInit = rndm->Rndm()*TMath::TwoPi();
        iX = radInit*TMath::Cos(phiInit);
        iY = radInit*TMath::Sin(phiInit); 
        // Z init between 0 and 50 km
        iZ = rndm->Rndm()*Conf()->GetNum("TestLightToEuso.ShowerTrack.fZInitMax")*km;

    } else if (sInitPos=="zfixed") {
        // Z coordinate of the interacion point is kept fixed
        Double_t rInitMax = Conf()->GetNum("TestLightToEuso.ShowerTrack.fRInitMax")*km;
        Double_t radInit = TMath::Sqrt(rndm->Rndm())*rInitMax;
        Double_t phiInit = rndm->Rndm()*TMath::TwoPi();
        Msg(EsafMsg::Info) << rInitMax << " " << phiInit << " " << radInit << MsgDispatch;
        iX = radInit*TMath::Cos(phiInit);
        iY = radInit*TMath::Sin(phiInit); 
        iZ = Conf()->GetNum("TestLightToEuso.ShowerTrack.fZInit")*km;

    } else {
        Msg(EsafMsg::Info) << "TestLightToEuso.cfg:"
             << " wrong value in TestLightToEuso.ShowerTrack.fPosInit:"
             << sInitPos << MsgDispatch;
        Msg(EsafMsg::Panic) << "ESAF Error: Wrong value in TestLightToEuso.cfg" << MsgDispatch;
    }
 
    Int_t nPhotons = (Int_t)Conf()->GetNum("TestLightToEuso.ShowerTrack.fPhotons");
    const TVector3& detPos = geo->GetPos();
    Double_t detRad = geo->GetRadius();
    
    TVector3 initPos(iX, iY, iZ);
    TVector3 diffPos;
    diffPos.SetMagThetaPhi(tLength,tTheta,tPhi);
    TVector3 endPos = diffPos + initPos;
    TVector3 dir = (endPos-initPos).Unit();
    TVector3 dL = (endPos-initPos)*(1./nPhotons);
    
    fTruth = new MCTruth;
    fTruth->SetFirstInt(initPos-detPos, 0);
    fTruth->SetThetaPhi(dir.Theta(), dir.Phi());
    fTruth->SetShowerMax(0.5*(endPos-initPos)+initPos-detPos, 0);

    if ( EEvent::GetCurrent() ) {
        EEventTruthAdder a(fTruth);
    
        EEvent::GetCurrent()->Fill( a );
    }

    EsafSpectrum *shSpec = new EsafSpectrum("nitro10km");

    Double_t wl;
    Double_t t;
    Double_t x[3] = {0.,0.,0.};
    Double_t y[3] = {0.,0.,0.};
    Double_t l;
    
    for(Int_t i(0); i<nPhotons; i++) {
        // get a random number, gaussian distributed between 0 and 1 sigma 
        switch ( tProfileCode) {
            case 1:
                for(l = rndm->Gaus(0.5,1/4.); l < 0 || l > 1;l = rndm->Gaus(0.5,1/4.) );
                break;
            case 0:
            default:
                l = rndm->Rndm();
        }

        TVector3 shPos = initPos+diffPos*l; 
        TVector3 inDir = -((shPos-detPos).Unit());
        Double_t radius = Sqrt(rndm->Rndm())*detRad;
        Double_t phi  = rndm->Rndm()*TMath::TwoPi();

        x[0] = (shPos-detPos)[0];
        x[1] = (shPos-detPos)[1];
        x[2] = (shPos-detPos)[2];
        
        y[0] = radius*Cos(phi)*mm;
        y[1] = radius*Sin(phi)*mm;
        y[2] = 0*mm;

        wl = shSpec->GetLambda();

        t = (l*tLength+(shPos-detPos).Mag())/Clight();
        ParentPhoton *p = new ParentPhoton(i,x,y,inDir.Theta(),inDir.Phi(),wl,t,Fluorescence);
        fPhOnPupil->Add( p );
    }
    
    
    return fPhOnPupil;
}
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