CXRadReader.cpp 82.9 KB
Newer Older
1 2 3 4 5
#include "CXRadReader.h"

#include <sys/stat.h>

#include "Riostream.h"
Jérémie Dudouet's avatar
Jérémie Dudouet committed
6 7
#include "TH1D.h"
#include "TH2D.h"
8 9 10 11 12 13 14 15 16 17 18 19 20 21
#include "TF1.h"

#include "CXMainWindow.h"

CXRadReader::CXRadReader()
{
    xxgd.numchs = 1;
    xxgd.many_cubes = 0;
    xxgd.cubefact[0] = xxgd.cubefact[1] = xxgd.cubefact[2] = xxgd.cubefact[3] = xxgd.cubefact[4] = 0.0;
    setcubenum(0);
}

CXRadReader::~CXRadReader()
{
Jérémie Dudouet's avatar
Jérémie Dudouet committed
22 23 24 25 26 27
    if(hLUTHist) delete hLUTHist;
    if(hTotalProj) delete hTotalProj;
    if(hBackground) delete hBackground;
    if(h2DProj) delete h2DProj;
    if(fECalibFunc) delete fECalibFunc;
    if(fEffFunc) delete fEffFunc;
28 29 30 31
}

Int_t CXRadReader::ReadCube(TString FileName)
{
Jérémie Dudouet's avatar
Jérémie Dudouet committed
32 33 34 35
    INFO_MESS << "Reading Cube file: " << FileName << ENDL;

    fName = FileName.Copy().ReplaceAll(".cub","");

36 37 38 39 40 41 42 43 44 45 46 47
    strcpy(xxgd.cubnam,FileName.Data());
    if (open3dfile(xxgd.cubnam, &xxgd.numchs)) {
        ERR_MESS<<xxgd.cubnam<<" not readable"<<ENDL;
        return 1;
    }
    if (xxgd.numchs > MAXCHS){
        ERR_MESS<<"Cube file has too many channels (" << xxgd.numchs << "), MAXCHS is "<< MAXCHS <<"  ==> EXIT"<<ENDL;
        return 1;
    }

    INFO_MESS<<xxgd.cubnam<<" loaded (" << xxgd.numchs << " channels )"<<ENDL;

Jérémie Dudouet's avatar
Jérémie Dudouet committed
48
    fCubeFileName = FileName;
49 50 51 52 53 54
    return 0;
}

Int_t CXRadReader::ReadTabFile(TString FileName)
{
    if(!gSystem->IsFileInIncludePath(FileName)){
Jérémie Dudouet's avatar
Jérémie Dudouet committed
55
        WARN_MESS<<"No LUT file detected, default values assumed (1keV/bin)"<<ENDL;
56 57 58 59 60 61
        xxgd.nclook  = xxgd.numchs;
        xxgd.lookmin = 1;
        xxgd.lookmax = xxgd.numchs;
        for(int i=0 ; i<16384 ; i++) xxgd.looktab[i] = i+1;
    }
    else {
Jérémie Dudouet's avatar
Jérémie Dudouet committed
62
        INFO_MESS << "Reading ADC-to-Cube lookup table file: " << FileName << ENDL;
63

Jérémie Dudouet's avatar
Jérémie Dudouet committed
64
        char  filnam[800];
65 66 67 68 69 70 71 72 73 74 75 76 77 78
        strcpy(filnam,FileName.Data());

        if (read_tab_file(filnam, &xxgd.nclook, &xxgd.lookmin, &xxgd.lookmax, xxgd.looktab, 16384)) {
            ERR_MESS<<FileName<<" not readable"<<ENDL;
            return 1;
        }
        if (xxgd.lookmax != xxgd.numchs) {
            ERR_MESS<<"No. of values in lookup table (" << xxgd.lookmax << ") is not coherent with cube dimension (" << xxgd.numchs << ")" <<ENDL;
            return 1;
        }
    }

    INFO_MESS<<Form("NCLook, LookMin, LookMax = %d, %d, %d",xxgd.nclook, xxgd.lookmin, xxgd.lookmax) << ENDL;

Jérémie Dudouet's avatar
Jérémie Dudouet committed
79
    fLUTFileName = FileName;
80 81 82 83 84 85 86
    return 0;
}

Int_t CXRadReader::ReadTotProj(TString FileName)
{
    // xxgd.bspec[0] ==> Total proj

Jérémie Dudouet's avatar
Jérémie Dudouet committed
87 88
    if(!gSystem->IsFileInIncludePath(FileName)) {
        WARN_MESS<<"No Total projection file detected, a default one will be built from the Cube"<<ENDL;
89

Jérémie Dudouet's avatar
Jérémie Dudouet committed
90 91 92 93 94 95 96
        if(AutoBuildProj(fCubeFileName,1) == 0) {
            ReadTotProj(fCubeFileName.Copy().ReplaceAll(".cub",".spe"));
        }
        else {
            ERR_MESS << "Error occured while building the total projection" << ENDL;
            return 1;
        }
97
    }
Jérémie Dudouet's avatar
Jérémie Dudouet committed
98 99 100
    else {

        INFO_MESS << "Reading Total projection file: " << FileName << ENDL;
101

Jérémie Dudouet's avatar
Jérémie Dudouet committed
102 103 104
        int   nch;
        char  filnam[800], namesp[16];
        xxgd.le2pro2d = 0;
105

Jérémie Dudouet's avatar
Jérémie Dudouet committed
106 107 108 109 110 111 112 113 114 115 116 117 118 119
        strcpy(filnam,FileName.Data());
        if (read_spe_file(filnam, &xxgd.bspec[0][0], namesp, &nch, MAXCHS)) {
            ERR_MESS<<filnam<<" not readable ==> EXIT"<<ENDL;
            return 1;
        }
        if (nch != xxgd.numchs){
            ERR_MESS<< filnam << " spectrum lenght (" << nch  << ") has wrong length ( must be " << xxgd.numchs << ") ==> EXIT"<<ENDL;
            return 1;
        }

        INFO_MESS<< "Total projection file loaded: " << filnam <<ENDL;
    }

    fTotalProjFileName = FileName;
120 121 122
    return 0;
}

Jérémie Dudouet's avatar
Jérémie Dudouet committed
123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
TH1 *CXRadReader::GetHistFromSpeFile(TString filename)
{
    if(!gSystem->IsFileInIncludePath(filename)) {
        WARN_MESS<<"File "<< filename << " not found "<<ENDL;
        return nullptr;
    }
    else {
        INFO_MESS << "Exporting radware spectrum: " << filename << " in ROOT histogram" << ENDL;

        TObjArray *arr = filename.Tokenize("/");
        TString name = arr->Last()->GetName();
        delete arr;

        float spec[MAXCHS];
        char  filnam[800],namesp[800];
        int   nch;

        strcpy(filnam,filename.Data());
        if (read_spe_file(filnam, &spec[0], namesp, &nch, MAXCHS)) {
            ERR_MESS<<filnam<<" not readable ==> EXIT"<<ENDL;
            return nullptr;
        }

        TH1F *h = new TH1F(name,name,nch,0,nch);
        for(int i=0 ; i<nch ; i++) {
            h->SetBinContent(i+1,spec[i]);
        }

        return h;
    }

    return nullptr;
}

void CXRadReader::SaveBackground(TString filename)
{
    char  filnam[800];
    strcpy(filnam,filename.Data());

    wspec(filnam, xxgd.background, xxgd.numchs);
}

