events_tool.cpp

Go to the documentation of this file.

// events_tool.cpp : Defines the entry point for the console application.
//


#include "stdafx.h"
#include "../newicp/icputils.h"
#include "../newicp/event_store.h"
#include "../newicp/TCPEventStreamDefs.h"
#include "../newicp/isisraw.h"
#include "../nexus/NeXusFile.hpp"

class EventsToolApp : public Poco::Util::Application, public ISIS::Base
{
public:
        EventsToolApp() : Poco::Util::Application() { }
    virtual int main(const std::vector<std::string>& args);
protected:
        void initialize(Poco::Util::Application& self)
        {
                loadConfiguration();
                Application::initialize(self);
                setLoggerName("EventsToolApp");
        }

        void uninitialize()
        {
                Poco::Util::Application::uninitialize();
        }

        void defineOptions(Poco::Util::OptionSet& options);
        void displayHelp();
        void handleHelp(const std::string& name, const std::string& value);
        int convertRawToEvents(const std::string& raw_file, int run_number);
        int processNeXusFile(const std::string& nexus_file);
        int liveData(const std::string& host);
};


int _tmain(int argc, _TCHAR* argv[])
{
        EventsToolApp app;
        app.init(argc, argv);
        int ret = app.run();
        return ret;
}

void EventsToolApp::defineOptions(Poco::Util::OptionSet& options)
{
        Application::defineOptions(options);

        options.addOption(Poco::Util::Option("help", "h", "display help information")
                .required(false)
                .repeatable(false)
                .callback(Poco::Util::OptionCallback<EventsToolApp>(this, &EventsToolApp::handleHelp))
        );

        options.addOption(Poco::Util::Option("run-number", "r", "specify run number")
                .required(false)
                .repeatable(false)
                .argument("value")
                .binding("events_tool.run_number")
        );

        options.addOption(Poco::Util::Option("file", "f", "specify raw file")
                .required(false)
                .repeatable(false)
                .argument("file")
                .binding("events_tool.raw_file")
        );

        options.addOption(Poco::Util::Option("event-dir", "d", "specify directory for event files")
                .required(false)
                .repeatable(false)
                .argument("value")
                .binding("events_tool.event_dir")
        );

        options.addOption(Poco::Util::Option("nxfile", "n", "specify nexus file")
                .required(false)
                .repeatable(false)
                .argument("file")
                .binding("events_tool.nexus_file")
        );

        options.addOption(Poco::Util::Option("live", "l", "live data host")
                .required(false)
                .repeatable(false)
                .argument("value")
                .binding("events_tool.live_host")
        );
        
//      options.addOption(Option("config-file", "f", "load configuration data from a file")
//              .required(false)
//              .repeatable(true)
//              .argument("file")
//              .callback(OptionCallback<InstCtrlApp>(this, &EventsToolApp::handleConfig)));

//      options.addOption(Option("bind", "b", "bind option value to test.property")
//              .required(false)
//              .argument("value")
//              .validator(new IntValidator(0, 100))
//              .binding("test.pro
}

void EventsToolApp::displayHelp()
{
        Poco::Util::HelpFormatter helpFormatter(options());
        helpFormatter.setCommand(commandName());
        helpFormatter.setUsage("OPTIONS");
        helpFormatter.setHeader(
                "A sample application that demonstrates some of the features "
                "of the Poco::Util::Application class.");
        helpFormatter.format(std::cout);
}

void EventsToolApp::handleHelp(const std::string& name, const std::string& value)
{
//      _helpRequested = true;
        displayHelp();
        stopOptionsProcessing();
}

static uint64_t all_events_c = 0;
static uint64_t good_events_c = 0;

static std::map<int,uint64_t> raw_protons;
static std::map<int,uint64_t> good_protons;
static std::map<int,uint64_t> raw_events;
static std::map<int,uint64_t> good_events;
static std::map<int,int> first_data_frame;

static Poco::Mutex the_mutex;

static void all_event_callback(const DAEEventHeader* head, const DetectorEvent32* det_ev, int n, int event_source_id, const int* mapping)
{
        Poco::Mutex::ScopedLock _lock(the_mutex);
        all_events_c += n;
}

//static std::map<int,int> spec_map;

