Commit fc5e0200 authored by legoc's avatar legoc

Corrected using time base for calculating rates. Added the lstrates128

application.
parent b4c64db9
......@@ -17,8 +17,8 @@
*/
#include "../../algorithms/nomad/Coincidence.h"
#include <iomanip>
#include <cmath>
using namespace std;
......@@ -211,7 +211,7 @@ void Coincidence::processEvent(Event const & event) {
int32_t index = _channelOffset[event.crate][event.board] + event.channel;
_channelCounters[index]++;
_channelRates[index] = 1.0E8 * (double)_channelCounters[index] / (double)event.time();
_channelRates[index] = pow(10, listModeContext.timeBase) * (double)_channelCounters[index] / (double)event.time();
}
void Coincidence::processEvent(EventBlockArrayIterator const & e) {
......@@ -221,7 +221,7 @@ void Coincidence::processEvent(EventBlockArrayIterator const & e) {
// Update the counter and rates
_channelCounters[index]++;
_channelRates[index] = 1.0E8 * (double)_channelCounters[index] / (double)(*e).time();
_channelRates[index] = pow(10, listModeContext.timeBase) * (double)_channelCounters[index] / (double)(*e).time();
}
void Coincidence::processEventNotInserted(Event const & event) {
......@@ -356,7 +356,7 @@ void Coincidence::updateRates() {
//cout << "time = " << _maxTime << " (" << (double)_maxTime * 1.0E-8 << ")" << endl;
for (int i = 0; i < 4 * 7; ++i) {
_coincidenceRates[i] = 1.0E8 * (double)_counters[i] / (double)_maxTime;
_coincidenceRates[i] = pow(10, listModeContext.timeBase) * (double)_counters[i] / (double)_maxTime;
}
//cout << "rate global = " << _rates[0] << endl;
......
......@@ -20,6 +20,7 @@
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;
......@@ -176,7 +177,7 @@ void Histogram::processEvent(EventBlockArrayIterator const & e) {
// Update the counter and rates
_channelCounters[index]++;
_channelRates[index] = 1.0E8 * (double)_channelCounters[index] / (double)(*e).time();
_channelRates[index] = pow(10, listModeContext.timeBase) * (double)_channelCounters[index] / (double)(*e).time();
// Check if e is a trigger event
if (isTrigger(*e)) {
......
......@@ -3,7 +3,9 @@ liblstdppnomad128_la_SOURCES = \
Coincidence.cpp \
Coincidence.h \
Histogram.cpp \
Histogram.h
Histogram.h \
Rates.cpp \
Rates.h
liblstdppnomad128_la_CPPFLAGS = $(LST_CXXFLAGS)
liblstdppnomad128_la_LDFLAGS = -version-number $(LIBRARY_VERSION) -export-dynamic $(LST_LDFLAGS)
......
/*
* Nomad Instrument Control Software
*
* Copyright 2011 Institut Laue-Langevin
*
* Licensed under the EUPL, Version 1.1 only (the "License");
* You may not use this work except in compliance with the Licence.
* You may obtain a copy of the Licence at:
*
* http://www.osor.eu/eupl
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the Licence is distributed on an "AS IS" basis,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the Licence for the specific language governing permissions and
* limitations under the Licence.
*/
#include "../../algorithms/nomad/Rates.h"
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;
namespace lstdpp128 {
Rates::Rates() :
_totalNbChannels(0),
_channelCounters(0),
_channelRates(0),
_maxBlockSize(0),
_outMemblock(0),
_numberOfGoodEvents(0),
_numberOfBadEvents(0),
_numberOfEventsNotInserted(0),
_numberOfEventsInConflict(0) {
}
void Rates::init(int32_t* headerMemoryBlock, uint32_t version, uint32_t timeBase, uint32_t numberOfBoards, int32_t maxBlockSize, std::string const & outFileName) {
// Generate the channel offset and the total size of the histograms
BoardIterator b;
int32_t currentNbChannels = (*b).nbChannels;
int32_t currentRatesSize = (*b).nbChannels * (1 << BoardDataResolution[(*b).boardType]);
++b;
while (!b.end()) {
_channelOffset[b.crate()][b.board()] = currentNbChannels;
cout << "channel offset " << (int)b.crate() << ", " << (int)b.board() << " = " << _channelOffset[b.crate()][b.board()] << endl;
currentNbChannels += (*b).nbChannels;
currentRatesSize += (*b).nbChannels * (1 << BoardDataResolution[(*b).boardType]);
++b;
}
// The total number of channels is the current number of channels
_totalNbChannels = currentNbChannels;
// Allocate the channel counters and rates
_channelCounters = new uint64_t[_totalNbChannels];
for (int32_t i = 0; i < _totalNbChannels; ++i) {
_channelCounters[i] = 0;
}
_channelRates = new double[_totalNbChannels];
for (int32_t i = 0; i < _totalNbChannels; ++i) {
_channelRates[i] = 0.0;
}
_numberOfGoodEvents = 0;
_numberOfBadEvents = 0;
_numberOfEventsNotInserted = 0;
_numberOfEventsInConflict = 0;
// Manage output file
_outFileName = outFileName;
if (_outFileName != "") {
string tempOutFileName = outFileName + ".temp";
_outDataFile.open(tempOutFileName.c_str());
uint32_t nbEvents = 0;
// Write the header.
_outDataFile.write((char *)&version, sizeof(uint32_t));
_outDataFile.write((char *)&timeBase, sizeof(uint32_t));
_outDataFile.write((char *)&nbEvents, sizeof(uint32_t));
_outDataFile.write((char *)&numberOfBoards, sizeof(uint32_t));
_outDataFile.write((char *)headerMemoryBlock, numberOfBoards * sizeof(uint32_t));
}
_maxBlockSize = maxBlockSize;
_outMemblock = new int32_t[maxBlockSize * 4];
}
Rates::~Rates() {
delete [] _channelCounters;
delete [] _channelRates;
delete [] _outMemblock;
}
int32_t Rates::numberOfGoodEvents() const {
return _numberOfGoodEvents;
}
int32_t Rates::numberOfBadEvents() const {
return _numberOfBadEvents;
}
int32_t Rates::numberOfEventsNotInserted() const {
return _numberOfEventsNotInserted;
}
int32_t Rates::numberOfEventsInConflict() const {
return _numberOfEventsInConflict;
}
bool Rates::checkEvent(Event const & event) {
return true;
}
int32_t Rates::bitshift(Event const & event) const {
return MAX_BITS - BoardDataResolution[event.boardType()];
}
void Rates::processEvent(EventBlockArrayIterator const & e) {
_numberOfGoodEvents++;
int32_t index = _channelOffset[(*e).crate][(*e).board] + (*e).channel;
// Update the counter and rates
_channelCounters[index]++;
// Convert to Hz using time base.
_channelRates[index] = pow(10, listModeContext.timeBase) * (double)_channelCounters[index] / (double)(*e).time();
}
void Rates::processBadEvent(Event const & event) {
_numberOfBadEvents++;
}
void Rates::processEventNotInserted(Event const & event) {
_numberOfGoodEvents++;
_numberOfEventsNotInserted++;
}
void Rates::finish() {
if (_outFileName != "") {
// Write the number of events.
// The pointer is 2 uint32 after the beginning of the file.
_outDataFile.seekp(2 * sizeof(uint32_t), ios::beg);
_outDataFile.write((char *)&_numberOfGoodEvents, sizeof(uint32_t));
_outDataFile.close();
string tempOutFileName = _outFileName + ".temp";
rename(tempOutFileName.c_str(), _outFileName.c_str());
}
}
void Rates::printRates() {
cout << "channel rates" << endl;
BoardIterator b;
int32_t channel = 0;
while (!b.end()) {
cout << "crate " << (int)b.crate() << " board " << (int)b.board() << endl;
for (int32_t c = 0; c < (*b).nbChannels; ++c) {
cout << "\t" << setw(5) << c << setw(10) << _channelRates[channel] << endl;
++channel;
}
cout << endl;
++b;
}
}
int32_t Rates::getTotalNumberOfChannels() const {
return _totalNbChannels;
}
double * Rates::getChannelRates() const {
return _channelRates;
}
}
/*
* Nomad Instrument Control Software
*
* Copyright 2011 Institut Laue-Langevin
*
* Licensed under the EUPL, Version 1.1 only (the "License");
* You may not use this work except in compliance with the Licence.
* You may obtain a copy of the Licence at:
*
* http://www.osor.eu/eupl
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the Licence is distributed on an "AS IS" basis,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the Licence for the specific language governing permissions and
* limitations under the Licence.
*/
#ifndef LSTDPP_RATES_H
#define LSTDPP_RATES_H
#include <cstring>
#include <algorithm>
#include <fstream>
#include <vector>
#include <time.h>
#include <map>
#include "../../Event.h"
#include "../../EventBlockArray.h"
namespace lstdpp128 {
class Rates {
public:
Rates();
~Rates();
void init(int32_t* headerMemoryBlock, uint32_t version, uint32_t timeBase, uint32_t numberOfBoards, int32_t maxBlockSize, std::string const & outFileName);
int32_t numberOfGoodEvents() const;
int32_t numberOfBadEvents() const;
int32_t numberOfEventsNotInserted() const;
int32_t numberOfEventsInConflict() const;
/**
* Functions to implement.
*/
bool checkEvent(Event const & event);
void processEvent(EventBlockArrayIterator const & e);
void processBadEvent(Event const & event);
void processEventNotInserted(Event const & event);
template<typename Iterator>
Iterator processBlock(Iterator first, Iterator last, bool lastBlock);
void finish();
void printRates();
int32_t getTotalNumberOfChannels() const;
double * getChannelRates() const;
private:
void writeBlock(EventBlock const & block);
int32_t bitshift(Event const & event) const;
// crate, board - channel offset
std::map<int32_t, std::map<int32_t, int32_t> > _channelOffset;
int32_t _totalNbChannels;
uint64_t * _channelCounters;
double * _channelRates;
int32_t _maxBlockSize;
int32_t * _outMemblock;
std::string _outFileName;
std::ofstream _outDataFile;
int32_t _numberOfGoodEvents;
int32_t _numberOfBadEvents;
int32_t _numberOfEventsNotInserted;
int32_t _numberOfEventsInConflict;
static const int16_t MAX_BITS = 15;
};
template<typename Iterator>
Iterator Rates::processBlock(Iterator first, Iterator last, bool lastBlock) {
if (_outFileName != "") {
int32_t * buffer = _outMemblock;
int32_t index = 0;
// iterating events
for (Iterator e = first; e != last; ++e) {
writeEvent(*e, buffer + index);
index += 4;
}
if (index <= _maxBlockSize * 4) {
_outDataFile.write((char *)_outMemblock, index * sizeof(int32_t));
} else {
std::cerr << "cannot write block with size " << index << std::endl;
}
}
return last;
}
}
#endif
......@@ -18,6 +18,7 @@
#include "CoincidenceOffline.h"
#include <iomanip>
#include <cmath>
using namespace std;
......@@ -228,7 +229,7 @@ void CoincidenceOffline::processEvent(Event & event) {
int32_t index = _channelOffset[event.crate][event.board] + event.channel;
_channelCounters[index]++;
_channelRates[index] = 1.0E8 * (double)_channelCounters[index] / (double)event.time();
_channelRates[index] = pow(10, listModeContext.timeBase) * (double)_channelCounters[index] / (double)event.time();
// Root
_runQA->_htimeVsEvt->SetBinContent(_numberOfGoodEvents,event.time());
......
/*
* Nomad Instrument Control Software
*
* Copyright 2011 Institut Laue-Langevin
*
* Licensed under the EUPL, Version 1.1 only (the "License");
* You may not use this work except in compliance with the Licence.
* You may obtain a copy of the Licence at:
*
* http://www.osor.eu/eupl
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the Licence is distributed on an "AS IS" basis,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the Licence for the specific language governing permissions and
* limitations under the Licence.
*/
#include <cstdlib>
#include <fstream>
#include <sstream>
#include <cstring>
#include <algorithm>
#include <time.h>
#include "../../algorithms/nomad/Rates.h"
#include "../../Event.h"
#include "../../EventBlockArray.h"
#include "../../EventBlockArrayAlgorithm.h"
#include "../../Reader.h"
using namespace std;
using namespace lstdpp128;
int32_t MAX_BLOCK_SIZE = 0;
int32_t BUFFER_SIZE = 1 << 12;
int percentage(int a, int b) {
return 100.0f * (float)a / (float)b;
}
int main(int argc, char * argv[]) {
if (argc < 2) {
cerr << "usage : lstrates <file> [maxBlocks] [blockSize]" << endl;
return EXIT_FAILURE;
}
string fileName = argv[1];
int32_t maxBlocks = -1;
if (argc >= 3) {
istringstream im(argv[2]);
im >> maxBlocks;
}
if (argc >= 4) {
istringstream is(argv[3]);
is >> BUFFER_SIZE;
}
MAX_BLOCK_SIZE = BUFFER_SIZE * 2;
cout << "reading " << fileName << ", max blocks " << maxBlocks << ", block size " << BUFFER_SIZE << endl;
string outFileName = fileName.substr(0, fileName.find_first_of('.')) + ".lsto";
Reader reader(BUFFER_SIZE);
bool success = reader.open(fileName);
if (!success) {
return EXIT_FAILURE;
}
// algorithm
EventBlockArrayAlgorithm<Rates> algorithm;
algorithm.processor().init(reader.boardHeader(), reader.version(), reader.timeBase(), reader.totalNumberOfBoards(), MAX_BLOCK_SIZE, outFileName);
// creating the events
EventBlockArray blocks(MAX_BLOCK_SIZE, maxBlocks);
int currentPercentage = 0;
time_t start = time(0);
time_t t1 = start;
// reading list mode
while (true) {
// reading next block
int32_t size = reader.read();
// exiting if no data
if (size == 0) {
break;
}
// processing block
algorithm.process(reader.buffer(), size, blocks);
}
time_t t2 = time(0);
algorithm.finish(blocks);
cout << "#events read = " << algorithm.processor().numberOfGoodEvents() << " (" << ((float)algorithm.processor().numberOfGoodEvents() / (float)(t2 - t1)) << " event/s)" << endl;
cout << "#events bad = " << algorithm.processor().numberOfBadEvents() << " (" << percentage(algorithm.processor().numberOfBadEvents(), algorithm.processor().numberOfGoodEvents()) << "%)" << endl;
cout << "#events not inserted = " << algorithm.processor().numberOfEventsNotInserted() << " (" << percentage(algorithm.processor().numberOfEventsNotInserted(), algorithm.processor().numberOfGoodEvents()) << "%)" << endl;
cout << "#events in conflict = " << algorithm.processor().numberOfEventsInConflict() << " (" << percentage(algorithm.processor().numberOfEventsInConflict(), algorithm.processor().numberOfGoodEvents()) << "%)" << endl;
algorithm.processor().printRates();
return EXIT_SUCCESS;
}
bin_PROGRAMS = \
lsthisto128 \
lstcoincidence128
lstcoincidence128 \
lstrates128
libs = ../../liblstdpp128.la ../../algorithms/nomad/liblstdppnomad128.la
......@@ -18,3 +19,12 @@ lstcoincidence128_SOURCES = \
lstcoincidence128_CPPFLAGS = $(LST_CXXFLAGS)
lstcoincidence128_LDFLAGS = $(LST_LDFLAGS)
lstcoincidence128_LDADD = $(libs) $(LST_LDLIBS)
lstrates128_SOURCES = \
IterativeRates.cpp
lstrates128_CPPFLAGS = $(LST_CXXFLAGS)
lstrates128_LDFLAGS = $(LST_LDFLAGS)
lstrates128_LDADD = $(libs) $(LST_LDLIBS)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment