RdBeamFormer.cc

Go to the documentation of this file.

#include "RdBeamFormer.h"

#include <evt/Event.h>
#include <evt/ShowerRecData.h>
#include <evt/ShowerRRecData.h>
#include <evt/ShowerSimData.h>
#include <evt/BeamQuantities.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/Channel.h>
#include <revt/StationRecData.h>

#include <det/Detector.h>
#include <rdet/RDetector.h>

#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/ErrorLogger.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <utl/Triple.h>
#include <utl/FFTDataContainerAlgorithm.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <utl/UTMPoint.h>
#include <utl/CoordinateSystem.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/ReferenceEllipsoid.h>
#include <utl/RadioGeometryUtilities.h>

/*
This headers should be uncommented if the footprint estimation fuctionality is used

#include <atm/AerosolDB.h>
#include <atm/AerosolZone.h>
#include <atm/ProfileResult.h>
#include <atm/USStdADBProfileModel.h>
#include <atm/MonthlyAvgDBProfileModel.h>
#include "ReadMolecularHL.h"
*/

#include <fwk/CentralConfig.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/MagneticFieldModel.h>

#include <sstream>
#include <vector>
#include <algorithm>
#include <cmath>

using namespace std;
using namespace revt;
using namespace utl;
using namespace fwk;
using namespace det;

/*
Should be uncommented if the footprint estimation fuctionality is used

using namespace atm;
using namespace ReadMolecularHLNS;
*/

//Just an abbreviation of repeated code from UserModule.cc skeleton
#define PRINT(...) {stringstream msg; msg << __VA_ARGS__; INFO(msg.str());}

#define SSTR( x ) dynamic_cast< std::ostringstream & >(( std::ostringstream() << std::dec << x ) ).str()

namespace RdBeamFormer {

  RdBeamFormer::RdBeamFormer() : fNumCalls(0)
  {
  }

  RdBeamFormer::~RdBeamFormer()
  {
  }

  VModule::ResultFlag
  RdBeamFormer::Init()
  {
    INFO("RdBeamFormer::Init()");

    //read in the configurations of the xml file
    Branch conf = CentralConfig::GetInstance()->GetTopBranch("RdBeamFormer");

    conf.GetChild("startBin").GetData(startbin);
    conf.GetChild("stopBin").GetData(stopbin);
    conf.GetChild("offsetFactor").GetData(offsetfactor);
    
    conf.GetChild("waveModel").GetData(waveModelConf);
    
    if (waveModelConf == "ePlane")
       waveModel = new RdPlaneWaveModel;
    else if (waveModelConf == "eSpherical")
       waveModel = new RdSphericalWaveModel;
    else if (waveModelConf == "eConical")
       waveModel = new RdConicalWaveModel;
    else if (waveModelConf == "eHyperbolic")
       waveModel = new RdHyperbolicWaveModel;
       
    std::string confstring;
    const Branch txtOut = conf.GetChild("asciiOutput");
    if (txtOut) {
      txtOut.GetChild("outputFileName").GetData(outfile);
      txtOut.GetChild("outputGridType").GetData(confstring);
      if (confstring == "ePolar")
         gridtype = polar;
      else if (confstring == "eCarthesian")
         gridtype = carthesian;
      else
         gridtype = none;      
    } else {
      gridtype = none; outfile = "";
    }

    return eSuccess;
  }