static void good_event_callback(const DAEEventHeader* head, const DetectorEvent32* det_ev, int n, int event_source_id, const int* mapping)
{
        Poco::Mutex::ScopedLock _lock(the_mutex);
        int code;
        for(int i=0; i<n; ++i)
        {
                code = 1000 * det_ev[i].spectrum + event_source_id;
                if ( first_data_frame.find(code) == first_data_frame.end() )
                {
                        first_data_frame[code] = head->frame_number;
                }
        }
        good_events_c += n;
}


static void convert_good_event_callback(const DAEEventHeader* head, const DetectorEvent32* det_ev, int n, int event_source_id, const int* mapping)
{
        Poco::Mutex::ScopedLock _lock(the_mutex);
        good_events_c += n;
//      if (event_source_id == 2)
//      {
//              for(int i=0; i<n; ++i)
//              {
//                      ++(spec_map[det_ev[i].spectrum]);
//              }
//      }
}

static void all_frame_callback(const DAEEventHeader* head, int event_source_id)
{
        Poco::Mutex::ScopedLock _lock(the_mutex);
//      std::cerr << "Frame " << head->frame_number << " time " << DAEEventList::DAETimeAsString(head->time) << std::endl;
        raw_protons[head->frame_number] = head->protons;
        raw_events[head->frame_number] += head->num_events;
}

static void good_frame_callback(const DAEEventHeader* head, int event_source_id)
{
        Poco::Mutex::ScopedLock _lock(the_mutex);
        good_protons[head->frame_number] = head->protons;
        good_events[head->frame_number] += head->num_events;
}

#define PPP_TO_UAMPH    1.738E-6

int EventsToolApp::main(const std::vector<std::string>& args)
{
//      Poco::Logger::root().setLevel(Poco::Message::PRIO_TRACE);  // default for new loggers
        Poco::Logger::root().setLevel(Poco::Message::PRIO_INFORMATION);  // default for new loggers
        logger().setLevel(Poco::Message::PRIO_TRACE);   // sets for current application logger

        if (loadAppConfig("isisicp") != 0)
        {
                return Poco::Util::Application::EXIT_CONFIG;
        }
        if (loadAppConfig("events_tool") != 0)
        {
                return Poco::Util::Application::EXIT_CONFIG;
        }

        int run_number = Poco::Util::Application::instance().config().getInt("events_tool.run_number", -1);
        std::string raw_file = Poco::Util::Application::instance().config().getString("events_tool.raw_file","");
        std::string nexus_file = Poco::Util::Application::instance().config().getString("events_tool.nexus_file","");
        std::string live_host = Poco::Util::Application::instance().config().getString("events_tool.live_host","");

        if (live_host.size() > 0)
        {
                return liveData(live_host);
        }
        if (raw_file.size() > 0)
        {
                return convertRawToEvents(raw_file, run_number);
        }
        if (nexus_file.size() > 0)
        {
                return processNeXusFile(nexus_file);
        }
        EventStore es;
        EventStore::FileEventSourceInfo info;
        es.addFileInputSource(run_number, info);
        LOGSTR_INFORMATION("Info file says " << info.raw_frames << " raw frames and " << info.raw_events << " raw events");
        es.addAllEventCallback(all_event_callback);
        es.addGoodEventCallback(good_event_callback);
        es.addAllFrameCallback(all_frame_callback);
        es.addGoodFrameCallback(good_frame_callback);
        es.start();
        while (!es.endHeaderSeen())
        {
                Sleep(500);
        }
        es.stop(false);
        LOGSTR_INFORMATION("Found " << raw_protons.size() << " raw frames and " << good_protons.size() << " good frames");
        LOGSTR_INFORMATION("Found " << all_events_c << " and " << sum_map_values(raw_events) << " raw events");
        LOGSTR_INFORMATION("Found " << good_events_c << " and " << sum_map_values(good_events) << " good events");
        uint64_t raw_ppp = sum_map_values(raw_protons);
        uint64_t good_ppp = sum_map_values(good_protons);
        LOGSTR_INFORMATION("Raw uAmpH " << raw_ppp * PPP_TO_UAMPH << " good uAmpH " << good_ppp * PPP_TO_UAMPH);
        for(std::map<int,int>::const_iterator it = first_data_frame.begin(); it != first_data_frame.end(); ++it)
        {
                int es = it->first % 1000;
                int spec = it->first / 1000;
                LOGSTR_INFORMATION("Source id " << es << " spectrum " << spec << " first data frame " << it->second);
        }
        return 0;
}

