Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

SlastLightToEuso - source file

// ESAF : Euso Simulation and Analysis Framework
// $Id: SlastLightToEuso.cc,v 1.7 2005/10/29 17:41:50 naumov Exp $
//  created Oct, 28 2002

#include "SlastLightToEuso.hh"
#include "EEventTruthAdder.hh"
#include "Config.hh"
#include "ShowerTrack.hh"
#include "EShower.hh"
#include "EShowerFiller.hh"
#include "Atmosphere.hh"
#include "EConst.hh"
#include "DetectorGeometry.hh"
#include <TVector3.h>
using namespace sou;

ClassImp(SlastLightToEuso)

// interface to fortran functions
extern "C" int slastinit_();
extern "C" int slastsimulation_();
extern "C" int slastdetectorgeometry_(Float_t*,Float_t*,Float_t*,Float_t*);
extern "C" {
    struct slast2esaf_header {
        Float_t Energy;
        Float_t Theta;
        Float_t Phi;
        Float_t X1;
        Float_t Rint[3];  
        Float_t AgeEarth;
        Float_t Rimpact[3];
        Float_t Rmax[3];
        Float_t Xmax;
        Float_t ShowerHitGround;
    };
    struct slast2esaf_shower {
        int sn;
        Float_t ShowerXi[10000];
        Float_t ShowerXf[10000];
        Float_t ShowerRxi[10000];
        Float_t ShowerRyi[10000];
        Float_t ShowerRzi[10000];
        Float_t ShowerRxf[10000];
        Float_t ShowerRyf[10000];
        Float_t ShowerRzf[10000];
        Float_t ShowerTimei[10000];
        Float_t ShowerTimef[10000];
        Float_t ShowerAgei[10000];
        Float_t ShowerAgef[10000];
        Float_t ShowerNe[10000];
    }; 
    struct slast2esaf_fluoresence {
        int fn;
        Float_t fxshower[500000];
        Float_t fyshower[500000];
        Float_t fzshower[500000];
        Float_t fxfocal[500000];
        Float_t fyfocal[500000];
        Float_t fzfocal[500000];
        Float_t ftheta1[500000];
        Float_t fphi1[500000];
        Float_t flambda[500000];
        Float_t ftime[500000];
    };
    struct slast2esaf_cherenkov {
        int cn;
        Float_t cxshower[500000];
        Float_t cyshower[500000];
        Float_t czshower[500000];
        Float_t cxfocal[500000];
        Float_t cyfocal[500000];
        Float_t czfocal[500000];
        Float_t ctheta1[500000];
        Float_t cphi1[500000];
        Float_t clambda[500000];
        Float_t ctime[500000];
    };
    extern struct slast2esaf_header      esafheader_;
    extern struct slast2esaf_shower      esafshower_;
    extern struct slast2esaf_fluoresence esaffluor_;
    extern struct slast2esaf_cherenkov   esafcherenkov_;
}

// constructor
//_________________________________________________________________________________________ 
 SlastLightToEuso::SlastLightToEuso() : LightToEuso("SLAST"),fPhotons(0),fTruth(0) {
    Init();
}

//_________________________________________________________________________________________ 
 void SlastLightToEuso::Init() {
    // Enabling the atmosphere

    Msg(EsafMsg::Info) << "Atmosphere enabled " << Atmosphere::Get()->GetType() << MsgDispatch;
    PrepareDataCards();
    slastinit_();
    fPhotons = new ListPhotonsOnPupil(new vector<ParentPhoton*>);

}


// destructor
//_________________________________________________________________________________________ 
 SlastLightToEuso::~SlastLightToEuso() {
    if ( fTruth )
        delete fTruth;
}

