// ESAF : Euso Simulation and Analysis Framework
// $Id: UnisimFileShowerGenerator.cc,v 1.20 2006/03/20 21:51:35 thea Exp $
// sergio bottai created Jan, 22 2004

#include "UnisimFileShowerGenerator.hh"
#include "MCTruth.hh"
#include "ShowerTrack.hh"
#include <zlib.h>
#include "EarthVector.hh"
#include "Atmosphere.hh"
#include <TMath.h>
#include <TSystem.h>

using namespace sou;
using TMath::RadToDeg;

ClassImp(UnisimFileShowerGenerator)

//______________________________________________________________________________
 UnisimFileShowerGenerator::UnisimFileShowerGenerator() : 
    EventGenerator("Unisimfile"), fTrack(0),fFile(NULL), fTruth(0) {
    //
    // Construtor
    //
    
    fFileName = "";
    fFirstEvent = 0;
    fCurrentEvent = 0;
}

//______________________________________________________________________________
 UnisimFileShowerGenerator::~UnisimFileShowerGenerator() {
    //
    // Destructor
    //

    if ( fTrack ) {
        delete fTrack;
        fTrack=NULL; 
    }
    if ( fTruth ) {
        delete fTruth;  
        fTruth=NULL;
    }
    
    
}

//______________________________________________________________________________
 void UnisimFileShowerGenerator::ReplaceInputFile(const char* fn) {
    //
    //
    //
    
    if ( fFile )
        throw runtime_error("UnisimFileShowerGenerator: Cannot change name when file is already open!");
    if (!fn) {
        Msg(EsafMsg::Error) << "UnisimFileShowerGenerator: NULL file name ignored" << MsgDispatch;
        return;
    }
    Conf()->ReplaceStr("UnisimFileShowerGenerator.FileName",fn);
    Msg(EsafMsg::Info) << "Input file name changed to " << fn << MsgDispatch;
}



//______________________________________________________________________________
 void UnisimFileShowerGenerator::Open() {
    //
    // Open file
    //

    if ( Conf()->IsDefined( "UnisimFileShowerGenerator.FileName" ) ) {
        fFileName = Conf()->GetStrPath("UnisimFileShowerGenerator.FileName");
        fFirstEvent = (Int_t)Conf()->GetNum("UnisimFileShowerGenerator.FirstEvent") ;
        if (fFirstEvent < 0 ) {
            Msg(EsafMsg::Warning) << "First event must be >= 0. Now forced fFirstEvent=0" << MsgDispatch;
            fFirstEvent = 0;
        }
    }
    else
        Msg(EsafMsg::Panic) << "Unknown input file name in UnisimFileShowerGenerator" << MsgDispatch;

    fFile = (FILE*)gzopen( fFileName.c_str(), "rb" );

    if( fFile == NULL ) 
        Msg(EsafMsg::Panic) << "Unknown input file name in UnisimFileShowerGenerator" << MsgDispatch;
    else
        Msg(EsafMsg::Info) << "File: " << fFileName << " successfully opened " << MsgDispatch;

}

//______________________________________________________________________________
 void UnisimFileShowerGenerator::Close() {
    //
    // Close file
    //

    if( fFile  ) gzclose( (void*)fFile ) ;
    
}