  double RdBeamFormer::shiftTraces(std::vector<revt::StationFFTDataContainer>& stationData, int& nStat, double& binning, 
                                      int& start, int& stop, evt::Event& event) {
    const evt::ShowerRRecData& rrec= event.GetRecShower().GetRRecShower();
    const utl::CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    coreCS = fwk::LocalCoordinateSystem::Create(rrec.GetCorePosition());
    waveModel->setOrigin(rrec.GetCorePosition());        
    Point skyPos=rrec.GetCorePosition()+rrec.GetAxis();
    waveModel->setSkyPos (skyPos);

    static FFTDataContainerAlgorithm fft;

    std::vector<int> start_bins;
    std::vector<int> stop_bins;
    std::vector<double> shifts;

    const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();

    REvent& rEvent = event.GetREvent();  

    REvent::CandidateStationIterator it = rEvent.CandidateStationsBegin(); 
    binning = it->GetConstStationTimeSeries().GetBinning();
    double offset = (stop - start) * binning * offsetfactor;

    //desired precision for binning comparison
    double eps = 1e-14;        
     
    //loop over the stations which were not rejected and calculate the shifts
    for (it = rEvent.CandidateStationsBegin(); it != rEvent.CandidateStationsEnd(); ++it) {
        Station& station = *it;
        if (station.IsSaturated()) //skip saturated stations
           continue;        

        static double old_binning = it->GetConstStationTimeSeries().GetBinning();
        double new_binning = it->GetConstStationTimeSeries().GetBinning();     
        if (abs(new_binning - old_binning) >= eps) {
           ERROR ("Traces have different binning!");
           return eFailure;           
        }
 
        //will be used to ensure that all traces have the same start and stop bins
        start_bins.push_back(it->GetConstStationTimeSeries().GetStart());       
        stop_bins.push_back(it->GetConstStationTimeSeries().GetStop());           

        const StationRecData& sRec = it->GetRecData();
   
        double traceStartTime = 0.;
        if (sRec.HasParameter(eTraceStartTime))
           traceStartTime = sRec.GetParameter(eTraceStartTime);
        else {
           ERROR ("Parameter eTraceStartTime is not set!");
           return eFailure;
        }              

        double signalTime = 0.;
        if (sRec.HasParameter(eSignalTime))
           signalTime = sRec.GetParameter(eSignalTime);
        else {
           ERROR ("Parameter eSignalTime is not set!");
           return eFailure;
        }    

        const Point pos = rDetector.GetStation(*it).GetPosition();
        double shift = waveModel->delay(pos, rrec);
        shift=shift-traceStartTime;       
        shifts.push_back(-shift + offset);       
    }

    //now we will find minimum and maximum shifts to be used later
    std::vector<double>::iterator min_shift_iter;
    std::vector<double>::iterator max_shift_iter;

    min_shift_iter = std::min_element(shifts.begin(), shifts.end());
    max_shift_iter = std::max_element(shifts.begin(), shifts.end());

    int min_pos = min_shift_iter - shifts.begin();
    int max_pos = max_shift_iter - shifts.begin();

    double min_shift = *min_shift_iter;
    double max_shift = *max_shift_iter;
    //cout<<"min_shift= "<<min_shift<<" max_shift= "<<max_shift<<"\n";

    //ensure that all traces have the same start and stop bins
    int max_start = *std::max_element(start_bins.begin(), start_bins.end());
    start = max_start;
    int min_stop = *std::min_element(stop_bins.begin(), stop_bins.end());
    stop = min_stop;
    //cout<<"max_start= "<<max_start<<" min_stop= "<<min_stop<<"\n";               

    /*
    Uncomment this part if you want to write to a file traces before and after shifting

    std::string trace_before = "trace_before_shift";
    static std::string trace_name_before = trace_before + "_";
    std::string trace_after = "trace_after_shift";
    static std::string trace_name_after = trace_after + "_"; 
    
    int eventId = event.GetREvent().GetHeader().GetId();
    std::string stringId = SSTR(eventId);
    //cout<<"eventId= "<<eventId<<"\n";
    static int Id = -1;
    //cout<<"oldId= "<<oldId<<"\n";
    if (eventId != Id) {
       //cout<<"eventId= "<<eventId<<"\n";    
       Id = eventId;
       //cout<<"oldId= "<<oldId<<"\n";
       std::string buffer = SSTR(Id) + "_";          
       trace_name_before.replace(trace_before.size()+1,buffer.size(),buffer);  
       trace_name_after.replace(trace_after.size()+1,buffer.size(),buffer);  
    }     
    
    if (waveModelConf == "eSpherical") {       
       const utl::Triple polarCoords = rrec.GetAxis().GetSphericalCoordinates(coreCS);
       double r = polarCoords.get<0>()/km;
       static double R = -1;
       //cout<<"oldR= "<<oldR<<" oldR-r= "<<oldR - r<<"\n";       
       if (abs(r - R) >= eps) {
          //cout<<"r= "<<r<<"\n";
          R = r;
          //cout<<"oldR= "<<oldR<<"\n";          
          std::string buffer = SSTR(R);          
          trace_name_before.replace(trace_before.size()+stringId.size()+2,buffer.size()+10,buffer); 
          trace_name_after.replace(trace_after.size()+stringId.size()+2,buffer.size()+10,buffer);            
       }
    }
    else if (waveModelConf == "eConical") {
       double Rho = rrec.GetParameter(eRho)/degree;
       //cout<<"eRho= "<<Rho<<"\n";
       static double old_Rho = -1;   
       //cout<<"oldRho= "<<oldRho<<"\n"; 
       if (abs(Rho - old_Rho) >= eps) {
          //cout<<"eRho= "<<Rho<<"\n";
          old_Rho = Rho;
          //cout<<"oldRho= "<<oldRho<<"\n";                    
          std::string buffer = SSTR(old_Rho);                           
          trace_name_before.replace(trace_before.size()+stringId.size()+2,buffer.size()+10,buffer); 
          trace_name_after.replace(trace_after.size()+stringId.size()+2,buffer.size()+10,buffer);                 
       }
    }
    else if (waveModelConf == "eHyperbolic") {
       double Rho = rrec.GetParameter(eRho)/degree;
       double b = rrec.GetParameter(eB)/ns;
       //double b = 3*ns;
       //cout<<"eRho= "<<Rho<<"\n";
       static double old_Rho = -1; 
       static double old_b = -1;  
       //cout<<"oldRho= "<<oldRho<<"\n"; 
       if (abs(Rho - old_Rho) >= eps || abs(b - old_b) >= eps) {
          //cout<<"eRho= "<<Rho<<"\n";
          old_Rho = Rho;
          old_b = b;
          //cout<<"oldRho= "<<oldRho<<"\n";                    
          std::string buffer1 = SSTR(old_Rho);       
          std::string buffer2 = SSTR(old_b);  
          std::string buffer = buffer1 + "_" + buffer2;                                             
          trace_name_before.replace(trace_before.size()+stringId.size()+2,buffer.size()+10,buffer); 
          trace_name_after.replace(trace_after.size()+stringId.size()+2,buffer.size()+10,buffer);                 
       }
    }

    std::string trace1 = trace_name_before + ".txt";
    std::string trace2 = trace_name_after + ".txt";

    ofstream fout1;
    fout1.open(trace1.c_str(), ios::out|ios::app);

    ofstream fout2;
    fout2.open(trace2.c_str(), ios::out|ios::app);
    */

    //loop over the stations which were not rejected and perform the shifts
    int i = 0;
    for (it = rEvent.CandidateStationsBegin(); it != rEvent.CandidateStationsEnd(); ++it) {
        Station& station = *it;
        if (station.IsSaturated()) //skip saturated stations
           continue;    

        stationData[i] = it->GetFFTDataContainer();   

        const StationRecData& sRec = it->GetRecData();
   
        double traceStartTime = sRec.GetParameter(eTraceStartTime);        

        /*
        Uncomment this part if you want to write to a file traces before and after shifting

        const StationTimeSeries& ts1 = stationData[i].GetTimeSeries();
        StationTimeSeries::ConstIterator iter1 = ts1.Begin() + start;
        for (int j = 0; j < (stop-start); ++j) {
            fout1<<(j+start)*binning<<" "<<iter1->GetMag()<<"\n";
            iter1++;
        }
      
        fout1<<"\n";
        fout1<<"\n"; 
        */
 
        //if all shifts have the same sign then we can shift only relatively to the smallest shift, which stays untouched
        if (min_shift >= 0 && max_shift >=0) {
           fft.ShiftTimeSeries(stationData[i], shifts[i] - min_shift);  
           //cout<<"shifts[i]= "<<shifts[i]<<" shifts[i]+offset= "<<shifts[i] + offset
               //<<" shifts[i]-min_shift= "<<shifts[i] - min_shift<<"\n";
        }
        else if (min_shift < 0 && max_shift < 0) {
           fft.ShiftTimeSeries(stationData[i], shifts[i]- max_shift); 
           //cout<<"shifts[i]= "<<shifts[i]<<" shifts[i]+offset= "<<shifts[i] + offset
               //<<" shifts[i]-max_shift= "<<shifts[i] - max_shift<<"\n";
        }   
        else if (min_shift <= 0 && max_shift >= 0) {
           fft.ShiftTimeSeries(stationData[i], shifts[i]);
           //cout<<"shifts[i]= "<<shifts[i]<<" shifts[i]+offset= "<<shifts[i] + offset<<"\n";
        }  

        //}
        ++i;
    }    

    //if after shifting there is no minimum required overlap (in our case it is 100 bins) for all stations, 
    //then we need to drop the most distant station. This is done iteratively untill the condition is fulfilled.
    if (min_shift >= 0 && max_shift >=0) {
       //cout<<"start= "<<start<<" stop= "<<stop<<"\n";
       start = (max_shift - min_shift)/binning + 1;
       //cout<<"start= "<<start<<" stop= "<<stop<<"\n";
       if (start < max_start) 
          start = max_start;
       //cout<<"start= "<<start<<" stop= "<<stop<<"\n";
       while (stop - start < 100) {
          shifts.erase(shifts.begin() + max_pos);
          stationData.erase(stationData.begin() + max_pos);

          min_shift_iter = std::min_element(shifts.begin(), shifts.end());
          max_shift_iter = std::max_element(shifts.begin(), shifts.end());

          min_pos = min_shift_iter - shifts.begin();
          max_pos = max_shift_iter - shifts.begin();

          min_shift = *min_shift_iter;
          max_shift = *max_shift_iter;

          //cout<<"min_shift= "<<min_shift<<" max_shift= "<<max_shift<<"\n";

          start = (max_shift - min_shift)/binning + 1;
          //cout<<"start= "<<start<<" stop= "<<stop<<"\n";

          if (start < max_start) 
             start = max_start;
          //cout<<"start= "<<start<<" stop= "<<stop<<"\n";             

          --nStat;
       }
    }
    else if (min_shift < 0 && max_shift < 0) {
       //cout<<"start= "<<start<<" stop= "<<stop<<"\n";
       stop = min_stop - (max_shift - min_shift)/binning;
       //cout<<"start= "<<start<<" stop= "<<stop<<"\n";
       if (stop > min_stop)
          stop = min_stop;
       //cout<<"start= "<<start<<" stop= "<<stop<<"\n";
       while (stop - start < 100) {
          shifts.erase(shifts.begin() + min_pos);
          stationData.erase(stationData.begin() + min_pos);

          min_shift_iter = std::min_element(shifts.begin(), shifts.end());
          max_shift_iter = std::max_element(shifts.begin(), shifts.end());

          min_pos = min_shift_iter - shifts.begin();
          max_pos = max_shift_iter - shifts.begin();

          min_shift = *min_shift_iter;
          max_shift = *max_shift_iter;

          //cout<<"min_shift= "<<min_shift<<" max_shift= "<<max_shift<<"\n";

          stop = min_stop - (max_shift - min_shift)/binning;
          //cout<<"start= "<<start<<" stop= "<<stop<<"\n";

          if (stop > min_stop)
              stop = min_stop;
          //cout<<"start= "<<start<<" stop= "<<stop<<"\n";

          --nStat;   
       }
    }
    else if (min_shift <= 0 && max_shift >= 0) {
       //cout<<"start= "<<start<<" stop= "<<stop<<"\n";
       start = max_shift/binning + 1;
       stop = min_stop - abs(min_shift)/binning;
       //cout<<"start= "<<start<<" stop= "<<stop<<"\n";
       while (stop - start < 100) {
          if (abs(max_shift) - abs(min_shift) > 0) {
             shifts.erase(shifts.begin() + max_pos);
             stationData.erase(stationData.begin() + max_pos);

             min_shift_iter = std::min_element(shifts.begin(), shifts.end());
             max_shift_iter = std::max_element(shifts.begin(), shifts.end());

             min_pos = min_shift_iter - shifts.begin();
             max_pos = max_shift_iter - shifts.begin();

             min_shift = *min_shift_iter;
             max_shift = *max_shift_iter;

             //cout<<"min_shift= "<<min_shift<<" max_shift= "<<max_shift<<"\n";

             start = max_shift/binning + 1;
             stop = min_stop - abs(min_shift)/binning;
             //cout<<"start= "<<start<<" stop= "<<stop<<"\n";

             if (start < max_start) 
                start = max_start;
             if (stop > min_stop)
                stop = min_stop;
             //cout<<"start= "<<start<<" stop= "<<stop<<"\n";             

             --nStat;
          }
          else {
             shifts.erase(shifts.begin() + min_pos);
             stationData.erase(stationData.begin() + min_pos);

             min_shift_iter = std::min_element(shifts.begin(), shifts.end());
             max_shift_iter = std::max_element(shifts.begin(), shifts.end());

             min_pos = min_shift_iter - shifts.begin();
             max_pos = max_shift_iter - shifts.begin();

             min_shift = *min_shift_iter;
             max_shift = *max_shift_iter;

             //cout<<"min_shift= "<<min_shift<<" max_shift= "<<max_shift<<"\n";

             start = max_shift/binning + 1;
             stop = min_stop - abs(min_shift)/binning;
             //cout<<"start= "<<start<<" stop= "<<stop<<"\n";

             if (start < max_start) 
                start = max_start;
             if (stop > min_stop)
                stop = min_stop;
             //cout<<"start= "<<start<<" stop= "<<stop<<"\n";

             --nStat;
          }
       }
    }    

    /*
    Uncomment this part if you want to write to a file traces before and after shifting

    for (int i = 0; i < nStat; ++i) { 

        const StationTimeSeries& ts2 = stationData[i].GetTimeSeries();
        StationTimeSeries::ConstIterator iter2 = ts2.Begin() + start;
        for (int k = 0; k < (stop-start); ++k) {
            fout2<<(k+start)*binning<<" "<<iter2->GetMag()<<"\n";
            iter2++;
        }
      
        fout2<<"\n";
        fout2<<"\n";
    }

    fout1.close();
    fout2.close();
    */

    return offset - start*binning;
  }

