tabulatedata.cc

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <string>
#include <TCanvas.h>
#include <TH2F.h>
#include <TMath.h>
#include <TFile.h>
#include <TGraph.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TColor.h>
#include <TTimeStamp.h>
#define AUGER_RADIO_ENABLED
#include "RecEventFile.h"
#include "RecEvent.h"
#include "RdEvent.h"
#include "DetectorGeometry.h"
#include "SDEvent.h"
#include "DetectorGeometry.h"
#include "SdRecStation.h"
#include "ShowerRRecDataQuantities.h"
#include "StationRRecDataQuantities.h"
#include "PolarHist.h"

//Added
#include <TVector3.h>
#include <TProfile.h>

#include <RdTrace.h>

const int    kGPS_UTC_Offset =  315957587; // To be updated

using namespace std;

//Nice Colours
void set_plot_style() {
    const Int_t NRGBs = 5;
    const Int_t NCont = 255;
    Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
    Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
    Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
    Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
    TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
    gStyle->SetNumberContours(NCont);
}

struct antenna {
    int id;
    float dist_axis;
};

bool compDist (antenna i, antenna j) {
    return (i.dist_axis<j.dist_axis);
};

int main() {
    set_plot_style();

    vector<string> v_filename;
    ifstream fileList ("./dataADSTlist.txt");
    char filename[10000] ;
    while (fileList >> filename) {
        v_filename.push_back(filename);
    }

    const TVector3 vmag(-884.62,-19699.87,+14154.43);
    RecEventFile *file = new RecEventFile(v_filename);
    RecEvent* theevent= new RecEvent();
    int nevents=file->GetNEvents();
    cout << nevents << endl;
    file->SetBuffers(&theevent);
    DetectorGeometry *geo =new DetectorGeometry();
    file->ReadDetectorGeometry(*geo);

    ofstream File_data("Amplitude_data.txt");

    float distance=0;
    float distance_before=0;
    int distance_stadtion_id=0;

    //Loop over events
    while(file->ReadNextEvent()==RecEventFile::eSuccess) {
        const RdEvent& theradio= theevent->GetRdEvent();
        const RdRecShower& rshow=theradio.GetRdRecShower();

        if (!theradio.HasRdStations()) continue ;
        const std::vector<RdRecStation>& vecRstat=theradio.GetRdStationVector();
        const SDEvent& sdev=theevent->GetSDEvent();
        const SdRecShower&  sshow=  sdev.GetSdRecShower();

        const TVector3& score=sshow.GetCoreSiteCS();
        const TVector3& rcore=rshow.GetCoreSiteCS();
        const TVector3& saxis=sshow.GetAxisSiteCS();

        double rtheta=TMath::RadToDeg()*rshow.GetZenith();
        double rphi=TMath::RadToDeg()*rshow.GetAzimuth();
        double stheta=TMath::RadToDeg()*sshow.GetZenith();
        double sphi=TMath::RadToDeg()*sshow.GetAzimuth();
        TVector3 raxis=rshow.GetAxisSiteCS();
        float scal=saxis.Dot(raxis);
        float snrm=saxis.Mag();

        scal=saxis.Dot(vmag);
        float mnrm=vmag.Mag();

        //SD quality check
        const int minRecLevel=3; // has reconstructed energy
        const int minCandidateStations=4;
        const int T4Trigger=2; // has T4 and SD triggered independently of FD
        const int T5Trigger=3; //has T5 Posterior, which is the T5 HAS in this case

        const float minZenithSD=62;
        const float maxZenithSD=80;

        if (sdev.GetRecLevel()<minRecLevel) continue;// alle events
        if (sdev.GetNumberOfCandidates()<minCandidateStations) continue; //alle events
        if (sdev.GetT4Trigger()<T4Trigger) continue;
        if (sdev.GetT5Trigger()<T5Trigger) continue;
        if (stheta>maxZenithSD) continue;
        if (stheta<minZenithSD) continue;

        //Calculate closest station

        antenna arr[160] ;
        for (int i=0; i<160; i++) {
            arr[i].dist_axis=100000000;
            arr[i].id=101+i;
        }
        for (int i=0; i<160; i++) {
            if(theradio.HasRdStation(101+i)) {
                arr[i].dist_axis=geo->GetRdStationAxisDistance(101+i,raxis,rcore);
            }
        }
        sort(arr, arr + 160, compDist);

        int how_many_stations_participate=0;
        for (int i=0; i<160; i++) {
            if(arr[i].dist_axis < 10000000) {
                how_many_stations_participate++;
            }
        }

        distance_stadtion_id=arr[0].id;
        distance_before=arr[0].dist_axis;

        //Average Amplitude
        float average_amp=0;
        float max_amp=0;

        double minDist2 = -1;
        unsigned int closestRdStationId = 0;

        int contained=-1000;

        int dutch_stations[46]= {102, 104, 105, 108, 109, 112, 134, 135, 136, 143, 144, 149, 158, 159, 160, 169, 170, 183, 184, 185, 188, 192, 193, 194, 200, 201, 202, 203, 210, 211, 212, 213, 222, 223, 224, 225, 226, 235, 236, 237, 244, 245, 246, 251, 252, 257};

        for (vector<RdRecStation>::const_iterator stiter=vecRstat.begin(); stiter!=vecRstat.end(); ++stiter) { //loop stations

            for(int i=0; i<46; i++) {
                if(((stiter->GetId())-100)==dutch_stations[i]) {
                    continue;
                }
            }

            TVector3 rd_station_coordinate=geo->GetRdStationPosition(stiter->GetId());
            double dist2 = (rd_station_coordinate.X()-score.X())*(rd_station_coordinate.X()-score.X()) + (rd_station_coordinate.Y()-score.Y())*(rd_station_coordinate.Y()-score.Y());

            if (closestRdStationId == 0 || sqrt(dist2) < minDist2) {
                closestRdStationId = stiter->GetId();
                minDist2 = sqrt(dist2);
            }

            double distance_to_closest_rd_station=minDist2;

            if(distance_to_closest_rd_station<= 500.0) {
                contained = 1;
            } else {
                contained = 0;
            }
        }

        for (vector<RdRecStation>::const_iterator stiter=vecRstat.begin(); stiter!=vecRstat.end(); ++stiter) { //loop stations

            if (stiter->IsRejected()) {
                continue; // Do not draw stations which are rejected
            }

            float axis_distance_to_sd=geo->GetStationAxisDistance(geo->GetRdStationPosition(stiter->GetId()), saxis , score);

            if(stiter->GetParameter(revt::eSignalToNoise)>10.0) {
                File_data << theradio.GetRdRunNumber() << " "
                          << theradio.GetRdEventId() << " "
                          << stiter->GetId() << " "
                          << stiter->GetParameter(revt::ePeakAmplitudeMag)*1e6 << " "
                          << stheta << " " << TMath::Log10(sshow.GetEnergy()) << " "
                          << contained << " " << axis_distance_to_sd << endl;
            }
        }
    }//loop over events
}