From 056baf97d087ca88d23147761f5b2a6bac4286a3 Mon Sep 17 00:00:00 2001 From: tlatorre Date: Wed, 19 Jun 2019 09:55:17 -0500 Subject: add FTP, RSP, and FTK results to the output file --- src/fit.c | 32 ++++++ src/zdab_utils.c | 316 +++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/zdab_utils.h | 215 +++++++++++++++++++++++++++++++++++++ 3 files changed, 563 insertions(+) diff --git a/src/fit.c b/src/fit.c index adaaaef..4e04c0f 100644 --- a/src/fit.c +++ b/src/fit.c @@ -5872,6 +5872,9 @@ int main(int argc, char **argv) zebraBank bmast, bmc, bmcgn, mctk, b; int rv; EVBank bev; + FTPVBank bftpv; + FTXKBank bftxk; + RSPBank bftxr; MCTKBank bmctk; MCVXBank bmcvx; event ev = {0}; @@ -6192,6 +6195,35 @@ skip_mc: fprintf(fout, " dc: 0x%08x\n", get_dc_word(&ev, f, &bmast, &b)); } + if (fout) { + if (get_ftpv(f,&b,&bftpv)) { + fprintf(stderr, "%s\n", zdab_err); + } else { + fprintf(fout, " ftp:\n"); + fprintf(fout, " x: %.2f\n", bftpv.x); + fprintf(fout, " y: %.2f\n", bftpv.y); + fprintf(fout, " z: %.2f\n", bftpv.z); + } + } + + if (fout) { + if (get_ftxk(f,&b,&bftxk)) { + fprintf(stderr, "%s\n", zdab_err); + } else { + fprintf(fout, " ftk:\n"); + fprintf(fout, " energy: %.2f\n", bftxk.energy); + } + } + + if (fout) { + if (get_rsp(f,&b,&bftxr)) { + fprintf(stderr, "%s\n", zdab_err); + } else { + fprintf(fout, " rsp:\n"); + fprintf(fout, " energy: %.2f\n", bftxr.ene); + } + } + if (nhit < min_nhit) goto skip_event; if (fout) fprintf(fout, " fit:\n"); diff --git a/src/zdab_utils.c b/src/zdab_utils.c index 20c33aa..ca0e3dc 100644 --- a/src/zdab_utils.c +++ b/src/zdab_utils.c @@ -24,6 +24,8 @@ #include "zebra.h" #include "misc.h" +char zdab_err[256]; + size_t get_nhit(event *ev) { /* Returns the number of PMT hits in event `ev`. @@ -264,6 +266,320 @@ void swap_int16(int16_t *val_pt, int count) return; } +/* Get and unpack the FTPV bank. + * + * Returns 0 if successful, or -1 on error. If there was an error, the error + * string is printed to zdab_err. */ +int get_ftpv(zebraFile *f, zebraBank *ev, FTPVBank *bftpv) +{ + int rv; + zebraBank ft, ftp, ftpv; + + if (ev->links[KEV_FT-1] == 0) { + sprintf(zdab_err, "FT link is zero!"); + goto err; + } + + rv = zebra_get_bank(f,&ft,ev->links[KEV_FT-1]); + + if (rv) { + sprintf(zdab_err, "error getting FT bank: %s", zebra_err); + goto err; + } + + if (ft.links[KFT_FTP-1] == 0) { + sprintf(zdab_err, "FTP link is zero!"); + goto err; + } + + rv = zebra_get_bank(f,&ftp,ft.links[KFT_FTP-1]); + + if (rv) { + sprintf(zdab_err, "error getting FTP bank: %s", zebra_err); + goto err; + } + + if (ftp.links[KFTX_FTXV-1] == 0) { + sprintf(zdab_err, "FTXV link is zero!"); + goto err; + } + + rv = zebra_get_bank(f,&ftpv,ftp.links[KFTX_FTXV-1]); + + if (rv) { + sprintf(zdab_err, "error getting FTXV bank: %s", zebra_err); + goto err; + } + + unpack_ftpv(ftpv.data,bftpv); + + return 0; + +err: + + return -1; +} + + +/* Get and unpack the FTXK bank from the FTP fitter. + * + * Returns 0 if successful, or -1 on error. If there was an error, the error + * string is printed to zdab_err. */ +int get_ftxk(zebraFile *f, zebraBank *ev, FTXKBank *bftxk) +{ + int rv; + zebraBank ft, ftp, ftxa, ftxk; + + if (ev->links[KEV_FT-1] == 0) { + sprintf(zdab_err, "FT link is zero!"); + goto err; + } + + rv = zebra_get_bank(f,&ft,ev->links[KEV_FT-1]); + + if (rv) { + sprintf(zdab_err, "error getting FT bank: %s", zebra_err); + goto err; + } + + if (ft.links[KFT_FTP-1] == 0) { + sprintf(zdab_err, "FTP link is zero!"); + goto err; + } + + rv = zebra_get_bank(f,&ftp,ft.links[KFT_FTP-1]); + + if (rv) { + sprintf(zdab_err, "error getting FTP bank: %s", zebra_err); + goto err; + } + + if (ftp.links[KFTX_FTXA-1] == 0) { + sprintf(zdab_err, "FTXA link is zero!"); + goto err; + } + + rv = zebra_get_bank(f,&ftxa,ftp.links[KFTX_FTXA-1]); + + if (rv) { + sprintf(zdab_err, "error getting FTXA bank: %s", zebra_err); + goto err; + } + + if (ftxa.links[KFTXA_FTXK-1] == 0) { + sprintf(zdab_err, "FTXK link is zero!"); + goto err; + } + + rv = zebra_get_bank(f,&ftxk,ftxa.links[KFTXA_FTXK-1]); + + if (rv) { + sprintf(zdab_err, "error getting FTXK bank: %s", zebra_err); + goto err; + } + + unpack_ftxk(ftxk.data, bftxk); + + return 0; + +err: + + return -1; +} + +/* Get and unpack the RSP bank from the FTXR bank from the FTP fitter. + * + * Returns 0 if successful, or -1 on error. If there was an error, the error + * string is printed to zdab_err. */ +int get_rsp(zebraFile *f, zebraBank *ev, RSPBank *brsp) +{ + int rv; + zebraBank ft, ftp, ftxa, ftxr; + + if (ev->links[KEV_FT-1] == 0) { + sprintf(zdab_err, "FT link is zero!"); + goto err; + } + + rv = zebra_get_bank(f,&ft,ev->links[KEV_FT-1]); + + if (rv) { + sprintf(zdab_err, "error getting FT bank: %s", zebra_err); + goto err; + } + + if (ft.links[KFT_FTP-1] == 0) { + sprintf(zdab_err, "FTP link is zero!"); + goto err; + } + + rv = zebra_get_bank(f,&ftp,ft.links[KFT_FTP-1]); + + if (rv) { + sprintf(zdab_err, "error getting FTP bank: %s", zebra_err); + goto err; + } + + if (ftp.links[KFTX_FTXA-1] == 0) { + sprintf(zdab_err, "FTXA link is zero!"); + goto err; + } + + rv = zebra_get_bank(f,&ftxa,ftp.links[KFTX_FTXA-1]); + + if (rv) { + sprintf(zdab_err, "error getting FTXA bank: %s", zebra_err); + goto err; + } + + if (ftxa.links[KFTXA_FTXR-1] == 0) { + sprintf(zdab_err, "FTXR link is zero!"); + goto err; + } + + rv = zebra_get_bank(f,&ftxr,ftxa.links[KFTXA_FTXR-1]); + + if (rv) { + sprintf(zdab_err, "error getting FTXR bank: %s", zebra_err); + goto err; + } + + unpack_rsp(ftxr.data, brsp); + + return 0; + +err: + + return -1; +} + +void unpack_rsp(uint32_t *data, RSPBank *b) +{ + unpack((uint8_t *) data, "f", &b->optical_response); + unpack((uint8_t *) (data+1), "f", &b->nwin); + unpack((uint8_t *) (data+2), "f", &b->nwin2); + unpack((uint8_t *) (data+3), "f", &b->ndark); + unpack((uint8_t *) (data+4), "f", &b->neff); + unpack((uint8_t *) (data+5), "f", &b->ncor); + unpack((uint8_t *) (data+6), "f", &b->ncormc); + unpack((uint8_t *) (data+7), "f", &b->nonline); + unpack((uint8_t *) (data+8), "f", &b->ncal); + unpack((uint8_t *) (data+9), "f", &b->nefficient); + unpack((uint8_t *) (data+10), "f", &b->nworking); + unpack((uint8_t *) (data+11), "f", &b->ene); + unpack((uint8_t *) (data+12), "f", &b->uncertainty); + unpack((uint8_t *) (data+13), "f", &b->quality); + unpack((uint8_t *) (data+14), "f", &b->r_d2o); + unpack((uint8_t *) (data+15), "f", &b->r_acr); + unpack((uint8_t *) (data+16), "f", &b->r_h2o); + unpack((uint8_t *) (data+17), "f", &b->r_fresnel); + unpack((uint8_t *) (data+18), "f", &b->r_mpe); + unpack((uint8_t *) (data+19), "f", &b->r_pmtr); + unpack((uint8_t *) (data+20), "f", &b->r_eff); + unpack((uint8_t *) (data+21), "f", &b->drift); + unpack((uint8_t *) (data+22), "f", &b->nhits); + unpack((uint8_t *) (data+23), "l", &b->fit_idx); + unpack((uint8_t *) (data+24), "f", &b->nwin_allq); + unpack((uint8_t *) (data+25), "f", &b->nhits_allq); + unpack((uint8_t *) (data+26), "f", &b->nhits_dqxx); + unpack((uint8_t *) (data+27), "f", &b->nwin_pt); + unpack((uint8_t *) (data+28), "f", &b->tshift); + unpack((uint8_t *) (data+29), "f", &b->pmt_response); + unpack((uint8_t *) (data+30), "f", &b->alt_energy); + unpack((uint8_t *) (data+31), "f", &b->nckv); + unpack((uint8_t *) (data+32), "f", &b->resolution); + unpack((uint8_t *) (data+33), "f", &b->fom); + unpack((uint8_t *) (data+34), "f", &b->ncd_shad_cor); + unpack((uint8_t *) (data+35), "f", &b->rlambda); + unpack((uint8_t *) (data+36), "f", &b->omega); + unpack((uint8_t *) (data+37), "f", &b->ckvprob); + unpack((uint8_t *) (data+38), "f", &b->chaneff); + unpack((uint8_t *) (data+39), "f", &b->pmteff); + unpack((uint8_t *) (data+40), "f", &b->mpe); + unpack((uint8_t *) (data+41), "f", &b->spare1); + unpack((uint8_t *) (data+42), "f", &b->spare2); + unpack((uint8_t *) (data+43), "f", &b->spare3); + unpack((uint8_t *) (data+44), "f", &b->spare4); + unpack((uint8_t *) (data+45), "f", &b->spare5); + unpack((uint8_t *) (data+46), "f", &b->spare6); + unpack((uint8_t *) (data+47), "f", &b->spare7); + unpack((uint8_t *) (data+48), "f", &b->spare8); + unpack((uint8_t *) (data+49), "f", &b->spare9); +} + +void unpack_ftpt(uint32_t *data, FTPTBank *b) +{ + unpack((uint8_t *) data, "f",&b->u); + unpack((uint8_t *) (data+1), "f",&b->v); + unpack((uint8_t *) (data+2), "f",&b->w); + unpack((uint8_t *) (data+3), "f",&b->du); + unpack((uint8_t *) (data+4), "f",&b->dv); + unpack((uint8_t *) (data+5), "f",&b->dw); +} + +void unpack_ftpv(uint32_t *data, FTPVBank *b) +{ + unpack((uint8_t *) data, "f",&b->x); + unpack((uint8_t *) (data+1), "f",&b->y); + unpack((uint8_t *) (data+2), "f",&b->z); + unpack((uint8_t *) (data+3), "f",&b->t); + unpack((uint8_t *) (data+4), "f",&b->dx); + unpack((uint8_t *) (data+5), "f",&b->dy); + unpack((uint8_t *) (data+6), "f",&b->dz); + unpack((uint8_t *) (data+7), "f",&b->dt); +} + +void unpack_ftxk(uint32_t *data, FTXKBank *b) +{ + unpack((uint8_t *) data, "f",&b->prob); + unpack((uint8_t *) (data+1), "f",&b->energy); + unpack((uint8_t *) (data+2), "f",&b->ene_merr); + unpack((uint8_t *) (data+3), "f",&b->ene_perr); + unpack((uint8_t *) (data+4), "f",&b->neff); + unpack((uint8_t *) (data+5), "f",&b->dir_scale); + unpack((uint8_t *) (data+6), "f",&b->dir_scale_sq); + unpack((uint8_t *) (data+7), "f",&b->scat_scale); + unpack((uint8_t *) (data+8), "f",&b->refl_scale); + unpack((uint8_t *) (data+9), "f",&b->refl_av1_scale); + unpack((uint8_t *) (data+10), "f",&b->refl_av2_scale); + unpack((uint8_t *) (data+11), "f",&b->refl_ncd_scale); + unpack((uint8_t *) (data+12), "f",&b->spare1); + unpack((uint8_t *) (data+13), "f",&b->spare2); + unpack((uint8_t *) (data+14), "f",&b->spare3); + unpack((uint8_t *) (data+15), "f",&b->spare4); + unpack((uint8_t *) (data+16), "f",&b->spare5); +} + +void unpack_ftk(uint32_t *data, FTKBank *b) +{ + unpack((uint8_t *) data, "l",&b->method); + unpack((uint8_t *) (data+1), "l",&b->retc); + unpack((uint8_t *) (data+2), "l",&b->method); + unpack((uint8_t *) (data+3), "l",&b->retc); + unpack((uint8_t *) (data+4), "l",&b->pmt_avail); + unpack((uint8_t *) (data+5), "l",&b->pmt_used); + unpack((uint8_t *) (data+6), "l",&b->early); + unpack((uint8_t *) (data+7), "l",&b->late); + unpack((uint8_t *) (data+8), "l",&b->iter); + unpack((uint8_t *) (data+9), "f",&b->prob); + unpack((uint8_t *) (data+10), "f",&b->energy); + unpack((uint8_t *) (data+11), "f",&b->ene_merr); + unpack((uint8_t *) (data+12), "f",&b->ene_perr); + unpack((uint8_t *) (data+13), "f",&b->neff); + unpack((uint8_t *) (data+14), "f",&b->dir_scale); + unpack((uint8_t *) (data+15), "f",&b->dir_scale_sq); + unpack((uint8_t *) (data+16), "f",&b->scat_scale); + unpack((uint8_t *) (data+17), "f",&b->refl_scale); + unpack((uint8_t *) (data+18), "f",&b->refl_av1_scale); + unpack((uint8_t *) (data+19), "f",&b->refl_av2_scale); + unpack((uint8_t *) (data+20), "f",&b->refl_ncd_scale); + unpack((uint8_t *) (data+21), "f",&b->spare1); + unpack((uint8_t *) (data+22), "f",&b->spare2); + unpack((uint8_t *) (data+23), "f",&b->spare3); + unpack((uint8_t *) (data+24), "f",&b->spare4); + unpack((uint8_t *) (data+25), "f",&b->spare5); +} + void unpack_mcgn(uint32_t *data, MCGNBank *b) { unpack((uint8_t *) data, "l",&b->id); diff --git a/src/zdab_utils.h b/src/zdab_utils.h index 5b5fe96..d4411d9 100644 --- a/src/zdab_utils.h +++ b/src/zdab_utils.h @@ -23,6 +23,8 @@ #include "event.h" #include "zebra.h" +extern char zdab_err[256]; + // the builder won't put out events with NHIT > 10000 // (note that these are possible due to hardware problems) // but XSNOED can write an event with up to 10240 channels @@ -172,6 +174,210 @@ /* Links for the MCTK bank. */ #define KMCTK_MCVX 1 +/* Structural links for the FT bank. */ +#define KFT_FTN 19 /* Not used; FTN (QPDF fitter for the NCD Phase) does not produce a ZEBRA bank. */ +#define KFT_FTA 18 /* Attenuation Fitter */ +#define KFT_FTC 17 /* Stopping Muon */ +#define KFT_FTZ 16 /* Time-only version of FTI */ +#define KFT_FTY 15 /* Charge-only version of FTI */ +#define KFT_FTH 14 /* Hough Fitter (not yet implemented) */ +#define KFT_FTI 13 /* Muon "Impact parameter" fitter */ +#define KFT_FTK 12 /* Energy fitter */ +#define KFT_FTR 11 /* TRack fitter */ +#define KFT_FTU 10 /* QPDF fitter */ +#define KFT_FTP 9 /* Path fitter */ +#define KFT_FTM 8 /* Muon fitter */ +#define KFT_FTG 7 /* Grid fitter */ +#define KFT_FTQ 6 /* Quad fitter */ +#define KFT_FTT 5 /* Time fitter */ +#define KFT_FTE 4 /* Elastic fitter */ +#define KFT_FT3 3 /* User fitter 3 */ +#define KFT_FT2 2 /* User fitter 2 */ +#define KFT_FT1 1 /* User fitter 1 */ + +/* Structural links for FTX bank. */ +#define KFTX_FTXV 1 /* Link to first fit vertex. */ +#define KFTX_FTXA 2 /* Link to fit analysis bank. */ + +/* Structural links for FTXA bank. */ +#define KFTXA_FTXC 1 /* Link to Classifier Parameters bank */ +#define KFTXA_FTXR 2 /* Link to FTXR bank */ +#define KFTXA_FTXK 3 /* Link to FTXK bank */ +#define KFTXA_FTXP 4 /* Link to FTXP bank */ + +/* Structural links for FTXV bank. */ +#define KFTXV_FTXT 1 /* Link to first fit track. */ + +typedef struct RSPBank { + /* Optical response for this event. */ + float optical_response; + float nwin; /* Number of in-time hits. */ + float nwin2; /* Total number of in-time hits (MC only). */ + float ndark; /* Average number of noise hits in window. */ + float neff; /* Effective hits, see note 2. */ + float ncor; /* Corrected hits, see note 3. */ + float ncormc; /* Scaled corrected hits, see note 4. */ + float nonline; /* Number of on-line tubes. */ + float ncal; /* Number of calibratable tubes, see note 5. */ + float nefficient; /* Number of tubes with efficiency > 0, see note 5. */ + float nworking; /* Number of working tubes, see note 5. */ + float ene; /* Estimated energy in MeV. */ + float uncertainty; /* Uncertainty on above. */ + float quality; /* Quality of energy calibration. */ + float r_d2o; /* D2O contribution to response. */ + float r_acr; /* Acrylic contribution to response. */ + float r_h2o; /* H2O contribution to response. */ + float r_fresnel; /* Fresnel scattering contribution to response, see note 6. */ + float r_mpe; /* MPE contribution to response. */ + float r_pmtr; /* PMT response contribution to response. */ + float r_eff; /* PMT efficiency contribution to response. */ + float drift; /* Scale drift correction */ + float nhits; /* Total number of working inward-looking PMTs which have fired in this event. */ + uint32_t fit_idx; /* Index of vertex fitter used. */ + float nwin_allq; /* Number of working inward looking tubes in the prompt time window which have fired, without the software charge threshold. See note 3. */ + float nhits_allq; /* Number of working inward looking tubes which have fired without the software charge threshold. See note 3. */ + float nhits_dqxx; /* Number of online inword looking tubes which have fired. See note 3. */ + float nwin_pt; /* In-time hits with walk-correction and charge cut. */ + float tshift; /* Time shift between walk corrected and non-corrected times. */ + float pmt_response; /* RSP PMT response. */ + float alt_energy; /* Alternate RSP energy see note 8. */ + float nckv; /* Alternate RSP Ncor equivalent. */ + float resolution; /* RSP resolution estimate, not yet defined. */ + float fom; /* RSP figure of merit, not yet defined. */ + float ncd_shad_cor; /* NCD shadow correction, not yet defined. */ + float rlambda; /* New energy optical response */ + float omega; /* New energy solid angle response */ + float ckvprob; /* New energy cerenkov distribution response */ + float chaneff; /* New energy channel efficiency response */ + float pmteff; /* New energy relative PMT efficiency response */ + float mpe; /* New energy MPE response */ + float spare1; + float spare2; + float spare3; + float spare4; + float spare5; + float spare6; + float spare7; + float spare8; + float spare9; +} RSPBank; + +typedef struct FTPTBank { + /* Fitted X direction cosine. */ + float u; + /* Fitted Y direction cosine. */ + float v; + /* Fitted Z direction cosine. */ + float w; + /* Error on fitted U. */ + float du; + /* Error on fitted V. */ + float dv; + /* Error on fitted W. */ + float dw; +} FTPTBank; + +typedef struct FTPVBank { + /* Fitted X. */ + float x; + /* Fitted Y. */ + float y; + /* Fitted Z. */ + float z; + /* Fitted time. */ + float t; + /* Error on fitted X. */ + float dx; + /* Error on fitted Y. */ + float dy; + /* Error on fitted Z. */ + float dz; + /* Error on fitted time. */ + float dt; +} FTPVBank; + +typedef struct FTXKBank { + /* LH value at the best fit point. */ + float prob; + /* Total energy of the event. */ + float energy; + /* Lower uncertainty on the energy (MeV). */ + float ene_merr; + /* Upper uncertainty on the energy (MeV). */ + float ene_perr; + /* Number of hits inside the FTK time cut. */ + float neff; + /* Probability of a direct photon being detected. */ + float dir_scale; + /* MPE correction term at the mean of the cerenkov photons distribution. */ + float dir_scale_sq; + /* Probability of a scattered photon being detected. */ + float scat_scale; + /* Probability of detecting a PMT reflection photon. */ + float refl_scale; + /* Probability of detecting a AV reflection photon which reflects off the + * first AV boundary. */ + float refl_av1_scale; + /* Probability of detecting a AV reflection photon which reflects off the + * second AV boundary. */ + float refl_av2_scale; + /* Probability of detecting a NCD reflection photon. */ + float refl_ncd_scale; + float spare1; + float spare2; + float spare3; + float spare4; + float spare5; +} FTXKBank; + +typedef struct FTKBank { + /* Fitter method. */ + uint32_t method; + /* Return code. */ + uint32_t retc; + /* Number of PMT hits available. */ + uint32_t pmt_avail; + /* Number of PMT hits used. */ + uint32_t pmt_used; + /* Number of early PMT hits rejected. */ + uint32_t early; + /* Number of late PMT hits rejected. */ + uint32_t late; + /* Number of fit iterations. */ + uint32_t iter; + /* Fit probability. */ + float prob; + /* Total energy of the event. */ + float energy; + /* Lower uncertainty on the energy (MeV). */ + float ene_merr; + /* Upper uncertainty on the energy (MeV). */ + float ene_perr; + /* Number of hits inside the FTK time cut. */ + float neff; + /* Probability of a direct photon being detected. */ + float dir_scale; + /* MPE correction term at the mean of the cerenkov photons distribution. */ + float dir_scale_sq; + /* Probability of a scattered photon being detected. */ + float scat_scale; + /* Probability of detecting a PMT reflection photon. */ + float refl_scale; + /* Probability of detecting a AV reflection photon which reflects off the + * first AV boundary. */ + float refl_av1_scale; + /* Probability of detecting a AV reflection photon which reflects off the + * second AV boundary. */ + float refl_av2_scale; + /* Probability of detecting a NCD reflection photon. */ + float refl_ncd_scale; + float spare1; + float spare2; + float spare3; + float spare4; + float spare5; +} FTKBank; + typedef struct MCGNBank { /* Particle id. */ uint32_t id; @@ -423,6 +629,15 @@ typedef struct PMTBank { float qrc; } PMTBank; +int get_ftpv(zebraFile *f, zebraBank *ev, FTPVBank *bftpv); +int get_ftxk(zebraFile *f, zebraBank *ev, FTXKBank *bftxk); +int get_rsp(zebraFile *f, zebraBank *ev, RSPBank *brsp); + +void unpack_rsp(uint32_t *data, RSPBank *b); +void unpack_ftpt(uint32_t *data, FTPTBank *b); +void unpack_ftpv(uint32_t *data, FTPVBank *b); +void unpack_ftxk(uint32_t *data, FTXKBank *b); +void unpack_ftk(uint32_t *data, FTKBank *b); void unpack_mcgn(uint32_t *data, MCGNBank *b); void unpack_mcvx(uint32_t *data, MCVXBank *b); void unpack_mctk(uint32_t *data, MCTKBank *b); -- cgit