  inline void RdBeamFormer::efieldproduct (const StationTimeSeries& in1, const StationTimeSeries& in2, 
                                              TraceD& out, int start, int stop) {
    StationTimeSeries::ConstIterator it1 = in1.Begin() + start;
    StationTimeSeries::ConstIterator it2 = in2.Begin() + start;
    int n = it1->GetSize(); //should be 3
    //cout<<n<<"\n";

    for (int i = 0; i < (stop-start); ++i) {
      double prod = 0.0;
      for (int j = 0; j < n; ++j) {
        prod += it1->At(j) * it2->At(j);
      }
      out[i] = prod;
      ++it1;
      ++it2;
    }
  }

  TraceD RdBeamFormer::crosscorr (std::vector<revt::StationFFTDataContainer>& stationData, int nStat, 
                                     double binning, int start, int stop) {
    TraceD sum(stop-start, binning);
    TraceD addend(stop-start, binning);
    TraceD res(stop-start, binning);

    double pair_of_Stat = nStat*(nStat-1)/2;
    
    for (int i=0; i < nStat; ++i) {
      for (int j=i+1; j < nStat; ++j) {
        efieldproduct (stationData[i].GetTimeSeries(),stationData[j].GetTimeSeries(), addend, start, stop);
        sum += addend;        
      }
    }

    /*
    In case you want to have a square root from the trace
    for (int i = 0; i < (stop-start); ++i) {
      if (sum[i]>0) 
         res[i]=sqrt(sum[i]);
      else 
         res[i]=-sqrt(abs(sum[i]));   
    }
    return res*(1.0/pair_of_Stat);
    */

    return sum*(1.0/pair_of_Stat);
  }  

  inline void RdBeamFormer::efieldtopower (const StationTimeSeries& in, TraceD& out, int start, int stop) {
    StationTimeSeries::ConstIterator it = in.Begin() + start;
    for (int i = 0; i < (stop-start); ++i) {
        out[i] = it->GetMag2();
        ++it;
    }
  }

  TraceD RdBeamFormer::powertrace (std::vector<revt::StationFFTDataContainer>& stationData, 
                                      int nStat, double binning, int start, int stop) {
    TraceD sum(stop-start, binning);
    TraceD addend(stop-start, binning);
    TraceD res(stop-start, binning);

    for (int i = 0; i<nStat; ++i) {
        const StationTimeSeries& ts = stationData[i].GetTimeSeries();
        efieldtopower (ts, addend, start, stop);
        sum += addend;      
    }

    /*
    In case you want to have a square root from the trace

    for (int i = 0; i < (stop-start); ++i) 
        res[i]=sqrt((sum[i]));
    return res*(1.0/nStat);
    */

    return sum*(1.0/nStat); 
  }

  TraceD RdBeamFormer::xtrace (TraceD ccTrace, TraceD ccNormal, TraceD pwNormal) {
    TraceD result(ccTrace);    
    TraceD::ConstIterator itCC = ccNormal.Begin();
    TraceD::ConstIterator itPW = pwNormal.Begin();
    for (TraceD::Iterator it = result.Begin(); it != result.End(); ++it) {
        (*it) *= abs ( (*itCC) / (*itPW) );
        ++itCC;
        ++itPW;
    }
    return result;
  }

  evt::BeamPeak RdBeamFormer::findPeak (const TraceD& trace, const double timeOffset) {
    TraceD::ConstIterator max = max_element(trace.Begin(), trace.End());
    const double maxval = *max;
    const int maxi = (max - trace.Begin());
    const double maxtime = maxi * trace.GetBinning();

    //cout<<"maxval= "<<maxval<<" maxi= "<<maxi<<" maxtime= "<<maxtime<<"\n";

    int start;
    int stop;

    //devide the trace in two parts and check which of them contains the maximum value
    if (maxi > (trace.GetStart() + trace.GetStop())/2) {
       start = trace.GetStart() + 1;
       stop = (trace.GetStart() + trace.GetStop())/2 - (trace.GetStop() - trace.GetStart())*0.05;
       //cout<<start<<" "<<stop<<"\n";
    } 
    else {
       start = (trace.GetStart() + trace.GetStop())/2 + (trace.GetStop() - trace.GetStart())*0.05;
       stop = trace.GetStop() - 1;
       //cout<<start<<" "<<stop<<"\n";
    }
 
    //calculate RMS and mean value for the part of the trace which doesn`t contain the maximum value
    static TraceAlgorithm algo;
    const double rms = algo.RootMeanSquare ( trace, start, stop);
    const double mean = algo.Mean ( trace, start, stop);
    const double snr = maxval/rms;

    //cout<<"mean= "<<mean*V/m<<" rms= "<<rms*V/m<<" SNR= "<<snr*V/m<<"\n";

    return evt::BeamPeak ( maxval, maxtime, rms, snr);
  }