//______________________________________________________________________________
 PhysicsData* UnisimFileShowerGenerator::Get() {
    //
    // Get a new event
    //


    // build and return a track

    if ( !fFile ) Open();  
    Int_t event=0;

    while (1) {

        if ( fTrack ) {
            delete fTrack;
            fTrack=NULL; 
        }
        if ( fTruth ) {
            delete fTruth;  
            fTruth=NULL;
        }

        if (fCurrentEvent >= fFirstEvent) {
            if ( LoadTrack(event) ){

                MsgForm(EsafMsg::Info,"UnisimFile event %d, header id %d", fCurrentEvent, (Int_t)fHeader["EVTNUMBR"]);
                MsgForm(EsafMsg::Info,
                        "ShowerTrack produced : Energy = %.2E MeV, Theta= %.2f deg, Phi= %.2f deg",
                        fTrack->GetEnergy(), fTrack->GetTheta()*RadToDeg(),fTrack->GetPhi()*RadToDeg());
                MsgForm(EsafMsg::Info,
                        "         First interaction X1 = %.3f at : (%.2f,%.2f,%.2f) km gives a shower with %d steps",
                        fTrack->fX1*cm2/g,fTrack->GetInitPos().X()/km,fTrack->GetInitPos().Y()/km,
                        fTrack->GetInitPos().Z()/km, fTrack->GetNumStep());
            } else {
                Msg(EsafMsg::Warning) << "UnisimFileShowerGenerator: End of file "
                    << fFileName << " reached? " << MsgDispatch;
                fTrack = (ShowerTrack*)NULL;
            }
            fCurrentEvent++;
            break;
        } else {
            if ( !LoadTrack(event) ){
                Msg(EsafMsg::Warning) << "UnisimFileShowerGenerator: End of file "
                    << fFileName << " reached? " << MsgDispatch;
                fTrack = (ShowerTrack*)NULL;
                fCurrentEvent++;
                break;
            }
        }

        fCurrentEvent++;

    }
    return fTrack;
}

//______________________________________________________________________________
 MCTruth* UnisimFileShowerGenerator::GetTruth() {
    // 
    // get pointer to MonteCarlo Truth
    //

    if ( fTrack ) {
        return fTruth;
    }
        else { Msg(EsafMsg::Warning) << "Event to be loaded with LoadTrack before GetTruth" << MsgDispatch;
        return (MCTruth*)0;
    }
    

}

//______________________________________________________________________________
 bool UnisimFileShowerGenerator::LoadTrack( Int_t event ) {
    //
    // Load one event from file
    //
    
    if ( !LoadHeader() ) return false;

    char s[5010]; 

    fTruth= new MCTruth();
     
    LoadTruth();

    fTrack= new ShowerTrack();
    fTrack->fEthrEl=0.;
    
    fTrack->fDirVers[0]=fHeader["ELDIRCOS"];
    fTrack->fDirVers[1]=fHeader["EMDIRCOS"];
    fTrack->fDirVers[2]=fHeader["ENDIRCOS"];
    
    fTrack->fInitPos[0]=fHeader["EX0COMES"]*m;
    fTrack->fInitPos[1]=fHeader["EY0COMES"]*m;
    fTrack->fInitPos[2]=fHeader["EZ0COMES"]*m;
    
    if(fHeader["HITGROUN"]==0) fTrack->fHitGround=false ;
    else fTrack->fHitGround=true ;
    
    fTrack->fEnergy    = fHeader["EPRMENRG"]*GeV;
    fTrack->fTheta     = fHeader["EPOLANGD"]*deg; 
    fTrack->fPhi       = fHeader["EAZIANGD"]*deg;
    fTrack->fX1        = fHeader["EPRMDPTH"]*g/cm2;
    fTrack->fEarthImpact.SetXYZ(fHeader["EXDMPPOS"]*m,fHeader["EYDMPPOS"]*m,fHeader["EZDMPPOS"]*m);
   
    while( 1 ) {
        char* ans = gzgets((void*)fFile,s,5000);

        if ( ans == Z_NULL ) {
            Msg(EsafMsg::Error) << "Error reading gzgets in UnisimFileShowerGenerator::Load()n" << MsgDispatch;
            delete fTrack;
            fTrack=NULL;
            return false;
        }
        if (strncmp(s," BEVTBEND EEND",14) == 0) return true;
        if (strncmp(s," STEPNUMB",9) == 0) {

            char key;
            int n;
            if ( sscanf(s,"%s %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",
                        &key,&n,&fStep.fXi,&fStep.fXf,&fStep.fXYZi[0],&fStep.fXYZi[1],&fStep.fXYZi[2],
                        &fStep.fXYZf[0],&fStep.fXYZf[1],&fStep.fXYZf[2],
                        &fStep.fTimei,&fStep.fTimef,&fStep.fAgei,&fStep.fAgef,&fStep.fNelectrons,
                        &fStep.fNcharged,&fStep.fEloss,&fStep.fNcherenkov) != 18 ) {
                Msg(EsafMsg::Error) << "Data format error in UnisimFileShowerGenerator::LoadHeader()" << MsgDispatch;
                Msg(EsafMsg::Error) << s << MsgDispatch;
                return false;
            }

            fStep.fXYZi[0]=fStep.fXYZi[0]*m;
            fStep.fXYZi[1]=fStep.fXYZi[1]*m;
            fStep.fXYZi[2]=fStep.fXYZi[2]*m;

            fStep.fXYZf[0]=fStep.fXYZf[0]*m;
            fStep.fXYZf[1]=fStep.fXYZf[1]*m;
            fStep.fXYZf[2]=fStep.fXYZf[2]*m;

            fStep.fTimei=fStep.fTimei*microsecond;	
            fStep.fTimef=fStep.fTimef*microsecond;	

            fStep.fXi=fStep.fXi*gram/cm2;
            fStep.fXf=fStep.fXf*gram/cm2;

            fTrack->Add(fStep);

        }
    }
       
    return true;
}

