Universita' di GenovaINFN Sezione di Genova  
AIRWATCH / EUSO Genova

UnisimPhotonsOnPupil - source file

// $Id: UnisimPhotonsOnPupil.cc,v 1.16 2005/10/04 02:55:48 thea Exp $
// Author: M. Pallavicini created Jul,  2 2002

/*****************************************************************************
 * ESAF: Euso Simulation and Analysis Framework                              *
 *                                                                           *
 *  UnisimId: PhotonsOnPupil                                                 *
 *  Package: RadiativeTransfer                                               *
 *  Coordinator: Sylvain.Moreggia                                            *
 *                                                                           *
 *****************************************************************************/                                                                                     
#include "UnisimPhotonsOnPupil.hh"
#include "utils.hh"
#include "EVector.hh"
#include "Config.hh"
#include "EsafRefFrame.hh"
#include <zlib.h>

using namespace sou;

ClassImp(UnisimPhotonsOnPupil)

//_________________________________________________________________________________________ 
 UnisimPhotonsOnPupil::UnisimPhotonsOnPupil() : fPhotonOndisk(kFALSE), fTheTruth(0) {
    //
    // Constructor
    //

    fPhotons = 0;
    // local buffer for photons
    fBufferSize = 200000;
    NewPhotonBuffer(fBufferSize); 

    /*  Unisim photons want to know where they must be put. 
        Now such parameters are set in FileUnisimLightToEuso.
        In the future they should be retrieved from a detector
        configuration set of informations 

     */
//    ConfigFileParser *cfp=Config::Get()->GetCF("LightToEuso","FileUnisimLightToEuso") ;

    fEvCounter = 0;

    fFrame = new EsafRefFrame();
}

//_________________________________________________________________________________________ 
 UnisimPhotonsOnPupil::~UnisimPhotonsOnPupil() {
    //
    // Destructor
    //
    DeletePhotonBuffer();
    if ( fTheTruth )
        delete fTheTruth;

    SafeDelete(fFrame);
}

//_________________________________________________________________________________________ 
 void UnisimPhotonsOnPupil::NewPhotonBuffer(Int_t s) {
    //
    // Create a larger buffer and copy what is already in the old one
    // 
    ParentPhoton** newbuff = (ParentPhoton**) new char[4*s];
    for(Int_t i=0; i<s; i++)
        newbuff[i]=0;
    if ( fPhotons ) {
        for(Int_t i=0; i<s; i++) {
            if (i<fBufferSize)
                newbuff[i] = fPhotons[i];
        }
        delete[] fPhotons;
    }
    fPhotons = newbuff;
}

//_________________________________________________________________________________________ 
 void UnisimPhotonsOnPupil::DeletePhotonBuffer() {
    // 
    // Clean the buffer and release memory
    // 
    if ( fPhotons ) {
        for(Int_t i=0; i<fBufferSize; i++) {
            if (fPhotons[i])
                delete fPhotons[i];
            fPhotons[i]=0;
        }
        delete[] fPhotons;
    }
}
//_________________________________________________________________________________________ 
 Bool_t UnisimPhotonsOnPupil::SkipEvent( FILE* fp ) {
    //
    // Skip one event 
    //

    Clear();
    char s[5010];  

    fEvCounter++;

    Bool_t firstLine = kTRUE;
    while( kTRUE ) {
        char* ans = gzgets((void*)fp,s,5000);
        if ( ans == Z_NULL ) {
            MsgForm(EsafMsg::Warning,"Error reading gzgets in UnisimPhotonsOnPupil::Load()");
            return kFALSE;
        }
        if (!firstLine && strncmp(s," RUNHEADR RUNH",14) == 0) 
            return kTRUE;

        if ( firstLine ) firstLine = kFALSE;
    }
}
    