165 166
Int_t CXRadReader::Read2DProj(TString FileName)
{
Jérémie Dudouet's avatar
Jérémie Dudouet committed
167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205
    if(!gSystem->IsFileInIncludePath(FileName)) {
        WARN_MESS<<"No 2D projection file detected, a default one will be built from the Cube"<<ENDL;

        if(AutoBuildProj(fCubeFileName,2) == 0) {
            Read2DProj(fCubeFileName.Copy().ReplaceAll(".cub",".2dp"));
        }
        else {
            ERR_MESS << "Error occured while building the 2D projection" << ENDL;
            return 1;
        }
    }
    else {

        int   reclen, ich;
        char  filnam[800];

        strcpy(filnam,FileName.Data());
        if (!inq_file(filnam, &reclen)) {
            ERR_MESS<<Form("File %s does not exist",filnam) <<ENDL;
            return 1;
        }
        FILE *file = fopen(filnam, "r");

        xxgd.luch[0] = 0;
        for (int iy = 1; iy < xxgd.numchs; ++iy) {
            xxgd.luch[iy] = xxgd.luch[iy-1] + xxgd.numchs - iy;
        }
        xxgd.matchs = xxgd.luch[xxgd.numchs - 1] + xxgd.numchs;

        reclen *= 4;
        for (int iy = 0; iy < xxgd.numchs; ++iy) {
            fseek(file, iy*reclen, SEEK_SET);
            fread(xxgd.spec[0], reclen, 1, file);
            for (int ix = 0; ix <= iy; ++ix) {
                ich = xxgd.luch[ix] + iy;
                xxgd.pro2d[ich] = xxgd.spec[0][ix];
            }
        }
        fclose(file);
206

Jérémie Dudouet's avatar
Jérémie Dudouet committed
207 208 209 210 211 212 213 214 215 216 217 218 219
        INFO_MESS<<"2D projection file loaded: "<<FileName<<ENDL;

        f2DProjFileName = FileName;

    }

    return 0;
}

Int_t CXRadReader::ReadGG(TH2 *h)
{
    if (h == nullptr) {
        ERR_MESS<<"input is null"<<ENDL;
220 221
        return 1;
    }
Jérémie Dudouet's avatar
Jérémie Dudouet committed
222 223 224 225 226 227 228 229

    xxgd.numchs = h->GetNbinsX();

    if(h->GetNbinsX() != h->GetNbinsY()) {
        ERR_MESS << " NBinsX(" << h->GetNbinsX() << ") and NBinsY(" << h->GetNbinsY() << ") are different ! IGNORED" << ENDL;
    }

    //Init
230 231 232 233 234 235 236 237

    xxgd.luch[0] = 0;
    for (int iy = 1; iy < xxgd.numchs; ++iy) {
        xxgd.luch[iy] = xxgd.luch[iy-1] + xxgd.numchs - iy;
    }
    xxgd.matchs = xxgd.luch[xxgd.numchs - 1] + xxgd.numchs;

    for (int iy = 0; iy < xxgd.numchs; ++iy) {
Jérémie Dudouet's avatar
Jérémie Dudouet committed
238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259
        for (int ii = 0; ii < 6; ii++) {
            xxgd.bspec[ii][iy] = 0;
        }
    }

    // LUT

    for (int i= 0; i < xxgd.numchs; ++i) {
        xxgd.elo_sp[i] = h->GetXaxis()->GetBinLowEdge(i+1);
        xxgd.ehi_sp[i] = h->GetXaxis()->GetBinUpEdge(i+1);

        double_t eg = (xxgd.elo_sp[i] + xxgd.ehi_sp[i]) / 2.0f;
        xxgd.eff_sp[i] = 1;
        xxgd.energy_sp[i] = eg;
        xxgd.ewid_sp[i] = xxgd.ehi_sp[i] - xxgd.elo_sp[i];
    }

    // 2D Proj
    memset(xxgd.pro2d,0,sizeof(xxgd.pro2d));

    // Fill 2D histos
    for (int iy = 0; iy < xxgd.numchs; ++iy) {
260
        for (int ix = 0; ix <= iy; ++ix) {
Jérémie Dudouet's avatar
Jérémie Dudouet committed
261 262
            int ich = xxgd.luch[ix] + iy;
            xxgd.pro2d[ich] = h->GetBinContent(ix+1,iy+1) * ((xxgd.ewid_sp[iy]*xxgd.ewid_sp[ix]) / (xxgd.ewid_sp[0]*xxgd.ewid_sp[0]));
263 264 265
        }
    }

Jérémie Dudouet's avatar
Jérémie Dudouet committed
266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318
    // Total Proj

    if(hTotalProj) delete hTotalProj;
    hTotalProj = h->ProjectionX();
    hTotalProj->Reset();
    hTotalProj->SetNameTitle("TotalProj","TotalProj");

    for (int y=0; y<xxgd.numchs; y++) {
        for (int x=0; x<=y; x++) {
            int ich = xxgd.luch[x] + y;
            xxgd.bspec[0][y] += xxgd.pro2d[ich];
            xxgd.bspec[0][x] += xxgd.pro2d[ich];
        }
    }

    for (int ich = 0; ich < xxgd.numchs; ++ich) {
        hTotalProj->SetBinContent(ich+1,xxgd.bspec[0][ich] / (xxgd.ewid_sp[ich]/xxgd.ewid_sp[0]));
    }

    // Background

    if(hBackground) delete hBackground;
    hBackground = (TH1*) hTotalProj->Clone();
    hBackground->Reset();
    hBackground->SetNameTitle("Background","Background");

    autobkgnd();

    for (int i = 0; i < xxgd.numchs; ++i){
        hBackground->SetBinContent(i+1,xxgd.background[i] / (xxgd.ewid_sp[i]/xxgd.ewid_sp[0]));
    }

    // Init gates spectra

    fGatedSpectrum = (TH1*) hTotalProj->Clone();
    fGatedSpectrum->Reset();
    fGatedSpectrum->SetNameTitle("GatedSpectrum","GatedSpectrum");
    fGatedSpectrum->SetXTitle("Energy (kev)");
    fGatedSpectrum->SetYTitle(Form("Counts (%g kev/bin)",fGatedSpectrum->GetBinWidth(1)));
    fGatedSpectrum->GetXaxis()->SetTitleOffset(1.2);
    fGatedSpectrum->GetYaxis()->SetTitleOffset(0.8);
    fGatedSpectrum->SetStats();

    fGatedSpectrumNoBG = (TH1*) hTotalProj->Clone();
    fGatedSpectrumNoBG->Reset();
    fGatedSpectrumNoBG->SetNameTitle("GatedSpectrumBG","GatedSpectrumBG");
    fGatedSpectrumNoBG->SetXTitle("Energy (kev)");
    fGatedSpectrumNoBG->SetYTitle(Form("Counts (%g kev/bin)",fGatedSpectrum->GetBinWidth(1)));
    fGatedSpectrum->GetXaxis()->SetTitleOffset(1.2);
    fGatedSpectrum->GetYaxis()->SetTitleOffset(0.8);
    fGatedSpectrum->SetStats();

    INFO_MESS<<"2D matrix file loaded: "<<h->GetName()<<" "<<h->GetTitle()<<ENDL;
319 320 321 322

    return 0;
}

Jérémie Dudouet's avatar
Jérémie Dudouet committed
323
void CXRadReader::ReEvalBackground(Float_t Xmin, Float_t XMax, Float_t LenghtFactor, Float_t Width, Float_t FWHM_0, Float_t FWHM_1, Float_t FWHM_2)
324
{
Jérémie Dudouet's avatar
Jérémie Dudouet committed
325 326 327 328
    if(hBackground == nullptr) {
        ERR_MESS<<"No existing background"<<ENDL;
        return;
    }
329

Jérémie Dudouet's avatar
Jérémie Dudouet committed
330 331 332 333 334
    fBGLenghtFactor = LenghtFactor;
    fDefaultWidth = Width;
    fFWHMPars[0] = FWHM_0;
    fFWHMPars[1] = FWHM_1;
    fFWHMPars[2] = FWHM_2;
335

Jérémie Dudouet's avatar
Jérémie Dudouet committed
336
    hBackground->Reset();
337

Jérémie Dudouet's avatar
Jérémie Dudouet committed
338 339 340 341
    autobkgnd(Xmin, XMax);

    for (int i = 0; i < xxgd.numchs; ++i){
        hBackground->SetBinContent(i+1,xxgd.background[i] / (xxgd.ewid_sp[i]/xxgd.ewid_sp[0]));
342
    }
Jérémie Dudouet's avatar
Jérémie Dudouet committed
343
}
344

Jérémie Dudouet's avatar
Jérémie Dudouet committed
345 346 347 348 349 350 351 352

Int_t CXRadReader::ReadBackground(TString FileName)
{
    // xxgd.bspec[1] ==> Background

    if(!gSystem->IsFileInIncludePath(FileName)) {
        WARN_MESS<<"No Background file detected, a default one will be built from the total projection"<<ENDL;
        autobkgnd();
353
    }
Jérémie Dudouet's avatar
Jérémie Dudouet committed
354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393
    else {

        INFO_MESS << "Reading Total projection file: " << FileName << ENDL;

        double fsum1;
        double bspec2[4][MAXCHS];
        int   nch;
        char  filnam[800], namesp[16];
        xxgd.le2pro2d = 0;

        strcpy(filnam,FileName.Data());
        if (read_spe_file(filnam, &xxgd.bspec[2][0], namesp, &nch, MAXCHS)) {
            ERR_MESS<<filnam<<" not readable ==> EXIT"<<ENDL;
            return 1;
        }
        if (nch != xxgd.numchs){
            ERR_MESS<< filnam << " spectrum lenght (" << nch  << ") has wrong length ( must be " << xxgd.numchs << ") ==> EXIT"<<ENDL;
            return 1;
        }

        fsum1 = 0.0f;
        for (int i = 0; i < xxgd.numchs; ++i) {
            fsum1 += xxgd.bspec[0][i];
            bspec2[0][i] = xxgd.bspec[0][i] - xxgd.bspec[2][i];
        }
        for (int j = 0; j < xxgd.numchs; ++j) {
            xxgd.background[j] = xxgd.bspec[2][j];
            bspec2[1][j] = bspec2[2][j] = bspec2[3][j] = 0.0f;
            xxgd.bspec[2][j] /= fsum1;
        }

        for (int i = 0; i < xxgd.numchs; ++i) {
            xxgd.bspec[1][i] = xxgd.bspec[2][i];
            xxgd.bspec[2][i] = bspec2[0][i];
            xxgd.bspec[3][i] = bspec2[1][i];
            xxgd.bspec[4][i] = bspec2[2][i];
            xxgd.bspec[5][i] = bspec2[3][i];
        }

        INFO_MESS<< "Background file loaded: " << filnam <<ENDL;
394

Jérémie Dudouet's avatar
Jérémie Dudouet committed
395 396
        fBackgroundFileName =FileName;
    }
397 398 399 400

    return 0;
}

Jérémie Dudouet's avatar
Jérémie Dudouet committed
401
Int_t CXRadReader::ReadCalibs(TString EnergyFileName, TString EfficiencyFileName, Float_t CompFactor)
402 403
{
    double f1, f2, x1, x2, x3;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
404
    double  x, eg, eff;
405
    int    i, j, jj, iorder, nterms;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
406
    char   title[800], filnam[800];
407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462

    /* get energy calibration */

    bool readcal = false;
    if(gSystem->IsFileInIncludePath(EnergyFileName)){
        INFO_MESS<<"Calibration file: " << EnergyFileName  << " found"<<ENDL;
        readcal = true;
    }
    else
        INFO_MESS<<"No Calibration file set or found, no energy calibration applied"<<ENDL;

    if(readcal){
        strcpy(filnam,EnergyFileName.Data());
        if (read_cal_file(filnam, title, &iorder, fGains)) {
            ERR_MESS<<EnergyFileName<<" not readable"<<ENDL;
            return 1;
        }
    }
    else {
        iorder = 1;
        nterms = iorder + 1;
        memset(fGains,0,sizeof(fGains));
        fGains[1] = 1.;
    }

    fCompFactor = CompFactor;
    INFO_MESS<<"Contraction factor: "<< fCompFactor <<ENDL;

    /* get efficiency calibration */

    readcal = false;
    if(gSystem->IsFileInIncludePath(EfficiencyFileName)){
        INFO_MESS<<"Efficiency file: " << EfficiencyFileName  << " found"<<ENDL;
        readcal = true;
    }
    else
        INFO_MESS<<"No Efficiency file set or found, no efficiency value applied"<<ENDL;

    if(readcal){
        strcpy(filnam,EfficiencyFileName.Data());
        if (read_eff_file(filnam, title, feff_pars)) {
            ERR_MESS<<EfficiencyFileName<<" not readable ==> EXIT"<<ENDL;
            return 1;
        }
    }
    else {
        memset(feff_pars,0,sizeof(feff_pars));
        feff_pars[7] = 100;
        feff_pars[8] = 1000;
    }

    /* calculate energy and efficiency spectra */
    jj = 0;

    for (i = 0; i < xxgd.nclook; ++i) {
        if (xxgd.looktab[i] != jj) {
Jérémie Dudouet's avatar
Jérémie Dudouet committed
463 464
//            x = ((double) i - 0.5f);
            x = ((double) i);
465 466 467 468
            eg = fGains[nterms - 1];
            for (j = nterms - 1; j >= 1; --j) {
                eg = fGains[j - 1] + eg * x;
            }
Jérémie Dudouet's avatar
Jérémie Dudouet committed
469
            if (jj != 0) xxgd.ehi_sp[jj - 1] = eg / fCompFactor;
470
            jj = xxgd.looktab[i];
Jérémie Dudouet's avatar
Jérémie Dudouet committed
471
            if (jj != 0) xxgd.elo_sp[jj - 1] = eg / fCompFactor;
472 473
        }
    }
Jérémie Dudouet's avatar
Jérémie Dudouet committed
474

475
    if (xxgd.looktab[xxgd.nclook - 1] != 0) {
Jérémie Dudouet's avatar
Jérémie Dudouet committed
476
        x = ((double) xxgd.nclook - 1.0f);
477 478 479 480
        eg = fGains[nterms - 1];
        for (j = nterms - 1; j >= 1; --j) {
            eg = fGains[j - 1] + eg * x;
        }
Jérémie Dudouet's avatar
Jérémie Dudouet committed
481
        xxgd.ehi_sp[xxgd.looktab[xxgd.nclook - 1] - 1] = eg / fCompFactor;
482
    }
Jérémie Dudouet's avatar
Jérémie Dudouet committed
483

484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507
    for (i = 0; i < xxgd.numchs; ++i) {
        eg = (xxgd.elo_sp[i] + xxgd.ehi_sp[i]) / 2.0f;
        x1 = log(eg / feff_pars[7]);
        x2 = log(eg / feff_pars[8]);
        f1 = feff_pars[0] + feff_pars[1] * x1 + feff_pars[2] * x1 * x1;
        f2 = feff_pars[3] + feff_pars[4] * x2 + feff_pars[5] * x2 * x2;
        if (f1 <= 0. || f2 <= 0.) {
            eff = 1.0f;
        } else {
            x3 = exp(-feff_pars[6] * log(f1)) + exp(-feff_pars[6] * log(f2));
            if (x3 <= 0.) {
                eff = 1.0f;
            } else {
                eff = exp(exp(-log(x3) / feff_pars[6]));
            }
        }
        xxgd.eff_sp[i] = eff;
        xxgd.energy_sp[i] = eg;

        xxgd.ewid_sp[i] = xxgd.ehi_sp[i] - xxgd.elo_sp[i];
    }

    // define edges for non constant binning
    for(int i=0 ; i<xxgd.numchs ; i++) {
Jérémie Dudouet's avatar
Jérémie Dudouet committed
508
        edges[i] = xxgd.elo_sp[i];
509
    }
Jérémie Dudouet's avatar
Jérémie Dudouet committed
510
    edges[xxgd.numchs] = xxgd.ehi_sp[xxgd.numchs-1];
511

Jérémie Dudouet's avatar
Jérémie Dudouet committed
512 513
    fECalFileName = EnergyFileName;
    fEffFileName = EfficiencyFileName;
514 515 516 517 518
    return 0;
}

