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

Starting Adding radware utilities

parent 5cecf7ef
Pipeline #32647 failed with stage
in 2 minutes and 37 seconds
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
This diff is collapsed.
#define RW_MAXFILES 9999 /* Max files per tape to read */
#define RW_MAXRECORDS 1000000000 /* Max number of records per tape to read */
#define RW_DATFILE "incub8r.inc" /* Defaults file for hypercube */
#define RW_CHNK_S 32 /* size of update chunk, in units of 1/4 MB */
/*
RW_LB2 = ceil((float)(ceil((float)((RW_MAXCH/8)*((RW_MAXCH/8)+1)*((RW_MAXCH/8)+2)/6)/128))/RW_CHNK_S)
RW_LB1 = RW_LB2*RW_CHNK_S
*/
//for 1400 channels
//#define RW_MAXCH 1400 /* Maximum channels on cube axis */
// /* (must be multiple of 8) */
//#define RW_LB2 222 /* length required for buf2, depends on RW_LB1 */
// /* RW_LB2 >= RW_LB1/RW_CHNK_S */
//#define RW_LB1 7104 /* length required for buf1, depends on RW_MAXCH */
// /* RW_LB1 > x(x+1)(x+2)/(6*128), */
// /* where x = RW_MAXCH/8 */
// /* must also be a multiple of RW_CHNK_S */
//for 4096channels
//#define RW_MAXCH 4096
//#define RW_LB1 175808
//#define RW_LB2 5494
//for 8196 channels
#define RW_MAXCH 8192
#define RW_LB1 1402208
#define RW_LB2 43819
#define RW_DB1 90 /* depth of buf1 */
#define RW_DB2 180 /* depth of buf2 not to exceed RW_SCR_RECL */
/* byte records in scr file */
#define RW_SCR_RECL 32768 /* Record length of scr file */
/* must be >= RW_DB2*2*(RW_DB1+1) + 4 */
#define RW_MAXMULT 40
/************************************************
* 3D CUBE FILE FORMAT
* 1024byte file header
* 4096byte data records
*
* each data record contains:
* variable number of bit-compressed 8x8x4 mini-cubes
* 8bytes header, 4088 bytes of data
*/
/* 3d cube file header */
typedef struct {
char id[16]; /* "Incub8r3/Pro4d " */
int numch; /* number of channels on axis */
int bpc; /* bytes per channel, = 4 */
int cps; /* 1/cps symmetry compression, = 6 */
int numrecs; /* number of 4kB data records in the file */
char resv[992]; /* FUTURE flags */
} FHead3D;
/* 3Drecord header */
typedef struct {
int minmc; /* start minicube number, starts at 0 */
short nummc; /* number of minicubes stored in here */
short offset; /* offset in bytes to first full minicube */
} RHead3D;
typedef struct {
RHead3D h;
unsigned char d[4088]; /* the bit compressed data */
} Record3D; /* see the compression alg for details */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int
compress3d (unsigned int *in, unsigned char *out)
{
static int sum_limits[] = { 16,32,64,128,
256,512,1024,2048,
4096,8192,17000,35000,
80000,160000,320000,640000 };
static int bitmask[] = { 0,1,3,7,15,31,63,127,255,511,1023,2047,
4095,8191,16383,32767 };
int cube_length;
int sum=0, i, j;
int nplanes, noverflows, stmp, nbits;
unsigned char *bitptr, *ovrptr;
int test1, test2, test3;
/* guess at depth of bit field */
for (i=0;i<256;i+=16) {
sum += in[i];
}
i = 7;
j = 8;
while (j>1) {
j = j / 2;
if (sum > sum_limits[i])
i += j;
else
i -= j;
}
if (sum > sum_limits[i+1])
nplanes = i + 1;
else
nplanes = i;
while (1) {
test1 = test2 = test3 = 0;
for (i=0;i<256;i++) {
test1 += in[i] >> nplanes;
test2 += in[i] >> (nplanes+1);
test3 += in[i] >> (nplanes+2);
}
if (test2 > 31) {
if (test3 > 31) {
nplanes += 3;
continue;
}
nplanes += 2;
noverflows = test3;
break;
}
else if (test1 > 31) {
nplanes += 1;
noverflows = test2;
break;
}
noverflows = test1;
break;
}
if (nplanes > 30)
fprintf(stderr,"Expecting core dump...\n");
/* insert length of compressed cube */
if (nplanes < 7) {
out[0] = 32*nplanes + noverflows;
bitptr = out+1;
cube_length = 1;
}
else {
out[0] = 224 + noverflows;
out[1] = nplanes-7;
bitptr = out+2;
cube_length = 2;
}
ovrptr = bitptr + nplanes*32;
cube_length += nplanes*32 + noverflows;
/* now, compress */
/* prepare bit planes and insert overflow bits... */
while (nplanes>=8) {
for (i=0; i<256; i++) {
*bitptr++ = in[i]&bitmask[8];
in[i] = in[i]>>8;
}
nplanes -= 8;
}
if (nplanes > 0) {
stmp = 0;
nbits = 0;
for (i=0; i<256; i++) {
/* insert nplanes number of bits */
stmp = (stmp << nplanes) + (in[i] & bitmask[nplanes]);
nbits += nplanes;
if (nbits > 7) {
*bitptr++ = stmp >> (nbits - 8);
nbits -= 8;
stmp &= bitmask[nbits];
}
/* append overflows */
noverflows = in[i] >> nplanes;
for(j=0; j<noverflows; j++)
*ovrptr++ = i;
}
}
else { /* just do overflows */
for (i=0; i<256; i++) {
for(j=0; j<(int)in[i]; j++) {
*ovrptr++ = i;
}
}
}
return cube_length;
}
void
decompress3d (unsigned char in[1024], unsigned int out[256])
{
static int bitmask[] = { 0,1,3,7,15,31,63,127,255,511,1023,2047,
4095,8191,16383,32767 };
int nplanes, noverflows, nbits, savenplanes;
unsigned char *bitptr;
int i, j, t;
nplanes = in[0] >> 5;
noverflows = in[0] & 31;
if (nplanes == 7) {
nplanes += in[1];
bitptr = in+2;
}
else {
bitptr = in+1;
}
/* printf("%d %d %d\n",nplanes,noverflows,*bitptr); */
/* extract bit planes */
savenplanes = nplanes;
for (i=0; i<256; i++)
out[i] = 0;
j = 0;
while (nplanes>=8) {
for (i=0; i<256; i++)
out[i] += ((unsigned int)*bitptr++)<<j;
nplanes -= 8;
j += 8;
}
if (nplanes > 0) {
nbits = 0;
for (i=0; i<256; i++) {
if (nbits+nplanes < 9) {
out[i] += ((*bitptr >> (8-nbits-nplanes)) & bitmask[nplanes])<<j;
nbits += nplanes;
if (nbits > 7) {
bitptr++;
nbits -= 8;
}
}
else {
t = nplanes-8+nbits;
out[i] += (((*bitptr & bitmask[8-nbits]) << t)
+ (*(bitptr+1) >> (8-t)))<<j;
bitptr++;
nbits = t;
}
}
}
/* extract overflows */
for (i=0; i<noverflows; i++)
out[*bitptr++] += 1 << savenplanes;
}
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <string.h>
#include "util.h"
#include "minig.h"
#include "gls.h"
#include "esclev.h"
extern FILE *infile, *cffile; /* command file flag and file descriptor*/
extern int cflog; /* command logging flag */
/* This file contains subroutines that are common to all versions */
/* of escl8r and levit8r. */
/* D.C. Radford Sept 1999 */
/* function defined in drawstring.c */
int drawstring(char *, int, char, float, float, float, float);
/* functions defined in escl8ra.c or lev4d.c */
int combine(char *ans, int nc); /* returns 1 on error */
int disp_dat(int iflag);
int energy(float x, float dx, float *eg, float *deg);
int examine(char *ans, int nc); /* returns 1 on error */
int num_gate(char *ans, int nc); /* returns 1 on error */
char listnam[] =
"ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz[]|";
/* ======================================================================= */
float get_energy(void)
{
float eg, x, y, dx, fj1, deg;
int nc;
char ans[80];
eg = -1.0f;
printf("Enter gamma energy using cursor on spectrum window.\n"
" ...Hit any character;\n"
" T to type energy, X or button 3 to exit.\n");
retic(&x, &y, ans);
if (*ans == 'X' || *ans == 'x') return 0;
if (*ans == 'T' || *ans == 't') {
nc = ask(ans, 80, "Gamma energy = ? (rtn to quit)");
if (nc > 0) ffin(ans, nc, &eg, &fj1, &fj1);
} else {
dx = 0.0f;
energy(x - 0.5f, dx, &eg, &deg);
}
return eg;
} /* get_energy */
/* ======================================================================= */
int get_fit_gam(int *jng, short *jgam)
{
float x1, y1;
int igam, ilev;
int i, k, iband, askit, nc, ki;
char li[80], ans[80], new_ans[80];
/* get gamma indices to be fitted in intensity and/or energy */
/* ask for gammas to be fitted */
askit = 0;
START:
if (*jng <= 0 || askit ||
(*ans != 'L' && *ans != 'l' && *ans != 'B' && * ans != 'b')) {
printf(" Gammas to be fitted = ?\n"
" ...X to end,\n"
" ...button 2 or L to specify initial level(s),\n"
" ...button 3 or B to specify initial band(s),\n"
" ...A for all gammas in level scheme.\n");
retic2(&x1, &y1, ans);
if (!strcmp(ans, "X") || !strcmp(ans, "x")) {
*jng = 0;
return 0;
}
}
askit = 1;
*jng = 0;
if (*ans == 'A' || *ans == 'a') {
for (i = 0; i < glsgd.ngammas; ++i) {
jgam[(*jng)++] = (short) i;
}
} else if (*ans == 'L' || *ans == 'l') {
while (1) {
if (!(nc = ask(new_ans, 80,
" ...Name of initial level = ? (rtn to end)"))) goto START;
if (nc >= 3 && nc <= 10) {
strcpy(li, " ");
strcpy(li + 10 - nc, new_ans);
for (ki = 0; ki < glsgd.nlevels; ++ki) {
if (!strcmp(glsgd.lev[ki].name, li)) break;
}
if (ki < glsgd.nlevels) break;
}
printf("*** Bad level name. ***\n");
}
for (k = 0; k < levemit.n[ki]; ++k) {
jgam[(*jng)++] = (short) levemit.l[ki][k];
}
} else if (*ans == 'B' || *ans == 'b') {
while (1) {
nc = ask(new_ans, 80, " ...Name of initial band = ? (rtn to end)");
if (nc == 0) goto START;
if (nc < 9) break;
printf("*** Bad band name. ***\n");
}
strcpy(li, " ");
strcpy(li + 8 - nc, new_ans);
for (ki = 0; ki < glsgd.nlevels; ++ki) {
if (!strncmp(glsgd.lev[ki].name, li, 8)) {
for (k = 0; k < levemit.n[ki]; ++k) {
jgam[(*jng)++] = (short) levemit.l[ki][k];
}
}
}
} else if (*ans == 'G') {
igam = nearest_gamma(x1, y1);
if (igam < 0) {
printf("No gamma selected, try again...\n");
goto START;
}
*jng = 1;
jgam[0] = (short) igam;
} else if (!strncmp(ans, "X2", 2)) {
ilev = nearest_level(x1, y1);
if (ilev < 0) {
printf("No level selected, try again...\n");
goto START;
}
for (k = 0; k < levemit.n[ilev]; ++k) {
jgam[(*jng)++] = (short) levemit.l[ilev][k];
}
} else if (!strncmp(ans, "X3", 2)) {
iband = nearest_band(x1, y1);
if (iband < 0 || iband >= glsgd.nbands) {
printf("No band selected, try again...\n");
goto START;
}
for (ilev = 0; ilev < glsgd.nlevels; ++ilev) {
if (glsgd.lev[ilev].band == iband) {
for (k = 0; k < levemit.n[ilev]; ++k) {
jgam[(*jng)++] = (short) levemit.l[ilev][k];
}
}
}
} else {
printf(" *** Bad entry, try again... ***\n");
goto START;
}
return 0;
} /* get_fit_gam */
/* ======================================================================= */
int listgam(void)
{
/* defined in esclev.h:
char listnam[55] = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz[]|" */
float x, y, eg, dx, fj1, fj2, deg;
int igam, iwin, graphics = 0, i, k, ni, nl;
char reply, ans[80];
/* ....define/edit list of gammas */
while (1) {
if (!ask(ans, 80, "List name = ? (A to Z, a to z)")) return 0;
reply = *ans;
for (nl = 0; nl < 52; ++nl) {
if (reply == listnam[nl]) break;
}
if (nl < 52) break;
}
ni = elgd.nlist[nl];
i = 0;
while (i < 40) {
if (graphics) {
/* get list energies from level scheme or spectrum window */
if (i < ni) {