//_________________________________________________________________________________________ 
 Bool_t UnisimPhotonsOnPupil::Load( FILE* fp ) {
    // 
    // Load one event from file
    //

    Clear();
    fEvCounter++;

    Int_t nstart=0;
    if ( !LoadHeader( fp ) )
        return false;

    //
    // set the detector's altitude
    // 

    fFrame->SetPos(0.,0.,fHeader["REUSHEIG"]*m);
    
    char s[5010];  
    while( 1 ) {
        char* ans = gzgets((void*)fp,s,5000);
        if ( ans == Z_NULL ) {
            MsgForm(EsafMsg::Warning,"Error reading gzgets in UnisimPhotonsOnPupil::Load()");
            return false;
        }
        if (strncmp(s," RUNHEADR RUNH",14) == 0) 
            return true;
        Int_t n;
        double x[3],tm,w,th,ph,r,ph2;
        Int_t t;
        if ( sscanf(s,"%d %lf %lf %lf %lf %lf %lf %lf %lf %lf %d",
                    &n,&x[0],&x[1],&x[2],&tm,&w,&th,&ph,&r,&ph2,&t) != 11 ) {
            MsgForm(EsafMsg::Warning,"Error reading file in UnisimPhotonsOnPupil");
            return false;
        }


        nstart++;
        if ( n != nstart )
            MsgForm(EsafMsg::Warning,"WARNING: Photon counter not consequtive");

        x[0]*=m;
        x[1]*=m;
        x[2]*=m;
        tm *= microsecond;

        // angles are all degrees in the file
        th *= deg;
        ph *= deg;
        ph2 *= deg;

        // Wavelength is forced to 357 nm later
        w *= nm;

        // compute position on the pupil 
        /*  unisim photon positions (r, ph2) are given on a disk normal to photon direction
            in local frame. Now we have to put them on a disk normal to nadir direction
         */


        TVector3 phdir ;
        phdir.SetX( sin(th)*sin(ph) ) ;
        phdir.SetY( sin(th)*cos(ph) ) ;
        phdir.SetZ( cos(th) ) ;

        phdir = phdir.Unit() ;


        TVector3 rdisk ;
        r *= m;
        // rdisk in local frame    
        rdisk.SetX( r*cos(ph2) ) ;
        rdisk.SetY( r*sin(ph2) ) ;
        rdisk.SetZ( 0 ) ;


        // now transform rdsik to global frame      
        rdisk.RotateUz(phdir) ;

        if ( !fPhotonOndisk ){
            double y[3];
            rdisk.GetXYZ(y);
            if (!fPhotons[fNphotons]) 

                fPhotons[fNphotons] = new ParentPhoton();
            fPhotons[fNphotons++]->Build(n,x,y,th,ph,w,tm,(ParentPhotonType)t);
        } else {

            // calculate the impact point of the photon on a disk perpendicular to nadir
            Double_t tprop ;
            tprop = - rdisk.Z()/phdir.Z() ;

            TVector3 ydiskunisim ;
            ydiskunisim = rdisk + tprop*phdir ;


            // put the photons on the EUSO-ESAF entrance disk at Z coordinate = zdisk


            Double_t rdiskactual ;
            rdiskactual = sqrt( ydiskunisim.X()*ydiskunisim.X() + ydiskunisim.Y()*ydiskunisim.Y() ) ;

            if ( fRmaxDisk>=  rdiskactual ) {

                double y[3];
                y[0] = ydiskunisim.X() ;
                y[1] = ydiskunisim.Y() ;
                y[2] = fZdisk;


                if (!fPhotons[fNphotons]) 

                    fPhotons[fNphotons] = new ParentPhoton();
                fPhotons[fNphotons++]->Build(n,x,y,th,ph,w,tm,(ParentPhotonType)t);

            }
        }

    }

    return true;
}