void CXRadReader::BuildHistos()
{
Jérémie Dudouet's avatar
Jérémie Dudouet committed
519 520
    gErrorIgnoreLevel = kError;

521 522
    // LUT
    if(hLUTHist) delete hLUTHist;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
523
    hLUTHist = new TH1D("LUT","LUT",xxgd.numchs,edges);
524 525 526 527 528
    hLUTHist->SetXTitle("Energy (kev)");
    hLUTHist->SetYTitle("BinWidth (kev)");

    // Total Proj
    if(hTotalProj) delete hTotalProj;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
529
    hTotalProj = new TH1D("TotalProj","TotalProj",xxgd.numchs,edges);
530
    hTotalProj->SetXTitle("Energy (kev)");
Jérémie Dudouet's avatar
Jérémie Dudouet committed
531
    hTotalProj->SetYTitle(Form("Counts (%g kev/bin)",hTotalProj->GetBinWidth(1)));
532 533 534

    // Background
    if(hBackground) delete hBackground;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
535
    hBackground = new TH1D("Background","Background",xxgd.numchs,edges);
536
    hBackground->SetXTitle("Energy (kev)");
Jérémie Dudouet's avatar
Jérémie Dudouet committed
537
    hBackground->SetYTitle(Form("Counts (%g kev/bin)",hBackground->GetBinWidth(1)));
538 539 540 541

    // 2DProj

    if(h2DProj) delete h2DProj;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
542
    h2DProj = new TH2D("2DProj","2DProj",xxgd.numchs,edges,xxgd.numchs,edges);
543 544 545 546 547 548 549 550 551 552 553 554 555 556
    h2DProj->SetXTitle("E#gamma_{1}");
    h2DProj->SetYTitle("E#gamma_{2}");

    // Calib functions

    if(fECalibFunc) delete fECalibFunc;
    fECalibFunc = new TF1("ECalib","pol5",0,xxgd.nclook);
    fECalibFunc->SetParameters(fGains);
    fECalibFunc->GetXaxis()->SetTitle("Energy (kev)");
    fECalibFunc->SetNpx(1000);

    if(fEffFunc) delete fEffFunc;
    fEffFunc = new TF1("EfficiencyFunc",this,&CXRadReader::EffFuncFormula,0,xxgd.nclook,10,"CXRadReader","EffFuncFormula");
    fEffFunc->SetParameters(feff_pars);
Jérémie Dudouet's avatar
Jérémie Dudouet committed
557
    //    fEffFunc->SetParameters(7.04,0.7,0.,5.273,-0.863,0.01,11,100.,1000.); //default from radware's web site
558 559 560 561 562 563

    fEffFunc->GetXaxis()->SetTitle("Energy (kev)");
    fEffFunc->SetNpx(1000);

    // Fill 1D histos
    for (int i = 0; i < xxgd.numchs; ++i){
Jérémie Dudouet's avatar
Jérémie Dudouet committed
564 565 566
        hLUTHist->SetBinContent(i+1,xxgd.ewid_sp[i]);
        hTotalProj->SetBinContent(i+1,xxgd.bspec[0][i] / (xxgd.ewid_sp[i]/xxgd.ewid_sp[0]));
        hBackground->SetBinContent(i+1,xxgd.background[i] / (xxgd.ewid_sp[i]/xxgd.ewid_sp[0]));
567 568 569 570 571 572
    }

    // Fill 2D histos
    for (int iy = 0; iy < xxgd.numchs; ++iy) {
        for (int ix = 0; ix <= iy; ++ix) {
            int ich = xxgd.luch[ix] + iy;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
573 574
            h2DProj->SetBinContent(ix+1,iy+1,xxgd.pro2d[ich] / ((xxgd.ewid_sp[iy]*xxgd.ewid_sp[ix])/(xxgd.ewid_sp[0]*xxgd.ewid_sp[0])));
            h2DProj->SetBinContent(iy+1,ix+1,xxgd.pro2d[ich] / ((xxgd.ewid_sp[iy]*xxgd.ewid_sp[ix])/(xxgd.ewid_sp[0]*xxgd.ewid_sp[0])));
575 576
        }
    }
Jérémie Dudouet's avatar
Jérémie Dudouet committed
577 578 579 580 581 582 583 584 585 586 587 588

    fGatedSpectrum = new TH1D("GatedSpectrum","GatedSpectrum",xxgd.numchs,edges);
    fGatedSpectrum->SetXTitle("Energy (kev)");
    fGatedSpectrum->SetYTitle(Form("Counts (%g kev/bin)",fGatedSpectrum->GetBinWidth(1)));
    fGatedSpectrum->SetStats();

    fGatedSpectrumNoBG = new TH1D("GatedSpectrumBG","GatedSpectrumBG",xxgd.numchs,edges);
    fGatedSpectrumNoBG->SetXTitle("Energy (kev)");
    fGatedSpectrumNoBG->SetYTitle(Form("Counts (%g kev/bin)",fGatedSpectrum->GetBinWidth(1)));
    fGatedSpectrumNoBG->SetStats();

    gErrorIgnoreLevel = kPrint;
589 590 591 592 593 594 595
}

double CXRadReader::EffFuncFormula(double *x, double *p)
{
    double EG = x[0];
    double eff=0;

596
    // Paramters from Radware's web site: https://radware.phy.ornl.gov/gf3/#5.3.
597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643
    // A=7.04 B=0.7 C=0.
    // D=5.273 E=-0.863 F=0.01
    // G=11

    //low energy part
    // Experience shows that, unless there is a good reason to do otherwise, the parameter C should usually be left fixed to zero.
    // In addition, many gamma-ray sources do not provide enough data points at low energy to unambiguously determine both B and G,
    // so that at least one of these parameters may also have to be fixed before beginning the fit.
    // If more data that better defines the turnover region is added later, B and/or G may then be freed.
    // Typical values for B and G, for coaxial Ge detectors, are of the order 1 and 20, respectively.

    double A=p[0];
    double B=p[1];
    double C=p[2];

    //high energy
    double D=p[3];
    double E=p[4];
    double F=p[5];

    // interaction parameter between the two regions
    // the larger G is, the sharper will be the turnover at the top, between the two curves.
    // If the efficiency turns over gently, G will be small.
    double G=p[6];


    double E1=p[7]; //100 keV
    double E2=p[8]; //1000 keV

    double x1 = log(EG / E1);
    double x2 = log(EG / E2);

    double f1 = A + B*x1 + C*x1*x1;
    double f2 = D + E*x2 + F*x2*x2;

    if (f1 <= 0. || f2 <= 0.)
        eff = 1.0f;
    else {
        double x3 = exp(-G * log(f1)) + exp(-G * log(f2));

        if (x3 <= 0.) eff = 1.0f;
        else eff = exp(exp(-log(x3) / G));
    }

    return eff;
}

Jérémie Dudouet's avatar
Jérémie Dudouet committed
644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003
TH1 *CXRadReader::Project(vector< pair<double, double> > gatesX, vector< pair<double, double> > gatesY, Bool_t BGSubtract)
{
    Int_t Status;

    if(gatesX.size()==0 && gatesY.size()!=0)
        Status = Project2D(gatesY);
    else if(gatesY.size()==0 && gatesX.size()!=0)
        Status = Project2D(gatesX);
    else if(gatesX.size()!=0 && gatesY.size()!=0)
        Status = Project3D(gatesX,gatesY);
    else {
        WARN_MESS<<"Empty gates, projection ignored"<<ENDL;
        return nullptr;
    }
    if(Status>0){
        ERR_MESS<<"Error in the Projection, ignored"<<ENDL;
        return nullptr;
    }

    if(BGSubtract)
        return fGatedSpectrum;
    else
        return fGatedSpectrumNoBG;
}

