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

Update Fits

Update CubePlayer
Update CubeBuilder
parent 78f02a99
Pipeline #35722 failed with stage
in 2 minutes and 4 seconds
......@@ -7,24 +7,18 @@ 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
obj_GenBinning = src/GenBinning.o src/util.o
ALL = CubeBuilder GenBinning levit8ra
ALL = CubeBuilder GenBinning
.PHONY: all
all: $(ALL)
levit8ra: $(obj_levit8ra)
$(CXX) $^ -o $@ $(LDFLAGS) $(CXXFLAGS)
CubeBuilder: $(obj_CubeBuilder)
$(CXX) $^ -o $@ $(LDFLAGS) $(CXXFLAGS)
GenBinning: $(obj_gen_binning)
GenBinning: $(obj_GenBinning)
$(CXX) $^ -o $@ $(LDFLAGS) $(CXXFLAGS)
.PHONY: clean
......
#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
......@@ -97,6 +97,7 @@ FILE *open_scr (char *name, int length, FILE *mesg);
TString fCubeFileName;
TString fTabFile;
Int_t fScratchSize;
Float_t fCompressionFactor=1;
TString fTreeFileName="";
TString fTreeName="";
......@@ -766,6 +767,12 @@ int ReadConfFile(TString ConfFileName){
delete loa;
cout<<"\e[1;92mRaw Scratch-file size (MB): "<<fScratchSize<<"\e[0m"<<endl;
}
else if(Buffer.BeginsWith("CompressionFactor")){
TObjArray *loa=Buffer.Tokenize(" ");
fCompressionFactor = ((TString)loa->At(1)->GetName()).Atoi();
delete loa;
cout<<"\e[1;92mRaw Compression factor: "<<fCompressionFactor<<"\e[0m"<<endl;
}
else if(Buffer.BeginsWith("TreeFile")){
TObjArray *loa=Buffer.Tokenize(" ");
fTreeFileName = ((TString)loa->At(1)->GetName());
......@@ -886,7 +893,7 @@ int FillFromTree(int *elist){
}
for(int i=0 ; i<fGammaMult ; i++) {
elist[i] = TMath::Nint(fEGamma[i]);
elist[i] = TMath::Nint(fEGamma[i]*fCompressionFactor);
// cout<<i<<"/"<<fGammaMult<<": "<<elist[i]<<" keV"<<endl;
}
// cout<<endl;
......
......@@ -129,7 +129,7 @@ int main(int argnc, char **argv)
if(loa->GetEntries() != 2) {
cout<<"Bad entry at line "<<linenb<<" ; skipped"<<endl; continue;
}
NumberOfChannelPerFWHM = ((TString)loa->At(1)->GetName()).Atoi();
NumberOfChannelPerFWHM = ((TString)loa->At(1)->GetName()).Atof();
}
if(Buffer.BeginsWith("CompressionFactor")) {
if(loa->GetEntries() != 2) {
......@@ -137,12 +137,12 @@ int main(int argnc, char **argv)
}
cfactor = ((TString)loa->At(1)->GetName()).Atof();
}
if(Buffer.BeginsWith("Range")) {
if(loa->GetEntries() != 3) {
if(Buffer.BeginsWith("RangeMax")) {
if(loa->GetEntries() != 2) {
cout<<"Bad entry at line "<<linenb<<" ; skipped"<<endl; continue;
}
Range.first = ((TString)loa->At(1)->GetName()).Atoi();
Range.second = ((TString)loa->At(2)->GetName()).Atoi();
Range.first = 0;
Range.second = ((TString)loa->At(1)->GetName()).Atoi();
}
if(Buffer.BeginsWith("FixedBinning")) {
fFixedBinning = true;
......@@ -157,7 +157,7 @@ int main(int argnc, char **argv)
chs = NumberOfChannelPerFWHM;
iclo = Range.first;
ichi = Range.second;
ichi = Range.second*cfactor;
cout<<" *-* "<<chs<<" channels per FWHM"<<endl;
cout<<" *-* Compression factor: "<<cfactor<<endl;
cout<<" *-* Channels Range: "<<iclo<<" -- "<<ichi<<endl;
......@@ -175,8 +175,8 @@ int main(int argnc, char **argv)
ich2 = 0;
while (lo <= ichi) {
++ich2;
x = ((float) lo / cfactor + 0.5f) / 1e3f;
w = sqrt(wa2 + wb2 * x + wc2 * x * x) * cfactor;
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;
......@@ -191,14 +191,14 @@ int main(int argnc, char **argv)
}
lo = iclo;
for (i = 1; i < ich2+1; ++i) {
x = ((float) lo / cfactor + 0.5f) / 1e3f;
w = sqrt(wa2 + wb2 * x + wc2 * x * x) * cfactor;
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);
wid_spec[i - 1] = (float) (hi - lo)/cfactor;
if ((lo = hi) == nclook) break;
}
......@@ -218,10 +218,10 @@ int main(int argnc, char **argv)
}
fclose(tabfile);
tell("The \"Width Spectrum\" has %d channels, and contains\n"
" the width (in ADC-chs) for each cube channel.\n", ich2);
wspec((char*)TabName.ReplaceAll(".tab","").Data(), wid_spec, ich2);
tell(" Width Spectrum written as file %s\n", TabName.ReplaceAll(".tab","").Append(".spe").Data());
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 <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) {
printf("Gamma %c%2.2i? (P for %.1f, T to type energies, X to end)\n",
reply, i+1, elgd.list[nl][i]);
} else {
printf("Gamma %c%2.2i? (T to type energies, X to end)\n",
reply, i+1);
}
retic3(&x, &y, ans, &iwin);
if (*ans == 'X' || *ans == 'x') {
elgd.nlist[nl] = i;
return 0;
}
if (*ans == 'T' || *ans == 't') {
graphics = 0;
} else if (i >= ni || (*ans != 'P' && *ans != 'p')) {
if (iwin == 2) {
if ((igam = nearest_gamma(x, y)) < 0) {
printf(" ... no gamma selected, try again...\n");
continue;
}
eg = glsgd.gam[igam].e;
} else {
dx = 0.0f;
energy(x - 0.5f, dx, &eg, &deg);
}
elgd.list[nl][i++] = eg;
} else {
i++;
}
} else {
/* get list energies from text window */
if (i < ni) {
k = ask(ans, 80, "Gamma %c%2.2i = ? (1 for %.1f,\n"
" G to enter energies from graphics, rtn to end)",
reply, i+1, elgd.list[nl][i]);
} else {
k = ask(ans, 80, "Gamma %c%2.2i = ?\n"
" (G to enter energies from graphics, rtn to end)",
reply, i+1);
}
if (k < 1) {
elgd.nlist[nl] = i;
return 0;
}
if (*ans == 'G' || *ans == 'g') {
graphics = 1;
} else if (ffin(ans, k, &eg, &fj1, &fj2)) {
printf(" ... bad entry, try again...\n");
} else if (i >= ni || eg != 1.0f) {
elgd.list[nl][i++] = eg;
} else {
i++;
}
}
}
elgd.nlist[nl] = 40;
return 0;
} /* listgam */
/* ======================================================================= */
int examine_gamma(void)
{
float x1, y1, egamma, dx, deg;
int igam, iwin, nc;
char ans[80], command[80];
printf("Button 1: set gate on gamma\n"
"Button 2: add gate on gamma\n"
"Button 3: show gamma\n"
" - : subtract gamma\n"
" . : logical-type combine gamma\n"
" b/f : back/forward in gate history\n"
" ...hit X on keyboard to exit...\n");
START:
igam = 0;
retic3(&x1, &y1, ans, &iwin);
if (!strcmp(ans, "X") || !strcmp(ans, "x")) {
return 0;
} else if (!strcmp(ans, "b") || !strcmp(ans, "B") ||
!strcmp(ans, "u") || !strcmp(ans, "U")) {
undo_esclev(-1, 4);
goto START;
} else if (!strcmp(ans, "f") || !strcmp(ans, "F") ||
!strcmp(ans, "r") || !strcmp(ans, "R")) {
undo_esclev(1, 4);
goto START;
}
if (iwin == 2) {
igam = nearest_gamma(x1, y1);
if (igam < 0) {
printf(" ... no gamma selected, try again...\n");
goto START;
}
egamma = glsgd.gam[igam].e;
printf("\n*** Egamma = %7.2f +- %4.2f Igamma = %6.2f +- %4.2f\n",
glsgd.gam[igam].e, glsgd.gam[igam].de,
glsgd.gam[igam].i, glsgd.gam[igam].di);
} else if (iwin == 1) {
dx = 0.0f;
energy(x1 - 0.5f, dx, &egamma, &deg);
} else {
printf("Click in spectrum window for gates on peaks,\n"
" or in level scheme window for gates on gammas.\n");
goto START;
}
sprintf(command, " %.2f", egamma);
nc = strlen(command);
if (ans[1] == '1') {
printf("%s\n", command);
hilite(-1);
if (num_gate(command, nc)) goto ERR;
disp_dat(0);
} else if (ans[1] == '2') {
*command = '+';
printf("%s\n", command);
if (combine(command, nc)) goto ERR;
disp_dat(0);
} else if (ans[1] == '3') {
*command = 'X';
printf("%s\n", command);
if (examine(command, nc)) goto ERR;
} else if (!strcmp(ans, "minus") || !strncmp(ans, "-", 1)) {
*command = '-';
printf("%s\n", command);
if (combine(command, nc)) goto ERR;
disp_dat(0);
} else if (!strcmp(ans, "period") || !strncmp(ans, ".", 1)) {
*command = '.';
printf("%s\n", command);
if (combine(command, nc)) goto ERR;
disp_dat(0);
}
goto START;
ERR:
printf("Bad energy...\n");
goto START;
} /* examine_gamma */
/* ======================================================================= */
int comfil(char *filnam, int nc)
{
static int cfopen = 0;
strncpy(filnam, " ", 2);
if (nc < 3) {
/* ask for command file name */
GETFILNAM:
infile = 0;
nc = askfn(filnam, 80, "", ".cmd", "Command file name = ?");
if (nc == 0) return 0;
}
setext(filnam, ".cmd", 80);
if (!strncmp(filnam, "END", 3) || !strncmp(filnam, "end", 3)) {
/* close command file, lu IR = console */
if (cfopen || cflog) {
fclose(cffile);
} else {
printf("Command file not open.\n");
}
cfopen = 0;
cflog = 0;
infile = 0;
} else if (!strncmp(filnam, "CON", 3) || !strncmp(filnam, "con", 3)) {
if (cfopen) {
infile = cffile;
} else {
printf("Command file not open.\n");
}
} else if (!strncmp(filnam, "CHK", 3) || !strncmp(filnam, "chk", 3)) {
if (cfopen) {
infile = 0;
if (!askyn("Proceed? (Y/N)")) return 0;
infile = cffile;
} else if (!cflog) {
printf("Command file not open.\n");
}
} else if (!strncmp(filnam, "ERA", 3) || !strncmp(filnam, "era", 3)) {
erase();
} else if (!strncmp(filnam, "LOG", 3) || !strncmp(filnam, "log", 3)) {
if (cfopen || cflog) fclose(cffile);
while (1) {
cfopen = 0;
cflog = 0;
infile = 0;
nc = askfn(filnam, 80, "", ".cmd", "File name for command logging = ?");
if (nc == 0) return 0;
if (strncmp(filnam, "END", 3) && strncmp(filnam, "end", 3) &&
strncmp(filnam, "CON", 3) && strncmp(filnam, "con", 3) &&
strncmp(filnam, "