  VModule::ResultFlag
  RdBeamFormer::Run(evt::Event& event)
  {
    DEBUGLOG("RdBeamFormer::Run()");

    //check if there is a radio event at all
    if (!event.HasREvent()) {
      WARNING("No radio event found!");
      return eContinueLoop;
    }

    REvent& rEvent = event.GetREvent();    
    //PRINT ("Radio event found with " << rEvent.GetNumberOfStations() << " stations!");

    const Header& rHeader = rEvent.GetHeader(); 
    PRINT ("Header ID: " << rHeader.GetId() << " and timestamp: " << rHeader.GetTime());
 
    if (!event.HasRecShower() || !event.GetRecShower().HasRRecShower()) {
      WARNING("No RecShower or RRecShower found!");
      return eContinueLoop;
    }    

    //calculate number of not rejected stations to be used for beam-forming
    int n = 0;    
    for (REvent::CandidateStationIterator it = rEvent.CandidateStationsBegin(); it != rEvent.CandidateStationsEnd(); ++it) {
        Station& station = *it;
        if (station.IsSaturated()) //skip saturated stations
           continue;        
        ++n;
    }

    //cout<<"Number of candidate stations is: "<<n<<"\n";

    //check if there is enough stations to perform beam-forming (minimum number of stations can be discussed and changed of course...)
    if (n < 2) {
       WARNING ("Not enough stations found to perform beam-forming!");
       return eContinueLoop;
    }    
    
    /*
    This part contains an attempt to estimate the shower footprint based on our analitical models for the radio emission
    from the extensive air-showers. This means estimating the signal power for every antenna station position and rejecting 
    that stations which are below the estimated noise power for this position. All calculations are based on these documents: 
    http://dx.doi.org/10.1016/j.astropartphys.2014.05.001
    http://dx.doi.org/10.1088/1475-7516/2015/05/018
    GAP-2014-073

    evt::ShowerRRecData& rrec= event.GetRecShower().GetRRecShower();   
    const utl::Point corePos = rrec.GetCorePosition();
    //get the core position in UTM coordinates
    utl::UTMPoint corePosUTM (corePos, utl::ReferenceEllipsoid::eWGS84);
    coreCS = fwk::LocalCoordinateSystem::Create(corePos);
    //get x and y components of the core position in the core coordinate system and z component in UTM coordinate system
    double corePosX = corePos.GetX(coreCS)/m;
    double corePosY = corePos.GetY(coreCS)/m;
    double corePosZ = corePosUTM.GetHeight()/m;  

    //cout<<"corePosX= "<<corePosX<<" corePosY= "<<corePosY<<" corePosZ= "<<corePosZ<<endl;

    const utl::Triple polarCoords = rrec.GetAxis().GetSphericalCoordinates(coreCS);
    double theta = polarCoords.get<1>();
    double cosTheta = cos(theta);

    //cout<<"theta= "<<polarCoords.get<1>()/degree<<" cosTheta= "<<cosTheta<<endl;

    evt::ShowerSimData& MCShower = event.GetSimShower();
    double mcEnergy = MCShower.GetEnergy()/eV;
    double xMax = MCShower.GetGHParameters().GetXMax()/(g/cm2);
    double distanceToXMax = MCShower.GetDistanceOfShowerMaximum();     
    const Atmosphere& theAtm = det::Detector::GetInstance().GetAtmosphere();
    ProfileResult depthVsHeight = theAtm.EvaluateDepthVsHeight();
    //calculate the vertical column density of the whole atmosphere by converting it from the height
    double xMaxAtm = depthVsHeight.Y(corePosZ)/(g/cm2);
    //calculate the distance to the shower maximum from the observer in g/cm2 not taking into account the curvature of the Earth
    double distanceXMax = xMaxAtm/cosTheta - xMax;

    //cout<<"mcEnergy= "<<mcEnergy<<" xMax= "<<xMax<<" distanceToXMax= "<<distanceToXMax/km<<" xMaxAtm= "<<xMaxAtm<<" distanceXMax= "<<distanceXMax<<endl;
 
    //calculate the magnetic field vector at the site of the shower
    Vector showerAxis (rrec.GetAxis());
    Vector magneticField = fwk::MagneticFieldModel::GetMagneticFieldVector(corePos, rHeader.GetTime()); 
    utl::RadioGeometryUtilities transformation = utl::RadioGeometryUtilities(showerAxis, coreCS, magneticField);

    //loop over the not rejected stations, calculate the expected signal power and reject that stations which are below the expected noise power 
    for (REvent::CandidateStationIterator it = rEvent.CandidateStationsBegin(); it != rEvent.CandidateStationsEnd(); ++it) {
        Station& station = *it;
        if (station.IsSaturated())   // skip saturated stations
           continue;     
        StationRecData& sRec = station.GetRecData();
        const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
        const Point stationPos = rDetector.GetStation(*it).GetPosition();
        double stationPosX = stationPos.GetX(coreCS)/m;
        double stationPosY = stationPos.GetY(coreCS)/m;
        double stationPosZ = stationPos.GetZ(coreCS)/m;

        //transform staion positions to the shower plane
        double stationPosX_VxB = 0.;
        double stationPosY_VxB = 0.;
        double stationPosZ_VxB = 0.;
        transformation.GetVectorInShowerPlaneVxB(stationPosX_VxB, stationPosY_VxB, stationPosZ_VxB, stationPos);
        stationPosX_VxB = stationPosX_VxB/m;
        stationPosY_VxB = stationPosY_VxB/m;
        stationPosZ_VxB = stationPosZ_VxB/m;

        //cout<<"Station "<<it->GetId()<<" stationPosX= "<<stationPosX<<" stationPosY= "<<stationPosY<<" stationPosZ= "<<stationPosZ<<endl;
        //cout<<"Station "<<it->GetId()<<" stationPosX_VxB= "<<stationPosX_VxB<<" stationPosY_VxB= "<<stationPosY_VxB<<" stationPosZ_VxB= "<<stationPosZ_VxB<<endl;

        double c0 = 0.;
        double c1 = 0.;
        double c2 = 0.;
        if (theta/degree <= 55) {
           c0 = 0.41;           
        } 
        else if (theta/degree > 55 && theta/degree <= 60) {
           c0 = 0.71;
        } 
        else if (theta/degree > 60) {
           c0 = 0.71;
        }

        if (theta/degree <= 10) {
           c1 = 8.0;
           c2 = -21.2;
        } 
        else if (theta/degree > 10 && theta/degree <= 20) {
           c1 = 10.0;
           c2 = -23.1;
        } 
        else if (theta/degree > 20 && theta/degree <= 30) {
           c1 = 12.0;
           c2 = -25.5;
        }
        else if (theta/degree > 30 && theta/degree <= 40) {
           c1 = 20.0;
           c2 = -32.0;
        }
        else if (theta/degree > 40 && theta/degree <= 50) {
           c1 = 25.1;
           c2 = -34.5;
        }
        else if (theta/degree > 50 && theta/degree <= 60) {
           c1 = 27.3;
           c2 = -9.8;
        }
        else if (theta/degree > 60) {
           c1 = 27.3;
           c2 = -9.8;
        }
        
        long double function1 = pow(10,-52.8)*pow(mcEnergy,2);
        long double function2 = 37.7 - 0.006*distanceXMax + 5.5*pow(10,-4)*pow(distanceXMax,2) - 3.7*pow(10,-7)*pow(distanceXMax,3);
        long double function3 = 26.9 - 0.041*distanceXMax + 2*pow(10,-4)*pow(distanceXMax,2);

        long double totalPower = function1*(exp(-(pow((stationPosX_VxB - corePosX - c1),2) + pow((stationPosY_VxB - corePosY),2))/pow(function2,2)) - 
                                    c0*exp(-(pow((stationPosX_VxB - corePosX - c2),2) + pow((stationPosY_VxB - corePosY),2))/pow(function3,2)))*(joule/m2);
        double noisePower = sRec.GetParameter(revt::eNoiseEnergyDensityMag)/(eV/m2);
        
        //cout<<"function1= "<<function1<<" function2= "<<function2<<" "<<pow((stationPosX_VxB - corePosX - c1),2) + pow((stationPosY_VxB - corePosY),2)<<" function3= "<<function3<<" "<<pow((stationPosX_VxB - corePosX - c2),2) + pow((stationPosY_VxB - corePosY),2)<<" totalPower= "<<totalPower/(eV/m2)<<" noisePower= "<<noisePower<<endl;   
        
        if (totalPower < noisePower) {  
           station.SetRejected(revt::eManuallyRejected);
           --n;
        }
    }

    //cout<<"Number of station is: "<<n<<endl;

    //check if there is still enough stations to perform beam-forming
    if (n < 2) {
       WARNING ("Not enough stations found to perform beam-forming!");
       return eContinueLoop;
    }
    */

    std::vector<StationFFTDataContainer> stationData(n);
    double binning;
    int start, stop;
    double timeOffset = shiftTraces(stationData, n, binning, start, stop, event);

    //cout<<"Number of station is: "<<n<<endl;

    //check if there is still enough stations to perform beam-forming
    if (n < 2) {
       WARNING ("Not enough stations found to perform beam-forming!");
       return eContinueLoop;
    }

    TraceD cc = crosscorr (stationData, n, binning, start, stop);
    TraceD power = powertrace (stationData, n, binning, start, stop);

    evt::BeamPeak ccPeak = findPeak (cc, timeOffset);
    evt::BeamPeak pwPeak = findPeak (power, timeOffset);

    TraceD x = xtrace(cc, cc-ccPeak.GetOffset(), power-pwPeak.GetOffset());
    evt::BeamPeak xPeak = findPeak (x, timeOffset);     

    WriteASCII (cc, power, x, ccPeak, pwPeak, xPeak, event);   

    /*
    Uncomment this part if you want to write to the files maximum beam-forming values of the scanned region 
    and other respective information of the set of events you want to analyze. There are two options: 
    finding the maximum based on the absolut values of the beam-forming traces or alternatively on their SNRs

    //get the core position
    evt::ShowerRRecData& rrec= event.GetRecShower().GetRRecShower();   
    const utl::Point corePos = rrec.GetCorePosition();

    //calculate the magnetic field vector at the site of the shower
    Vector showerAxis (rrec.GetAxis());
    Vector magneticField = fwk::MagneticFieldModel::GetMagneticFieldVector(corePos, rHeader.GetTime()); 

    static double MCEnergy = MCShower.GetEnergy();
    static double distanceShowerMax = MCShower.GetDistanceOfShowerMaximum();
    static double XMax = MCShower.GetGHParameters().GetXMax();
    static double cosGeomagneticAngle = (showerAxis*magneticField)/(sqrt(showerAxis*showerAxis)*sqrt(magneticField*magneticField)); 
    static double sinGeomagneticAngle = sqrt(1 - cosGeomagneticAngle*cosGeomagneticAngle);
    //cout<<"MCEnergy= "<<MCEnergy<<" sinGeomagneticAngle= "<<sinGeomagneticAngle
        //<<" MCEnergy/sinGeomagneticAngle= "<<MCEnergy*sinGeomagneticAngle
        //<<" distanceShowerMax= "<<distanceShowerMax/km<<"\n";      

    static std::string cc_max_values = "cc_max_values.txt";
    static std::string pw_max_values = "pw_max_values.txt";
    static std::string x_max_values = "x_max_values.txt";

    static std::string cc_max_snr_values = "cc_max_snr_values.txt";
    static std::string pw_max_snr_values = "pw_max_snr_values.txt";
    static std::string x_max_snr_values = "x_max_snr_values.txt";

    ofstream fout3;
    fout3.open(cc_max_values.c_str(), ios::out|ios::app);

    ofstream fout4;
    fout4.open(pw_max_values.c_str(), ios::out|ios::app);

    ofstream fout5;
    fout5.open(x_max_values.c_str(), ios::out|ios::app);

    ofstream fout6;
    fout6.open(cc_max_snr_values.c_str(), ios::out|ios::app);

    ofstream fout7;
    fout7.open(pw_max_snr_values.c_str(), ios::out|ios::app);

    ofstream fout8;
    fout8.open(x_max_snr_values.c_str(), ios::out|ios::app);


    static std::vector<double> max_cc;
    static std::vector<double> max_pw;
    static std::vector<double> max_x;
    static std::vector<double> max_th;
    static std::vector<double> max_phi;
    static std::vector<double> max_r;
    static std::vector<double> max_rho;
    static std::vector<double> max_b;
    static std::vector<double> max_cc_snr;
    static std::vector<double> max_pw_snr;
    static std::vector<double> max_x_snr;

    if (waveModelConf == "eSpherical") {
       int eventId = event.GetREvent().GetHeader().GetId();
       static int previousId = eventId;
       if (eventId == previousId) {
          max_cc.push_back(ccPeak.GetHeight());
          max_pw.push_back(pwPeak.GetHeight());
          max_x.push_back(xPeak.GetHeight());
          max_cc_snr.push_back(ccPeak.GetOffset());
          max_pw_snr.push_back(pwPeak.GetOffset()); 
          max_x_snr.push_back(xPeak.GetOffset()); 
          max_th.push_back(polarCoords.get<1>()/degree);
          max_phi.push_back(polarCoords.get<2>()/degree);
          max_r.push_back(polarCoords.get<0>()/km);                  
       } 
       else {          
          //information about current event is written to the files when the analysis switches to the next event. 
          //To get the last event you are interested in recordered to the file you need to include one more event after it to your set of events.
          vector <double>::iterator cc_max_iter;
          cc_max_iter = std::max_element(max_cc.begin(), max_cc.end());
          int pos = cc_max_iter - max_cc.begin(); 

          vector <double>::iterator cc_max_snr_iter;
          cc_max_snr_iter = std::max_element(max_cc_snr.begin(), max_cc_snr.end());
          int pos_snr = cc_max_snr_iter - max_cc_snr.begin(); 
         
          fout3<<previousId<<" "<<max_cc[pos]<<" "<<max_cc_snr[pos]<<" "<<max_th[pos]<<" "<<max_phi[pos]<<" "<<MCEnergy<<" "
               <<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "<<XMax/(g/cm2)<<" "<<max_r[pos]<<"\n";
          fout4<<previousId<<" "<<max_pw[pos]<<" "<<max_pw_snr[pos]<<" "<<max_th[pos]<<" "<<max_phi[pos]<<" "<<MCEnergy<<" "
               <<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "<<XMax/(g/cm2)<<" "<<max_r[pos]<<"\n";
          fout5<<previousId<<" "<<max_x[pos]<<" "<<max_x_snr[pos]<<" "<<max_th[pos]<<" "<<max_phi[pos]<<" "<<MCEnergy<<" "
               <<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "<<XMax/(g/cm2)<<" "<<max_r[pos]<<"\n";

          fout6<<previousId<<" "<<max_cc[pos_snr]<<" "<<max_cc_snr[pos_snr]<<" "<<max_th[pos_snr]<<" "
               <<max_phi[pos_snr]<<" "<<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "
               <<XMax/(g/cm2)<<" "<<max_r[pos_snr]<<"\n";
          fout7<<previousId<<" "<<max_pw[pos_snr]<<" "<<max_pw_snr[pos_snr]<<" "<<max_th[pos_snr]<<" "
               <<max_phi[pos_snr]<<" "<<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "
               <<XMax/(g/cm2)<<" "<<max_r[pos_snr]<<"\n";
          fout8<<previousId<<" "<<max_x[pos_snr]<<" "<<max_x_snr[pos_snr]<<" "<<max_th[pos_snr]<<" "
               <<max_phi[pos_snr]<<" "<<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "
               <<XMax/(g/cm2)<<" "<<max_r[pos_snr]<<"\n";

          previousId = eventId;     
          MCEnergy = MCShower.GetEnergy();
          distanceShowerMax = MCShower.GetDistanceOfShowerMaximum(); 
          XMax = MCShower.GetGHParameters().GetXMax();                   
          cosGeomagneticAngle = (showerAxis*magneticField)/(sqrt(showerAxis*showerAxis)*sqrt(magneticField*magneticField)); 
          sinGeomagneticAngle = sqrt(1 - cosGeomagneticAngle*cosGeomagneticAngle);          

          max_cc.clear();
          max_pw.clear();
          max_x.clear();
          max_cc_snr.clear();
          max_pw_snr.clear();
          max_x_snr.clear();
          max_th.clear();
          max_phi.clear();
          max_r.clear();          

          max_cc.push_back(ccPeak.GetHeight());
          max_pw.push_back(pwPeak.GetHeight());
          max_x.push_back(xPeak.GetHeight());
          max_cc_snr.push_back(ccPeak.GetOffset());
          max_pw_snr.push_back(pwPeak.GetOffset()); 
          max_x_snr.push_back(xPeak.GetOffset()); 
          max_th.push_back(polarCoords.get<1>()/degree);
          max_phi.push_back(polarCoords.get<2>()/degree);
          max_r.push_back(polarCoords.get<0>()/km);
       }
    }
    else if (waveModelConf == "eConical") {
       int eventId = event.GetREvent().GetHeader().GetId();
       static int previousId = eventId;
       if (eventId == previousId) {
          max_cc.push_back(ccPeak.GetHeight());
          max_pw.push_back(pwPeak.GetHeight());
          max_x.push_back(xPeak.GetHeight());
          max_cc_snr.push_back(ccPeak.GetOffset());
          max_pw_snr.push_back(pwPeak.GetOffset()); 
          max_x_snr.push_back(xPeak.GetOffset()); 
          max_th.push_back(polarCoords.get<1>()/degree);
          max_phi.push_back(polarCoords.get<2>()/degree);
          max_rho.push_back(rrec.GetParameter(eRho)/degree);
       } 
       else {
          //information about current event is written to the files when the analysis switches to the next event. 
          //To get the last event you are interested in recordered to the file you need to include one more event after it to your set of events.          
          vector <double>::iterator cc_max_iter;
          cc_max_iter = std::max_element(max_cc.begin(), max_cc.end());
          int pos = cc_max_iter - max_cc.begin();     

          vector <double>::iterator cc_max_snr_iter;
          cc_max_snr_iter = std::max_element(max_cc_snr.begin(), max_cc_snr.end());
          int pos_snr = cc_max_snr_iter - max_cc_snr.begin();      
         
          fout3<<previousId<<" "<<max_cc[pos]<<" "<<max_cc_snr[pos]<<" "<<max_th[pos]<<" "<<max_phi[pos]<<" "<<MCEnergy<<" "
               <<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "<<XMax/(g/cm2)<<" "<<max_rho[pos]<<"\n";
          fout4<<previousId<<" "<<max_pw[pos]<<" "<<max_pw_snr[pos]<<" "<<max_th[pos]<<" "<<max_phi[pos]<<" "<<MCEnergy<<" "
               <<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "<<XMax/(g/cm2)<<" "<<max_rho[pos]<<"\n";
          fout5<<previousId<<" "<<max_x[pos]<<" "<<max_x_snr[pos]<<" "<<max_th[pos]<<" "<<max_phi[pos]<<" "<<MCEnergy<<" "
               <<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "<<XMax/(g/cm2)<<" "<<max_rho[pos]<<"\n";

          fout6<<previousId<<" "<<max_cc[pos_snr]<<" "<<max_cc_snr[pos_snr]<<" "<<max_th[pos_snr]<<" "
               <<max_phi[pos_snr]<<" "<<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "
               <<XMax/(g/cm2)<<" "<<max_rho[pos_snr]<<"\n";
          fout7<<previousId<<" "<<max_pw[pos_snr]<<" "<<max_pw_snr[pos_snr]<<" "<<max_th[pos_snr]<<" "
               <<max_phi[pos_snr]<<" "<<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "
               <<XMax/(g/cm2)<<" "<<max_rho[pos_snr]<<"\n";
          fout8<<previousId<<" "<<max_x[pos_snr]<<" "<<max_x_snr[pos_snr]<<" "<<max_th[pos_snr]<<" "
               <<max_phi[pos_snr]<<" "<<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "
               <<XMax/(g/cm2)<<" "<<max_rho[pos_snr]<<"\n";

          previousId = eventId;     
          MCEnergy = MCShower.GetEnergy();
          distanceShowerMax = MCShower.GetDistanceOfShowerMaximum();  
          XMax = MCShower.GetGHParameters().GetXMax();                  
          cosGeomagneticAngle = (showerAxis*magneticField)/(sqrt(showerAxis*showerAxis)*sqrt(magneticField*magneticField)); 
          sinGeomagneticAngle = sqrt(1 - cosGeomagneticAngle*cosGeomagneticAngle);  

          max_cc.clear();
          max_pw.clear();
          max_x.clear();
          max_cc_snr.clear();
          max_pw_snr.clear();
          max_x_snr.clear();
          max_th.clear();
          max_phi.clear();
          max_rho.clear();

          max_cc.push_back(ccPeak.GetHeight());
          max_pw.push_back(pwPeak.GetHeight());
          max_x.push_back(xPeak.GetHeight());
          max_cc_snr.push_back(ccPeak.GetOffset());
          max_pw_snr.push_back(pwPeak.GetOffset()); 
          max_x_snr.push_back(xPeak.GetOffset()); 
          max_th.push_back(polarCoords.get<1>()/degree);
          max_phi.push_back(polarCoords.get<2>()/degree);
          max_rho.push_back(rrec.GetParameter(eRho)/degree);
       }
    }
    else if (waveModelConf == "eHyperbolic") {
       int eventId = event.GetREvent().GetHeader().GetId();
       static int previousId = eventId;
       if (eventId == previousId) {
          max_cc.push_back(ccPeak.GetHeight());
          max_pw.push_back(pwPeak.GetHeight());
          max_x.push_back(xPeak.GetHeight());
          max_cc_snr.push_back(ccPeak.GetOffset());
          max_pw_snr.push_back(pwPeak.GetOffset()); 
          max_x_snr.push_back(xPeak.GetOffset()); 
          max_th.push_back(polarCoords.get<1>()/degree);
          max_phi.push_back(polarCoords.get<2>()/degree);
          max_rho.push_back(rrec.GetParameter(eRho)/degree);
          max_b.push_back(rrec.GetParameter(eB)/ns);
       } 
       else {
          //information about current event is written to the files when the analysis switches to the next event. 
          //To get the last event you are interested in recordered to the file you need to include one more event after it to your set of events.
          vector <double>::iterator cc_max_iter;
          cc_max_iter = std::max_element(max_cc.begin(), max_cc.end());
          int pos = cc_max_iter - max_cc.begin();     

          vector <double>::iterator cc_max_snr_iter;
          cc_max_snr_iter = std::max_element(max_cc_snr.begin(), max_cc_snr.end());
          int pos_snr = cc_max_snr_iter - max_cc_snr.begin();      
         
          fout3<<previousId<<" "<<max_cc[pos]<<" "<<max_cc_snr[pos]<<" "<<max_th[pos]<<" "<<max_phi[pos]<<" "
               <<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "<<XMax/(g/cm2)<<" "
               <<max_rho[pos]<<" "<<max_b[pos]<<"\n";
          fout4<<previousId<<" "<<max_pw[pos]<<" "<<max_pw_snr[pos]<<" "<<max_th[pos]<<" "<<max_phi[pos]<<" "
               <<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "<<XMax/(g/cm2)<<" "
               <<max_rho[pos]<<" "<<max_b[pos]<<"\n";
          fout5<<previousId<<" "<<max_x[pos]<<" "<<max_x_snr[pos]<<" "<<max_th[pos]<<" "<<max_phi[pos]<<" "
               <<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "<<XMax/(g/cm2)<<" "
               <<max_rho[pos]<<" "<<max_b[pos]<<"\n";

          fout6<<previousId<<" "<<max_cc[pos_snr]<<" "<<max_cc_snr[pos_snr]<<" "<<max_th[pos_snr]<<" "
               <<max_phi[pos_snr]<<" "<<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "
               <<XMax/(g/cm2)<<" "<<max_rho[pos_snr]<<" "<<max_b[pos_snr]<<"\n";
          fout7<<previousId<<" "<<max_pw[pos_snr]<<" "<<max_pw_snr[pos_snr]<<" "<<max_th[pos_snr]<<" "
               <<max_phi[pos_snr]<<" "<<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "
               <<XMax/(g/cm2)<<" "<<max_rho[pos_snr]<<" "<<max_b[pos_snr]<<"\n";
          fout8<<previousId<<" "<<max_x[pos_snr]<<" "<<max_x_snr[pos_snr]<<" "<<max_th[pos_snr]<<" "
               <<max_phi[pos_snr]<<" "<<MCEnergy<<" "<<MCEnergy*sinGeomagneticAngle<<" "<<distanceShowerMax/km<<" "
               <<XMax/(g/cm2)<<" "<<max_rho[pos_snr]<<" "<<max_b[pos_snr]<<"\n";

          previousId = eventId;     
          MCEnergy = MCShower.GetEnergy();
          distanceShowerMax = MCShower.GetDistanceOfShowerMaximum();
          XMax = MCShower.GetGHParameters().GetXMax();                    
          cosGeomagneticAngle = (showerAxis*magneticField)/(sqrt(showerAxis*showerAxis)*sqrt(magneticField*magneticField)); 
          sinGeomagneticAngle = sqrt(1 - cosGeomagneticAngle*cosGeomagneticAngle);  

          max_cc.clear();
          max_pw.clear();
          max_x.clear();
          max_cc_snr.clear();
          max_pw_snr.clear();
          max_x_snr.clear();
          max_th.clear();
          max_phi.clear();
          max_rho.clear();
          max_b.clear();

          max_cc.push_back(ccPeak.GetHeight());
          max_pw.push_back(pwPeak.GetHeight());
          max_x.push_back(xPeak.GetHeight());
          max_cc_snr.push_back(ccPeak.GetOffset());
          max_pw_snr.push_back(pwPeak.GetOffset()); 
          max_x_snr.push_back(xPeak.GetOffset()); 
          max_th.push_back(polarCoords.get<1>()/degree);
          max_phi.push_back(polarCoords.get<2>()/degree);
          max_rho.push_back(rrec.GetParameter(eRho)/degree);
          max_b.push_back(rrec.GetParameter(eB)/ns);
       }
    }

    fout3.close();
    fout4.close();
    fout5.close();
    fout6.close();
    fout7.close();
    fout8.close();
    */  

    /*
    This functionality is not finished yet...

    if (rrec.HasCurrentBeamMap()) {
      DEBUGLOG ("filling BeamQuantities into map");

      evt::BeamQuantities quants (rrec.GetAxis().GetTheta(coreCS) ,
                                  rrec.GetAxis().GetPhi(coreCS) ,
                                  waveModel->curvature(rrec.GetAxis(),coreCS),
                                  waveModel->coneAngle(rrec.GetAxis(),coreCS),
                                  ccPeak, pwPeak, xPeak);

      evt::BeamMap map = rrec.GetCurrentBeamMap();
      map.AddQuantities(quants);
    }
    if (!rrec.HasBeamResultQuantities()
        || (rrec.GetBeamResultQuantities().GetCCPeak().GetHeight() / rrec.GetBeamResultQuantities().GetCCPeak().GetRMS()
            < ccPeak.GetHeight() / ccPeak.GetRMS() )) {
      evt::BeamQuantities resultQuants (showerAxis.GetTheta(coreCS) ,
                                        showerAxis.GetPhi(coreCS) ,
                                        waveModel->curvature(rrec.GetAxis(),coreCS),
                                        waveModel->coneAngle(rrec.GetAxis(),coreCS),
                                        ccPeak, pwPeak, xPeak);
      rrec.SetBeamResultQuantities(resultQuants);
    }
    */

   return eSuccess;
  }