int CXRadReader::Project2D(vector< pair<double, double> > gates)
{
    /* read gate of energy elo to ehi and calculate expected spectrum */
    /* spec[0][] = background-subtracted gate */
    /* spec[1][] = calculated spectrum from gammas outside gate */
    /* spec[2][] = calculated spectrum from all gammas */
    /* spec[3][] = non-background-subtracted gate */
    /* spec[4][] = square of statistical uncertainty */
    /* spec[5][] = square of statistical plus systematic uncertainties */

    strcpy(xxgd.old_name_gat, xxgd.name_gat);

    /* copy data from 2D projection */
    xxgd.bf1 = 0.0f;
    xxgd.bf2 = 0.0f;
    xxgd.bf4 = 0.0f;
    xxgd.bf5 = 0.0f;

    for (int j = 0; j < 6; ++j) {
        for (int ix = 0; ix < xxgd.numchs; ++ix) {
            xxgd.old_spec[j][ix] = xxgd.spec[j][ix];
            xxgd.spec[j][ix] = 0.0f;
        }
    }

    TString Name = "Rad2D_";
    if(fName != "") Name = "RadCube_";

    for(int i=0 ; i<gates.size() ; i++){

        double elo = gates[i].first;
        double ehi = gates[i].second;

        Name += Form("%s%.1f(%.1f)", (i==0) ? "" : "+", (ehi+elo)*0.5,(ehi-elo)*0.5);

        add_gate(elo, ehi);
    }

    sprintf(xxgd.name_gat, "%s", Name.Data());

    /* subtract background */
    for (int ix = 0; ix < xxgd.numchs; ++ix) {
        xxgd.spec[0][ix] =  xxgd.spec[0][ix]- xxgd.bf1 * xxgd.bspec[0][ix]
                - xxgd.bf2 * xxgd.bspec[1][ix]
                - xxgd.bf4 * xxgd.bspec[3][ix]
                - xxgd.bf5 * xxgd.bspec[4][ix];
    }

    /* correct for width of channel in keV */
    for (int i = 0; i < xxgd.numchs; ++i) {
        for (int j = 0; j < 4; ++j) {
            xxgd.spec[j][i] /= xxgd.ewid_sp[i];
        }
    }

    fGatedSpectrum->Reset();
    fGatedSpectrum->SetNameTitle(Name,Name);

    fGatedSpectrumNoBG->Reset();
    fGatedSpectrumNoBG->SetNameTitle(Form("%s_NoBG",Name.Data()),Form("%s_NoBG",Name.Data()));

    for (int i = 0; i < xxgd.numchs; ++i) {
        fGatedSpectrum->SetBinContent(i+1,xxgd.spec[0][i] * xxgd.ewid_sp[0]);
        fGatedSpectrumNoBG->SetBinContent(i+1,xxgd.spec[3][i] * xxgd.ewid_sp[0]);
    }

    return 0;
}


int CXRadReader::Project3D(vector< pair<double, double> > gatesX, vector< pair<double, double> > gatesY)
{
    strcpy(xxgd.old_name_gat, xxgd.name_gat);

    /* copy data from 2D projection */
    xxgd.bf1 = 0.0f;  xxgd.bf2 = 0.0f;  xxgd.bf4 = 0.0f;  xxgd.bf5 = 0.0f;

    for (int j = 0; j < 6; ++j) {
        for (int ix = 0; ix < xxgd.numchs; ++ix) {
            xxgd.old_spec[j][ix] = xxgd.spec[j][ix];
            xxgd.spec[j][ix] = 0.0f;
        }
    }

    int nsum = 0;

    TString NameX="";
    TString NameY="";

    for(int i=0 ; i<gatesX.size() ; i++){

        double eloX = gatesX[i].first;
        double ehiX = gatesX[i].second;

        NameX += Form("%s%.1f(%.1f)", (i==0) ? "" : "+",(eloX+ehiX)*0.5,TMath::Abs(eloX-ehiX)*0.5);

        for(int j=0 ; j<gatesY.size() ; j++) {

            double eloY = gatesY[j].first;
            double ehiY = gatesY[j].second;

            NameY += Form("%s%.1f(%.1f)", (j==0) ? "" : "+",(eloY+ehiY)*0.5,TMath::Abs(eloY-ehiY)*0.5);

            ++nsum;
            add_trip_gate(eloX, ehiX, eloY, ehiY);
        }
    }

    sprintf(xxgd.name_gat, "%s/%s", NameX.Data(), NameY.Data());

    if (nsum > 1)
        cout<<Form("Sum of %d double-gates...", nsum)<<endl;

    /* subtract background */
    for (int ix = 0; ix < xxgd.numchs; ++ix) {
        xxgd.spec[0][ix] =  xxgd.spec[0][ix]- xxgd.bf1 * xxgd.bspec[0][ix]
                - xxgd.bf2 * xxgd.bspec[1][ix]
                - xxgd.bf4 * xxgd.bspec[3][ix]
                - xxgd.bf5 * xxgd.bspec[4][ix];
    }

    /* correct for width of channel in keV */
    for (int i = 0; i < xxgd.numchs; ++i) {
        for (int j = 0; j < 4; ++j) xxgd.spec[j][i] /= xxgd.ewid_sp[i];
        if (!xxgd.many_cubes) xxgd.spec[4][i] = xxgd.spec[3][i] / xxgd.ewid_sp[i];
        xxgd.spec[5][i] = xxgd.spec[4][i];
        if (xxgd.spec[5][i] < 1.0f) xxgd.spec[5][i] = 1.0f;
    }

    TString Name = Form("RadCube_%s/%s", NameX.Data(), NameY.Data());

    fGatedSpectrum->Reset();
    fGatedSpectrum->SetNameTitle(Name,Name);

    fGatedSpectrumNoBG->Reset();
    fGatedSpectrumNoBG->SetNameTitle(Form("%s_NoBG",Name.Data()),Form("%s_NoBG",Name.Data()));

    for (int i = 0; i < xxgd.numchs; ++i) {
        fGatedSpectrum->SetBinContent(i+1,xxgd.spec[0][i] * xxgd.ewid_sp[0]);
        fGatedSpectrumNoBG->SetBinContent(i+1,xxgd.spec[3][i] * xxgd.ewid_sp[0]);
    }

    return 0;
}

