Example 2 of Using Scaler Classes to Analyze ROC10/11.
Robert Michaels, rom@jlab.org, Jefferson
Lab Hall A, May 2003
This file:
hallaweb.jlab.org/equipment/daq/tscalroc11_main.html
The following example shows how to use the
hana_scaler package to analyze a CODA file and analyze the
scaler data that comes with ROC10 and ROC11. These data are
not any better or worse than the event type 140 scalers, except
for two issues: 1) The ROC10/11 data may come more frequently,
and 2) The ROC10/11 data are well synchronized to the event stream.
In contrast, the event type 140 scalers are inserted by the ET
system and are synchronized to approx. only 1 second. If synchronizing
is an issue, one would need to use clocks in the events and in
the scalers. However, for ROC10/11 one is guaranteed that the
physics trigger before/after really came before/after.
//--------------------------------------------------------
// tscalroc11_main.C
//
// Test of the ROC10/11 scalers which are read in the datastream
// typically every 100 events. As of Aug 2002 this is a new
// readout. (Actually there was also a similar ROC11 since
// Aug 2001, but the data format has changed due to G0 mode.)
// The focus is on the ring buffer data.
// R. Michaels, Oct 2002
//--------------------------------------------------------
// The ROC to find scalers (10/11 are R/L).
#define MYROC 11
// To printout (1) or not (0)
#define PRINTOUT 0
#define NBCM 6
#define MAXRING 20
#include
#include "THaCodaFile.h"
#include "THaEvData.h"
#ifndef __CINT__
#include "TROOT.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TNtuple.h"
#include "TRandom.h"
#endif
using namespace std;
int loadHelicity();
int ranBit(unsigned int& seed);
unsigned int getSeed();
// Global variables --------------------------------------
#define NBIT 24
int hbits[NBIT]; // The NBIT shift register
int present_reading; // present reading of helicity
int predicted_reading; // prediction of present reading (should = present_reading)
int present_helicity; // present helicity (using prediction)
#define NDELAY 2 // number of quartets between
// present reading and present helicity
// (i.e. the delay in reporting the helicity)
int recovery_flag = 0; // flag to determine if we need to recover
// from an error (1) or not (0)
unsigned int iseed; // value of iseed for present_helicity
unsigned int iseed_earlier; // iseed for predicted_reading
int main(int argc, char* argv[]) {
int helicity, qrt, gate, timestamp;
int len, data, status, nscaler, header;
int numread, badread, i, found, index;
int ring_clock, ring_qrt, ring_helicity;
int ring_trig, ring_bcm, ring_l1a, ring_v2fh;
int sum_clock, sum_trig, sum_bcm, sum_l1a;
int inquad, nrread, q1_helicity;
int ring_data[MAXRING], rloc;
if (argc < 2) {
cout << "You made a mistake... bye bye !\n" << endl;
cout << "Usage: 'tscalroc" << MYROC << " filename'" << endl;
cout << " where 'filename' is the CODA file." << endl;
exit(0);
}
// Setup
// Pedestals. Left, Right Arms. u1,u3,u10,d1,d3,d10
Float_t bcmpedL[NBCM] = { 188.2, 146.2, 271.6, 37.8, 94.2, 260.2 };
Float_t bcmpedR[NBCM] = { 53.0, 41.8, 104.1, 0., 101.6, 254.6 };
Float_t bcmped[NBCM];
if (MYROC == 11) {
cout << "Using Left Arm BCM pedestals" << endl;
for (i = 0; i < NBCM; i++) {
bcmped[i] = bcmpedL[i];
}
} else {
cout << "Using Right Arm BCM pedestals" << endl;
for (i = 0; i < NBCM; i++) {
bcmped[i] = bcmpedR[i];
}
}
// Initialize root and output.
TROOT scalana("scalroot","Hall A scaler analysis");
TFile hfile("scaler.root","RECREATE","Scaler data in Hall A");
// Define the ntuple here
// 0 1 2 3 4 5 6
char rawrates[]="evt:clk:trig:bcm:l1a:mhel:phel:";
// Asymmetries delayed various amounts, using RING BUFFER
// 7 8 9 10 11 12 13 14 15 16
char asydelay[]="a3d1:a3d2:a3d3:a3d4:a3d5:a3d6:a3d7:a3d8:a3d9:a3d10";
int nlen = strlen(rawrates) + strlen(asydelay);
char *string_ntup = new char[nlen+1];
strcpy(string_ntup,rawrates);
strcat(string_ntup,asydelay);
TNtuple *ntup = new TNtuple("ascal","Scaler Rates",string_ntup);
Float_t farray_ntup[17]; // Dimension = size of string_ntup (i.e. nlen+1)
THaCodaFile *coda = new THaCodaFile(TString(argv[1]));
THaEvData evdata;
inquad = 0;
q1_helicity = 0;
rloc = 0;
status = 0;
sum_clock = 0; sum_trig = 0; sum_bcm = 0; sum_l1a = 0;
nrread = 0;
memset(ring_data, 0, MAXRING*sizeof(int));
while (status == 0) {
status = coda->codaRead();
if (status != 0) break;
evdata.LoadEvent(coda->getEvBuffer());
len = evdata.GetRocLength(MYROC);
if (len <= 4) continue;
data = evdata.GetRawData(MYROC,3);
helicity = (data & 0x10) >> 4;
qrt = (data & 0x20) >> 5;
gate = (data & 0x40) >> 6;
timestamp = evdata.GetRawData(MYROC,4);
found = 0;
index = 5;
while ( !found ) {
data = evdata.GetRawData(MYROC,index++);
if ( (data & 0xffff0000) == 0xfb0b0000) found = 1;
if (index >= len) break;
}
if (!found) break;
nscaler = data & 0x7;
if (PRINTOUT) {
cout << hex << "helicity " << helicity << " qrt " << qrt;
cout << " gate " << gate << " time stamp " << timestamp << endl;
cout << "nscaler in this event " << nscaler << endl;
}
if (nscaler <= 0) continue;
if (nscaler > 2) nscaler = 2; // shouldn't be necessary
// 32 channels of scaler data for two helicities.
if (PRINTOUT) cout << "Synch event ----> " << endl;
for (int ihel = 0; ihel < nscaler; ihel++) {
header = evdata.GetRawData(MYROC,index++);
if (PRINTOUT) {
cout << "Scaler for helicity = " << dec << ihel;
cout << " unique header = " << hex << header << endl;
}
for (int ichan = 0; ichan < 32; ichan++) {
data = evdata.GetRawData(MYROC,index++);
if (PRINTOUT) {
cout << "channel # " << dec << ichan+1;
cout << " (hex) data = " << hex << data << endl;
}
}
}
numread = evdata.GetRawData(MYROC,index++);
badread = evdata.GetRawData(MYROC,index++);
if (PRINTOUT) cout << "FIFO num of last good read " << dec << numread << endl;
if (badread != 0) {
cout << "DISASTER: There are bad readings " << endl;
cout << "FIFO num of last bad read " << endl;
}
// Subset of scaler channels stored in a 30 Hz ring buffer.
int nring = 0;
while (index < len && nring == 0) {
header = evdata.GetRawData(MYROC,index++);
if ((header & 0xffff0000) == 0xfb1b0000) {
nring = header & 0x3ff;
}
}
if (PRINTOUT) cout << "Num in ring buffer = " << dec << nring << endl;
// The following assumes three are 6 pieces of data per 'iring'
// This was true after Jan 10, 2003.
for (int iring = 0; iring < nring; iring++) {
ring_clock = evdata.GetRawData(MYROC,index++);
data = evdata.GetRawData(MYROC,index++);
ring_qrt = (data & 0x10) >> 4;
ring_helicity = (data & 0x1);
present_reading = ring_helicity;
if (ring_qrt) {
inquad = 1;
if (loadHelicity()) {
if (present_reading != predicted_reading) {
cout << "DISASTER: The helicity is wrong !!"< 12) {
farray_ntup[0] = (float)nrread;
farray_ntup[1] = ring_clock;
farray_ntup[2] = ring_trig;
farray_ntup[3] = ring_bcm;
farray_ntup[4] = ring_l1a;
farray_ntup[5] = (float)present_reading;
farray_ntup[6] = (float)present_helicity;
for (int ishift = 3; ishift <= 12; ishift++) {
Float_t pdat = ring_data[(rloc-ishift)%MAXRING];
Float_t mdat = ring_data[(rloc-ishift+1)%MAXRING];
Float_t sum = pdat + mdat;
Float_t asy = 99999;
if (sum != 0) {
if (present_reading == 1) {
asy = (pdat - mdat)/sum;
} else {
asy = (mdat - pdat)/sum;
}
}
farray_ntup[4+ishift] = asy;
}
ntup->Fill(farray_ntup);
}
rloc++;
if (PRINTOUT) {
cout << "buff [" << dec << iring << "] ";
cout << " clock " << ring_clock << " qrt " << ring_qrt;
cout << " helicity " << ring_helicity;
cout << " trigger " << ring_trig << " bcm " << ring_bcm;
cout << " L1a " << ring_l1a << endl;
cout << " ring_v2fh (helicity in v2f) " << ring_v2fh << endl;
cout << " inquad " << inquad << endl;
cout << " sums: " << endl;
cout << " clock " << sum_clock << " trig " << sum_trig << endl;
cout << " bcm " << sum_bcm << " l1a " << sum_l1a << endl;
}
}
}
cout << "All done" << endl;
hfile.Write();
hfile.Close();
return 0;
}
// *************************************************************
// loadHelicity()
// Loads the helicity to determine the seed.
// After loading (nb==NBIT), predicted helicities are available.
// *************************************************************
int loadHelicity() {
int i;
static int nb;
if (recovery_flag) nb = 0;
recovery_flag = 0;
if (nb < NBIT) {
hbits[nb] = present_reading;
nb++;
return 0;
} else if (nb == NBIT) { // Have finished loading
iseed_earlier = getSeed();
for (i = 0; i < NBIT+1; i++)
predicted_reading = ranBit(iseed_earlier);
iseed = iseed_earlier;
for (i = 0; i < NDELAY; i++)
present_helicity = ranBit(iseed);
nb++;
return 1;
} else {
predicted_reading = ranBit(iseed_earlier);
present_helicity = ranBit(iseed);
if (PRINTOUT) {
cout << "helicities " << predicted_reading << " " <<
present_reading << " " << present_helicity << endl;
}
return 1;
}
}
// *************************************************************
// This is the random bit generator according to the G0
// algorithm described in "G0 Helicity Digital Controls" by
// E. Stangland, R. Flood, H. Dong, July 2002.
// Argument:
// ranseed = seed value for random number.
// This value gets modified.
// Return value:
// helicity (0 or 1)
// *************************************************************
int ranBit(unsigned int& ranseed) {
static int IB1 = 1; // Bit 1
static int IB3 = 4; // Bit 3
static int IB4 = 8; // Bit 4
static int IB24 = 8388608; // Bit 24
static int MASK = IB1+IB3+IB4+IB24;
if(ranseed & IB24) {
ranseed = ((ranseed^MASK) << 1) | IB1;
return 1;
} else {
ranseed <<= 1;
return 0;
}
};
// *************************************************************
// getSeed
// Obtain the seed value from a collection of NBIT bits.
// This code is the inverse of ranBit.
// Input:
// int hbits[NBIT] -- global array of bits of shift register
// Return:
// seed value
// *************************************************************
unsigned int getSeed() {
int seedbits[NBIT];
unsigned int ranseed = 0;
if (NBIT != 24) {
cout << "ERROR: NBIT is not 24. This is unexpected." << endl;
cout << "Code failure..." << endl; // admittedly awkward, but
return 0; // you probably won't care.
}
for (int i = 0; i < 20; i++) seedbits[23-i] = hbits[i];
seedbits[3] = hbits[20]^seedbits[23];
seedbits[2] = hbits[21]^seedbits[22]^seedbits[23];
seedbits[1] = hbits[22]^seedbits[21]^seedbits[22];
seedbits[0] = hbits[23]^seedbits[20]^seedbits[21]^seedbits[23];
for (int i=NBIT-1; i >= 0; i--) ranseed = ranseed << 1|(seedbits[i]&1);
ranseed = ranseed&0xFFFFFF;
return ranseed;
}