Commit ce136134 authored by Jérémie Dudouet's avatar Jérémie Dudouet
Browse files

Starting Adding radware utilities

parent 5037e0aa
ROOTCFLAGS := $(shell root-config --cflags)
ROOTLIBS := $(shell root-config --libs)
CXX = g++
CXXFLAGS = -O3 -I . -w
CXXFLAGS = -O3 -I ./src -w -g $(ROOTCFLAGS) -Wno-write-strings -Wno-unused-result
LDFLAGS = ${ROOTLIBS}
#src/levit8ra.o src/util.o src/esclev.o src/lev4d.o
obj_levit8ra = src/levit8ra.o src/util.o
obj_CubeBuilder = src/CubeBuilder.o src/progressbar.o
obj_gen_binning = src/gen_binning.o src/util.o
ALL = CubeBuilder GenBinning levit8ra
.PHONY: all
all: $(ALL)
levit8ra: $(obj_levit8ra)
$(CXX) $^ -o $@ $(LDFLAGS) $(CXXFLAGS)
CubeBuilder: $(obj_CubeBuilder)
$(CXX) $^ -o $@ $(LDFLAGS) $(CXXFLAGS)
GenBinning: $(obj_gen_binning)
$(CXX) $^ -o $@ $(LDFLAGS) $(CXXFLAGS)
.PHONY: clean
clean:
rm -f $(ALL) src/*.o
ReadMe CubeBuilder:
1 - Generate the adaptative binning file:
./GenBinning ConfBinning.txt
This will produce:
- The binning file: Binning.tab
- A width spectrum file: Binning.spe
2 - Generate the Cube
./CubeBuilder ConfBuildCube.conf
This will produce:
- The cube file: MyCube.cub
3 - Generate the total projection
(radware prog) pro3d => MyCube.cub
This will produce:
- The total 2D projection spectrum: MyCube.2dp
- The total 1D projection spectrum: MyCube.spe
4 - Generate the background (in gf3)
in gf3:
4-1: Open the 1D projection: MyCube.spe
4-2: Divide by the width spectrum: DV, then Binning.spe
4-3: Calculate background: BG, then ok, then click limits on the graph
4-4: Multiply by the width spectrum: MS, then Binning.spe
4-5: Save the background spectrum: WS, then file name
#Name of the tab file
TabFile test.tab
#To fit the FWHM vs Energy function:
#Fit function: sqrt([0]*[0] + [1]*[1]*x/1000. + [2]*[2]*x/1000.*x/1000.)
#For each new point: FWHM Energy (keV) Value (keV)
FWHM 122. 1.35
FWHM 1333. 2.35
FWHM 2500. 3.00
#Compression factor (to obtain binning <1 keV, one needs to scale the energy)
#For example, to obtain a minimal 0.5 KeV/Channel, we fill the Cube on 8k channel for 4keV range ==> Compression factor = 2
CompressionFactor 1
#Number of bin per FWHM
Precision 2
#Range
Range 0 8191
#######################
### Cube Properties ###
#######################
# Name of the 3D cube file
CubeFileName MyCube.cub
# name of the look-up table file to be used to map ADC chs to cube chs
TabFile test.tab
# Size to be used for the "scratch" disk file, in MB
ScratchSize 512
################################
## Input Root Tree Properties ##
################################
#Name of the input TTree (regular expressions can be used as for building TChain)
TreeFile MyTree.root
#Name of the TTree
TreeName TreeName
#Events to process (0 => All)
NEvents 1000
#Name of the branch corresponding to the gamma-ray energy
EGamma EGamma
#Name of the branch corresponding to the gamma-ray muliplicity
GammaMult GammaMult
#Cut to be apllied (see TTreeFormulaExpression)
Cut None
### Cube file
CubeFileName MyCube.cub
### 2D projection file
2DProj MyCube.2dp
### LUT file
LUT test.tab
### Contraction factor (default 1)
CompressFact 1
### Calibration file (not mandatory)
#CalFile MyCube.aca
### Efficiency file (not mandatory)
#EffFile MyCube.aef
/***************************************************************************
incub8r3: Replay program for compressed triple coincidence cubes.
Allows use of non-linear gains through look-up tables
to convert gain-matched ADC outputs to cube channels.
Uses disk buffering to .scr "scratch" file.
This version for 1/6 cubes.
D.C Radford Sept 1997
****************************************************************************/
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <sys/ioctl.h>
#ifdef VMS
#include <unixio.h>
#else
/* #include <sys/mtio.h> */
#endif
#include <math.h>
#include <string.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <fcntl.h>
#include <time.h>
#include <signal.h>
#include "CubeBuilder.h"
#include "rwutils.cpp"
#include "progressbar.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TTree.h"
#include "TFile.h"
#include "TString.h"
#include "TMath.h"
#include "TChain.h"
#include "TTreeFormula.h"
//#define FillRandom
#define GetEventFromTree
using namespace std;
/* Some global data */
struct GlobalData {
FILE *CubeFile;
FILE *ScrFile;
FILE *LogFile;
char CubeFileName[256];
int ScrSize;
char TapeDev[256];
int SwapTapeBytes;
int length; /* length of axis on cube */
int nummc; /* number of minicubes */
int numscr; /* number of records written to scrfile */
int adcmax;
int numrecs; /* current number of records in cube file */
int gotsignal; /* a signal was caught by signal handler */
} gd;
int dbinfo[10];
/* lookup table stuff */
char *tbuf; /* tape buffer */
short luch[16384]; /* look-up table, maps ADC channels to cube channels */
int lumx[RW_MAXCH]; /* look-up table, maps 3d ch to linear minicube */
int lumy[RW_MAXCH];
int lumz[RW_MAXCH];
/* look up minicube addr */
#define LUMC(x,y,z) lumx[x] + lumy[y] + lumz[z]
/* look up addr in minicube */
#define LUIMC(x,y,z) (x&7) + ((y&7)<<3) + ((z&3)<<6)
/* increment buffers */
unsigned short buf1[RW_LB1][RW_DB1];
unsigned short buf2[RW_LB2][RW_DB2][RW_DB1+1];
int nbuf1[RW_LB1]; /* number in each row of buf1 */
int nbuf2[RW_LB2]; /* number in each row of buf2 */
int scrptr[RW_LB2]; /* last record # written to scr file */
/* for each row of buf2 */
#define MCLEN(mcptr) (mcptr[0] + (mcptr[0] >= 7*32 ? mcptr[1]*32+2 : 1))
FILE *open_3cube (char *name, int length, FILE *mesg);
FILE *open_scr (char *name, int length, FILE *mesg);
TString fCubeFileName;
TString fTabFile;
Int_t fScratchSize;
TString fTreeFileName="";
TString fTreeName="";
TString fEGammaName="";
TString fGammaMultName="";
TString fCut="";
ULong64_t fNEvents=0;
ULong64_t fCurrentEntry = 0;
ULong64_t fNumberOfFilledEvents=0;
TChain *fChain = nullptr;
TBranch *b_EGamma = nullptr;
TBranch *b_MultGamma = nullptr;
Int_t fGammaMult;
Float_t fEGamma[500];
TTreeFormula *fTreeFormula = nullptr;
/***************************************************************************
setup_replay: reads in data file that describes cube, opens and
initializes cube, scr and log
*/
int setup_replay(void)
{
char cbuf[256], title[256], tnam[256], snam[256], cnam[256];
char lnam[256], scrs[256];
int nclook, lmin, lmax, i, ovrtprompt = 1;
cout<<"This is a replay program for compressed-format RadWare triple-coincidence cubes."<<endl;
cout<<"Allows use of non-linear gains through look-up tables"<<endl;
cout<<"to convert gain-matched ADC outputs to cube channels."<<endl;
cout<<"This version is for 1/6 cubes."<<endl;
cout<<"Author D.C Radford Sept 1997"<<endl;
cout<<"Sub-author J.Dudouet Jan 2019"<<endl;
strcpy(title,fCubeFileName.Data());
strcpy(snam,fCubeFileName.Copy().Append(".scr").Data());
strcpy(cnam,fCubeFileName.Copy().Append(".cub").Data());
strcpy(lnam,fTabFile.Data());
strcpy(scrs,TString::Itoa(fScratchSize,10).Data());
gd.SwapTapeBytes = 0;
gd.ScrSize = fScratchSize;
printf("\n Scratch file name = %s\n"
" Cube file name = %s\n"
" Look-up table file name = %s\n"
" Scratch File Size (MBytes) = %s\n",
snam,cnam,lnam,scrs);
/* test values */
if (gd.ScrSize < 1) {
printf("\007 ERROR -- scratch file size too small.\n");
return 1;
}
strcpy(tnam,cnam);
setfnamext(tnam,".log",1);
if ((gd.LogFile=fopen(tnam,"r"))) {
fclose(gd.LogFile);
if (ovrtprompt) {
printf(" WARNING: Log file (%s) already exists\n",tnam);
printf(" Enter new name or <Return> to overwrite: ");
if (fgets(cbuf,256,stdin) && cbuf[0] != '\n') {
cbuf[strlen(cbuf)-1]=0;
strcpy(tnam,cbuf);
setfnamext(tnam,".log",0);
return 10;
}
}
else
printf(" Overwriting Log file (%s).\n",tnam);
}
if (!(gd.LogFile=fopen(tnam,"w"))) {
printf("\007 ERROR -- Could not open %s for writing.\n",tnam);
printf(" Enter new log file name: ");
if (fgets(cbuf,256,stdin) && cbuf[0] != '\n') {
strncpy(tnam,cbuf,strlen(cbuf)-1);
setfnamext(tnam,".log",0);
return 11;
}
}
/* read in lookup table */
if (read_tab(lnam,&nclook,&lmin,&lmax,luch)) {
fclose(gd.LogFile);
cout<<"error reading lookup table"<<endl;
return 2;
}
gd.adcmax = nclook;
printf("\nNo. of values in lookup table = %d\n",lmax);
if (lmax > RW_MAXCH) {
printf("\007 ERROR -- number of values in lookup file is too"
"large (max = %d).\n",RW_MAXCH);
fclose(gd.LogFile);
return 3;
}
gd.length = lmax;
/* Inform the log file */
fprintf(gd.LogFile,"Replay log for program Incub8r3, %s\n\n",datim());
fprintf(gd.LogFile,"\n Scratch file name = %s\n"
" Cube file name = %s\n"
" Look-up table file name = %s\n"
" Scratch File Size (MBytes) = %s\n",
snam,cnam,lnam,scrs);
fprintf(gd.LogFile," No. of values in lookup table = %d\n\n",lmax);
fprintf(gd.LogFile,"-------------------------------------------\n\n");
/* calculate look-up tables to convert (x,y,z) cube address
to linearized minicube address. x<y<z
lum{x,y,z}[ch] returns the number of subcubes to skip over
to get to the subcube with ch in it.
*/
for (i=0;i<8;i++) {
lumx[i] = 0;
lumy[i] = 0;
lumz[i] = i/4;
}
for (i=8;i<gd.length;i++) {
lumx[i] = (i/8)*2;
lumy[i] = lumy[i-8] + lumx[i];
lumz[i] = lumz[i-8] + lumy[i];
}
/* open cube and scratch */
if (!(gd.CubeFile=open_3cube(cnam, gd.length, gd.LogFile))) {
fclose(gd.LogFile);
return 4;
}
if (!(gd.ScrFile = open_scr(snam, gd.ScrSize, gd.LogFile))) {
fclose(gd.CubeFile);
fclose(gd.LogFile);
return 5;
}
/* clear buf1 and buf2 */
for(i=0;i<RW_LB1;i++)
nbuf1[i] = 0;
for(i=0;i<RW_LB2;i++) {
nbuf2[i] = 0;
scrptr[i] = 0;
}
fflush(gd.LogFile);
}
/***************************************************************************
open_3cube: returns file descriptor on success or NULL on failure
name - file name on disk
length - number of channels on an axis
mesg - stream to print error messages (in addition to stdout)
*/
FILE *open_3cube (char *name, int length, FILE *mesg)
{
FILE *file;
int i, nummc;
FHead3D head;
Record3D rec;
if (!name)
return NULL;
if (!(file=fopen(name,"r+"))) {
if (!(file=fopen(name,"w+"))) {
printf("\007 ERROR -- Cannot open file %s.\n", name);
return NULL;
}
}
rewind(file);
nummc = (length + 7)/8;
nummc = nummc*(nummc+1)*(nummc+2)/3; /* minicubes in cube */
gd.nummc = nummc;
switch (fread(&head,1,1024,file)) {
default:
fprintf(stderr," WARNING -- partial header on cube. Junking it.\n");
case 0:
gd.numrecs = (nummc + 4087)/4088; /* note, an empty minicube is 1 byte */
strncpy(head.id, "Incub8r3/Pro4d ", 16);
head.numch = length;
head.bpc = 4;
head.cps = 6;
head.numrecs = gd.numrecs;
memset (head.resv, 0, 992);
fprintf(mesg," Creating new cube: %d channels per side...\n",length);
printf(" Creating new cube: %d channels per side...\n",length);
rewind(file);
if (!fwrite(&head,1024,1,file)) {
fprintf(stderr, "\007 ERROR -- Could not write cube header... "
"aborting.\n");
fclose(file);
return NULL;
}
/* initialize cube */
for (i=0;i<4088;i++)
rec.d[i] = 0;
for (i=0; i<gd.numrecs; i++) { /* loop through records */
rec.h.minmc = i*4088;
rec.h.nummc = 4088;
rec.h.offset = 8;
if(!fwrite(&rec,4096,1,file)) {
fprintf(stderr, "\007 ERROR -- Could not write record %d... "
"aborting.\n",i);
fclose(file);
return NULL;
}
}
fflush(file);
printf(" ...Done creating new cube.\n");
break;
case 1024:
fprintf(mesg," Checking existing cube file...\n");
printf(" Checking existing cube file...\n");
if (strncmp(head.id,"Incub8r3/Pro4d ",16) ||
head.bpc != 4 ||
head.cps != 6) {
fprintf(stderr,"\007 ERROR -- Invalid header... aborting.\n");
fclose(file);
return NULL;
}
else if (head.numch != length) {
fprintf(stderr,"\007 ERROR -- Different axis length in cube (%d)..\n",
head.numch);
fclose(file);
return NULL;
}
gd.numrecs = head.numrecs;
printf(" ...Okay, %d records.\n", gd.numrecs);
}
fprintf(mesg,
"Axis length of 3d cube is %d.\n"
"3d cube has %d minicubes and %d records.\n",
length, nummc, head.numrecs);
printf("Axis length of 3d cube is %d.\n"
"3d cube has %d minicubes and %d records.\n",
length, nummc, head.numrecs);
strcpy(gd.CubeFileName, name);
return file;
}
/***************************************************************************
open_scr: returns file descriptor on success or NULL on failure
name - file name on disk
length - size on Megabytes of scr file
mesg - stream to print error messages (in addition to stdout)
*/
FILE *open_scr (char *name, int length, FILE *mesg)
{
FILE *file;
int b[256], i;
if(!name || length < 1 || !(file=fopen(name,"w+")))
return NULL;
rewind(file);
fprintf(mesg," Making sure we have enough scratch space...\n");
printf(" Making sure we have enough scratch space...\n");
for (i=0;i<length*1024;i++) /* writing to make sure we have space */
if (!fwrite(b, 1024, 1, file)) {
fprintf(stderr,"\007 ERROR -- Could not allocate scratch space on disk"
"(%dM).\n",length);
fclose(file);
return NULL;
}
gd.numscr = 0;
rewind(file);
fprintf(mesg," ...Done.\n");
printf(" ...Done.\n");
return file;
}
#include "compress3d.cpp"
void dscr(void)
{
int numchunks, chunknum;
Record3D recIn, recOut;
int j, k, mc, mcl;
unsigned char *mcptrIn, *mcptrOut, cbuf[2048];
int minmc, maxmc;
int recnumIn, recnumOut, nbytes;
int addr13b, addr8b, addr21b;
int nmcperc=RW_CHNK_S*256;
unsigned short inbuf[RW_DB2+1][RW_DB1+1];
FHead3D filehead;
unsigned int *chunk;
FILE *OutFile;
char filname[512];
while (!(chunk=(unsigned int*)malloc(RW_CHNK_S*256*1024))) {
printf("Ooops, chunk malloc failed...\n"
" ... please free up some memory and press return.\n");
fgets(filname,256,stdin);
}
/* open new output file for incremented copy of cube */
strcpy(filname, gd.CubeFileName);
strcat(filname, ".tmp-increment-copy");
if (!(OutFile=fopen(filname,"w+"))) {
printf("\007 ERROR -- Cannot open new output cube file.\n");
exit (-1);
}
fseek (gd.CubeFile, 0, SEEK_SET);
if (!fread (&filehead, 1024, 1, gd.CubeFile)) {
printf("\007 ERROR -- Cannot read file header, aborting...\n");
exit(-1);
}
fseek (OutFile, 0, SEEK_SET);
if (!fwrite (&filehead, 1024, 1, OutFile)) {
printf("\007 ERROR -- Cannot write file header, aborting...\n");
exit(-1);
}
/* a chunk is RW_CHNK_S/16 Mchannels, nmcperc minicubes */
numchunks = (gd.nummc+nmcperc-1)/nmcperc;
fprintf(gd.LogFile," ...Updating cube from scratch file: %s\n", datim());
fflush(gd.LogFile);
printf( " ...Updating cube from scratch file: %s\n", datim());
printf("There are %d chunks to increment...\n",numchunks);
/* read in first record */
recnumIn = 0;
fseek(gd.CubeFile, 1024, SEEK_SET);
if (!fread(&recIn, 4096, 1, gd.CubeFile)) {
printf("\007 ERROR -- Corrupted cube, aborting...\n");
exit(-1);
}
mcptrIn = recIn.d;
/* init the first output record */
recnumOut = 0;
recOut.h.minmc = recIn.h.minmc;
recOut.h.offset = 8;
mcptrOut = recOut.d;
/* loop through all the chunks in the file */
for (chunknum=0; chunknum<numchunks; chunknum++) {
minmc = chunknum*nmcperc;
maxmc = minmc+nmcperc-1;
if (maxmc>gd.nummc-1)
maxmc = gd.nummc - 1;
dbinfo[0]=11;