Int_t CXRadReader::AutoBuildProj(TString FileName, Int_t Mode)
{
    int length = gd3d->length;
    char outnam[256];
    unsigned int mc[512], *data2d;
    unsigned short mc2[512];
    int i, j, k, x, y, z, ix, iy, iz;
    int ax[512], ay[512], az[512];
    float data1d[RW_MAXCH];

    unsigned char lobyte[RW_MAXCH+6];
    unsigned short overflows[RW_MAXCH][2], numoverflows;
    int overflowaddr;

    /* malloc the 2D projection space */
    if (!(data2d=(uint *)malloc(4*(length+8)*length))) {
        ERR_MESS << "Ooops, data2d malloc failed... please free up some memory" << ENDL;
        return 1;
    }

    /* calculate the mappings from 3d channel space to 2d channel space */
    i = 0;
    for (z=0; z<8; z++) {
        for (y=0; y<8; y++) {
            for (x=0; x<8; x++) {
                ax[i] = x;
                ay[i] = y;
                az[i] = z;
                i++;
            }
        }
    }

    memset(data2d, 0, 4*(length+8)*length);

    if(Mode==1)
        printf("Projecting cube to 1D...\n");
    else if(Mode==2)
        printf("Projecting cube to 2D...\n");
    else
        return 1;

    if (gd3d->cubetype==0) {  /* compressed 1/6 cube */
        for (iz=0; iz<length; iz+=8) {
            for (iy=0; iy<=iz; iy+=8) {
                for (ix=0; ix<=iy; ix+=8) {
                    j = LUMC3D(ix,iy,iz);
                    getmc3d(j,mc);
                    getmc3d(j+1,&mc[256]);
                    if ((j&511)==0) {
                        cout << Form("\rx y z: %4d %4d %4d",ix,iy,iz);
                    }
                    for (i=0; i<512; i++) {
                        *(data2d + (iy+ay[i])*length + ix+ax[i]) += mc[i];
                        *(data2d + (iz+az[i])*length + ix+ax[i]) += mc[i];
                        *(data2d + (iz+az[i])*length + iy+ay[i]) += mc[i];
                    }
                }
            }
        }
    }

    else if (gd3d->cubetype==1) {  /* compressed 1/2 cube */
        j = 0;
        for (iz=0; iz<length; iz+=8) {
            for (iy=0; iy<=iz; iy+=8) {
                for (ix=0; ix<length; ix+=8) {
                    getmc3d(j++,mc);
                    if ((ix+4) < length) {
                        getmc3d(j++,&mc[256]);
                    } else {
                        memset(&mc[256], 0, 1024);
                    }
                    if ((j&511)==0) {
                        cout << Form("\rx y z: %4d %4d %4d",ix,iy,iz);
                    }
                    for (i=0; i<512; i++)
                        *(data2d + (iz+ay[i])*length + iy+ax[i]) += mc[i];
                }
            }
        }
        for (iz=0; iz<length; iz++)
            *(data2d + iz*length + iz) = *(data2d + iz*length + iz)/2;
    }

    else if (gd3d->cubetype==2) {  /* uncompressed 1/6 cube */
        j = 0;
        for (iz=0; iz<length; iz+=8) {
            for (iy=0; iy<=iz; iy+=8) {
                for (ix=0; ix<=iy; ix+=8) {
                    if (!fread(mc2,1024,1,gd3d->CubeFile3d)) {
                        ERR_MESS << Form("Could not read record %d... aborting.", j) << ENDL;
                        return 1;
                    }
                    if (gd3d->swapbytes) {
                        for (i=0; i<512; i++)
                            swap2(&mc2[i]);
                    }
                    if (((j++)&255)==0) {
                        cout << Form("\rx y z: %4d %4d %4d",ix,iy,iz);
                    }
                    for (i=0; i<512; i++) {
                        *(data2d + (iy+ay[i])*length + ix+ax[i]) += mc2[i];
                        *(data2d + (iz+az[i])*length + ix+ax[i]) += mc2[i];
                        *(data2d + (iz+az[i])*length + iy+ay[i]) += mc2[i];
                    }
                }
            }
        }
    }

    else if (gd3d->cubetype==3) {  /* uncompressed 1/2 cube */
        j = 1024;
        k = 1 + ((length*(length+1)/2 * (length+6) + 1023) / 1024);
        k = k * 1024;
        for (iz=0; iz<length; iz++) {
            cout << Form("\r z: %4d ", iz);
            for (iy=0; iy<=iz; iy++) {
                fseek(gd3d->CubeFile3d, ((long)j), SEEK_SET);
                if (!fread(lobyte,length+6,1,gd3d->CubeFile3d)) {
                    ERR_MESS << "Could not read cube file... aborting" << ENDL;
                    return 1;
                }
                j += length+6;
                for (ix=0; ix<length; ix++)
                    *(data2d + iz*length + iy) += lobyte[ix];

                memcpy(&numoverflows, &lobyte[length], 2);
                if (numoverflows == 0)
                    continue;
                memcpy(&overflowaddr, &lobyte[length+2], 4);
                if (gd3d->swapbytes) {
                    swap2(&numoverflows);
                    swap4(&overflowaddr);
                }
                fseek(gd3d->CubeFile3d, (long)k+(long)(overflowaddr*4), SEEK_SET);
                if (!fread(overflows,4*numoverflows,1,gd3d->CubeFile3d)) {
                    ERR_MESS << "Could not read cube file... aborting" << ENDL;
                    return 1;
                }
                for (ix=0; ix<numoverflows; ix++) {
                    if (gd3d->swapbytes) swap2(&overflows[ix][1]);
                    *(data2d + iz*length + iy) += (int)overflows[ix][1];
                }
            }
            *(data2d + iz*length + iz) = *(data2d + iz*length + iz)/2;
        }
    }

    FILE *file;

    if(Mode == 2) {
        TString nameout = FileName.ReplaceAll(".cub",".2dp");
        strcpy(outnam,nameout.Data());

        file=fopen(outnam,"w+");

        for (y=0; y<length; y++) {
            for (x=0; x<=y; x++) {
                data1d[x] = (float)(*(data2d + y*length + x));
            }
            data1d[y] = 2.0*data1d[y];
            for (x=y+1; x<length; x++) {
                data1d[x] = (float)(*(data2d + x*length + y));
            }
            if (!(fwrite(data1d, 4*length, 1, file))) {
                ERR_MESS << "Cannot write 2d file, aborting..." << ENDL;
                return 1;
            }
        }
        fclose(file);
    }
    if(Mode == 1) {
        TString nameout = FileName.ReplaceAll(".cub",".spe");
        strcpy(outnam,nameout.Data());

        for (y=0; y<length; y++) {
            data1d[y] = 0.0;
            for (x=0; x<=y; x++) {
                data1d[x] += (float)(*(data2d + y*length + x));
                data1d[y] += (float)(*(data2d + y*length + x));
            }
        }

        wspec(outnam, data1d, length);
    }
    cout<<endl;
    return 0;
}

1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015
//******************************* RW util.c Stuffs ***********************************//

/* ======================================================================= */
int CXRadReader::read_tab_file(char *filnam, int *nclook, int *lookmin, int *lookmax, short *looktab, int dimlook)
{
    /* read lookup table from .tab file = filnam
     into array looktab of dimension dimlook
     Nclook = number of channels read
     lookmin, lookmax = min/max values in lookup table
     returns 1 for open/read error
     default file extension = .tab */

Jérémie Dudouet's avatar
Jérémie Dudouet committed
1016
    char cbuf[800];
1017 1018 1019
    int  j, rl, swap = -1;
    FILE *file;

Jérémie Dudouet's avatar
Jérémie Dudouet committed
1020
    setext(filnam, ".tab", 800);
1021 1022 1023 1024 1025 1026

    /* OPEN (ILU,FILE=FILNAM,FORM='UNFORMATTED',STATUS='OLD')
     READ (ILU) NCLOOK,LOOKMIN,LOOKMAX
     READ (ILU) (LOOKTAB(I),I=1,NCLOOK) */

    if (!(file = open_readonly(filnam))) return 1;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
1027
    rl = get_file_rec(file, cbuf, 800, 0);
1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367
    if (rl != 12  && rl != -12) {
        file_error("read lookup-table from", filnam);
        fclose(file);
        return 1;
    }
    if (rl < 0) {
        /* file has wrong byte order - swap bytes */
        cout<<"*** Swapping bytes read from file " << filnam <<endl;
        swap = 1;
        for (j = 0; j < 3; ++j) {
            swapb4(cbuf + 4*j);
        }
    }
    memcpy(nclook,  cbuf,     4);
    memcpy(lookmin, cbuf + 4, 4);
    memcpy(lookmax, cbuf + 8, 4);

    if (*nclook < 2 || *nclook > dimlook) {
        file_error("read lookup-table from", filnam);
        fclose(file);
        return 1;
    }

    rl = get_file_rec(file, looktab, 2*dimlook, swap);
    fclose(file);
    if (rl != 2*(*nclook) && rl != -2*(*nclook)) {
        file_error("read lookup-table from", filnam);
        return 1;
    }
    if (rl < 0) {
        /* file has wrong byte order - swap bytes */
        for (j = 0; j < *nclook; ++j) {
            swapb2((char *) (looktab + j));
        }
    }

    return 0;
} /* read_tab_file */


/* ====================================================================== */
int CXRadReader::setext(char *filnam, const char *cext, int filnam_len)
{
    /* set default extension of filename filnam to cext
     leading spaces are first removed from filnam
     if extension is present, it is left unchanged
     if no extension is present, cext is used
     returned value pointer to the dot of the .ext
     cext should include the dot plus a three-letter extension */

    int nc, iext;

    /* remove leading spaces from filnam */
    nc = strlen(filnam);
    if (nc > filnam_len) nc = filnam_len;
    while (nc > 0 && filnam[0] == ' ') {
        memmove(filnam, filnam+1, nc--);
        filnam[nc] = '\0';
    }
    /* remove trailing spaces from filnam */
    while (nc > 0 && filnam[nc-1] == ' ') {
        filnam[--nc] = '\0';
    }
    /* look for file extension in filnam
     if there is none, put it to cext */
    iext = 0;
    if (nc > 0) {
        for (iext = nc-1;
             (iext > 0 &&
              filnam[iext] != ']' &&
              filnam[iext] != '/' &&
              filnam[iext] != ':');
             iext--) {
            if (filnam[iext] == '.') return iext;
        }
        iext = nc;
    }
    strncpy(&filnam[iext], cext, filnam_len - iext - 1);
    return iext;
} /* setext */