// returns the appropriate PhotonsOnPupil object for one event
//_________________________________________________________________________________________ 
 PhotonsOnPupil *SlastLightToEuso::Get(const DetectorGeometry* dg) {

    // Slast banner
    Msg(EsafMsg::Info) << "SLAST.f77 called" << MsgDispatch;
    // call SLAST fortran code
    Float_t ISSX = (dg->GetPos()).X()/km; 
    Float_t ISSY = (dg->GetPos()).Y()/km; 
    Float_t ISSZ = (dg->GetPos()).Z()/km; 
    Float_t radius = dg->GetRadius()*mm/m;
    slastdetectorgeometry_(&radius,&ISSX,&ISSY,&ISSZ);
    slastsimulation_();
    Msg(EsafMsg::Info) << "SLAST.f77 completed succesfully" << MsgDispatch;
    // convert data to ListPhotonsOnPupil and return it
    fPhotons->Clear();
    // loop on photons coming from slast
    double x[3];
    double y[3];
    double th=0.,ph=0.,w=0.,t=0.;
    int fn = esaffluor_.fn;
    if( esaffluor_.fn == 0 && esafcherenkov_.cn == 0) return fPhotons;
    if (fn > 500000) fn = 500000;
    // Filling fluorescence information from SLAST
    for(int i=0; i<fn; i++) 
    {
        x[0] = esaffluor_.fxshower[i]*m; // in m
        x[1] = esaffluor_.fyshower[i]*m; // in m
        x[2] = esaffluor_.fzshower[i]*m; // in m
        y[0] = esaffluor_.fxfocal[i]*m;
        y[1] = esaffluor_.fyfocal[i]*m;
        y[2] = esaffluor_.fzfocal[i]*m;
        th   = esaffluor_.ftheta1[i]*deg;
        ph   = esaffluor_.fphi1[i]*deg;
        w    = esaffluor_.flambda[i]*nm ;
        t    = esaffluor_.ftime[i]*microsecond;
        ParentPhoton *p = new ParentPhoton(i,x,y,th,ph,w,t,(ParentPhotonType)1);
        fPhotons->Add( p );
    }
    int cn = esafcherenkov_.cn;
    if (cn > 500000) cn = 500000;
    // Filling cherenkov information from SLAST
    for(int i=0; i<cn; i++) 
    {
        x[0] = esafcherenkov_.cxshower[i]*m;
        x[1] = esafcherenkov_.cyshower[i]*m;
        x[2] = esafcherenkov_.czshower[i]*m;
        y[0] = esafcherenkov_.cxfocal[i]*m;
        y[1] = esafcherenkov_.cyfocal[i]*m;
        y[2] = esafcherenkov_.czfocal[i]*m;
        th   = esafcherenkov_.ctheta1[i]*deg;
        ph   = esafcherenkov_.cphi1[i]*deg;
        w    = esafcherenkov_.clambda[i]*nm ;
        t    = esafcherenkov_.ctime[i]*microsecond;
        ParentPhoton *p = new ParentPhoton(i+fn,x,y,th,ph,w,t,(ParentPhotonType)2);
        fPhotons->Add( p );
    }

    MsgForm(EsafMsg::Info,"Fluorescent %d Cherenkov %d Total %d photons generated",fn,cn,fn+cn);
    // set smontecarlo truth for the shower (clearing previous event first)
    SafeDelete(fTruth);

    EEvent* event = EEvent::GetCurrent();

    if ( event ) {
        EEventTruthAdder a( GetTruth());
        event->Fill( a );
        // Fill the Shower Track Object

        ShowerTrack *track = new ShowerTrack();
        ShowerStep s;
        int sn = esafshower_.sn;
        if (sn > 10000) sn = 10000;
        for(int i=0;i<sn;i++) {
            s.fXi = esafshower_.ShowerXi[i]*g/cm2;
            s.fXf = esafshower_.ShowerXf[i]*g/cm2;
            s.fXYZi.SetXYZ(esafshower_.ShowerRxi[i]*m,esafshower_.ShowerRyi[i]*m,esafshower_.ShowerRzi[i]*m);
            s.fXYZf.SetXYZ(esafshower_.ShowerRxf[i]*m,esafshower_.ShowerRyf[i]*m,esafshower_.ShowerRzf[i]*m);
            s.fTimei = esafshower_.ShowerTimei[i]*microsecond;
            s.fTimef = esafshower_.ShowerTimef[i]*microsecond;
            s.fAgei  = esafshower_.ShowerAgei[i];
            s.fAgef  = esafshower_.ShowerAgef[i];
            s.fNelectrons = esafshower_.ShowerNe[i];
            track->Add(s);
        }
        track->fEnergy    = esafheader_.Energy;
        track->fTheta     = esafheader_.Theta;
        track->fPhi       = esafheader_.Phi;
        track->fX1        = esafheader_.X1*g/cm2;
        track->fEthrEl    = 0;
        track->fDirVers.SetXYZ(TMath::Sin(esafheader_.Theta)*TMath::Cos(esafheader_.Phi),
                TMath::Sin(esafheader_.Theta)*TMath::Sin(esafheader_.Phi),
                TMath::Cos(esafheader_.Theta));
        track->fInitPos.SetXYZ(esafheader_.Rint[0]*m,esafheader_.Rint[1]*m,esafheader_.Rint[2]*m);
        track->fHitGround = (esafheader_.ShowerHitGround == 1 ? kTRUE: kFALSE);

        EShowerFiller f( track );
        event->Fill( f );
        delete track;
    }

    system("rm -rf TAPE*");
    return fPhotons;
}