//______________________________________________________________________________
 bool UnisimFileShowerGenerator::LoadHeader() {
    // 
    // Load header and store it into a map (key,value)
    //


    fHeader.clear();

    char s[5010];
    while(1) {
        char* ans = gzgets((void*)fFile,s,5000);

        if ( ans == Z_NULL ) {
            Msg(EsafMsg::Error) << "Error reading gzgets in UnisimFileShowerGenerator::LoadHeader() (2) " << MsgDispatch;
            return false;
        }
        if (strncmp(s," BEVTBUFF EVTB",14) == 0) break;

        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 ) {
                Msg(EsafMsg::Error) << "Data format error in UnisimFileShowerGenerator::LoadHeader()" << MsgDispatch;
                Msg(EsafMsg::Error) << s << MsgDispatch;
                return false;
            }
            fHeader[key]=val;
        }
    }
    return true;
}


//______________________________________________________________________________
 void UnisimFileShowerGenerator::LoadTruth() {
    //
    //
    //

    fTruth->SetEnergy(fHeader["EPRMENRG"]*GeV); // in ESAF units (MeV)
    fTruth->SetThetaPhi(fHeader["EPOLANGD"]*deg, fHeader["EAZIANGD"]*deg);
    EVector v;
    v[0] = fHeader["EX0COMES"]*m;
    v[1] = fHeader["EY0COMES"]*m;
    v[2] = fHeader["EZ0COMES"]*m;
    fTruth->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.
    fTruth->SetEarthImpact(v,2);     
    
    // recalculate TOAimpact
    EarthVector first_inter(fHeader["EX0COMES"]*m,fHeader["EY0COMES"]*m,fHeader["EZ0COMES"]*m);
    EarthVector Omega(1);
    EarthVector toaimpact(1);
    Omega.SetMagThetaPhi(1.,fHeader["EPOLANGD"]*deg,fHeader["EAZIANGD"]*deg);
    toaimpact = Atmosphere::Get()->ImpactAtTOA(first_inter,EarthVector(-Omega));
    fTruth->SetTOAImpact(toaimpact);
            
    // In case no value for ShowerMax is present  Puting arbitrary number 0.
    v[0] = 0;
    v[1] = 0;
    v[2] = 0;
    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;         


    fTruth->SetShowerMax(v,fHeader["ESHOWMAX"]*g/cm2);
    
    fTruth-> SetHclouds(0.);
    p = fHeader.find("ESTAHEIG") ;
    if ( p != fHeader.end() ) fTruth-> SetHclouds(fHeader["ESTAHEIG"]*km); // in ESAF units

}