Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • jlab/hallc/analyzer_software/hcana
  • whit/hcana
2 results
Show changes
Showing
with 97 additions and 2362 deletions
......@@ -105,6 +105,7 @@ THcShowerArray::~THcShowerArray()
delete [] fAdcTimeWindowMin; fAdcTimeWindowMin = 0;
delete [] fAdcTimeWindowMax; fAdcTimeWindowMax = 0;
delete [] fPedDefault; fPedDefault = 0;
}
//_____________________________________________________________________________
......@@ -280,6 +281,7 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
fAdcTimeWindowMin = new Double_t [fNelem];
fAdcTimeWindowMax = new Double_t [fNelem];
fPedDefault = new Int_t [fNelem];
DBRequest list1[]={
{"cal_arr_ped_limit", fPedLimit, kInt, static_cast<UInt_t>(fNelem),1},
......@@ -287,12 +289,14 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
{"cal_arr_gain_cor", cal_arr_gain_cor, kDouble, static_cast<UInt_t>(fNelem)},
{"cal_arr_AdcTimeWindowMin", fAdcTimeWindowMin, kDouble, static_cast<UInt_t>(fNelem),1},
{"cal_arr_AdcTimeWindowMax", fAdcTimeWindowMax, kDouble, static_cast<UInt_t>(fNelem),1},
{"cal_arr_PedDefault", fPedDefault, kInt, static_cast<UInt_t>(fNelem),1},
{0}
};
for(Int_t ip=0;ip<fNelem;ip++) {
fAdcTimeWindowMin[ip] = -1000.;
fAdcTimeWindowMax[ip] = 1000.;
fPedDefault[ip] = 0;
}
gHcParms->LoadParmValues((DBRequest*)&list1, prefix);
......@@ -873,6 +877,8 @@ void THcShowerArray::FillADC_DynamicPedestal()
{
Double_t StartTime = 0.0;
if( fglHod ) StartTime = fglHod->GetStartTime();
Double_t OffsetTime = 0.0;
if( fglHod ) OffsetTime = fglHod->GetOffsetTime();
for (Int_t ielem=0;ielem<frAdcPulseInt->GetEntries();ielem++) {
Int_t npad = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
......@@ -881,20 +887,15 @@ void THcShowerArray::FillADC_DynamicPedestal()
Double_t pulseInt = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetData();
Double_t pulseAmp = ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(ielem))->GetData();
Double_t pulseTime = ((THcSignalHit*) frAdcPulseTime->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime;
Bool_t errorflag = ((THcSignalHit*) frAdcErrorFlag->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
Bool_t pulseTimeCut = (adctdcdiffTime > fAdcTimeWindowMin[npad]) && (adctdcdiffTime < fAdcTimeWindowMax[npad]);
if (!errorflag)
{
fGoodAdcMult.at(npad) += 1;
}
if (!errorflag && pulseTimeCut) {
if (pulseTimeCut) {
fTotNumAdcHits++;
fGoodAdcPulseIntRaw.at(npad) = pulseIntRaw;
if(fGoodAdcPulseIntRaw.at(npad) > fThresh[npad]) {
if(fGoodAdcPulseIntRaw.at(npad) > fThresh[npad] && fGoodAdcPulseInt.at(npad)==0) {
fTotNumGoodAdcHits++;
fGoodAdcPulseInt.at(npad) = pulseInt;
fE.at(npad) = fGoodAdcPulseInt.at(npad)*fGain[npad];
......@@ -982,6 +983,24 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
} else {
((THcSignalHit*) frAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(padnum,1);
}
if (rawAdcHit.GetPulseAmpRaw(thit) <= 0) {
Double_t PeakPedRatio= rawAdcHit.GetF250_PeakPedestalRatio();
Int_t NPedSamples= rawAdcHit.GetF250_NPedestalSamples();
Double_t AdcToC = rawAdcHit.GetAdcTopC();
Double_t AdcToV = rawAdcHit.GetAdcTomV();
if (fPedDefault[padnum-1] !=0) {
Double_t tPulseInt = AdcToC*(rawAdcHit.GetPulseIntRaw(thit) - fPedDefault[padnum-1]*PeakPedRatio);
((THcSignalHit*) frAdcPulseInt->ConstructedAt(nrAdcHits))->Set(padnum, tPulseInt);
((THcSignalHit*) frAdcPedRaw->ConstructedAt(nrAdcHits))->Set(padnum, fPedDefault[padnum-1]);
((THcSignalHit*) frAdcPed->ConstructedAt(nrAdcHits))->Set(padnum, float(fPedDefault[padnum-1])/float(NPedSamples)*AdcToV);
}
((THcSignalHit*) frAdcPulseAmp->ConstructedAt(nrAdcHits))->Set(padnum, 0.);
}
++nrAdcHits;
}
ihit++;
......
......@@ -150,6 +150,7 @@ protected:
static const Int_t kADCSampIntDynPed=3;
Double_t *fAdcTimeWindowMin ;
Double_t *fAdcTimeWindowMax ;
Int_t *fPedDefault ;
Double_t fAdcThreshold ;
Double_t fAdcTdcOffset;
......
......@@ -7,6 +7,7 @@
#include <iterator>
#include <iostream>
#include <memory>
#include <vector>
#include "TMath.h"
using namespace std;
......
......@@ -612,6 +612,21 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
} else {
((THcSignalHit*) frPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(padnum,1);
}
if (rawPosAdcHit.GetPulseAmpRaw(thit) <= 0) {
Double_t PeakPedRatio= rawPosAdcHit.GetF250_PeakPedestalRatio();
Int_t NPedSamples= rawPosAdcHit.GetF250_NPedestalSamples();
Double_t AdcToC = rawPosAdcHit.GetAdcTopC();
Double_t AdcToV = rawPosAdcHit.GetAdcTomV();
Int_t PedDefaultTemp = static_cast<THcShower*>(fParent)->GetPedDefault(padnum-1,fLayerNum-1,0);
if (PedDefaultTemp !=0) {
Double_t tPulseInt = AdcToC*(rawPosAdcHit.GetPulseIntRaw(thit) - PedDefaultTemp*PeakPedRatio);
((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(nrPosAdcHits))->Set(padnum, tPulseInt);
((THcSignalHit*) frPosAdcPedRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, PedDefaultTemp);
((THcSignalHit*) frPosAdcPed->ConstructedAt(nrPosAdcHits))->Set(padnum, float(PedDefaultTemp)/float(NPedSamples)*AdcToV);
}
((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(nrPosAdcHits))->Set(padnum, 0.);
}
++nrPosAdcHits;
fTotNumAdcHits++;
fTotNumPosAdcHits++;
......@@ -637,6 +652,21 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
} else {
((THcSignalHit*) frNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(padnum,1);
}
if (rawNegAdcHit.GetPulseAmpRaw(thit) <= 0) {
Double_t PeakPedRatio= rawNegAdcHit.GetF250_PeakPedestalRatio();
Int_t NPedSamples= rawNegAdcHit.GetF250_NPedestalSamples();
Double_t AdcToC = rawNegAdcHit.GetAdcTopC();
Double_t AdcToV = rawNegAdcHit.GetAdcTomV();
Int_t PedDefaultTemp = static_cast<THcShower*>(fParent)->GetPedDefault(padnum-1,fLayerNum-1,1);
if (PedDefaultTemp !=0) {
Double_t tPulseInt = AdcToC*(rawNegAdcHit.GetPulseIntRaw(thit) - PedDefaultTemp*PeakPedRatio);
((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(nrNegAdcHits))->Set(padnum, tPulseInt);
((THcSignalHit*) frNegAdcPedRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, PedDefaultTemp);
((THcSignalHit*) frNegAdcPed->ConstructedAt(nrNegAdcHits))->Set(padnum, float(PedDefaultTemp)/float(NPedSamples)*AdcToV);
}
((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(nrNegAdcHits))->Set(padnum, 0.);
}
++nrNegAdcHits;
fTotNumAdcHits++;
fTotNumNegAdcHits++;
......@@ -744,6 +774,8 @@ void THcShowerPlane::FillADC_DynamicPedestal()
{
Double_t StartTime = 0.0;
if( fglHod ) StartTime = fglHod->GetStartTime();
Double_t OffsetTime = 0.0;
if( fglHod ) OffsetTime = fglHod->GetOffsetTime();
for (Int_t ielem=0;ielem<frNegAdcPulseInt->GetEntries();ielem++) {
Int_t npad = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
Double_t pulseInt = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetData();
......@@ -751,18 +783,11 @@ void THcShowerPlane::FillADC_DynamicPedestal()
Double_t pulseAmp = ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(ielem))->GetData();
Double_t pulseIntRaw = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
Double_t pulseTime = ((THcSignalHit*) frNegAdcPulseTime->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime;
Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
Double_t threshold = ((THcSignalHit*) frNegAdcThreshold->ConstructedAt(ielem))->GetData();
Bool_t errorflag = ((THcSignalHit*) frNegAdcErrorFlag->ConstructedAt(ielem))->GetData();
Bool_t pulseTimeCut = (adctdcdiffTime > static_cast<THcShower*>(fParent)->GetWindowMin(npad,fLayerNum-1,1)) && (adctdcdiffTime < static_cast<THcShower*>(fParent)->GetWindowMax(npad,fLayerNum-1,1) );
if (!errorflag)
{
fGoodNegAdcMult.at(npad) += 1;
}
if (!errorflag && pulseTimeCut) {
if (pulseTimeCut) {
fGoodNegAdcPulseIntRaw.at(npad) =pulseIntRaw;
if(fGoodNegAdcPulseIntRaw.at(npad) > threshold && fGoodNegAdcPulseInt.at(npad)==0) {
......@@ -793,18 +818,10 @@ void THcShowerPlane::FillADC_DynamicPedestal()
Double_t pulseInt = ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(ielem))->GetData();
Double_t pulseIntRaw = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
Double_t pulseTime = ((THcSignalHit*) frPosAdcPulseTime->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime;
Bool_t errorflag = ((THcSignalHit*) frPosAdcErrorFlag->ConstructedAt(ielem))->GetData();
Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
Bool_t pulseTimeCut = (adctdcdiffTime > static_cast<THcShower*>(fParent)->GetWindowMin(npad,fLayerNum-1,0)) && (adctdcdiffTime < static_cast<THcShower*>(fParent)->GetWindowMax(npad,fLayerNum-1,0) );
if (!errorflag)
{
fGoodPosAdcMult.at(npad) += 1;
}
if (!errorflag && pulseTimeCut) {
if (pulseTimeCut) {
fGoodPosAdcPulseIntRaw.at(npad) = pulseIntRaw;
if(fGoodPosAdcPulseIntRaw.at(npad) > threshold && fGoodPosAdcPulseInt.at(npad)==0) {
......
......@@ -29,7 +29,7 @@ public:
};
void SetXY(Double_t x, Double_t y) {fX = x; fY = y;};
void Clear(Option_t* opt="") {fNHits=0; fNCombos=0; fHits.clear();};
void Clear(Option_t* /* opt */ ="") {fNHits=0; fNCombos=0; fHits.clear();};
void AddHit(THcDCHit* hit) {
Hit newhit;
newhit.dchit = hit;
......
......@@ -211,7 +211,8 @@ void THcTrigDet::Clear(Option_t* opt) {
THaAnalysisObject::Clear(opt);
// Reset all data.
for (int i=0; i<fNumAdc; ++i) {
fTdcRefTime = kBig;
for (int i=0; i<fNumAdc; ++i) {
fAdcPedRaw[i] = 0;
fAdcPulseIntRaw[i] = 0;
fAdcPulseAmpRaw[i] = 0;
......@@ -255,7 +256,7 @@ Int_t THcTrigDet::Decode(const THaEvData& evData) {
UInt_t good_hit=999;
for (UInt_t thit=0; thit<rawAdcHit.GetNPulses(); ++thit) {
Int_t TestTime=rawAdcHit.GetPulseTimeRaw(thit);
if (TestTime>fAdcTimeWindowMin[cnt]&&TestTime<fAdcTimeWindowMax[cnt]&&good_hit==999) {
if (TestTime>=fAdcTimeWindowMin[cnt]&&TestTime<=fAdcTimeWindowMax[cnt]&&good_hit==999) {
good_hit=thit;
}
}
......@@ -272,14 +273,21 @@ Int_t THcTrigDet::Decode(const THaEvData& evData) {
}
else if (hit->fPlane == 2) {
THcRawTdcHit rawTdcHit = hit->GetRawTdcHit();
if (rawTdcHit.GetNHits() >0 && rawTdcHit.HasRefTime() && fTdcRefTime == kBig) fTdcRefTime=rawTdcHit.GetRefTime() ;
UInt_t good_hit=999;
UInt_t closest_hit=999;
Int_t TimeDiff=1000000;
for (UInt_t thit=0; thit<rawTdcHit.GetNHits(); ++thit) {
Int_t TestTime= rawTdcHit.GetTimeRaw(thit);
if (TestTime>fTdcTimeWindowMin[cnt]&&TestTime<fTdcTimeWindowMax[cnt]&&good_hit==999) {
if (abs(TestTime-fTdcTimeWindowMin[cnt]) < TimeDiff) {
closest_hit=thit;
TimeDiff=abs(TestTime-fTdcTimeWindowMin[cnt]);
}
if (TestTime>=fTdcTimeWindowMin[cnt]&&TestTime<=fTdcTimeWindowMax[cnt]&&good_hit==999) {
good_hit=thit;
}
}
if (good_hit == 999 and closest_hit != 999) good_hit=closest_hit;
if (good_hit<rawTdcHit.GetNHits()) {
fTdcTimeRaw[cnt] = rawTdcHit.GetTimeRaw(good_hit);
fTdcTime[cnt] = rawTdcHit.GetTime(good_hit)*fTdcChanperNS+fTdcOffset;
......@@ -334,11 +342,11 @@ Int_t THcTrigDet::ReadDatabase(const TDatime& date) {
fTdcTimeWindowMax = new Double_t [fNumTdc];
for(Int_t ip=0;ip<fNumAdc;ip++) { // Set a large default window
fAdcTimeWindowMin[ip]=0;
fAdcTimeWindowMax[ip]=10000;
fAdcTimeWindowMax[ip]=100000;
}
for(Int_t ip=0;ip<fNumTdc;ip++) { // Set a large default window
fTdcTimeWindowMin[ip]=0;
fTdcTimeWindowMax[ip]=10000;
fTdcTimeWindowMax[ip]=100000;
}
DBRequest list2[]={
{"_AdcTimeWindowMin", fAdcTimeWindowMin, kDouble, (UInt_t) fNumAdc, 1},
......@@ -396,7 +404,16 @@ Int_t THcTrigDet::DefineVariables(THaAnalysisObject::EMode mode) {
std::vector<TString> adcPulseIntTitle(fNumAdc), adcPulseIntVar(fNumAdc);
std::vector<TString> adcPulseAmpTitle(fNumAdc), adcPulseAmpVar(fNumAdc);
std::vector<TString> adcMultiplicityTitle(fNumAdc), adcMultiplicityVar(fNumAdc);
TString RefTimeTitle= "TdcRefTime";
TString RefTimeVar= "fTdcRefTime";
RVarDef entryRefTime {
RefTimeTitle.Data(),
RefTimeTitle.Data(),
RefTimeVar.Data()
};
vars.push_back(entryRefTime);
for (int i=0; i<fNumAdc; ++i) {
adcPedRawTitle.at(i) = fAdcNames.at(i) + "_adcPedRaw";
adcPedRawVar.at(i) = TString::Format("fAdcPedRaw[%d]", i);
......
......@@ -82,6 +82,7 @@ class THcTrigDet : public THaDetector, public THcHitList {
Int_t fTdcMultiplicity[fMaxTdcChannels];
Int_t fAdcMultiplicity[fMaxAdcChannels];
Double_t fTdcRefTime;
TString fSpectName;
std::vector<Int_t> eventtypes;
......
// Preamble to HallC_LinkDef.h file
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ nestedclass;
#pragma link C++ nestedtypedef;
......@@ -19,14 +19,12 @@
#pragma link C++ class hallc::data::DriftChamber+;
#pragma link C++ class hallc::data::PulseWaveForm+;
#pragma link C++ class std::vector<hallc::data::PulseWaveForm>+;
#pragma link C++ class std::vector < hallc::data::PulseWaveForm>+;
#pragma link C++ class podd2::HitLogging<TObject>+;
#pragma link C++ class podd2::HitLogging < TObject>+;
#pragma link C++ global gHcParms;
#pragma link C++ global gHcDetectorMap;
#pragma link C++ class Decoder::Scaler9001+;
#pragma link C++ class Decoder::Scaler9250+;
......@@ -55,6 +53,7 @@
#pragma link C++ class THcHodoEff+;
#pragma link C++ class THcHelicity+;
#pragma link C++ class THcHelicityReader+;
#pragma link C++ class THcHelicityScaler+;
#pragma link C++ class THcHodoHit+;
#pragma link C++ class THcHodoscope+;
#pragma link C++ class THcInterface+;
......@@ -93,6 +92,5 @@
#pragma link C++ class hcana::Scandalizer+;
#pragma link C++ class hcana::TrackingEfficiency+;
#pragma link C++ class hcana::OnlineMonitor+;
// Postamble for HallC_Linkdef.h file
#endif
......@@ -15,7 +15,7 @@ namespace hallc {
static const UInt_t MaxNSamples = 511;
// From THcRawAdcHit.h
PulseWaveForm() {}
PulseWaveForm(Int_t* buf, Int_t size = MaxNSamples) { std::copy_n(buf, MaxNSamples, std::begin(_buffer)); }
PulseWaveForm(Int_t* buf, Int_t size = MaxNSamples) { std::copy_n(buf, size, std::begin(_buffer)); }
virtual ~PulseWaveForm() {}
void ZeroBuffer() { std::fill(std::begin(_buffer), std::end(_buffer), 0); }
......
Hcana tests
===========
## Current tests:
- ep elastic tests
## Tests to add:
- Start time check
- Focal plane times test
- "Tracking efficiency" test
- Invariant mass check (jpsi)
#include <fmt/core.h>
#include <fmt/ostream.h>
#include <vector>
#include "TString.h"
R__LOAD_LIBRARY(libHallC.so)
#include "hcana/HallC_Data.h"
#include "DecData.h"
//R__LOAD_LIBRARY(libScandalizer.so)
//#include "monitor/DetectorDisplay.h"
//#include "monitor/DisplayPlots.h"
//#include "monitor/MonitoringDisplay.h"
//#include "scandalizer/PostProcessors.h"
//#include "scandalizer/ScriptHelpers.h"
//
//#include "THaPostProcess.h"
//#include "monitor/ExperimentMonitor.h"
//#include "scandalizer/PostProcessors.h"
#include "THcAnalyzer.h"
#include "THaCut.h"
#include "THcGlobals.h"
#include "THcHallCSpectrometer.h"
#include "THcDetectorMap.h"
#include "THcCherenkov.h"
#include "THcDC.h"
#include "THcHodoscope.h"
#include "THcParmList.h"
#include "THaGoldenTrack.h"
#include "THcHodoEff.h"
#include "THcScalerEvtHandler.h"
#include "THcShower.h"
#include "THcReactionPoint.h"
#include "THcExtTarCor.h"
#include "THcRasteredBeam.h"
#include "THcRun.h"
#include "THcCoinTime.h"
#include "THcConfigEvtHandler.h"
#include "THcTrigDet.h"
#include "THcTrigApp.h"
#include "THcSecondaryKine.h"
#include "THcAerogel.h"
#include "THcPrimaryKine.h"
#include "THaReactionPoint.h"
void elastic_coin_replay(Int_t RunNumber = 0, Int_t MaxEvent = -1) {
using namespace std;
// Get RunNumber and MaxEvent if not provided.
if( RunNumber<=0 ) {
std::exit(-1);
}
// Create file name patterns.
const char* RunFileNamePattern = "coin_all_%05d.dat";
vector<TString> pathList;
pathList.push_back(".");
pathList.push_back("./raw");
pathList.push_back("./raw/../raw.copiedtotape");
pathList.push_back("./cache");
//const char* RunFileNamePattern = "raw/coin_all_%05d.dat";
const char* ROOTFileNamePattern = "ROOTfiles/coin_replay_production_%d_%d.root";
// Load global parameters
gHcParms->Define("gen_run_number", "Run Number", RunNumber);
gHcParms->AddString("g_ctp_database_filename", "DBASE/COIN/standard.database");
gHcParms->Load(gHcParms->GetString("g_ctp_database_filename"), RunNumber);
gHcParms->Load(gHcParms->GetString("g_ctp_parm_filename"));
gHcParms->Load(gHcParms->GetString("g_ctp_kinematics_filename"), RunNumber);
// Load params for COIN trigger configuration
gHcParms->Load("PARAM/TRIG/tcoin.param");
// Load fadc debug parameters
gHcParms->Load("PARAM/HMS/GEN/h_fadc_debug.param");
gHcParms->Load("PARAM/SHMS/GEN/p_fadc_debug.param");
// const char* CurrentFileNamePattern = "low_curr_bcm/bcmcurrent_%d.param";
// gHcParms->Load(Form(CurrentFileNamePattern, RunNumber));
// Load the Hall C detector map
gHcDetectorMap = new THcDetectorMap();
gHcDetectorMap->Load("MAPS/COIN/DETEC/coin.map");
// Dec data
gHaApps->Add(new Podd::DecData("D","Decoder raw data"));
//=:=:=:=
// SHMS
//=:=:=:=
// Set up the equipment to be analyzed.
THcHallCSpectrometer* SHMS = new THcHallCSpectrometer("P", "SHMS");
SHMS->SetEvtType(1);
SHMS->AddEvtType(4);
SHMS->AddEvtType(5);
SHMS->AddEvtType(6);
SHMS->AddEvtType(7);
gHaApps->Add(SHMS);
// Add Noble Gas Cherenkov to SHMS apparatus
THcCherenkov* pngcer = new THcCherenkov("ngcer", "Noble Gas Cherenkov");
SHMS->AddDetector(pngcer);
// Add drift chambers to SHMS apparatus
THcDC* pdc = new THcDC("dc", "Drift Chambers");
SHMS->AddDetector(pdc);
// Add hodoscope to SHMS apparatus
THcHodoscope* phod = new THcHodoscope("hod", "Hodoscope");
SHMS->AddDetector(phod);
// Add Heavy Gas Cherenkov to SHMS apparatus
THcCherenkov* phgcer = new THcCherenkov("hgcer", "Heavy Gas Cherenkov");
SHMS->AddDetector(phgcer);
// Add Aerogel Cherenkov to SHMS apparatus
THcAerogel* paero = new THcAerogel("aero", "Aerogel");
SHMS->AddDetector(paero);
// Add calorimeter to SHMS apparatus
THcShower* pcal = new THcShower("cal", "Calorimeter");
SHMS->AddDetector(pcal);
// THcBCMCurrent* hbc = new THcBCMCurrent("H.bcm", "BCM current check");
// gHaPhysics->Add(hbc);
// Add rastered beam apparatus
THaApparatus* pbeam = new THcRasteredBeam("P.rb", "Rastered Beamline");
gHaApps->Add(pbeam);
// Add physics modules
// Calculate reaction point
THcReactionPoint* prp = new THcReactionPoint("P.react", "SHMS reaction point", "P", "P.rb");
gHaPhysics->Add(prp);
// Calculate extended target corrections
THcExtTarCor* pext = new THcExtTarCor("P.extcor", "HMS extended target corrections", "P", "P.react");
gHaPhysics->Add(pext);
// Calculate golden track quantites
THaGoldenTrack* pgtr = new THaGoldenTrack("P.gtr", "SHMS Golden Track", "P");
gHaPhysics->Add(pgtr);
// Calculate the hodoscope efficiencies
THcHodoEff* peff = new THcHodoEff("phodeff", "SHMS hodo efficiency", "P.hod");
gHaPhysics->Add(peff);
// Add event handler for scaler events
THcScalerEvtHandler* pscaler = new THcScalerEvtHandler("P", "Hall C scaler event type 1");
pscaler->AddEvtType(1);
pscaler->AddEvtType(4);
pscaler->AddEvtType(5);
pscaler->AddEvtType(6);
pscaler->AddEvtType(7);
pscaler->AddEvtType(129);
pscaler->SetDelayedType(129);
pscaler->SetUseFirstEvent(kTRUE);
gHaEvtHandlers->Add(pscaler);
//=:=:=
// HMS
//=:=:=
// Set up the equipment to be analyzed.
THcHallCSpectrometer* HMS = new THcHallCSpectrometer("H", "HMS");
HMS->SetEvtType(2);
HMS->AddEvtType(4);
HMS->AddEvtType(5);
HMS->AddEvtType(6);
HMS->AddEvtType(7);
gHaApps->Add(HMS);
// Add drift chambers to HMS apparatus
THcDC* hdc = new THcDC("dc", "Drift Chambers");
HMS->AddDetector(hdc);
// Add hodoscope to HMS apparatus
THcHodoscope* hhod = new THcHodoscope("hod", "Hodoscope");
HMS->AddDetector(hhod);
// Add Cherenkov to HMS apparatus
THcCherenkov* hcer = new THcCherenkov("cer", "Heavy Gas Cherenkov");
HMS->AddDetector(hcer);
// Add Aerogel Cherenkov to HMS apparatus
// THcAerogel* haero = new THcAerogel("aero", "Aerogel");
// HMS->AddDetector(haero);
// Add calorimeter to HMS apparatus
THcShower* hcal = new THcShower("cal", "Calorimeter");
HMS->AddDetector(hcal);
// Add rastered beam apparatus
THaApparatus* hbeam = new THcRasteredBeam("H.rb", "Rastered Beamline");
gHaApps->Add(hbeam);
// Add physics modules
// Calculate reaction point
THcReactionPoint* hrp = new THcReactionPoint("H.react", "HMS reaction point", "H", "H.rb");
gHaPhysics->Add(hrp);
// Calculate extended target corrections
THcExtTarCor* hext = new THcExtTarCor("H.extcor", "HMS extended target corrections", "H", "H.react");
gHaPhysics->Add(hext);
// Calculate golden track quantities
THaGoldenTrack* hgtr = new THaGoldenTrack("H.gtr", "HMS Golden Track", "H");
gHaPhysics->Add(hgtr);
// Calculate the hodoscope efficiencies
THcHodoEff* heff = new THcHodoEff("hhodeff", "HMS hodo efficiency", "H.hod");
gHaPhysics->Add(heff);
// Add event handler for scaler events
THcScalerEvtHandler *hscaler = new THcScalerEvtHandler("H", "Hall C scaler event type 4");
hscaler->AddEvtType(2);
hscaler->AddEvtType(4);
hscaler->AddEvtType(5);
hscaler->AddEvtType(6);
hscaler->AddEvtType(7);
hscaler->AddEvtType(129);
hscaler->SetDelayedType(129);
hscaler->SetUseFirstEvent(kTRUE);
gHaEvtHandlers->Add(hscaler);
//=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=
// Kinematics Modules
//=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=
// Add Physics Module to calculate primary (scattered electrons) beam kinematics
THcPrimaryKine* hkin_primary = new THcPrimaryKine("H.kin.primary", "HMS Single Arm Kinematics", "H", "H.rb");
gHaPhysics->Add(hkin_primary);
// Add Physics Module to calculate secondary (scattered hadrons) beam kinematics
THcSecondaryKine* pkin_secondary = new THcSecondaryKine("P.kin.secondary", "SHMS Single Arm Kinematics", "P", "H.kin.primary");
gHaPhysics->Add(pkin_secondary);
//=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=
// Global Objects & Event Handlers
//=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=
// Add trigger apparatus
THaApparatus* TRG = new THcTrigApp("T", "TRG");
gHaApps->Add(TRG);
// Add trigger detector to trigger apparatus
THcTrigDet* coin = new THcTrigDet("coin", "Coincidence Trigger Information");
// Suppress missing reference time warnings for these event types
coin->SetEvtType(1);
coin->AddEvtType(2);
TRG->AddDetector(coin);
//Add coin physics module THcCoinTime::THcCoinTime (const char *name, const char* description, const char* hadArmName,
// const char* elecArmName, const char* coinname) :
THcCoinTime* coinTime = new THcCoinTime("CTime", "Coincidende Time Determination", "P", "H", "T.coin");
gHaPhysics->Add(coinTime);
// Add event handler for prestart event 125.
THcConfigEvtHandler* ev125 = new THcConfigEvtHandler("HC", "Config Event type 125");
gHaEvtHandlers->Add(ev125);
// Add event handler for EPICS events
THaEpicsEvtHandler* hcepics = new THaEpicsEvtHandler("epics", "HC EPICS event type 180");
gHaEvtHandlers->Add(hcepics);
// Set up the analyzer - we use the standard one,
// but this could be an experiment-specific one as well.
// The Analyzer controls the reading of the data, executes
// tests/cuts, loops over Acpparatus's and PhysicsModules,
// and executes the output routines.
THcAnalyzer* analyzer = new THcAnalyzer;
// A simple event class to be output to the resulting tree.
// Creating your own descendant of THaEvent is one way of
// defining and controlling the output.
THaEvent* event = new THaEvent;
// Define the run(s) that we want to analyze.
// We just set up one, but this could be many.
THcRun* run = new THcRun( pathList, Form(RunFileNamePattern, RunNumber) );
// Set to read in Hall C run database parameters
run->SetRunParamClass("THcRunParameters");
// Eventually need to learn to skip over, or properly analyze the pedestal events
run->SetEventRange(1, MaxEvent); // Physics Event number, does not include scaler or control events.
run->SetNscan(1);
run->SetDataRequired(0x7);
run->Print();
// Define the analysis parameters
TString ROOTFileName = Form(ROOTFileNamePattern, RunNumber, MaxEvent);
analyzer->SetCountMode(2); // 0 = counter is # of physics triggers
// 1 = counter is # of all decode reads
// 2 = counter is event number
analyzer->SetEvent(event);
// Set EPICS event type
analyzer->SetEpicsEvtType(180);
// Define crate map
analyzer->SetCrateMapFileName("MAPS/db_cratemap.dat");
// Define output ROOT file
analyzer->SetOutFile(ROOTFileName.Data());
// Define DEF-file+
analyzer->SetOdefFile("DEF-files/COIN/PRODUCTION/coin_production_hElec_pProt.def");
// Define cuts file
analyzer->SetCutFile("DEF-files/COIN/PRODUCTION/CUTS/coin_production_cuts.def"); // optional
// File to record accounting information for cuts
analyzer->SetSummaryFile(Form("REPORT_OUTPUT/COIN/PRODUCTION/summary_production_%d_%d.report", RunNumber, MaxEvent)); // optional
// Start the actual analysis.
analyzer->Process(run);
// Create report file from template
analyzer->PrintReport("TEMPLATES/COIN/PRODUCTION/coin_production.template",
Form("REPORT_OUTPUT/COIN/PRODUCTION/replay_coin_production_%d_%d.report", RunNumber, MaxEvent)); // optional
}
#include "nlohmann/json.hpp"
#include <cmath>
#include <iostream>
#include "ROOT/RDataFrame.hxx"
#include "ROOT/RVec.hxx"
#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
#include "Math/VectorUtil.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TStyle.h"
#include "TSystem.h"
R__LOAD_LIBRARY(libMathMore.so)
R__LOAD_LIBRARY(libGenVector.so)
R__LOAD_LIBRARY(libfmt.so)
#include "fmt/core.h"
#include "THStack.h"
#ifdef __cpp_lib_filesystem
#include <filesystem>
namespace fs = std::filesystem;
#else
#include <experimental/filesystem>
namespace fs = std::experimental::filesystem;
#endif
using RDFNode = decltype(ROOT::RDataFrame{0}.Filter(""));
using Histo1DProxy =
decltype(ROOT::RDataFrame{0}.Histo1D(ROOT::RDF::TH1DModel{"", "", 128u, 0., 0.}, ""));
struct RDFInfo {
RDFNode& df;
const std::string title;
RDFInfo(RDFNode& df, std::string_view title) : df{df}, title{title} {}
};
constexpr const double M_P = .938272;
constexpr const double M_P2 = M_P * M_P;
constexpr const double M_pion = 0.139;
constexpr const double M_pion2 = M_pion * M_pion;
constexpr const double M_e = .000511;
// =================================================================================
// Cuts
// =================================================================================
std::string goodTrackSHMS = "P.gtr.dp > -10 && P.gtr.dp < 22";
std::string goodTrackHMS = "H.gtr.dp > -8 && H.gtr.dp < 8";
std::string piCutSHMS = "P.cal.etottracknorm<1.0";
//std::string piCutSHMS = "P.aero.npeSum > 1.0 && P.cal.eprtracknorm < 0.2 && P.cal.etottracknorm<1.0";
std::string eCutHMS = "H.cal.etottracknorm > 0.50 && H.cal.etottracknorm < 2. && "
"H.cer.npeSum > 1.";
std::string epiCut = "P.aero.npeSum > 1.0 && P.cal.eprtracknorm < 0.2 && "
"H.cer.npeSum > 1.0 && H.cal.etottracknorm > 0.6 && "
"H.cal.etottracknorm < 2.0 && P.cal.etottracknorm<1.0";
using Pvec3D = ROOT::Math::XYZVector;
using Pvec4D = ROOT::Math::PxPyPzMVector;
// =================================================================================
// reconstruction
// =================================================================================
auto p_pion = [](double px, double py, double pz) {
return Pvec4D{px * 0.996, py * 0.996, pz * 0.996, M_e};
};
auto p_electron = [](double px, double py, double pz) {
return Pvec4D{px * 0.994, py * 0.994, pz * 0.994, M_e};
};
auto t = [](const double Egamma, Pvec4D& jpsi) {
Pvec4D beam{0, 0, Egamma, 0};
return (beam - jpsi).M2();
};
bool root_file_exists(std::string rootfile) {
bool found_good_file = false;
if (!gSystem->AccessPathName(rootfile.c_str())) {
TFile file(rootfile.c_str());
if (file.IsZombie()) {
std::cout << rootfile << " is a zombie.\n";
std::cout
<< " Did your replay finish? Check that the it is done before running this script.\n";
return false;
// return;
} else {
std::cout << " using : " << rootfile << "\n";
return true;
}
}
return false;
}
void elastic_test(int RunNumber = 6012, int nevents = 50000, int prompt = 0, int update = 0,
int default_count_goal = 10000, int redo_timing = 0) {
// ===============================================================================================
// Initialization
// ===============================================================================================
using nlohmann::json;
json j;
{
std::ifstream json_input_file("db2/run_list_coin.json");
try {
json_input_file >> j;
} catch (json::parse_error) {
std::cerr << "error: json file, db2/run_list.json, is incomplete or has broken syntax.\n";
std::quick_exit(-127);
}
}
auto runnum_str = std::to_string(RunNumber);
if (j.find(runnum_str) == j.end()) {
std::cout << "Run " << RunNumber << " not found in db2/run_list_coin.json\n";
std::cout << "Check that run number and replay exists. \n";
std::cout << "If problem persists please contact Sylvester (217-848-0565)\n";
}
double P0_shms_setting = j[runnum_str]["spectrometers"]["shms_momentum"].get<double>();
double P0_shms = std::abs(P0_shms_setting);
bool found_good_file = false;
std::string rootfile =
fmt::format("full_online/coin_replay_production_{}_{}.root", RunNumber, nevents);
found_good_file = root_file_exists(rootfile.c_str());
if (!found_good_file) {
rootfile =
fmt::format("ROOTfiles_volatile/coin_replay_production_{}_{}.root", RunNumber, nevents);
found_good_file = root_file_exists(rootfile.c_str());
}
if (!found_good_file) {
rootfile = fmt::format("ROOTfiles_csv/coin_replay_production_{}_{}.root", RunNumber, nevents);
found_good_file = root_file_exists(rootfile.c_str());
}
if (!found_good_file) {
rootfile = fmt::format("ROOTfiles/coin_replay_production_{}_{}.root", RunNumber, nevents);
found_good_file = root_file_exists(rootfile.c_str());
}
if (!found_good_file) {
std::cout << " Error: suitable root file not found\n";
return;
}
// ===============================================================================================
// Dataframe
// ===============================================================================================
ROOT::EnableImplicitMT(24);
//---------------------------------------------------------------------------
// Detector tree
ROOT::RDataFrame d("T", rootfile);
//// SHMS Scaler tree
//ROOT::RDataFrame d_sh("TSP", rootfile);
//// int N_scaler_events = *(d_sh.Count());
auto d_coin = d.Filter("fEvtHdr.fEvtType == 4");
// Good track cuts
auto dHMSGoodTrack = d_coin.Filter(goodTrackHMS);
auto dSHMSGoodTrack = d_coin.Filter(goodTrackSHMS);
auto dCOINGoodTrack = dHMSGoodTrack.Filter(goodTrackSHMS)
.Define("p_electron", p_electron, {"H.gtr.px", "H.gtr.py", "H.gtr.pz"})
.Define("p_pion", p_pion, {"P.gtr.px", "P.gtr.py", "P.gtr.pz"});
// PID cuts
auto dHMSEl = dHMSGoodTrack.Filter(eCutHMS);
auto dSHMSEl = dSHMSGoodTrack.Filter(piCutSHMS);
auto dCOINEl = dCOINGoodTrack.Filter(eCutHMS + " && " + piCutSHMS);
//.Filter(
// [=](double npe, double dp) {
// double p_track = P0_shms * (100.0 + dp) / 100.0;
// // no cerenkov cut needed when momentum is below 2.8 GeV/c
// if (p_track < 2.8) {
// return true;
// }
// return npe > 1.0;
// },
// {"P.hgcer.npeSum", "P.gtr.dp"});
// Timing cuts
// Find the timing peak
// Find the coin peak
double coin_peak_center = 0;
if (redo_timing) {
auto h_coin_time =
dCOINEl.Histo1D({"coin_time", "coin_time", 8000, 0, 1000}, "CTime.ePositronCoinTime_ROC2");
h_coin_time->DrawClone();
int coin_peak_bin = h_coin_time->GetMaximumBin();
coin_peak_center = h_coin_time->GetBinCenter(coin_peak_bin);
std::cout << "COINCIDENCE time peak found at: " << coin_peak_center << std::endl;
} else {
//coin_peak_center = 43.0; // run 7240-7241
coin_peak_center = 23.0; // run 6012
std::cout << "COINCIDENCE time peak: using pre-calculated value at: " << coin_peak_center
<< std::endl;
;
}
// timing cut lambdas
// TODO: evaluate timing cut and offset for random background
auto timing_cut = [=](double coin_time) { return std::abs(coin_time - coin_peak_center) < 2.; };
// anti-timing set to 5x width of regular
auto anti_timing_cut = [=](double coin_time) {
return std::abs(coin_time - coin_peak_center - 28.) < 10.;
};
// timing counts
auto dHMSElInTime = dHMSEl.Filter(timing_cut, {"CTime.ePositronCoinTime_ROC2"});
auto dHMSElRandom = dHMSEl.Filter(anti_timing_cut, {"CTime.ePositronCoinTime_ROC2"});
auto dSHMSElInTime = dSHMSEl.Filter(timing_cut, {"CTime.ePositronCoinTime_ROC2"});
auto dSHMSElRandom = dSHMSEl.Filter(anti_timing_cut, {"CTime.ePositronCoinTime_ROC2"});
auto dCOINElInTime = dCOINEl.Filter(timing_cut, {"CTime.ePiCoinTime_ROC2"});
auto dCOINElRandom = dCOINEl.Filter(anti_timing_cut, {"CTime.ePiCoinTime_ROC2"});
// Output root file
//auto out_file =
// new TFile(fmt::format("monitoring/{}/good_csv_counter.root", RunNumber).c_str(), "UPDATE");
//out_file->cd();
// =========================================================================================
// Histograms
// =========================================================================================
// 2D correlations
auto hTheta2DNoCuts = d_coin.Histo2D(
{"theta2D", "No cuts;#theta_{SHMS};#theta_{HMS};#counts", 50, -.1, .1, 50, -.1, .1},
"P.gtr.th", "H.gtr.th");
auto hTheta2DTracking = dCOINGoodTrack.Histo2D(
{"theta2D", "Cuts: tracking;#theta_{SHMS};#theta_{HMS};#counts", 50, -.1, .1, 50, -.1, .1},
"P.gtr.th", "H.gtr.th");
auto hTheta2DPID =
dCOINEl.Histo2D({"theta2D", "Cuts: tracking+PID;#theta_{SHMS};#theta_{HMS};#counts", 50, -.1,
.1, 50, -.1, .1},
"P.gtr.th", "H.gtr.th");
auto hTheta2DTiming =
dCOINElInTime.Histo2D({"theta2D", "Cuts: tracking+PID;#theta_{SHMS};#theta_{HMS};#counts", 50,
-.1, .1, 50, -.1, .1},
"P.gtr.th", "H.gtr.th");
// timing
auto hCoinTimeNoCuts =
d_coin.Histo1D({"coin_time.NoCuts", "No Cuts;coin_time;counts", 8000, 0, 1000},
"CTime.ePositronCoinTime_ROC2");
auto hCoinTimeTracking = dCOINGoodTrack.Histo1D(
{"coin_time.Tracking", "Cuts: Tracking;coin_time;counts", 8000, 0, 1000},
"CTime.ePositronCoinTime_ROC2");
auto hCoinTimePID =
dCOINEl.Histo1D({"coin_time.PID", "Cuts: Tracking+PID;coin_time;counts", 8000, 0, 1000},
"CTime.ePositronCoinTime_ROC2");
auto hCoinTimeTiming = dCOINElInTime.Histo1D(
{"coin_time.Timing", "Cuts: Tracking+PID+Timing;coin_time;counts", 8000, 0, 1000},
"CTime.ePositronCoinTime_ROC2");
auto hRandCoinTimePID = dCOINElRandom.Histo1D(
{"rand_coin_time.PID", "Cuts: Tracking+PID;coin_time;counts", 8000, 0, 1000},
"CTime.ePositronCoinTime_ROC2");
// P.gtr.dp
auto hPdpNoCuts =
d_coin.Histo1D({"P.gtr.dp.NoCuts", "No Cuts;#deltap [%];counts", 200, -30, 40}, "P.gtr.dp");
auto hPdpTracking = dSHMSGoodTrack.Histo1D(
{"P.gtr.dp.Tracking", "Cuts: Tracking;#deltap [%];counts", 200, -30, 40}, "P.gtr.dp");
auto hPdpPID = dSHMSEl.Histo1D(
{"P.gtr.dp.PID", "Cuts: Tracking+PID;#deltap [%];counts", 200, -30, 40}, "P.gtr.dp");
auto hPdpTiming = dSHMSElInTime.Histo1D(
{"P.gtr.dp.Timing", "Cuts: Tracking+PID+Timing;#deltap [%];counts", 200, -30, 40},
"P.gtr.dp");
// P.gtr.th
auto hPthNoCuts = d_coin.Histo1D(
{"P.gtr.th.NoCuts", "No Cuts;#theta_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.th");
auto hPthTracking = dSHMSGoodTrack.Histo1D(
{"P.gtr.th.Tracking", "Cuts: Tracking;#theta_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.th");
auto hPthPID = dSHMSEl.Histo1D(
{"P.gtr.th.PID", "Cuts: Tracking+PID;#theta_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.th");
auto hPthTiming = dSHMSElInTime.Histo1D(
{"P.gtr.th.Timing", "Cuts: Tracking+PID+Timing;#theta_{SHMS};counts", 200, -0.1, 0.1},
"P.gtr.th");
// P.gtr.ph
auto hPphNoCuts =
d_coin.Histo1D({"P.gtr.ph.NoCuts", "No Cuts;#phi_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.ph");
auto hPphTracking = dSHMSGoodTrack.Histo1D(
{"P.gtr.ph.Tracking", "Cuts: Tracking;#phi_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.ph");
auto hPphPID = dSHMSEl.Histo1D(
{"P.gtr.ph.PID", "Cuts: Tracking+PID;#phi_{SHMS};counts", 200, -0.1, 0.1}, "P.gtr.ph");
auto hPphTiming = dSHMSElInTime.Histo1D(
{"P.gtr.ph.Timing", "Cuts: Tracking+PID+Timing;#phi_{SHMS};counts", 200, -0.1, 0.1},
"P.gtr.ph");
// P.gtr.y
auto hPyNoCuts =
d_coin.Histo1D({"P.gtr.y.NoCuts", "No Cuts;ytar;counts", 200, -10., 10.}, "P.gtr.y");
auto hPyTracking = dSHMSGoodTrack.Histo1D(
{"P.gtr.y.Tracking", "Cuts: Tracking;ytar;counts", 200, -10., 10.}, "P.gtr.y");
auto hPyPID =
dSHMSEl.Histo1D({"P.gtr.y.PID", "Cuts: Tracking+PID;ytar;counts", 200, -10., 10.}, "P.gtr.y");
auto hPyTiming = dSHMSElInTime.Histo1D(
{"P.gtr.y.Timing", "Cuts: Tracking+PID+Timing;ytar;counts", 200, -10., 10.}, "P.gtr.y");
// P.cal.etottracknorm
auto hPcalEPNoCuts =
d_coin.Histo1D({"P.cal.etottracknorm.NoCuts", "No Cuts;SHMS E/P;counts", 200, -.5, 1.5},
"P.cal.etottracknorm");
auto hPcalEPTracking = dSHMSGoodTrack.Histo1D(
{"P.cal.etottracknorm.Tracking", "Cuts: Tracking;SHMS E/P;counts", 200, -.5, 1.5},
"P.cal.etottracknorm");
auto hPcalEPPID = dSHMSEl.Histo1D(
{"P.cal.etottracknorm.PID", "Cuts: Tracking+PID;SHMS E/P;counts", 200, -.5, 1.5},
"P.cal.etottracknorm");
auto hPcalEPAll = dCOINElInTime.Histo1D(
{"P.cal.etottracknorm.All", "Cuts: Tracking+PID+Coincidence;SHMS E/P;counts", 200, -.5, 1.5},
"P.cal.etottracknorm");
// P.ngcer.npeSum
auto hPcerNpheNoCuts = d_coin.Histo1D(
{"P.ngcer.npeSum.NoCuts", "No Cuts;SHMS NGC #phe;counts", 200, -5, 76}, "P.ngcer.npeSum");
auto hPcerNpheTracking = dSHMSGoodTrack.Histo1D(
{"P.ngcer.npeSum.Tracking", "Cuts: Tracking;SHMS NGC #phe;counts", 200, -5, 76},
"P.ngcer.npeSum");
auto hPcerNphePID = dSHMSEl.Histo1D(
{"P.ngcer.npeSum.PID", "Cuts: Tracking+PID;SHMS NGC #phe;counts", 200, -5, 76},
"P.ngcer.npeSum");
auto hPcerNpheAll = dCOINElInTime.Histo1D(
{"P.ngcer.npeSum.All", "Cuts: Tracking+PID+Coincidence;SHMS NGC #phe;counts", 200, -5, 76},
"P.ngcer.npeSum");
// P.hgcer.npeSum
auto hPhgcerNpheNoCuts = d_coin.Histo1D(
{"P.hgcer.npeSum.NoCuts", "No Cuts;SHMS HGC #phe;counts", 200, -5, 76}, "P.hgcer.npeSum");
auto hPhgcerNpheTracking = dSHMSGoodTrack.Histo1D(
{"P.hgcer.npeSum.Tracking", "Cuts: Tracking;SHMS HGC #phe;counts", 200, -5, 76},
"P.hgcer.npeSum");
auto hPhgcerNphePID = dSHMSEl.Histo1D(
{"P.hgcer.npeSum.PID", "Cuts: Tracking+PID;SHMS HGC #phe;counts", 200, -5, 76},
"P.hgcer.npeSum");
auto hPhgcerNpheAll = dCOINElInTime.Histo1D(
{"P.hgcer.npeSum.All", "Cuts: Tracking+PID+Coincidence;SHMS HGC #phe;counts", 200, -5, 76},
"P.hgcer.npeSum");
// H.cal.etottracknorm
auto hHcalEPNoCuts =
d_coin.Histo1D({"H.cal.etottracknorm.NoCuts", "No Cuts;HMS E/P;counts", 200, -.5, 1.5},
"H.cal.etottracknorm");
auto hHcalEPTracking = dHMSGoodTrack.Histo1D(
{"H.cal.etottracknorm.Tracking", "Cuts: Tracking;HMS E/P;counts", 200, -.5, 1.5},
"H.cal.etottracknorm");
auto hHcalEPPID = dHMSEl.Histo1D(
{"H.cal.etottracknorm.PID", "Cuts: Tracking+PID;HMS E/P;counts", 200, -.5, 1.5},
"H.cal.etottracknorm");
auto hHcalEPAll = dCOINElInTime.Histo1D(
{"H.cal.etottracknorm.All", "Cuts: Tracking+PID+Coincidence;HMS E/P;counts", 200, -.5, 1.5},
"H.cal.etottracknorm");
// H.cer.npeSum
auto hHcerNpheNoCuts = d_coin.Histo1D(
{"H.cer.npeSum.NoCuts", "No Cuts;HMS Cer #phe;counts", 200, -1, 15}, "H.cer.npeSum");
auto hHcerNpheTracking = dHMSGoodTrack.Histo1D(
{"H.cer.npeSum.Tracking", "Cuts: Tracking;HMS Cer #phe;counts", 200, -1, 15}, "H.cer.npeSum");
auto hHcerNphePID = dHMSEl.Histo1D(
{"H.cer.npeSum.PID", "Cuts: Tracking+PID+Coincidence;HMS Cer #phe;counts", 200, -1, 15},
"H.cer.npeSum");
auto hHcerNpheAll = dCOINElInTime.Histo1D(
{"H.cer.npeSum.PID", "Cuts: Tracking+PID+Coincidence;HMS Cer #phe;counts", 200, -1, 15},
"H.cer.npeSum");
// H.gtr.dp
auto hHdpNoCuts =
d_coin.Histo1D({"H.gtr.dp.NoCuts", "No Cuts;#deltap [%];counts", 200, -30, 40}, "H.gtr.dp");
auto hHdpTracking = dHMSGoodTrack.Histo1D(
{"H.gtr.dp.Tracking", "Cuts: Tracking;#deltap [%];counts", 200, -30, 40}, "H.gtr.dp");
auto hHdpPID = dHMSEl.Histo1D(
{"H.gtr.dp.PID", "Cuts: Tracking+PID;#deltap [%];counts", 200, -30, 40}, "H.gtr.dp");
auto hHdpTiming = dHMSElInTime.Histo1D(
{"H.gtr.dp.Timing", "Cuts: Tracking+PID+Timing;#deltap [%];counts", 200, -30, 40},
"H.gtr.dp");
// H.gtr.th
auto hHthNoCuts = d_coin.Histo1D(
{"H.gtr.th.NoCuts", "No Cuts;#theta_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.th");
auto hHthTracking = dHMSGoodTrack.Histo1D(
{"H.gtr.th.Tracking", "Cuts: Tracking;#theta_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.th");
auto hHthPID = dHMSEl.Histo1D(
{"H.gtr.th.PID", "Cuts: Tracking+PID;#theta_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.th");
auto hHthTiming = dHMSElInTime.Histo1D(
{"H.gtr.th.Timing", "Cuts: Tracking+PID+Timing;#theta_{HMS};counts", 200, -0.1, 0.1},
"H.gtr.th");
// H.gtr.ph
auto hHphNoCuts =
d_coin.Histo1D({"H.gtr.ph.NoCuts", "No Cuts;#phi_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.ph");
auto hHphTracking = dHMSGoodTrack.Histo1D(
{"H.gtr.ph.Tracking", "Cuts: Tracking;#phi_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.ph");
auto hHphPID = dHMSEl.Histo1D(
{"H.gtr.ph.PID", "Cuts: Tracking+PID;#phi_{HMS};counts", 200, -0.1, 0.1}, "H.gtr.ph");
auto hHphTiming = dHMSElInTime.Histo1D(
{"H.gtr.ph.Timing", "Cuts: Tracking+PID+Timing;#phi_{HMS};counts", 200, -0.1, 0.1},
"H.gtr.ph");
// H.gtr.y
auto hHyNoCuts =
d_coin.Histo1D({"H.gtr.y.NoCuts", "No Cuts;ytar;counts", 200, -10., 10.}, "H.gtr.y");
auto hHyTracking = dHMSGoodTrack.Histo1D(
{"H.gtr.y.Tracking", "Cuts: Tracking;ytar;counts", 200, -10., 10.}, "H.gtr.y");
auto hHyPID =
dHMSEl.Histo1D({"H.gtr.y.PID", "Cuts: Tracking+PID;ytar;counts", 200, -10., 10.}, "H.gtr.y");
auto hHyTiming = dHMSElInTime.Histo1D(
{"H.gtr.y.Timing", "Cuts: Tracking+PID+Timing;ytar;counts", 200, -10., 10.}, "H.gtr.y");
// scalers
//auto total_charge = d_sh.Max("P.BCM4B.scalerChargeCut");
//auto shms_el_real_scaler = d_sh.Max("P.pEL_REAL.scaler");
//auto hms_el_real_scaler = d_sh.Max("P.hEL_REAL.scaler");
//auto time_1MHz = d_sh.Max("P.1MHz.scalerTime");
//auto time_1MHz_cut = d_sh.Max("P.1MHz.scalerTimeCut");
auto yield_all = d.Count();
// 5 timing cut widths worth of random backgrounds
auto yield_coin = d_coin.Count();
auto yield_HMSGoodTrack = dHMSGoodTrack.Count();
auto yield_SHMSGoodTrack = dSHMSGoodTrack.Count();
auto yield_COINGoodTrack = dCOINGoodTrack.Count();
auto yield_HMSEl = dHMSEl.Count();
auto yield_SHMSEl = dSHMSEl.Count();
auto yield_COINEl = dCOINEl.Count();
auto yield_HMSElInTime = dHMSElInTime.Count();
auto yield_HMSElRandom = dHMSElRandom.Count();
auto yield_SHMSElInTime = dSHMSElInTime.Count();
auto yield_SHMSElRandom = dSHMSElRandom.Count();
auto yield_COINElInTime = dCOINElInTime.Count();
auto yield_COINElRandom = dCOINElRandom.Count();
auto yield_coin_raw = dCOINElInTime.Count();
auto yield_coin_random = dCOINElRandom.Count();
// -------------------------------------
// End lazy eval
// -------------------------------------
auto n_coin_good = *yield_coin_raw - *yield_coin_random / 5.;
auto n_HMSElGood = *yield_HMSElInTime - *yield_HMSElRandom / 5;
auto n_SHMSElGood = *yield_SHMSElInTime - *yield_SHMSElRandom / 5;
auto n_COINElGood = *yield_COINElInTime - *yield_COINElRandom / 5;
hPdpNoCuts->DrawCopy();
//std::cout << " coin COUNTS : " << *(d_coin.Count()) << "\n";
//std::cout << " yield_HMSEl : " << *(yield_HMSEl) << "\n";
std::cout << " yield_COINEl : " << *(yield_COINEl) << "\n";
//std::cout << " ALL COUNTS : " << *yield_all << "\n";
//std::cout << " GOOD COUNTS : " << n_COINElGood << "\n";
//
if( 4 != (*(yield_COINEl)) ){
std::exit(-1);
}
std::exit(0);
}
#!/bin/bash
#echo "This is the elastic testing..."
#echo " "
#echo "There are currently 0 tests to run!"
#which hcana
#
#ls -lrth
#ls -lrth build
#
#git clone git@eicweb.phy.anl.gov:jlab/hallc/exp/CSV/hallc_replay_csv.git
#git clone git@eicweb.phy.anl.gov:jlab/hallc/exp/CSV/online_csv.git
#cd online_csv
#ln -s ../hallc_reaply_csv/PARAM
## and the reset
#
#mkdir raw
#pushd raw
# wget coin.dat
#popd
singularity help build/Singularity.hcana.simg
singularity exec build/Singularity.hcana.simg hcana tests/elastic_test.cxx
#!/bin/bash
echo "This is the elastic testing..."
echo " "
echo "There are currently 0 tests to run!"
which hcana
ls -lrth
ls -lrth build
git clone git@eicweb.phy.anl.gov:jlab/hallc/exp/CSV/hallc_replay_csv.git
git clone git@eicweb.phy.anl.gov:jlab/hallc/exp/CSV/online_csv.git
cd online_csv
ln -s ../hallc_reaply_csv/PARAM
# and the reset
mkdir raw
pushd raw
wget coin.dat
popd
singularity help build/Singularity.hcana.simg
singularity exec build/Singularity.hcana.simg which hcana
singularity exec build/Singularity.hcana.simg hcana tests/my_root_script.cxx
echo " WOOOO"
void my_root_script() {
std::cout << "Hello from my_root_script.cxx!\n";
std::cout << "This should be run with singularity\n";
double pi = 3.14;
auto pi_equals_3 = (3 == pi);
std::cout << " pi_equals_3 = " << pi_equals_3 << "\n";
if( pi_equals_3) {
std::cout << "what the hell?\n";
std::exit( 0 );
}
/* else */
std::exit( -1 );
}
#!/bin/bash
echo "This is the elastic testing..."
echo " "
echo "There are currently 0 tests to run!"
which hcana
ls -lrth
ls -lrth build
mkdir -p ROOTfiles
git clone https://eicweb.phy.anl.gov/whit/ci_test_data.git
git clone https://eicweb.phy.anl.gov/jlab/hallc/exp/CSV/hallc_replay_csv.git
git clone https://eicweb.phy.anl.gov/jlab/hallc/exp/CSV/online_csv.git
cd online_csv
mkdir -p logs
ln -s ../hallc_replay_csv/PARAM
ln -s ../hallc_replay_csv/DBASE
ln -s ../hallc_replay_csv/CALIBRATION
ln -s ../hallc_replay_csv/DEF-files
ln -s ../hallc_replay_csv/MAPS
ln -s ../hallc_replay_csv/SCRIPTS
ln -s ../hallc_replay_csv/DATFILES
ln -s ../ci_test_data/raw
ln -s ../ROOTfiles
# and the reset
ls -lrth
ls -lrth raw/
ls -lrth ROOTfiles/
pwd
# run replay script
df -h
singularity exec ../build/Singularity.hcana.simg hcana -b -q "../tests/elastic_coin_replay.cxx+(6012,50000)"
echo "hcana calls... the coin replay script and outputs blah.root"
file(GLOB APP_SOURCES RELATIVE ${CMAKE_CURRENT_LIST_DIR} "*.cpp")
foreach(exe_src ${APP_SOURCES})
string(REPLACE ".cpp" "" exe ${exe_src})
add_executable(${exe} ${exe_src})
target_include_directories(${exe}
PUBLIC #$<INSTALL_INTERFACE:include>
$<INSTALL_INTERFACE:include>
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include/hcana>
$<BUILD_INTERFACE:${Podd_DIR}/../../include/podd2>
$<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}> # for hc_compiledata.h
$<BUILD_INTERFACE:${SPDLOG_INCLUDE_DIR}>
$<INSTALL_INTERFACE:${SPDLOG_INCLUDE_DIR}>
$<BUILD_INTERFACE:${FMT_INCLUDE_DIR}>
$<BUILD_INTERFACE:${CODA_ET_INCLUDE_DIR}>
$<BUILD_INTERFACE:${EVIO_INCLUDE_DIR}>
)
target_link_libraries(${exe}
PUBLIC
${PROJECT_NAME}::HallC
Podd::Podd
Podd::Decode
coda_et::coda_et
EVIO::evioxx
)
install(TARGETS ${exe} DESTINATION ${CMAKE_INSTALL_BINDIR})
endforeach(exe_src ${APP_SOURCES})
/*----------------------------------------------------------------------------*
* Copyright (c) 1998 Southeastern Universities Research Association, *
* Thomas Jefferson National Accelerator Facility *
* *
* This software was developed under a United States Government license *
* described in the NOTICE file included as part of this distribution. *
* *
* Author: Carl Timmer *
* timmer@jlab.org Jefferson Lab, MS-12B3 *
* Phone: (757) 269-5130 12000 Jefferson Ave. *
* Fax: (757) 269-6248 Newport News, VA 23606 *
* *
*----------------------------------------------------------------------------*
*
* Description:
* ET system sample event consumer
*
*----------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <signal.h>
#include <sys/time.h>
#include <getopt.h>
#include <limits.h>
#include <time.h>
#include "et.h"
/* prototype */
static void *signal_thread (void *arg);
int main(int argc,char **argv) {
int i, j, c, i_tmp, status, numRead, locality;
int flowMode=ET_STATION_SERIAL, position=ET_END, pposition=ET_END;
int errflg=0, chunk=1, qSize=0, verbose=0, remote=0, blocking=1, dump=0, readData=0;
int multicast=0, broadcast=0, broadAndMulticast=0;
int con[ET_STATION_SELECT_INTS];
int sendBufSize=0, recvBufSize=0, noDelay=0;
int debugLevel = ET_DEBUG_ERROR;
unsigned short port=0;
char stationName[ET_STATNAME_LENGTH], et_name[ET_FILENAME_LENGTH], host[256], interface[16];
char localAddr[16];
int mcastAddrCount = 0, mcastAddrMax = 10;
char mcastAddr[mcastAddrMax][16];
pthread_t tid;
et_att_id attach1;
et_stat_id my_stat;
et_sys_id id;
et_statconfig sconfig;
et_event **pe;
et_openconfig openconfig;
sigset_t sigblock;
struct timespec timeout;
#if defined __APPLE__
struct timeval t1, t2;
#else
struct timespec t1, t2;
#endif
/* statistics variables */
double rate=0.0, avgRate=0.0;
int64_t count=0, totalCount=0, totalT=0, time, time1, time2, bytes=0, totalBytes=0;
/* 4 multiple character command-line options */
static struct option long_options[] =
{ {"host", 1, NULL, 1},
{"nb", 0, NULL, 2},
{"pos", 1, NULL, 3},
{"ppos", 1, NULL, 4},
{"rb", 1, NULL, 5},
{"sb", 1, NULL, 6},
{"nd", 0, NULL, 7},
{"dump", 0, NULL, 8},
{"read", 0, NULL, 9},
{0,0,0,0}};
memset(host, 0, 256);
memset(interface, 0, 16);
memset(mcastAddr, 0, (size_t) mcastAddrMax*16);
memset(et_name, 0, ET_FILENAME_LENGTH);
memset(stationName, 0, ET_STATNAME_LENGTH);
while ((c = getopt_long_only(argc, argv, "vbmhrn:s:p:f:c:q:a:i:", long_options, 0)) != EOF) {
if (c == -1)
break;
switch (c) {
case 'c':
i_tmp = atoi(optarg);
if (i_tmp > 0 && i_tmp < 1001) {
chunk = i_tmp;
printf("Setting chunk to %d\n", chunk);
} else {
printf("Invalid argument to -c. Must < 1001 & > 0.\n");
exit(-1);
}
break;
case 'q':
i_tmp = atoi(optarg);
if (i_tmp > 0) {
qSize = i_tmp;
printf("Setting queue size to %d\n", qSize);
} else {
printf("Invalid argument to -q. Must > 0.\n");
exit(-1);
}
break;
case 's':
if (strlen(optarg) >= ET_STATNAME_LENGTH) {
fprintf(stderr, "Station name is too long\n");
exit(-1);
}
strcpy(stationName, optarg);
break;
case 'p':
i_tmp = atoi(optarg);
if (i_tmp > 1023 && i_tmp < 65535) {
port = (unsigned short)i_tmp;
} else {
printf("Invalid argument to -p. Must be < 65535 & > 1023.\n");
exit(-1);
}
break;
case 'f':
if (strlen(optarg) >= ET_FILENAME_LENGTH) {
fprintf(stderr, "ET file name is too long\n");
exit(-1);
}
strcpy(et_name, optarg);
break;
case 'i':
if (strlen(optarg) > 15 || strlen(optarg) < 7) {
fprintf(stderr, "interface address is bad\n");
exit(-1);
}
strcpy(interface, optarg);
break;
case 'a':
if (strlen(optarg) >= 16) {
fprintf(stderr, "Multicast address is too long\n");
exit(-1);
}
if (mcastAddrCount >= mcastAddrMax) break;
strcpy(mcastAddr[mcastAddrCount++], optarg);
multicast = 1;
break;
/* case host: */
case 1:
if (strlen(optarg) >= 255) {
fprintf(stderr, "host name is too long\n");
exit(-1);
}
strcpy(host, optarg);
break;
/* case nb: */
case 2:
blocking = 0;
break;
/* case pos: */
case 3:
i_tmp = atoi(optarg);
if (i_tmp > 0) {
position = i_tmp;
} else {
printf("Invalid argument to -pos. Must be > 0.\n");
exit(-1);
}
break;
/* case ppos: */
case 4:
i_tmp = atoi(optarg);
if (i_tmp > -3 && i_tmp != 0) {
pposition = i_tmp;
flowMode=ET_STATION_PARALLEL;
} else {
printf("Invalid argument to -ppos. Must be > -3 and != 0.\n");
exit(-1);
}
break;
/* case rb */
case 5:
i_tmp = atoi(optarg);
if (i_tmp < 1) {
printf("Invalid argument to -rb. Recv buffer size must be > 0.\n");
exit(-1);
}
recvBufSize = i_tmp;
break;
/* case sb */
case 6:
i_tmp = atoi(optarg);
if (i_tmp < 1) {
printf("Invalid argument to -sb. Send buffer size must be > 0.\n");
exit(-1);
}
sendBufSize = i_tmp;
break;
/* case nd */
case 7:
noDelay = 1;
break;
/* case dump */
case 8:
dump = 1;
break;
/* case read */
case 9:
readData = 1;
break;
case 'v':
verbose = 1;
debugLevel = ET_DEBUG_INFO;
break;
case 'r':
remote = 1;
break;
case 'm':
multicast = 1;
break;
case 'b':
broadcast = 1;
break;
case ':':
case 'h':
case '?':
default:
errflg++;
}
}
if (!multicast && !broadcast) {
/* Default to local host if direct connection */
if (strlen(host) < 1) {
strcpy(host, ET_HOST_LOCAL);
}
}
if (optind < argc || errflg || strlen(et_name) < 1) {
fprintf(stderr,
"\nusage: %s %s\n%s\n%s\n%s\n%s\n%s\n%s\n\n",
argv[0], "-f <ET name> -s <station name>",
" [-h] [-v] [-nb] [-r] [-m] [-b] [-nd] [-read] [-dump]",
" [-host <ET host>] [-p <ET port>]",
" [-c <chunk size>] [-q <Q size>]",
" [-pos <station pos>] [-ppos <parallel station pos>]",
" [-i <interface address>] [-a <mcast addr>]",
" [-rb <buf size>] [-sb <buf size>]");
fprintf(stderr, " -f ET system's (memory-mapped file) name\n");
fprintf(stderr, " -host ET system's host if direct connection (default to local)\n");
fprintf(stderr, " -s create station of this name\n");
fprintf(stderr, " -h help\n\n");
fprintf(stderr, " -v verbose output (also prints data if reading with -read)\n");
fprintf(stderr, " -read read data (1 int for each event)\n");
fprintf(stderr, " -dump dump events back into ET (go directly to GC) instead of put\n");
fprintf(stderr, " -c number of events in one get/put array\n");
fprintf(stderr, " -r act as remote (TCP) client even if ET system is local\n");
fprintf(stderr, " -p port, TCP if direct, else UDP\n\n");
fprintf(stderr, " -nb make station non-blocking\n");
fprintf(stderr, " -q queue size if creating non-blocking station\n");
fprintf(stderr, " -pos position of station (1,2,...)\n");
fprintf(stderr, " -ppos position of within a group of parallel stations (-1=end, -2=head)\n\n");
fprintf(stderr, " -i outgoing network interface address (dot-decimal)\n");
fprintf(stderr, " -a multicast address(es) (dot-decimal), may use multiple times\n");
fprintf(stderr, " -m multicast to find ET (use default address if -a unused)\n");
fprintf(stderr, " -b broadcast to find ET\n\n");
fprintf(stderr, " -rb TCP receive buffer size (bytes)\n");
fprintf(stderr, " -sb TCP send buffer size (bytes)\n");
fprintf(stderr, " -nd use TCP_NODELAY option\n\n");
fprintf(stderr, " This consumer works by making a direct connection to the\n");
fprintf(stderr, " ET system's server port and host unless at least one multicast address\n");
fprintf(stderr, " is specified with -a, the -m option is used, or the -b option is used\n");
fprintf(stderr, " in which case multi/broadcasting used to find the ET system.\n");
fprintf(stderr, " If multi/broadcasting fails, look locally to find the ET system.\n");
fprintf(stderr, " This program gets all events from the given station and puts them back.\n\n");
exit(2);
}
timeout.tv_sec = 2;
timeout.tv_nsec = 0;
/* allocate some memory */
pe = (et_event **) calloc((size_t)chunk, sizeof(et_event *));
if (pe == NULL) {
printf("%s: out of memory\n", argv[0]);
exit(1);
}
/*************************/
/* setup signal handling */
/*************************/
/* block all signals */
sigfillset(&sigblock);
status = pthread_sigmask(SIG_BLOCK, &sigblock, NULL);
if (status != 0) {
printf("%s: pthread_sigmask failure\n", argv[0]);
exit(1);
}
/* spawn signal handling thread */
pthread_create(&tid, NULL, signal_thread, (void *)NULL);
/******************/
/* open ET system */
/******************/
et_open_config_init(&openconfig);
if (broadcast && multicast) {
broadAndMulticast = 1;
}
/* if multicasting to find ET */
if (multicast) {
if (mcastAddrCount < 1) {
/* Use default mcast address if not given on command line */
status = et_open_config_addmulticast(openconfig, ET_MULTICAST_ADDR);
}
else {
/* add multicast addresses to use */
for (j = 0; j < mcastAddrCount; j++) {
if (strlen(mcastAddr[j]) > 7) {
status = et_open_config_addmulticast(openconfig, mcastAddr[j]);
if (status != ET_OK) {
printf("%s: bad multicast address argument\n", argv[0]);
exit(1);
}
printf("%s: adding multicast address %s\n", argv[0], mcastAddr[j]);
}
}
}
}
if (broadAndMulticast) {
printf("Broad and Multicasting\n");
if (port == 0) {
port = ET_UDP_PORT;
}
et_open_config_setport(openconfig, port);
et_open_config_setcast(openconfig, ET_BROADANDMULTICAST);
et_open_config_sethost(openconfig, ET_HOST_ANYWHERE);
}
else if (multicast) {
printf("Multicasting\n");
if (port == 0) {
port = ET_UDP_PORT;
}
et_open_config_setport(openconfig, port);
et_open_config_setcast(openconfig, ET_MULTICAST);
et_open_config_sethost(openconfig, ET_HOST_ANYWHERE);
}
else if (broadcast) {
printf("Broadcasting\n");
if (port == 0) {
port = ET_UDP_PORT;
}
et_open_config_setport(openconfig, port);
et_open_config_setcast(openconfig, ET_BROADCAST);
et_open_config_sethost(openconfig, ET_HOST_ANYWHERE);
}
else {
if (port == 0) {
port = ET_SERVER_PORT;
}
et_open_config_setserverport(openconfig, port);
et_open_config_setcast(openconfig, ET_DIRECT);
if (strlen(host) > 0) {
et_open_config_sethost(openconfig, host);
}
et_open_config_gethost(openconfig, host);
printf("Direct connection to %s\n", host);
}
/* Defaults are to use operating system default buffer sizes and turn off TCP_NODELAY */
et_open_config_settcp(openconfig, recvBufSize, sendBufSize, noDelay);
if (strlen(interface) > 6) {
et_open_config_setinterface(openconfig, interface);
}
if (remote) {
printf("Set as remote\n");
et_open_config_setmode(openconfig, ET_HOST_AS_REMOTE);
}
/* If responses from different ET systems, return error. */
et_open_config_setpolicy(openconfig, ET_POLICY_ERROR);
/* debug level */
et_open_config_setdebugdefault(openconfig, debugLevel);
et_open_config_setwait(openconfig, ET_OPEN_WAIT);
if (et_open(&id, et_name, openconfig) != ET_OK) {
printf("%s: et_open problems\n", argv[0]);
exit(1);
}
et_open_config_destroy(openconfig);
/*-------------------------------------------------------*/
/* Find out if we have a remote connection to the ET system */
et_system_getlocality(id, &locality);
if (locality == ET_REMOTE) {
printf("ET is remote\n\n");
et_system_gethost(id, host);
et_system_getlocaladdress(id, localAddr);
printf("Connect to ET, from ip = %s to %s\n", localAddr, host);
}
else {
printf("ET is local\n\n");
}
/* set level of debug output (everything) */
et_system_setdebug(id, debugLevel);
/* define station to create */
et_station_config_init(&sconfig);
et_station_config_setflow(sconfig, flowMode);
if (!blocking) {
et_station_config_setblock(sconfig, ET_STATION_NONBLOCKING);
if (qSize > 0) {
et_station_config_setcue(sconfig, qSize);
}
}
if ((status = et_station_create_at(id, &my_stat, stationName, sconfig, position, pposition)) != ET_OK) {
if (status == ET_ERROR_EXISTS) {
/* my_stat contains pointer to existing station */
printf("%s: station already exists\n", argv[0]);
}
else if (status == ET_ERROR_TOOMANY) {
printf("%s: too many stations created\n", argv[0]);
goto error;
}
else {
printf("%s: error in station creation\n", argv[0]);
goto error;
}
}
et_station_config_destroy(sconfig);
if (et_station_attach(id, my_stat, &attach1) != ET_OK) {
printf("%s: error in station attach\n", argv[0]);
goto error;
}
/* read time for future statistics calculations */
#if defined __APPLE__
gettimeofday(&t1, NULL);
time1 = 1000L*t1.tv_sec + t1.tv_usec/1000L; /* milliseconds */
#else
clock_gettime(CLOCK_REALTIME, &t1);
time1 = 1000L*t1.tv_sec + t1.tv_nsec/1000000L; /* milliseconds */
#endif
while (1) {
/**************/
/* get events */
/**************/
/* example of single, timeout read */
//status = et_events_get(id, attach1, pe, ET_TIMED | ET_MODIFY, &timeout, chunk, &numRead);
/* example of reading array of up to "chunk" events */
status = et_events_get(id, attach1, pe, ET_SLEEP, NULL, chunk, &numRead);
if (status == ET_OK) {
;
}
else if (status == ET_ERROR_DEAD) {
printf("%s: ET system is dead\n", argv[0]);
goto error;
}
else if (status == ET_ERROR_TIMEOUT) {
printf("%s: got timeout\n", argv[0]);
goto end;
}
else if (status == ET_ERROR_EMPTY) {
printf("%s: no events\n", argv[0]);
goto end;
}
else if (status == ET_ERROR_BUSY) {
printf("%s: station is busy\n", argv[0]);
goto end;
}
else if (status == ET_ERROR_WAKEUP) {
printf("%s: someone told me to wake up\n", argv[0]);
goto error;
}
else if ((status == ET_ERROR_WRITE) || (status == ET_ERROR_READ)) {
printf("%s: socket communication error\n", argv[0]);
goto error;
}
else {
printf("%s: get error, status = %d\n", argv[0], status);
goto error;
}
/*******************/
/* read/print data */
/*******************/
if (readData) {
size_t len;
int *data, endian, swap;
for (j = 0; j < numRead; j++) {
et_event_getdata(pe[j], (void **) &data);
et_event_getlength(pe[j], &len);
et_event_getendian(pe[j], &endian);
et_event_needtoswap(pe[j], &swap);
bytes += len;
totalBytes += len;
if (verbose) {
printf("data byte order = %s\n", (endian == ET_ENDIAN_BIG ? "BIG" : "LITTLE"));
if (swap) {
printf(" data (len = %d) needs swapping, swapped int = %d\n", (int) len, ET_SWAP32(data[0]));
}
else {
printf(" data (len = %d) does NOT need swapping, int = %d\n", (int) len, data[0]);
}
et_event_getcontrol(pe[j], con);
printf("control array = {");
for (i = 0; i < ET_STATION_SELECT_INTS; i++) {
printf("%d ", con[i]);
}
printf("}\n");
}
}
}
else {
size_t len;
for (j = 0; j < numRead; j++) {
et_event_getlength(pe[j], &len);
/* include ET overhead by adding commented out portions */
bytes += len /* +52 */;
totalBytes += len /* +52 */;
}
/*
bytes += 20;
totalBytes += 20;
*/
}
/*******************/
/* put/dump events */
/*******************/
if (!dump) {
/* putting array of events */
status = et_events_put(id, attach1, pe, numRead);
}
else {
/* dumping array of events */
status = et_events_dump(id, attach1, pe, numRead);
}
if (status == ET_ERROR_DEAD) {
printf("%s: ET is dead\n", argv[0]);
goto error;
}
else if ((status == ET_ERROR_WRITE) || (status == ET_ERROR_READ)) {
printf("%s: socket communication error\n", argv[0]);
goto error;
}
else if (status != ET_OK) {
printf("%s: put error\n", argv[0]);
goto error;
}
count += numRead;
end:
/* statistics */
#if defined __APPLE__
gettimeofday(&t2, NULL);
time2 = 1000L*t2.tv_sec + t2.tv_usec/1000L; /* milliseconds */
#else
clock_gettime(CLOCK_REALTIME, &t2);
time2 = 1000L*t2.tv_sec + t2.tv_nsec/1000000L; /* milliseconds */
#endif
time = time2 - time1;
if (time > 5000) {
/* reset things if necessary */
if ( (totalCount >= (LONG_MAX - count)) ||
(totalT >= (LONG_MAX - time)) ) {
bytes = totalBytes = totalT = totalCount = count = 0;
time1 = time2;
continue;
}
rate = 1000.0 * ((double) count) / time;
totalCount += count;
totalT += time;
avgRate = 1000.0 * ((double) totalCount) / totalT;
/* Event rates */
printf("%s Events: %3.4g Hz, %3.4g Avg.\n", argv[0], rate, avgRate);
/* Data rates */
rate = ((double) bytes) / time;
avgRate = ((double) totalBytes) / totalT;
printf("%s Data: %3.4g kB/s, %3.4g Avg.\n\n", argv[0], rate, avgRate);
bytes = count = 0;
#if defined __APPLE__
gettimeofday(&t1, NULL);
time1 = 1000L*t1.tv_sec + t1.tv_usec/1000L;
#else
clock_gettime(CLOCK_REALTIME, &t1);
time1 = 1000L*t1.tv_sec + t1.tv_nsec/1000000L;
#endif
}
} /* while(1) */
error:
printf("%s: ERROR\n", argv[0]);
return 0;
}
/************************************************************/
/* separate thread to handle signals */
static void *signal_thread (void *arg) {
sigset_t signal_set;
int sig_number;
sigemptyset(&signal_set);
sigaddset(&signal_set, SIGINT);
/* Wait for Control-C */
sigwait(&signal_set, &sig_number);
printf("Got control-C, exiting\n");
exit(1);
}
#include "ConfigOption.h"
#include "PRadETChannel.h"
#include "et.h"
#include "evio/evioUtil.hxx"
#include "evio/evioFileChannel.hxx"
#include <csignal>
#include <thread>
#include <chrono>
#include <iostream>
#define PROGRESS_COUNT 100
using namespace std::chrono;
volatile std::sig_atomic_t gSignalStatus;
void signal_handler(int signal) {
gSignalStatus = signal;
}
int main(int argc, char* argv[]) try
{
// setup input arguments
ConfigOption conf_opt;
conf_opt.AddLongOpt(ConfigOption::help_message, "help");
conf_opt.AddOpt(ConfigOption::arg_require, 'h');
conf_opt.AddOpt(ConfigOption::arg_require, 'p');
conf_opt.AddOpt(ConfigOption::arg_require, 'f');
conf_opt.AddOpt(ConfigOption::arg_require, 'i');
conf_opt.SetDesc("usage: %0 <data_file>");
conf_opt.SetDesc('h', "host address of the ET system, default \"localhost\".");
conf_opt.SetDesc('p', "port to connect, default 11111.");
conf_opt.SetDesc('f', "memory mapped et file, default \"/tmp/et_feeder\".");
conf_opt.SetDesc('i', "interval in milliseconds to write data, default \"10\"");
if (!conf_opt.ParseArgs(argc, argv) || conf_opt.NbofArgs() != 1) {
std::cout << conf_opt.GetInstruction() << std::endl;
return -1;
}
std::string host = "localhost";
int port = 11111;
std::string etf = "/tmp/et_feeder";
int interval = 10;
for (auto &opt : conf_opt.GetOptions()) {
switch (opt.mark) {
case 'h':
host = opt.var.String();
break;
case 'c':
port = opt.var.Int();
break;
case 'f':
etf = opt.var.String();
break;
case 'i':
interval = opt.var.Int();
break;
default :
std::cout << conf_opt.GetInstruction() << std::endl;
return -1;
}
}
// attach to ET system
auto ch = new PRadETChannel();
ch->Open(host.c_str(), port, etf.c_str());
ch->NewStation("Data Feeder");
ch->AttachStation();
// evio file reader
evio::evioFileChannel *chan = new evio::evioFileChannel(conf_opt.GetArgument(0).c_str(), "r");
chan->open();
// install signal handler
std::signal(SIGINT, signal_handler);
int count = 0;
while (chan->read()) {
if (gSignalStatus == SIGINT) {
std::cout << "Received control-C, exiting..." << std::endl;
ch->ForceClose();
break;
}
system_clock::time_point start(system_clock::now());
system_clock::time_point next(start + std::chrono::milliseconds(interval));
if (++count % PROGRESS_COUNT == 0) {
std::cout << "Read and feed " << count << " events to ET, rate is 1 event per "
<< interval << " ms.\r" << std::flush;
}
ch->Write((void*) chan->getBuffer(), chan->getBufSize() * sizeof(uint32_t));
std::this_thread::sleep_until(next);
}
std::cout << "Read and feed " << count << " events to ET, rate is 1 event per "
<< interval << " ms." << std::endl;
chan->close();
return 0;
} catch (PRadException e) {
std::cerr << e.FailureType() << ": " << e.FailureDesc() << std::endl;
return -1;
} catch (evio::evioException e) {
std::cerr << e.toString() << endl;
} catch (...) {
std::cerr << "?unknown exception" << endl;
}
/*----------------------------------------------------------------------------*
* Copyright (c) 1998 Southeastern Universities Research Association, *
* Thomas Jefferson National Accelerator Facility *
* *
* This software was developed under a United States Government license *
* described in the NOTICE file included as part of this distribution. *
* *
* Author: Carl Timmer *
* timmer@jlab.org Jefferson Lab, MS-12H *
* Phone: (757) 269-5130 12000 Jefferson Ave. *
* Fax: (757) 269-5800 Newport News, VA 23606 *
* *
*----------------------------------------------------------------------------*
*
* Description:
* ET system sample event producer
*
*----------------------------------------------------------------------------*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <signal.h>
#include <sys/types.h>
#include <unistd.h>
#include <limits.h>
#include <getopt.h>
#include <time.h>
#include <sys/time.h>
#include <math.h>
#include "et.h"
/* prototype */
static void * signal_thread (void *arg);
int main(int argc,char **argv) {
int i, j, c, i_tmp, status, junkData, numRead, locality;
int startingVal=0, errflg=0, group=1, chunk=1, size=32, writeData=0, localEndian=1;
int delay=0, remote=0, multicast=0, broadcast=0, broadAndMulticast=0;
int sendBufSize=0, recvBufSize=0, noDelay=0, blast=0, noAllocFlag=0;
int debugLevel = ET_DEBUG_ERROR;
unsigned short port=0;
char et_name[ET_FILENAME_LENGTH], host[256], interface[16], localAddr[16];
void *fakeData;
int mcastAddrCount = 0, mcastAddrMax = 10;
char mcastAddr[mcastAddrMax][16];
et_att_id attach1;
et_sys_id id;
et_openconfig openconfig;
et_event **pe;
struct timespec timeout;
#if defined __APPLE__
struct timeval t1, t2;
#else
struct timespec t1, t2;
#endif
sigset_t sigblock;
pthread_t tid;
/* statistics variables */
double rate = 0.0, avgRate = 0.0;
int64_t count = 0, totalCount = 0, totalT = 0, time, time1, time2;
/* control int array for event header if writing junk data */
int control[] = {17, 8, -1, -1, 0, 0};
/* 4 multiple character command-line options */
static struct option long_options[] =
{{"host", 1, NULL, 1},
{"rb", 1, NULL, 2},
{"sb", 1, NULL, 3},
{"nd", 0, NULL, 4},
{"blast", 0, NULL, 5},
{0, 0, 0, 0}
};
memset(host, 0, 256);
memset(interface, 0, 16);
memset(mcastAddr, 0, mcastAddrMax*16);
memset(et_name, 0, ET_FILENAME_LENGTH);
while ((c = getopt_long_only(argc, argv, "vbmhrn:a:s:p:d:f:c:g:i:w:", long_options, 0)) != EOF) {
if (c == -1)
break;
switch (c) {
case 'c':
i_tmp = atoi(optarg);
if (i_tmp > 0 && i_tmp < 1001) {
chunk = i_tmp;
printf("Setting chunk to %d\n", chunk);
} else {
printf("Invalid argument to -c. Must < 1001 & > 0.\n");
exit(-1);
}
break;
case 's':
i_tmp = atoi(optarg);
if (i_tmp > -1) {
size = i_tmp;
} else {
printf("Invalid argument to -s. Must be a positive integer.\n");
exit(-1);
}
break;
case 'w':
writeData = 1;
i_tmp = atoi(optarg);
if (i_tmp != 1) {
localEndian = 0;
}
break;
case 'd':
delay = atoi(optarg);
if (delay < 0) {
printf("Delay must be >= 0\n");
exit(-1);
}
break;
case 'p':
i_tmp = atoi(optarg);
if (i_tmp > 1023 && i_tmp < 65535) {
port = i_tmp;
} else {
printf("Invalid argument to -p. Must be < 65535 & > 1023.\n");
exit(-1);
}
break;
case 'g':
i_tmp = atoi(optarg);
if (i_tmp > 0 && i_tmp < 501) {
group = i_tmp;
} else {
printf("Invalid argument to -g. Must be 501 > g > 0.\n");
exit(-1);
}
break;
case 'f':
if (strlen(optarg) >= ET_FILENAME_LENGTH) {
fprintf(stderr, "ET file name is too long\n");
exit(-1);
}
strcpy(et_name, optarg);
break;
case 'i':
if (strlen(optarg) > 15 || strlen(optarg) < 7) {
fprintf(stderr, "interface address is bad\n");
exit(-1);
}
strcpy(interface, optarg);
break;
case 'a':
if (strlen(optarg) >= 16) {
fprintf(stderr, "Multicast address is too long\n");
exit(-1);
}
if (mcastAddrCount >= mcastAddrMax) break;
strcpy(mcastAddr[mcastAddrCount++], optarg);
multicast = 1;
break;
case 1:
if (strlen(optarg) >= 255) {
fprintf(stderr, "host name is too long\n");
exit(-1);
}
strcpy(host, optarg);
break;
/* case rb */
case 2:
i_tmp = atoi(optarg);
if (i_tmp < 1) {
printf("Invalid argument to -rb. Recv buffer size must be > 0.\n");
exit(-1);
}
recvBufSize = i_tmp;
break;
/* case sb */
case 3:
i_tmp = atoi(optarg);
if (i_tmp < 1) {
printf("Invalid argument to -sb. Send buffer size must be > 0.\n");
exit(-1);
}
sendBufSize = i_tmp;
break;
/* case nd */
case 4:
noDelay = 1;
break;
/* case blast */
case 5:
blast = 1;
break;
case 'v':
debugLevel = ET_DEBUG_INFO;
break;
case 'r':
remote = 1;
break;
case 'm':
multicast = 1;
break;
case 'b':
broadcast = 1;
break;
case ':':
case 'h':
case '?':
default:
errflg++;
}
}
if (!multicast && !broadcast) {
/* Default to local host if direct connection */
if (strlen(host) < 1) {
strcpy(host, ET_HOST_LOCAL);
}
}
if (optind < argc || errflg || strlen(et_name) < 1) {
fprintf(stderr,
"\nusage: %s %s\n%s\n%s\n%s\n%s\n%s\n%s\n\n",
argv[0], "-f <ET name>",
" [-h] [-v] [-r] [-m] [-b] [-nd] [-blast]",
" [-host <ET host>] [-w <local endian? 0/1>]",
" [-s <event size>] [-c <chunk size>] [-g <group>]",
" [-d <delay>] [-p <ET port>]",
" [-i <interface address>] [-a <mcast addr>]",
" [-rb <buf size>] [-sb <buf size>]");
fprintf(stderr, " -f ET system's (memory-mapped file) name\n");
fprintf(stderr, " -host ET system's host if direct connection (default to local)\n");
fprintf(stderr, " -h help\n");
fprintf(stderr, " -v verbose output\n\n");
fprintf(stderr, " -s event size in bytes\n");
fprintf(stderr, " -c number of events in one get/put array\n");
fprintf(stderr, " -g group from which to get new events (1,2,...)\n");
fprintf(stderr, " -d delay in millisec between each round of getting and putting events\n\n");
fprintf(stderr, " -p ET port (TCP for direct, UDP for broad/multicast)\n");
fprintf(stderr, " -r act as remote (TCP) client even if ET system is local\n");
fprintf(stderr, " -w write data (1 sequential int per event), 1 local endian, 0 else\n");
fprintf(stderr, " -blast if remote, use external data buf (no mem allocation),\n");
fprintf(stderr, " do not write data (overrides -w)\n\n");
fprintf(stderr, " -i outgoing network interface address (dot-decimal)\n");
fprintf(stderr, " -a multicast address(es) (dot-decimal), may use multiple times\n");
fprintf(stderr, " -m multicast to find ET (use default address if -a unused)\n");
fprintf(stderr, " -b broadcast to find ET\n\n");
fprintf(stderr, " -rb TCP receive buffer size (bytes)\n");
fprintf(stderr, " -sb TCP send buffer size (bytes)\n");
fprintf(stderr, " -nd use TCP_NODELAY option\n\n");
fprintf(stderr, " This producer works by making a direct connection to the\n");
fprintf(stderr, " ET system's server port and host unless at least one multicast address\n");
fprintf(stderr, " is specified with -a, the -m option is used, or the -b option is used\n");
fprintf(stderr, " in which case multi/broadcasting used to find the ET system.\n");
fprintf(stderr, " If multi/broadcasting fails, look locally to find the ET system.\n");
fprintf(stderr, " This program gets new events from the system and puts them back.\n\n");
exit(2);
}
/* fake data for blasting */
fakeData = (void *) malloc(size);
if (fakeData == NULL) {
printf("%s: out of memory\n", argv[0]);
exit(1);
}
/* delay is in milliseconds */
if (delay > 0) {
timeout.tv_sec = delay / 1000;
timeout.tv_nsec = (delay - (delay / 1000) * 1000) * 1000000;
}
/* allocate some memory */
pe = (et_event **) calloc(chunk, sizeof(et_event *));
if (pe == NULL) {
printf("%s: out of memory\n", argv[0]);
exit(1);
}
/*************************/
/* setup signal handling */
/*************************/
/* block all signals */
sigfillset(&sigblock);
status = pthread_sigmask(SIG_BLOCK, &sigblock, NULL);
if (status != 0) {
printf("%s: pthread_sigmask failure\n", argv[0]);
exit(1);
}
/* spawn signal handling thread */
pthread_create(&tid, NULL, signal_thread, (void *) NULL);
/******************/
/* open ET system */
/******************/
et_open_config_init(&openconfig);
if (broadcast && multicast) {
broadAndMulticast = 1;
}
/* if multicasting to find ET */
if (multicast) {
if (mcastAddrCount < 1) {
/* Use default mcast address if not given on command line */
status = et_open_config_addmulticast(openconfig, ET_MULTICAST_ADDR);
}
else {
/* add multicast addresses to use */
for (j = 0; j < mcastAddrCount; j++) {
if (strlen(mcastAddr[j]) > 7) {
status = et_open_config_addmulticast(openconfig, mcastAddr[j]);
if (status != ET_OK) {
printf("%s: bad multicast address argument\n", argv[0]);
exit(1);
}
printf("%s: adding multicast address %s\n", argv[0], mcastAddr[j]);
}
}
}
}
if (broadAndMulticast) {
printf("Broad and Multicasting\n");
if (port == 0) {
port = ET_UDP_PORT;
}
et_open_config_setport(openconfig, port);
et_open_config_setcast(openconfig, ET_BROADANDMULTICAST);
et_open_config_sethost(openconfig, ET_HOST_ANYWHERE);
}
else if (multicast) {
printf("Multicasting\n");
if (port == 0) {
port = ET_UDP_PORT;
}
et_open_config_setport(openconfig, port);
et_open_config_setcast(openconfig, ET_MULTICAST);
et_open_config_sethost(openconfig, ET_HOST_ANYWHERE);
}
else if (broadcast) {
printf("Broadcasting\n");
if (port == 0) {
port = ET_UDP_PORT;
}
et_open_config_setport(openconfig, port);
et_open_config_setcast(openconfig, ET_BROADCAST);
et_open_config_sethost(openconfig, ET_HOST_ANYWHERE);
}
else {
if (port == 0) {
port = ET_SERVER_PORT;
}
et_open_config_setserverport(openconfig, port);
et_open_config_setcast(openconfig, ET_DIRECT);
if (strlen(host) > 0) {
et_open_config_sethost(openconfig, host);
}
et_open_config_gethost(openconfig, host);
printf("Direct connection to %s\n", host);
}
/*et_open_config_sethost(openconfig, ET_HOST_REMOTE);*/
/* Defaults are to use operating system default buffer sizes and turn off TCP_NODELAY */
et_open_config_settcp(openconfig, recvBufSize, sendBufSize, noDelay);
if (strlen(interface) > 6) {
et_open_config_setinterface(openconfig, interface);
}
if (remote) {
printf("Set as remote\n");
et_open_config_setmode(openconfig, ET_HOST_AS_REMOTE);
}
/* If responses from different ET systems, return error. */
et_open_config_setpolicy(openconfig, ET_POLICY_ERROR);
/* debug level */
et_open_config_setdebugdefault(openconfig, debugLevel);
et_open_config_setwait(openconfig, ET_OPEN_WAIT);
if (et_open(&id, et_name, openconfig) != ET_OK) {
printf("%s: et_open problems\n", argv[0]);
exit(1);
}
et_open_config_destroy(openconfig);
/*-------------------------------------------------------*/
/* Make things self-consistent by not taking time to write data if blasting.
* Blasting flag takes precedence. */
if (blast) {
writeData = 0;
}
/* Find out if we have a remote connection to the ET system
* so we know if we can use external data buffer for events
* for blasting - which is quite a bit faster. */
et_system_getlocality(id, &locality);
if (locality == ET_REMOTE) {
if (blast) {
noAllocFlag = ET_NOALLOC;
}
printf("ET is remote\n\n");
et_system_gethost(id, host);
et_system_getlocaladdress(id, localAddr);
printf("Connect to ET, from ip = %s to %s\n", localAddr, host);
}
else {
/* local blasting is just the same as local producing */
blast = 0;
printf("ET is local\n\n");
}
/* set level of debug output (everything) */
et_system_setdebug(id, debugLevel);
/* attach to grandcentral station */
if (et_station_attach(id, ET_GRANDCENTRAL, &attach1) < 0) {
printf("%s: error in et_station_attach\n", argv[0]);
exit(1);
}
/* read time for future statistics calculations */
#if defined __APPLE__
gettimeofday(&t1, NULL);
time1 = 1000L*t1.tv_sec + t1.tv_usec/1000L; /* milliseconds */
#else
clock_gettime(CLOCK_REALTIME, &t1);
time1 = 1000L*t1.tv_sec + t1.tv_nsec/1000000L; /* milliseconds */
#endif
while (1) {
status = et_events_new_group(id, attach1, pe, ET_SLEEP | noAllocFlag,
NULL, size, chunk, group, &numRead);
if (status == ET_OK) {
;
}
else if (status == ET_ERROR_DEAD) {
printf("%s: ET system is dead\n", argv[0]);
break;
}
else if (status == ET_ERROR_TIMEOUT) {
printf("%s: got timeout\n", argv[0]);
break;
}
else if (status == ET_ERROR_EMPTY) {
printf("%s: no events\n", argv[0]);
break;
}
else if (status == ET_ERROR_BUSY) {
printf("%s: grandcentral is busy\n", argv[0]);
break;
}
else if (status == ET_ERROR_WAKEUP) {
printf("%s: someone told me to wake up\n", argv[0]);
break;
}
else if ((status == ET_ERROR_WRITE) || (status == ET_ERROR_READ)) {
printf("%s: socket communication error\n", argv[0]);
goto error;
}
else {
printf("%s: request error\n", argv[0]);
goto error;
}
/* if blasting data (and remote), don't write anything, just use what's in buffer when allocated */
if (blast) {
for (i=0; i < numRead; i++) {
et_event_setlength(pe[i], size);
et_event_setdatabuffer(id, pe[i], fakeData);
}
}
/* write data, set control values here */
else if (writeData) {
char *pdata;
for (i=0; i < numRead; i++) {
junkData = i + startingVal;
if (!localEndian) {
junkData = ET_SWAP32(junkData);
et_event_setendian(pe[i], ET_ENDIAN_NOTLOCAL);
}
et_event_getdata(pe[i], (void **) &pdata);
memcpy((void *)pdata, (const void *) &junkData, sizeof(int));
/* Send all data even though we only wrote one int. */
et_event_setlength(pe[i], size);
et_event_setcontrol(pe[i], control, sizeof(control)/sizeof(int));
}
startingVal += numRead;
}
else {
for (i = 0; i < numRead; i++) {
et_event_setlength(pe[i], size);
}
}
/* put events back into the ET system */
status = et_events_put(id, attach1, pe, numRead);
if (status == ET_OK) {
;
}
else if (status == ET_ERROR_DEAD) {
printf("%s: ET is dead\n", argv[0]);
break;
}
else if ((status == ET_ERROR_WRITE) || (status == ET_ERROR_READ)) {
printf("%s: socket communication error\n", argv[0]);
goto error;
}
else {
printf("%s: put error, status = %d\n", argv[0], status);
goto error;
}
count += numRead;
if (delay > 0) {
nanosleep(&timeout, NULL);
}
/* statistics */
#if defined __APPLE__
gettimeofday(&t2, NULL);
time2 = 1000L*t2.tv_sec + t2.tv_usec/1000L; /* milliseconds */
#else
clock_gettime(CLOCK_REALTIME, &t2);
time2 = 1000L*t2.tv_sec + t2.tv_nsec/1000000L; /* milliseconds */
#endif
time = time2 - time1;
if (time > 5000) {
/* reset things if necessary */
if ( (totalCount >= (LONG_MAX - count)) ||
(totalT >= (LONG_MAX - time)) ) {
totalT = totalCount = count = 0;
time1 = time2;
continue;
}
rate = 1000.0 * ((double) count) / time;
totalCount += count;
totalT += time;
avgRate = 1000.0 * ((double) totalCount) / totalT;
/* Event rates */
printf("%s Events: %3.4g Hz, %3.4g Avg.\n", argv[0], rate, avgRate);
/* Data rates */
rate = ((double) count) * size / time;
avgRate = ((double) totalCount) * size / totalT;
printf("%s Data: %3.4g kB/s, %3.4g Avg.\n\n", argv[0], rate, avgRate);
/* If including msg overhead in data rates, need to do the following
avgRate = 1000.0 * (((double) totalCount) * (size + 52) + 20)/ totalT; */
count = 0;
#if defined __APPLE__
gettimeofday(&t1, NULL);
time1 = 1000L*t1.tv_sec + t1.tv_usec/1000L;
#else
clock_gettime(CLOCK_REALTIME, &t1);
time1 = 1000L*t1.tv_sec + t1.tv_nsec/1000000L;
#endif
}
} /* while(1) */
error:
printf("%s: ERROR\n", argv[0]);
exit(0);
}
/************************************************************/
/* separate thread to handle signals */
static void * signal_thread (void *arg) {
sigset_t signal_set;
int sig_number;
sigemptyset(&signal_set);
sigaddset(&signal_set, SIGINT);
/* Not necessary to clean up as ET system will do it */
sigwait(&signal_set, &sig_number);
printf("Got control-C, exiting\n");
exit(1);
}