int EventsToolApp::convertRawToEvents(const std::string& raw_file, int run_number)
{
        int i, ndet;
        DAEstatus status;
        ISISRAW *ir = new ISISRAW;
        ir->readFromFile(raw_file.c_str());
        ir->printInfo(std::cout);
        ndet = ir->i_det;
        EventStore es;
        std::vector<int> cards(ndet);
        std::copy(ir->crat, ir->crat + ndet, cards.begin());
        std::sort(cards.begin(), cards.end());
        std::vector<int>::iterator new_end = std::unique(cards.begin(), cards.end());
        cards.resize(new_end - cards.begin());
        time_t run_start = time(NULL);
        for(i=0; i<cards.size(); ++i)
        {
                es.addRawFileInputSource(ir, cards[i], run_start);
        }
        es.setOutputFileName(run_number, true, status);
        es.addGoodEventCallback(convert_good_event_callback);
        es.start();
        while (!es.endHeaderSeen())
        {
                Sleep(500);
        }
        es.stop(false);
        LOGSTR_INFORMATION("good events " << good_events_c);
//      LOGSTR_INFORMATION("spec_map[15] " << spec_map[15]);
        return 0;
}

int EventsToolApp::processNeXusFile(const std::string& nexus_file)
{
        NeXus::File nxfile(nexus_file);
        std::vector<uint64_t> event_index;
        nxfile.openPath("/raw_data_1/detector_1_events/event_index");
        nxfile.getData(event_index);
        int nframes = event_index.size(); // event_index.size();
        nframes = std::min(nframes, 2400);
        int nframes_stuff = std::min(nframes, 100);
        uint64_t iend = event_index[nframes-1];
        std::cerr << "Reading " << nframes << " frames and " << iend << " events from " << nexus_file << std::endl;
        std::vector<uint32_t> event_id(iend);
        nxfile.openPath("/raw_data_1/detector_1_events/event_id");
//      nxfile.getData(event_id);
        std::vector<int> istart, isize;
        istart.push_back(0);
        isize.push_back(iend);
        nxfile.getSlab(&(event_id[0]), istart, isize);

        std::vector<int> crat, spec;
        nxfile.openPath("/raw_data_1/isis_vms_compat/CRAT");
        nxfile.getData(crat);
        nxfile.openPath("/raw_data_1/isis_vms_compat/SPEC");
        nxfile.getData(spec);

        nxfile.close();

        std::map<int,int> crat_map;
        for(int i=0; i<spec.size(); ++i)
        {
                crat_map[spec[i]] = crat[i];
        }


        std::cerr << "Starting to process events" << std::endl;
          
        uint32_t maxspec = 0;
        for(uint32_t i=0; i < nframes-1; ++i)   // i is frame number
        {
                for(uint64_t j=event_index[i]; j < event_index[i+1]; ++j)   // j is index of events for frame i
                {
                        maxspec = std::max(maxspec, event_id[j]); 
                }
        }
        std::cout << " Frames " << nframes << " max spec number " << maxspec << std::endl;

        typedef std::vector< std::vector<uint32_t> > stuff_t;
        stuff_t stuff(nframes_stuff, 0);
        std::vector<int> first_counts(maxspec+1, -1), last_counts(maxspec+1, -1), counts_per_frame(nframes, 0);
        for(uint32_t i=0; i<nframes_stuff; ++i)
        {
                stuff[i].resize(maxspec+1, 0);
        }
        for(uint32_t i=0; i < nframes-1; ++i)   // i is frame number
        {
                for(uint64_t j=event_index[i]; j < event_index[i+1]; ++j)   // j is index of events for frame i
                {
                        ++(counts_per_frame[i]);
                        if ( first_counts[event_id[j]] == -1 )
                        {
                                first_counts[event_id[j]] = i;
                        }
                        last_counts[event_id[j]] = i;
                        if ( i < nframes_stuff )
                        {
                                ++(stuff[i][event_id[j]]);
                        }
                }
        }
        for(uint32_t i=0; i < nframes; ++i)  // frame
        {
//              std::cout << "Total counts in frame " << i << " = " << std::accumulate(stuff[i].begin(), stuff[i].end(), 0) << std::endl;
                std::cout << "Total counts in frame " << i << " = " << counts_per_frame[i] << std::endl;
        }
        for(uint32_t j=0; j <= maxspec; ++j)  // spectra
        {
                std::cout << "Spec " << j << " card " << crat_map[j] << " first counts frame " << first_counts[j] << " last counts frame " << last_counts[j] << std::endl;
//              if (last_counts[j] != -1)
//              {
//                      std::cout << "Spec " << j << " first counts frame " << first_counts[j] << " last counts frame " << last_counts[j] << std::endl;
//              }
        }
//      for(uint32_t j=77830; j <= maxspec; ++j)
//      {
//              std::cout << "* Per frame distribution for spec " << j << std::endl;
//              for(uint32_t i=0; i < nframes; ++i)  // frame
//              {
//                      std::cout << "Frame " << i << " counts " << stuff[i][j] << std::endl;
//              }
//      }
//                      if ( stuff[i].find(j) != stuff[i].end() )
//                      {
//                              std::cout << "Frame " << i << " spectra " << j << " counts " << stuff[i][j] << std::endl;
//                      }
        return 0;
}