/* ======================================================================= */
FILE *CXRadReader::open_readonly(char *filename)
{
    /* open old file for READONLY (if possible)
     filename  = file name
     provided for VMS compatibility */

    FILE *file;

    if (!(file = fopen(filename, "r"))) file_error("open", filename);
    return file;
} /* open_readonly */


/* ======================================================================= */
int CXRadReader::get_file_rec(FILE *fd, void *data, int maxbytes, int swap_bytes)
{
    /* read one fortran-unformatted style binary record into data */
    /* for unix systems, swap_bytes controls how get_file_rec deals with
     swapping the bytes of the record length tags at the start and end
     of the records.  Set swap_bytes to
       0 to try to automatically sense if the bytes need swapping
       1 to force the byte swap, or
      -1 to force no byte swap */
    /* returns number of bytes read in record,
     or number of bytes * -1 if bytes need swapping,
     or 0 for error */
    int  reclen, j1, j2;

    if (fread(&reclen, 4, 1, fd) != 1) return 0;
    if (reclen == 0) return 0;
    j1 = reclen;
    if ((swap_bytes == 1) ||
            (swap_bytes == 0 && reclen >= 65536)) swapb4((char *) &reclen);
    if (reclen > maxbytes) goto ERR1;
    if (fread(data, reclen, 1, fd) != 1 ||
            fread(&j2, 4, 1, fd) != 1) goto ERR2;
    /* if (j1 != j2) goto ERR2; */
    if (reclen == j1) return reclen;
    return (-reclen);

ERR1:
    ERR_MESS<<Form("ERROR: record is too big for get_file_rec\n"
                   "       max size = %d, record size = %d.\n",
                   maxbytes, reclen)<<ENDL;
    return 0;
ERR2:
    ERR_MESS<<Form("ERROR during read in get_file_rec.\n")<<ENDL;
    return 0;
} /* get_file_rec */


/* ======================================================================= */
int CXRadReader::file_error(const char *error_type, char *filename)
{
    /* write error message */
    /* cannot perform operation error_type on file filename */

    if (strlen(error_type) + strlen(filename) > 58) {
        ERR_MESS<<Form("ERROR - cannot %s file\n%s\n", error_type, filename)<<ENDL;
    } else {
        ERR_MESS<<Form("ERROR - cannot %s file %s\n", error_type, filename)<<ENDL;
    }
    return 0;
} /* file_error */

/* ======================================================================= */
void CXRadReader::swapb8(char *buf)
{
    char c;
    c = buf[7]; buf[7] = buf[0]; buf[0] = c;
    c = buf[6]; buf[6] = buf[1]; buf[1] = c;
    c = buf[5]; buf[5] = buf[2]; buf[2] = c;
    c = buf[4]; buf[4] = buf[3]; buf[3] = c;
} /* swapb8 */
void CXRadReader::swapb4(char *buf)
{
    char c;
    c = buf[3]; buf[3] = buf[0]; buf[0] = c;
    c = buf[2]; buf[2] = buf[1]; buf[1] = c;
} /* swapb4 */
void CXRadReader::swapb2(char *buf)
{
    char c;
    c = buf[1]; buf[1] = buf[0]; buf[0] = c;
} /* swapb2 */

int CXRadReader::open3dfile(char *name, int *numchs_ret)
{
    int length, i, nummc, sr;
    FHead3D head;
    Record3D rec;
    char cname[256], *tmp;

    strcpy(cname, name);
    if ((tmp = strchr(cname,' '))) *tmp = 0;

    if (!gd3d) gd3d = &gd3dn[0];

    if (!cname[0] || !(gd3d->CubeFile3d=fopen(cname,"r"))) {
        printf("\007  ERROR -- Cannot open file %s for reading.\n", cname);
        return 1;
    }

#define QUIT close3dfile(); return 1;
    if (!fread(&head,1024,1,gd3d->CubeFile3d)) {
        printf("\007  ERROR -- Cannot read file %s.\n",name);
        QUIT;
    }
    if (head.cps > 6) {
        gd3d->swapbytes = 1;
        swap4(&head.numch);
        swap4(&head.bpc);
        swap4(&head.cps);
        swap4(&head.numrecs);
    }
    else
        gd3d->swapbytes = 0;

    if (!(strncmp(head.id,"Incub8r3/Pro4d  ",16)) && head.cps == 6)
        gd3d->cubetype = 0;
    else if (!(strncmp(head.id,"Foldout         ",16)) && head.cps == 2)
        gd3d->cubetype = 1;
    else if (!(strncmp(head.id,"INCUB8R2        ",16)) && head.cps == 6)
        gd3d->cubetype = 2;
    else if (!(strncmp(head.id,"FOLDOUT         ",16)) && head.cps == 2)
        gd3d->cubetype = 3;
    else {
        printf("\007  ERROR -- Invalid header in file %s.\n",name);
        QUIT;
    }

    length = head.numch;
    *numchs_ret = length;
    gd3d->length = length;
    if (length > RW_MAXCH) {
        printf("ERROR: Number of channels in cube (%d)\n"
               "         is greater than RW_MAXCH (%d) in pro3d.h\n",
               length, RW_MAXCH);
        exit(-1);
    }
    gd3d->numrecs = head.numrecs;
    nummc = (length+7)/8;  /* minicubes per side */
    if (gd3d->cubetype==0 || gd3d->cubetype==2)
        nummc = nummc*(nummc+1)*(nummc+2)/3; /* minicubes in cube */
    else if (gd3d->cubetype==1 || gd3d->cubetype==3)
        nummc = (nummc*(nummc+1)/2) * ((length+3)/4); /* minicubes in cube */
    gd3d->nummc = nummc;

    /* calculate look-up tables to convert (x,y,z) 3d cube address
        to linearized 3d minicube address.  x<y<z
     lum{x,y,z}3d[ch] returns the number of subcubes to skip over
        to get to the subcube with ch in it.
     note that a 3d subcube on disk is 4x8x8 channels, but we always read them
        in pairs, i.e. we always read 8x8x8 channels.
  */
    for (i=0;i<8;i++) {
        lumx3d[i] = 0;
        lumy3d[i] = 0;
        lumz3d[i] = 0;
    }
    for (i=8;i<gd3d->length+16;i++) {
        lumx3d[i] = lumx3d[i-8] + 2;
        lumy3d[i] = lumy3d[i-8] + lumx3d[i];
        lumz3d[i] = lumz3d[i-8] + lumy3d[i];
    }

    printf("Axis length of cube is %d.\n", length);
    if (gd3d->cubetype>1) return 0;

    gd3d->numrecs = head.numrecs;
    printf("Cube has %d minicubes and %d records.\n", nummc, gd3d->numrecs);
    gd3d->numsr = (head.numrecs + 15)/16;

    /* read SR list from end of file */
    while (!(gd3d->mcnum = (int *)malloc(4*(1+gd3d->numsr))))
        printf("trying to malloc gd3d->mcnum...\n");
    if ((fseek(gd3d->CubeFile3d, (((long) head.numrecs)*4+1)*1024, SEEK_SET)) ||
            (!fread(gd3d->mcnum,4*(1+gd3d->numsr),1,gd3d->CubeFile3d))) {

        /* build SuperRecord physical # and minicube # lookup table */
        printf("  Building 3d SuperRecord list...\n");
        for (sr=0; sr<gd3d->numsr; sr++) {
            if ((fseek(gd3d->CubeFile3d, (((long) sr)*64+1)*1024, SEEK_SET)) ||
                    (!fread(&rec,4096,1,gd3d->CubeFile3d))) {
                printf("\007  ERROR -- Could not read Record %d... aborting.\n", sr*16);
                QUIT;
            }
            gd3d->mcnum[sr] = rec.h.minmc;
        }
        gd3d->mcnum[gd3d->numsr] = nummc;
        if (gd3d->swapbytes)
            swap4(&gd3d->mcnum[gd3d->numsr]);
        printf("  ...Done building SuperRecord list,\n"
               "       will now save it at end of cube file.\n");

        /* save SR list to end of file */
        fclose(gd3d->CubeFile3d);
        if (!(gd3d->CubeFile3d=fopen(cname,"a+"))) {
            printf("\007  ERROR -- Cannot open file %s for writing.\n", name);
            QUIT;

        }
        if (!fwrite(gd3d->mcnum,4*(1+gd3d->numsr),1,gd3d->CubeFile3d)) {
            printf("\007  ERROR -- Could not write SuperRecord list... aborting.\n");
            QUIT;
        }
        fflush(gd3d->CubeFile3d);
    }
    if (gd3d->swapbytes) {
        for (sr=0; sr<=gd3d->numsr; sr++)
            swap4(&gd3d->mcnum[sr]);
    }
#undef QUIT
    return 0;
}