  VModule::ResultFlag
  RdBeamFormer::Finish()
  {
    INFO("RdBeamFormer::Finish()");

    delete (waveModel);

    return eSuccess;
  }

  void RdBeamFormer::WriteASCII (const TraceD& cc, const TraceD& power, const TraceD& x,
                                    const evt::BeamPeak& ccPeak, const evt::BeamPeak& pwPeak, const evt::BeamPeak& xPeak,
                                    const evt::Event& event) const {
    fNumCalls++;
      
    double eps = 1e-14; //desired precision for variables comparision

    static std::string outfile_name = outfile + "_";
 
    //write to a new file if analysis switches to the next event 
    const evt::ShowerRRecData& rrec= event.GetRecShower().GetRRecShower();
    int eventId = event.GetREvent().GetHeader().GetId();
    std::string stringId = SSTR(eventId);
    //cout<<"eventId= "<<eventId<<"\n";
    static int oldId = -1;
    //cout<<"oldId= "<<oldId<<"\n";
    if (eventId != oldId) {
       //cout<<"eventId= "<<eventId<<"\n";    
       oldId = eventId;
       //cout<<"oldId= "<<oldId<<"\n";
       std::string buffer = SSTR(oldId) + "_";          
       outfile_name.replace(outfile.size()+1,buffer.size(),buffer);    
    }     
    
    //write to a new file if the respective wavefront parameter switches to the next step
    if (waveModelConf == "eSpherical") {       
       const utl::Triple polarCoords = rrec.GetAxis().GetSphericalCoordinates(coreCS);
       double r = polarCoords.get<0>()/km;
       static double oldR = -1;
       //cout<<"oldR= "<<oldR<<" oldR-r= "<<oldR - r<<"\n";       
       if (abs(r - oldR) >= eps) {
          //cout<<"r= "<<r<<"\n";
          oldR = r;
          //cout<<"oldR= "<<oldR<<"\n";          
          std::string buffer = SSTR(oldR);          
          outfile_name.replace(outfile.size()+stringId.size()+2,buffer.size()+10,buffer);         
       }
    }
    else if (waveModelConf == "eConical") {
       double Rho = rrec.GetParameter(eRho)/degree;
       //cout<<"eRho= "<<Rho<<"\n";
       static double oldRho = -1;   
       //cout<<"oldRho= "<<oldRho<<"\n"; 
       if (abs(Rho - oldRho) >= eps) {
          //cout<<"eRho= "<<Rho<<"\n";
          oldRho = Rho;
          //cout<<"oldRho= "<<oldRho<<"\n";                    
          std::string buffer = SSTR(oldRho);                           
          outfile_name.replace(outfile.size()+stringId.size()+2,buffer.size()+10,buffer);          
       }
    }
    else if (waveModelConf == "eHyperbolic") {
       double Rho = rrec.GetParameter(eRho)/degree;
       double b = rrec.GetParameter(eB)/ns;
       //double b = 3*ns;
       //cout<<"eRho= "<<Rho<<"\n";
       static double old_Rho = -1; 
       static double old_b = -1;  
       //cout<<"oldRho= "<<oldRho<<"\n"; 
       if (abs(b - old_b) >= eps) {   //  abs(Rho - old_Rho) >= eps || 
          //cout<<"eRho= "<<Rho<<"\n";
          old_Rho = Rho;
          old_b = b;
          //cout<<"oldRho= "<<oldRho<<"\n";                    
          std::string buffer1 = SSTR(old_Rho);       
          std::string buffer2 = SSTR(old_b);  
          std::string buffer = buffer1 + "_" + buffer2;                                             
          outfile_name.replace(outfile.size()+stringId.size()+2,buffer.size()+10,buffer); 
       }
    }   
 
    /*
    Uncomment this part if you want to write to the files the whole beam-forming traces as well

    WriteTrace(cc, ccPeak, "cc_"+outfile_name+".txt", rrec);
    WriteTrace(power, pwPeak, "pw_"+outfile_name+".txt", rrec);
    WriteTrace(x, xPeak, "x_"+outfile_name+".txt", rrec);
    */

    WritePeak(ccPeak, "cc_max_"+outfile_name+".txt", rrec);
    WritePeak(pwPeak, "pw_max_"+outfile_name+".txt", rrec);
    WritePeak(xPeak, "x_max_"+outfile_name+".txt", rrec);
  }
            