// re-configure SLAST
//_________________________________________________________________________________________ 
 void SlastLightToEuso::Configure() {
}

// returns Montecarlo shower truth informations for this event
//_________________________________________________________________________________________ 
 MCTruth* SlastLightToEuso::GetTruth() {
    if (!fTruth) {
        fTruth = new MCTruth();
        fTruth->SetLightToEuso( this );
        fTruth->SetEnergy( esafheader_.Energy*eV);
        fTruth->SetThetaPhi(esafheader_.Theta,esafheader_.Phi);
        EVector v;
        v[0] = esafheader_.Rint[0]*m;
        v[1] = esafheader_.Rint[1]*m;
        v[2] = esafheader_.Rint[2]*m;
        fTruth->SetFirstInt(v,esafheader_.X1*g/cm2);
        v[0] = esafheader_.Rimpact[0]*m;
        v[1] = esafheader_.Rimpact[1]*m;
        v[2] = esafheader_.Rimpact[2]*m;
        fTruth->SetEarthImpact(v,esafheader_.AgeEarth);
        v[0] = esafheader_.Rmax[0]*m;
        v[1] = esafheader_.Rmax[1]*m;
        v[2] = esafheader_.Rmax[2]*m;
        fTruth->SetShowerMax(v,esafheader_.Xmax*g/cm2);
    }
    return fTruth;
}

// returns simulated shower additional infos for this event
//_________________________________________________________________________________________ 
 PhysicsData* SlastLightToEuso::GetPhysics() {
    return NULL;
}