int CXRadReader::close3dfile(void)
{
    if (gd3d->CubeFile3d) {
        fclose(gd3d->CubeFile3d);
        free(gd3d->mcnum);
        gd3d->CubeFile3d = 0;
    }
    return 0;
}

void CXRadReader::swap4(int *in)
{
    char *c, tmp;

    c = (char *) in;
    tmp = *c;
    *c = *(c+3);
    *(c+3) = tmp;
    tmp = *(c+1);
    *(c+1) = *(c+2);
    *(c+2) = tmp;
}

void CXRadReader::swap2(unsigned short *in)
{
    char *c, tmp;

    c = (char *) in;
    tmp = *c;
    *c = *(c+1);
    *(c+1) = tmp;
}

int CXRadReader::read_spe_file(char *filnam, float *sp, char *namesp, int *numch, int idimsp)
{
    /* read spectrum from .spe file = filnam
     into array sp of dimension idimsp
     numch = number of channels read
     namesp = name of spectrum (character*8)
     returns 1 for open/read error
     default file extension = .spe */

Jérémie Dudouet's avatar
Jérémie Dudouet committed
1368
    char cbuf[800];
1369 1370 1371
    int  idim1, idim2, j, rl, swap = -1;
    FILE *file;

Jérémie Dudouet's avatar
Jérémie Dudouet committed
1372
    setext(filnam, ".spe", 800);
1373 1374 1375 1376 1377 1378 1379

    /* read spectrum in standard gf3 format
     OPEN(ILU,FILE=FILNAM,FORM='UNFORMATTED',STATUS='OLD')
     READ(ILU) NAMESP, IDIM1, IDIM2, IRED1, IRED2
     READ(ILU,ERR=90) (SP(I),I=1,NUMCH) */

    if (!(file = open_readonly(filnam))) return 1;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
1380
    rl = get_file_rec(file, cbuf, 800, 0);
1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424
    if (rl != 24  && rl != -24) {
        file_error("read", filnam);
        fclose(file);
        return 1;
    }
    if (rl < 0) {
        /* file has wrong byte order - swap bytes */
        cout<<Form("*** Swapping bytes read from file %s", filnam)<<endl;
        swap = 1;
        for (j = 2; j < 4; ++j) {
            swapb4(cbuf + 4*j);
        }
    }
    strncpy(namesp, cbuf, 8);
    memcpy(&idim1, cbuf + 8, 4);
    memcpy(&idim2, cbuf + 12, 4);
    *numch = idim1 * idim2;

    if (*numch > idimsp) {
        *numch = idimsp;
        cout<<Form("First %d chs only taken.", idimsp)<<endl;
    }
    rl = get_file_rec(file, sp, 4*idimsp, swap);
    fclose(file);
    if (rl != 4*(*numch) && rl != -4*(*numch)) {
        file_error("read spectrum from", filnam);
        return 1;
    }
    if (rl < 0) {
        /* file has wrong byte order - swap bytes */
        for (j = 0; j < *numch; ++j) {
            swapb4((char *) (sp + j));
        }
    }
    return 0;
} /* read_spe_file */

/* ======================================================================= */
int CXRadReader::inq_file(char *filename, int *reclen)
{
    /* inquire for file existence and record length in longwords
     returns 0 for file not exists, 1 for file exists */

    int  ext;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
1425
    char jfile[800];
1426 1427 1428 1429 1430 1431
    struct stat statbuf;

    *reclen = 0;
    if (stat(filename, &statbuf)) return 0;

    ext = 0;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
1432 1433
    strncpy(jfile, filename, 800);
    ext = setext(jfile, "    ", 800);
1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453
    if (!strcmp(&jfile[ext], ".mat") ||
            !strcmp(&jfile[ext], ".MAT") ||
            !strcmp(&jfile[ext], ".esc") ||
            !strcmp(&jfile[ext], ".ESC")) {
        *reclen = 2048;
    } else if (!strcmp(&jfile[ext], ".spn") ||
               !strcmp(&jfile[ext], ".SPN") ||
               !strcmp(&jfile[ext], ".m4b") ||
               !strcmp(&jfile[ext], ".M4B") ||
               !strcmp(&jfile[ext], ".e4k") ||
               !strcmp(&jfile[ext], ".E4K")) {
        *reclen = 4096;
    } else if (!strcmp(&jfile[ext], ".cub") ||
               !strcmp(&jfile[ext], ".CUB")) {
        *reclen = 256;
    } else if (!strcmp(&jfile[ext], ".2dp") ||
               !strcmp(&jfile[ext], ".2DP")) {
        if (statbuf.st_size <= 0) {
            *reclen = 0;
        } else {
Jérémie Dudouet's avatar
Jérémie Dudouet committed
1454
            *reclen = (int) (0.5 + sqrt((double) (statbuf.st_size/4)));
1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468
        }
    }
    return 1;
} /* inq_file */


/* ======================================================================= */
int CXRadReader::read_cal_file(char *filnam, char *title, int *iorder, double *gain)
{
    /* read energy calib parameters from .aca/.cal file = filnam
     into title, iorder, gain
     returns 1 for open/read error
     default file extension = .aca, .cal */

Jérémie Dudouet's avatar
Jérémie Dudouet committed
1469
    char cbuf[800];
1470
    int  j, rl;
Jérémie Dudouet's avatar
Jérémie Dudouet committed
1471
    char fn1[800], fn2[800], line[120], *c;
1472 1473 1474
    FILE *file;
    struct stat buf;

Jérémie Dudouet's avatar
Jérémie Dudouet committed
1475 1476 1477 1478
    strncpy(fn1, filnam, 800);
    strncpy(fn2, filnam, 800);
    j = setext(fn1, ".aca", 800);
    j = setext(fn2, ".cal", 800);
1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1496 1497 1498 1499
    if (strcmp(fn1+j, ".aca") && strcmp(fn2+j, ".cal")) {
        ERR_MESS<<"ERROR - filename .ext must be .aca or .cal!"<<ENDL;
        return 1;
    } else if (strcmp(fn1+j, ".aca") && stat(fn2, &buf)) {
        ERR_MESS<<Form("ERROR - file %s does not exist!", fn2)<<ENDL;
        return 1;
    } else if (strcmp(fn2+j, ".cal") && stat(fn1, &buf)) {
        ERR_MESS<<Form("ERROR - file %s does not exist!", fn1)<<ENDL;
        return 1;
    } else if (stat(fn1, &buf) && stat(fn2, &buf)) {
        ERR_MESS<<Form("ERROR - neither file %s nor file %s exists!", fn1, fn2)<<ENDL;
        return 1;
    }

#define ERR { file_error("read", fn1); fclose(file); return 1; }
    /* first try ascii cal file (.aca) */
    if (!strcmp(fn1+j, ".aca") &&
            (file = fopen(fn1, "r"))) {
        if (!(fgets(line, 120, file)) ||
                !strstr(line, "ENCAL OUTPUT FILE") ||
                !(fgets(line, 120, file)) ||
Jérémie Dudouet's avatar
Jérémie Dudouet committed
1500
                !strncpy(title, line, 800) ||
1501 1502 1503 1504 1505 1506