  void RdBeamFormer::WriteTrace(const TraceD& trace, const evt::BeamPeak &Peak,
                                   const string& filename, const evt::ShowerRRecData& rrec) const {
    ofstream fout;    
    fout.open(filename.c_str(),ios::out|ios::app);
    
    if (gridtype == polar) {
       if (waveModelConf == "ePlane" || waveModelConf == "eSpherical") { 
          const utl::Triple polarCoords = rrec.GetAxis().GetSphericalCoordinates(coreCS);     
          for (StationTimeSeries::SizeType i = 0; i < trace.GetSize(); ++i) {           
              fout<<polarCoords.get<0>()/km<<" "<<polarCoords.get<1>()/degree<<" "<<polarCoords.get<2>()/degree<<" "
                  <<trace[i]<<" "<<trace.GetBinning()*(i+startbin+1)/nanosecond<<" "
                  <<Peak.GetHeight()<<" "<<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";                 
          }
          fout<<endl;
          fout<<endl;
       }
       else if (waveModelConf == "eConical") {
          const utl::Triple polarCoords = rrec.GetAxis().GetSphericalCoordinates(coreCS);     
          for (StationTimeSeries::SizeType i = 0; i < trace.GetSize(); ++i) {           
              fout<<rrec.GetParameter(eRho)/degree<<" "<<polarCoords.get<0>()/km<<" "<<polarCoords.get<1>()/degree<<" "
                  <<polarCoords.get<2>()/degree<<" "<<trace[i]<<" "<<trace.GetBinning()*(i+startbin+1)/nanosecond<<" "
                  <<Peak.GetHeight()<<" "<<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";                 
          }
          fout<<endl;
          fout<<endl;
       }
       else if (waveModelConf == "eHyperbolic") {
          const utl::Triple polarCoords = rrec.GetAxis().GetSphericalCoordinates(coreCS);     
          for (StationTimeSeries::SizeType i = 0; i < trace.GetSize(); ++i) {           
              fout<<rrec.GetParameter(eB)/ns<<" "<<rrec.GetParameter(eRho)/degree<<" "
                  <<polarCoords.get<0>()/km<<" "<<polarCoords.get<1>()/degree<<" "<<polarCoords.get<2>()/degree<<" "
                  <<trace[i]<<" "<<trace.GetBinning()*(i+startbin+1)/nanosecond<<" "
                  <<Peak.GetHeight()<<" "<<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";                 
          }
          fout<<endl;
          fout<<endl;
       }
    }
    else if (gridtype == carthesian) {
       if (waveModelConf == "ePlane" || waveModelConf == "eSpherical") { 
          const utl::Triple xyz = rrec.GetAxis().GetCoordinates(coreCS);
          for (StationTimeSeries::SizeType i = 0; i < trace.GetSize(); ++i) {           
              fout<<xyz.get<0>()/km<<" "<<xyz.get<1>()/km<<" "<<xyz.get<2>()/km<<" "
                  <<trace[i]<<" "<<trace.GetBinning()*(i+startbin+1)/nanosecond<<" "
                  <<Peak.GetHeight()<<" "<<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";                 
          }
          fout<<endl;
          fout<<endl;
       }
       else if (waveModelConf == "eConical") {
          const utl::Triple xyz = rrec.GetAxis().GetCoordinates(coreCS);
          for (StationTimeSeries::SizeType i = 0; i < trace.GetSize(); ++i) {           
              fout<<rrec.GetParameter(eRho)/degree<<" "<<xyz.get<0>()/km<<" "<<xyz.get<1>()/km<<" "<<xyz.get<2>()/km<<" "
                  <<trace[i]<<" "<<trace.GetBinning()*(i+startbin+1)/nanosecond<<" "
                  <<Peak.GetHeight()<<" "<<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";                 
          }
          fout<<endl;
          fout<<endl;
       }
       else if (waveModelConf == "eHyperbolic") {
          const utl::Triple xyz = rrec.GetAxis().GetCoordinates(coreCS);
          for (StationTimeSeries::SizeType i = 0; i < trace.GetSize(); ++i) {           
              fout<<rrec.GetParameter(eB)/ns<<" "<<rrec.GetParameter(eRho)/degree<<" "
                  <<xyz.get<0>()/km<<" "<<xyz.get<1>()/km<<" "<<xyz.get<2>()/km<<" "
                  <<trace[i]<<" "<<trace.GetBinning()*(i+startbin+1)/nanosecond<<" "
                  <<Peak.GetHeight()<<" "<<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";                 
          }
          fout<<endl;
          fout<<endl;
       }
    }
    
    fout.close();
  }