// interface to fortran code
//_________________________________________________________________________________________ 
 void SlastLightToEuso::PrepareDataCards() {
    //
    // Prepare data card to read by slast77 engine
    //
    ConfigFileParser *pConfig = Config::Get()->GetCF("General","Euso");
    Float_t  EarthRadius    = EConst::EarthRadius()/km;
    
    // reading common to slast77 and slast++ options from GeneratorLightToEuso file    
    pConfig = Config::Get()->GetCF("LightToEuso","GeneratorLightToEuso");
    Float_t  FOV                = pConfig->GetNum("GeneratorLightToEuso.FoV"); 
    Float_t  InteractionVectorX = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorX");
    Float_t  InteractionVectorY = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorY");
    Float_t  InteractionVectorZ = pConfig->GetNum("GeneratorLightToEuso.InteractionVectorZ");
    Float_t  ThetaRangeMin  = pConfig->GetNum("GeneratorLightToEuso.ThetaRangeMin");
    Float_t  ThetaRangeMax  = pConfig->GetNum("GeneratorLightToEuso.ThetaRangeMax");
    Float_t  PhiRangeMin    = pConfig->GetNum("GeneratorLightToEuso.PhiRangeMin");
    Float_t  PhiRangeMax    = pConfig->GetNum("GeneratorLightToEuso.PhiRangeMax");
    Float_t  EnergyRangeMin = pConfig->GetNum("GeneratorLightToEuso.EnergyRangeMin");
    Float_t  EnergyRangeMax = pConfig->GetNum("GeneratorLightToEuso.EnergyRangeMax");
    Float_t  EnergySlope    = pConfig->GetNum("GeneratorLightToEuso.EnergySlope");
    string fType            = pConfig->GetStr("GeneratorLightToEuso.UhecrType");
    Float_t UhecrType = EConst::AtomicMass(fType);
    
    // reading specific to slast77 options from its own config
    Float_t  WaveRangeMin   =  Conf()->GetNum("SlastLightToEuso.WaveRangeMin");
    Float_t  WaveRangeMax   =  Conf()->GetNum("SlastLightToEuso.WaveRangeMax");
    string DoCherenkov      = Conf()->GetStr("SlastLightToEuso.DoCherenkov");
    string DoFluorescence   = Conf()->GetStr("SlastLightToEuso.DoFluorescence");
    string AtmType          = Conf()->GetStr("SlastLightToEuso.AtmosphericType");
    Float_t  AtmTemperature = Conf()->GetNum("SlastLightToEuso.AtmTemperature");
    Float_t  Albedo         = Conf()->GetNum("SlastLightToEuso.Albedo");
    Float_t  GTU            = Conf()->GetNum("SlastLightToEuso.GTU");
    string AtmCurvature     = Conf()->GetStr("SlastLightToEuso.AtmCurvature");
    string EnergyDistributionParametrization = Conf()->GetStr("SlastLightToEuso.EnergyDistributionParametrization");
    string ShowerParametrization = Conf()->GetStr("SlastLightToEuso.ShowerParametrization");
    
    FILE *datacard = fopen("auxilar/global.datacard","w");
    int doch=0,dofl=0,atm_type=2,curvature=2,energydistribution=1,showerparametrization=1;
    fprintf(datacard,"LISTn");
    fprintf(datacard,"%s%fn","ERAD ",EarthRadius);
    fprintf(datacard,"%s%fn","FOV  ",FOV);
    fprintf(datacard,"%s%f %fn","WAVE ",WaveRangeMin,WaveRangeMax);
    if(DoCherenkov == "yes")    doch=1;
    if(DoFluorescence == "yes") dofl=1;
    fprintf(datacard,"%s%dn","DOCH ",doch);
    fprintf(datacard,"%s%dn","DOFL ",dofl);
    if(AtmType == "Isothermic") atm_type=1;
    if(AtmType == "USStandard") atm_type=2;
    fprintf(datacard,"%s%dn","ATMO ",atm_type);
    fprintf(datacard,"%s%fn","TEMP ",AtmTemperature);
    fprintf(datacard,"%s%fn","ALBE ",Albedo);
    fprintf(datacard,"%s%fn","GTU  ",GTU);
    if(AtmCurvature == "Curved") curvature = 1;
    if(AtmCurvature == "Planar") curvature = 2;
    fprintf(datacard,"%s%dn","CURV ",curvature);
    fprintf(datacard,"%s%f %f %fn","RINT ",InteractionVectorX,InteractionVectorY,InteractionVectorZ);
    fprintf(datacard,"%s%f %fn","THET ",ThetaRangeMin,ThetaRangeMax);
    fprintf(datacard,"%s%f %fn","PHI  ",PhiRangeMin,PhiRangeMax);
    fprintf(datacard,"%s%e %en","EINT ",EnergyRangeMin,EnergyRangeMax);
    fprintf(datacard,"%s%fn","ERAN ",EnergySlope);
    fprintf(datacard,"%s%fn","TYPE ",UhecrType);
    if(EnergyDistributionParametrization == "Hillas") energydistribution=1;
    fprintf(datacard,"%s%dn","DIST ",energydistribution);
    if(ShowerParametrization == "GIL") showerparametrization=1;
    fprintf(datacard,"%s%dn","SHOW ",showerparametrization);
    fclose(datacard);
}
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