//_________________________________________________________________________________________ 
 Photon* UnisimPhotonsOnPupil::GetPhoton() {
    // 
    // returns a Photon object
    //
    //if ( (( fCurrent % 200 ) == 0) && (fCurrent) )
    //    cout << "Event: " << event << "  Photon: " << fCurrent << "n";
    if ( fCurrent == fNphotons )
        return (Photon*)0;
    return &(fPhotons[fCurrent++]->GetPhoton());
}

//_________________________________________________________________________________________ 
 void UnisimPhotonsOnPupil::Clear() {
    // 
    // Reset 
    // 
    fNphotons = 0;
    fCurrent = 0;
    fHeader.clear();
    fTheTruth = 0;
}

//_________________________________________________________________________________________ 
 Bool_t UnisimPhotonsOnPupil::LoadHeader( FILE* fp) {
    //
    // Load header and store it into a map (key,value)
    //
    char s[5010];
    while(1) {
        char* ans = gzgets((void*)fp,s,5000);
        if ( ans == Z_NULL ) {
            MsgForm(EsafMsg::Warning, "Error reading gzgets in UnisimPhotonsOnPupil::LoadHeader() (2)");
            return false;
        }
        if (strncmp(s," BEVTBUFF EVTB",14) != 0) {
            char key[200];
            double val=0.;
            if ((strncmp(s," EVTHEADR EVTH",14) != 0) && (strncmp(s," RUNHEADR RUNH",14) != 0)) {
                if ( sscanf(s,"%s %lf",key,&val) != 2 ) {
                    MsgForm(EsafMsg::Warning,"Data format error in UnisimPhotonsOnPupil::LoadHeader()");
                    Msg(EsafMsg::Warning) << s << MsgDispatch;
                    return false;
                }
                fHeader[key]=val;
            }
        } else {
            for(Int_t i=0;i<11;i++)  // skip 11 lines
                gzgets((void*)fp,s,5000);
            break;
        }
    }
    return true;
}
//_________________________________________________________________________________________ 
 MCTruth* UnisimPhotonsOnPupil::Truth() {
    if (!fTheTruth) {
        fTheTruth = new MCTruth();
        fTheTruth->SetEnergy(fHeader["EPRMENRG"]*GeV); // now in MeV
        fTheTruth->SetThetaPhi(fHeader["EPOLANGD"]*deg, fHeader["EAZIANGD"]*deg);
        EVector v;
        v[0] = fHeader["EX0COMES"]*m;
        v[1] = fHeader["EY0COMES"]*m;
        v[2] = fHeader["EZ0COMES"]*m;
        fTheTruth->SetFirstInt(v,fHeader["EPRMDPTH"]*g/cm2);
        v[0] = fHeader["EXDMPPOS"]*m;
        v[1] = fHeader["EYDMPPOS"]*m;
        v[2] = fHeader["EZDMPPOS"]*m;
        // FIXME, Age on the Earth is not stored in files. Puting arbitrary number 2.
        fTheTruth->SetEarthImpact(v,2);             
        // In case no value for ShowerMax is present  Puting arbitrary number -100.
        v[0] = -100;
        v[1] = -100;
        v[2] = -100;
        map<string,double>::iterator p;
        p = fHeader.find("EXMAXPOS") ;
        if ( p != fHeader.end() ) v[0] = fHeader["EXMAXPOS"]*m;
        p = fHeader.find("EYMAXPOS") ;
        if ( p != fHeader.end() ) v[1] = fHeader["EYMAXPOS"]*m;
        p = fHeader.find("EZMAXPOS") ;
        if ( p != fHeader.end() ) v[2] = fHeader["EZMAXPOS"]*m;         


        fTheTruth->SetShowerMax(v,fHeader["ESHOWMAX"]*g/cm2);
	
;
    
        fTheTruth->SetHclouds(0.);
        p = fHeader.find("ESTAHEIG") ;
        if ( p != fHeader.end() ) fTheTruth-> SetHclouds(fHeader["ESTAHEIG"]*km); // in ESAF units
	
    }
  return fTheTruth;
}
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