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

Merge branch 'Cubix' into 'preprod'

New Cubix version, including Radwares utilities

See merge request IPNL_GAMMA/gammaware!43
parents 71674e3a 58f2e6c7
Pipeline #38413 passed with stages
in 5 minutes and 58 seconds
......@@ -57,6 +57,21 @@ link_directories(${CMAKE_LIBRARY_OUTPUT_DIRECTORY})
#
# configure file
#
#TEST current configuration OS
if ( UNIX )
set(CMAKE_OS_TYPE 1)
endif()
if ( WIN32 )
set(CMAKE_OS_TYPE 2)
endif()
if ( CYGWIN )
set(CMAKE_OS_TYPE 3)
endif()
if ( APPLE )
set(CMAKE_OS_TYPE 4)
endif()
configure_file (
"${CMAKE_CURRENT_SOURCE_DIR}/GwConfig.h.cmake"
"${CMAKE_BINARY_DIR}/include/GwConfig.h"
......
......@@ -43,3 +43,14 @@
#define GW_MacrosPath "@GW_DATA_DIR@/macros";
#endif
//////////////////////////////
////// Define OS_TYPE ////////
//////////////////////////////
#define OS_LINUX 1
#define OS_WINDOWS 2
#define OS_CYGWIN 3 // rather useless ==> to be removed ??
#define OS_APPLE 4
#define OS_TYPE @CMAKE_OS_TYPE@
......@@ -38,7 +38,7 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -pthread -Wno-unused-paramete
# root
INCLUDE_DIRECTORIES( ${ROOT_INCLUDE_DIR} )
LINK_DIRECTORIES( ${ROOT_LIBRARY_DIR} )
SET(EXTRA_EXTERNAL_LIBRARIES ${EXTRA_EXTERNAL_LIBRARIES} Matrix)
SET(EXTRA_EXTERNAL_LIBRARIES ${EXTRA_EXTERNAL_LIBRARIES} Matrix Rint)
# gw
set(EXTRA_INTERNAL_LIBRARIES GWPHYSICS )
......
ENCAL OUTPUT FILE
TMEUHS.STO EU152.SOU TM 159 EU152 HIGH SPIN CALIBRATION
2 0.000000000000D+00 1.000000000000D+00 0.000000000000D+00
0.000000000000D+00 0.000000000000D+00 0.000000000000D+00
### Cube file
CubeFileName MyCube.cub
### 2D projection file
2DProj MyCube.2dp
### LUT file
LUT MyCube.tab
### Contraction factor (default 1)
CompressFact 1
### Total projection file
TotalProj MyCube.spe
### Backgroud file
BackgroudFile MyCubeBG.spe
### Calibration file (not mandatory)
#CalFile MyCube.aca
### Efficiency file (not mandatory)
#EffFile MyCube.aef
#######################
### 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
# 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
# Must be similar to the value in the GenBinning conf file
CompressionFactor 1
################################
## 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
EFFIT PARAMETER FILE 19-FEB-92 10:14:04
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+02 0.00000000E+03
#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 1
#RangeMax (keV)
RangeMax 8192
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}
obj_CubeBuilder = src/CubeBuilder.o src/progressbar.o
obj_GenBinning = src/GenBinning.o src/util.o
ALL = CubeBuilder GenBinning
.PHONY: all
all: $(ALL)
CubeBuilder: $(obj_CubeBuilder)
$(CXX) $^ -o $@ $(LDFLAGS) $(CXXFLAGS)
GenBinning: $(obj_GenBinning)
$(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
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 <unistd.h>
#include <math.h>
#include <string.h>
#include <fstream>
#include "util.h"
#include "TGraph.h"
#include "TF1.h"
#include "TString.h"
#include "TCanvas.h"
#include "TApplication.h"
#include "TSystem.h"
using namespace std;
float wid_spec[16384];
float wa = 3.0f, wb = 1.0f, wc = 4.0f, cfactor = 1.0f, chs = 2.0f;
float w, x, f1, f2, f3, wa2, wb2, wc2;
int iclo = 200, ichi = 5000;
int i1, i2, i3, hi, lo, nc, i, j, ich2;
int nclook, lookmin, lookmax, rl = 0;
short looktab[16384];
char buf[32];
TString TabName = "FWHM.tab";
Bool_t fFixedBinning = false;
bool CalcFWHMParameters(vector< pair<float,float> > calibpoints) {
TApplication *theApp = new TApplication("App", 0, 0);
TGraph *CalibGraph = new TGraph;
CalibGraph->SetNameTitle("CalibFWHM","CalibFWHM");
TF1 *CalibFunc = new TF1("CalibFunc","sqrt([0]*[0] + [1]*[1]*x/1000. + [2]*[2]*x/1000.*x/1000.)");
for(uint i=0 ; i<calibpoints.size() ; i++) {
CalibGraph->SetPoint(CalibGraph->GetN(),calibpoints.at(i).first,calibpoints.at(i).second);
}
TCanvas *c = new TCanvas;
CalibGraph->SetMarkerStyle(20);
CalibGraph->SetMarkerSize(2);
CalibGraph->SetMarkerColor(kBlue);
CalibGraph->Draw("APL");
CalibGraph->Fit(CalibFunc,"SQ");
wa2 = CalibFunc->GetParameter(0)*CalibFunc->GetParameter(0);
wb2 = CalibFunc->GetParameter(1)*CalibFunc->GetParameter(1);
wc2 = CalibFunc->GetParameter(2)*CalibFunc->GetParameter(2);
cout<<" -- Fit Function: FWHM = sqrt(F*F + G*G*x + H*H*x*x) (x = ch.no. /1000)"<<endl;
cout<<" --> F = "<< CalibFunc->GetParameter(0)<<endl;
cout<<" --> G = "<< CalibFunc->GetParameter(1)<<endl;
cout<<" --> H = "<< CalibFunc->GetParameter(2)<<endl;
c->Update();
c->Modified();
gSystem->Sleep(200);
gSystem->ProcessEvents();
gSystem->DispatchOneEvent(true);
bool ok = askyn("Agree with the fit ? (Y/N)");
delete theApp;
return ok;
}
int main(int argnc, char **argv)
{
FILE *tabfile = nullptr;
cout<<" -- This program is used to create lookup tables for incub8r and 4play replay tasks"<<endl;
cout<<" Incub8r and 4play allow non-linear gains"<<endl;
cout<<" This is an adaptation of the lufwhm program of the Radware package"<<endl;
if(argnc != 2)
{
cout<<" -- To generate a binning file, use the command: ./gen_binning Conf_binning.txt"<<endl;
cout<<" --> see an example of conf file in conf/Conf_binning.txt"<<endl;
return 1;
}
vector< pair<float,float> > calibpoints;
pair<int, int> Range(0,8191);
Float_t NumberOfChannelPerFWHM = 1;
string line;
TString Buffer;
ifstream FileConf;
FileConf.open(argv[1], std::ifstream::in);
if(!FileConf) {
cout<<"No binning conf file ==> Exit"<<endl;
return 1;
}
int linenb=0;
while(FileConf) {
getline(FileConf,line);
linenb++;
Buffer = line;
TObjArray *loa=Buffer.ReplaceAll("\t"," ").Tokenize(" ");
if(Buffer.BeginsWith("#") || Buffer.Length()==0 ) continue;
if(Buffer.BeginsWith("TabFile")) {
if(loa->GetEntries() != 2) {
cout<<"Bad entry at line "<<linenb<<" ; skipped"<<endl; continue;
}
TabName = loa->At(1)->GetName();
if(!TabName.EndsWith(".tab")) TabName.Append(".tab");
}
if(Buffer.BeginsWith("FWHM")) {
if(loa->GetEntries() != 3) {
cout<<"Bad entry at line "<<linenb<<" ; skipped"<<endl; continue;
}
pair<float,float> value (((TString)loa->At(1)->GetName()).Atof(),((TString)loa->At(2)->GetName()).Atof());
calibpoints.push_back(value);
}
if(Buffer.BeginsWith("Precision")) {
if(loa->GetEntries() != 2) {
cout<<"Bad entry at line "<<linenb<<" ; skipped"<<endl; continue;
}
NumberOfChannelPerFWHM = ((TString)loa->At(1)->GetName()).Atof();
}
if(Buffer.BeginsWith("CompressionFactor")) {
if(loa->GetEntries() != 2) {
cout<<"Bad entry at line "<<linenb<<" ; skipped"<<endl; continue;
}
cfactor = ((TString)loa->At(1)->GetName()).Atof();
}
if(Buffer.BeginsWith("RangeMax")) {
if(loa->GetEntries() != 2) {
cout<<"Bad entry at line "<<linenb<<" ; skipped"<<endl; continue;
}
Range.first = 0;
Range.second = ((TString)loa->At(1)->GetName()).Atoi();
}
if(Buffer.BeginsWith("FixedBinning")) {
fFixedBinning = true;
}
delete loa;
}
tabfile = open_new_file(((char*)TabName.Data()), 1);
cout<<"Output tab file: "<<TabName.Data()<<endl;
chs = NumberOfChannelPerFWHM;
iclo = Range.first;
ichi = Range.second*cfactor;
cout<<" *-* "<<chs<<" channels per FWHM"<<endl;
cout<<" *-* Compression factor: "<<cfactor<<endl;
cout<<" *-* Channels Range: "<<iclo<<" -- "<<ichi<<endl;
if(fFixedBinning) {
cout<<" *-* Fixed binning used"<<endl;
wa2=1.;
wb2=0.;
wc2=0.;
}
else if(CalcFWHMParameters(calibpoints) == false)
return 1;
lo = iclo;
ich2 = 0;
while (lo <= ichi) {
++ich2;
x = ((float) lo + 0.5f) / 1e3f;
w = sqrt(wa2 + wb2 * x/cfactor + wc2 * x/cfactor * x/cfactor) * cfactor;
int incr = (int) ((w / chs) + 0.5);
if(incr==0) incr = 1;
lo += incr;
}
if (!askyn("No. of chs. in cube = %d --- Okay? (Y/N)", ich2))
return 1;
nclook = ichi+1;
for (j = 0; j < nclook; ++j) {
looktab[j] = 0;
}
lo = iclo;
for (i = 1; i < ich2+1; ++i) {
x = ((float) lo + 0.5f) / 1e3f;
w = sqrt(wa2 + wb2 * x/cfactor + wc2 * x/cfactor * x/cfactor) * cfactor;
hi = lo + (int) ((w / chs) + 0.5);
if (hi > nclook) hi = nclook;
for (j = lo; j < hi; ++j) {
looktab[j] = i;
}
wid_spec[i - 1] = (float) (hi - lo)/cfactor;
if ((lo = hi) == nclook) break;
}
/* write out look-up table to disk file */
lookmin = 1;
lookmax = ich2;
rewind(tabfile);
#define W(a,b) { memcpy(buf + rl, a, b); rl += b; }
W(&nclook, 4);
W(&lookmin, 4);
W(&lookmax, 4);
#undef W
if (put_file_rec(tabfile, buf, rl) ||
put_file_rec(tabfile, looktab, 2*nclook)) {
file_error("write to", ((char*)TabName.Data()));
return 1;
}
fclose(tabfile);
TString SpectrumFileName = TabName.Copy().Append(".spe");
tell("The \"Width Spectrum\" has %d channels, and contains the width (in ADC-chs) for each cube channel.\n", ich2);
wspec((char*)SpectrumFileName.Data(), wid_spec, ich2);
tell("=> Width Spectrum written as file %s\n", SpectrumFileName.Data());
return 0;
}
#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;
}