int EventsToolApp::liveData(const std::string& host)
{
        static char* junk_buffer[10000];
        Poco::Net::StreamSocket s;
        unsigned port = 10000;
        Poco::Net::SocketAddress address(host, port);
        s.connect(address);
        TCPStreamEventDataSetup setup;
        while( s.available() < sizeof(setup) )
        {
                Poco::Thread::sleep(1000);
        }
        s.receiveBytes(&setup, sizeof(setup));
        if ( !setup.isValid() )
        {
                throw std::runtime_error("version wrong");
        }
        std::cerr << "run number " << setup.head_setup.run_number << std::endl; 
        TCPStreamEventDataNeutron events;
        while(true)
        {
                while(s.available() < sizeof(events.head))
                {
                        Poco::Thread::sleep(100);
                }
                s.receiveBytes(&events.head, sizeof(events.head));
                if ( !events.head.isValid() )
                {
                        throw std::runtime_error("corrupt stream - you should reconnect");
                }
                if ( !(events.head.type == TCPStreamEventHeader::Neutron) )
                {
                        throw std::runtime_error("corrupt stream - you should reconnect");
                }
                s.receiveBytes(junk_buffer, events.head.length - sizeof(events.head));
                while(s.available() < sizeof(events.head_n))
                {
                        Poco::Thread::sleep(100);
                }
                s.receiveBytes(&events.head_n, sizeof(events.head_n));
                if ( !events.head_n.isValid() )
                {
                        throw std::runtime_error("corrupt stream - you should reconnect");
                }
                s.receiveBytes(junk_buffer, events.head_n.length - sizeof(events.head_n));
                events.data.resize(events.head_n.nevents);
                int nread = 0;
                while( nread < events.head_n.nevents )
                {
                        int ntoread = s.available() / sizeof(TCPStreamEventNeutron);
                        if ( ntoread > (events.head_n.nevents - nread) )
                        {
                                ntoread = events.head_n.nevents - nread;
                        }
                        if (ntoread > 0)
                        {
                                s.receiveBytes(&(events.data[nread]), ntoread * sizeof(TCPStreamEventNeutron));
                                nread += ntoread;
                        }
                        else
                        {
                                Poco::Thread::sleep(100);
                        }
                }
                if (!events.isValid())
                {
                        throw std::runtime_error("corrupt stream - you should reconnect");
                }
                TCPStreamEventHeader& head = events.head;
                TCPStreamEventHeaderNeutron& head_n = events.head_n;
                std::cerr << "Read " << nread << " events for frame number " << head_n.frame_number << " time " << head_n.frame_time_zero << std::endl;
                for(int i=0; i<std::min(10,nread); ++i)
                {
                        std::cerr << events.data[i].time_of_flight << " " << events.data[i].spectrum << std::endl;
                }
        }
        s.close();
        return 0;
}