  void RdBeamFormer::WritePeak(const evt::BeamPeak& Peak, const string& filename, const evt::ShowerRRecData& rrec) const {
    int stepsPhi = 6; //have the same purpose as fNumcalls. You have to change it manually accordingly to the number of steps in phi set in .xml.
                       //It always equals stepsPhi+1 if stepsPhi is not equal to 360, otherwise it should be 360.
    ofstream fout;    
    fout.open(filename.c_str(),ios::out|ios::app);

    if (gridtype == polar) {
       if (waveModelConf == "ePlane" || waveModelConf == "eSpherical") { 
          const utl::Triple polarCoords = rrec.GetAxis().GetSphericalCoordinates(coreCS);
          fout<<polarCoords.get<0>()/km<<" "<<polarCoords.get<1>()/degree<<" "<<polarCoords.get<2>()/degree<<" "
              <<Peak.GetHeight()<<" "<<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";
          if ((fNumCalls % stepsPhi) == 0) 
             fout<<"\n";
       }
       else if (waveModelConf == "eConical") {
          const utl::Triple polarCoords = rrec.GetAxis().GetSphericalCoordinates(coreCS);
          fout<<rrec.GetParameter(eRho)/degree<<" "<<polarCoords.get<0>()/km<<" "<<polarCoords.get<1>()/degree<<" "
              <<polarCoords.get<2>()/degree<<" "<<Peak.GetHeight()<<" "<<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";
          if ((fNumCalls % stepsPhi) == 0) 
             fout<<"\n";
       }
       else if (waveModelConf == "eHyperbolic") {
          const utl::Triple polarCoords = rrec.GetAxis().GetSphericalCoordinates(coreCS);
          fout<<rrec.GetParameter(eB)/ns<<" "<<rrec.GetParameter(eRho)/degree<<" "<<polarCoords.get<0>()/km<<" "
              <<polarCoords.get<1>()/degree<<" "<<polarCoords.get<2>()/degree<<" "<<Peak.GetHeight()<<" "
              <<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";
          if ((fNumCalls % stepsPhi) == 0) 
             fout<<"\n";
       }
    }
    else if (gridtype == carthesian) {
       if (waveModelConf == "ePlane" || waveModelConf == "eSpherical") { 
          const utl::Triple xyz = rrec.GetAxis().GetCoordinates(coreCS);
          fout<<xyz.get<0>()/km<<" "<<xyz.get<1>()/km<<" "<<xyz.get<2>()/km<<" "<<Peak.GetHeight()<<" "
              <<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";
          if ((fNumCalls % stepsPhi) == 0) 
             fout<<"\n";
       }
       else if (waveModelConf == "eConical") {
          const utl::Triple xyz = rrec.GetAxis().GetCoordinates(coreCS);
          fout<<rrec.GetParameter(eRho)/degree<<" "<<xyz.get<0>()/km<<" "<<xyz.get<1>()/km<<" "<<xyz.get<2>()/km<<" "
              <<Peak.GetHeight()<<" "<<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";
          if ((fNumCalls % stepsPhi) == 0) 
             fout<<"\n";
       }
       else if (waveModelConf == "eHyperbolic") {
          const utl::Triple xyz = rrec.GetAxis().GetCoordinates(coreCS);
          fout<<rrec.GetParameter(eB)/ns<<" "<<rrec.GetParameter(eRho)/degree<<" "<<xyz.get<0>()/km<<" "<<xyz.get<1>()/km<<" "
              <<xyz.get<2>()/km<<" "<<Peak.GetHeight()<<" "<<Peak.GetTime()<<" "<<Peak.GetOffset()<<"\n";
          if ((fNumCalls % stepsPhi) == 0) 
             fout<<"\n";
       }
    }

    fout.close();
  }
}

#undef PRINT