From 5fc4906e894957761e6ef118b4cf3ec2f036e87a Mon Sep 17 00:00:00 2001 From: Bea Cobo Date: Thu, 18 Feb 2021 08:57:15 +0100 Subject: [PATCH] If weird oscillations in the record => numPulses=0 in that record and WARNING --- doc/SIRENAfunctions.rst | 29 +++--- doc/_build/html/SIRENAfunctions.html | 24 +++-- libsixt/tasksSIRENA.cpp | 129 +++++++++++++++------------ libsixt/tasksSIRENA.h | 2 +- libsixt/versionSIRENA.h | 4 +- 5 files changed, 110 insertions(+), 78 deletions(-) diff --git a/doc/SIRENAfunctions.rst b/doc/SIRENAfunctions.rst index cf03537..3084c82 100755 --- a/doc/SIRENAfunctions.rst +++ b/doc/SIRENAfunctions.rst @@ -5146,7 +5146,7 @@ Search functions by name at :ref:`genindex`. Status -.. cpp:function:: int procRecord (ReconstructInitSIRENA** reconstruct_init, double tstartRecord, double samprate, fitsfile *dtcObject, gsl_vector *record, gsl_vector *recordWithoutConvert2R, PulsesCollection *foundPulses) +.. cpp:function:: int procRecord (ReconstructInitSIRENA** reconstruct_init, double tstartRecord, double samprate, fitsfile *dtcObject, gsl_vector *record, gsl_vector *recordWithoutConvert2R, PulsesCollection *foundPulses, int oscillations) Located in file: *tasksSIRENA.cpp* @@ -5158,7 +5158,9 @@ Search functions by name at :ref:`genindex`. 3) (Low-pass filtering and) differentiation - 4) Find the events (pulses) in the record + 4) If there are weird oscillations in the record, it is not processed => numPulses = 0 + + 5) Find the events (pulses) in the record - If production mode (:option:`opmode` = 1): @@ -5170,19 +5172,19 @@ Search functions by name at :ref:`genindex`. - If calibration mode (:option:`opmode` = 0): 'findPulsesCAL' - 5) Calculate the end time of the found pulses and check if the pulse is saturated + 6) Calculate the end time of the found pulses and check if the pulse is saturated - 6) Calculate the baseline (mean and standard deviation) before a pulse (in general *before*) :math:`\Rightarrow` To be written in **BSLN** and **RMSBSLN** columns in the output FITS file + 7) Calculate the baseline (mean and standard deviation) before a pulse (in general *before*) :math:`\Rightarrow` To be written in **BSLN** and **RMSBSLN** columns in the output FITS file - 7) Obtain the approximate rise and fall times of each pulse + 8) Obtain the approximate rise and fall times of each pulse - 8) Load the found pulses data in the input/output *foundPulses* structure + 9) Load the found pulses data in the input/output *foundPulses* structure - 9) Write test info (if *reconstruct_init->intermediate* = 1) + 10) Write test info (if *reconstruct_init->intermediate* = 1) - 10) Write pulses info in intermediate output FITS file (if *reconstruct_init->intermediate* = 1) + 11) Write pulses info in intermediate output FITS file (if *reconstruct_init->intermediate* = 1) - 11) Free allocated GSL vectors + 12) Free allocated GSL vectors **Members/Variables** @@ -5226,6 +5228,10 @@ Search functions by name at :ref:`genindex`. Photon ID (from the input file) to be propagated + int oscillations + + 1 (there are weird oscillations in the record) or 0 (record without weird oscillations) + .. cpp:member:: ReconstructInitSIRENA** reconstruct_init Member of *ReconstructInitSIRENA* structure to initialize the reconstruction parameters (pointer and values) @@ -5266,6 +5272,9 @@ Search functions by name at :ref:`genindex`. Photon ID (from the input file) to be propagated + .. cpp:member:: int oscillations + + 1 (there are weird oscillations in the record) or 0 (record without weird oscillations) .. cpp:function:: int pulseGrading(ReconstructInitSIRENA *reconstruct_init, int grade1, int grade2, int OFlength_strategy, int *pulseGrade, long *OFlength) @@ -5644,7 +5653,7 @@ Search functions by name at :ref:`genindex`. 4) Store the input record in *invector* (:cpp:func:`loadRecord`) - 5) Detect weird oscillations in some GSFC records + 5) Detect weird oscillations in some GSFC records providing a warning (no pulses detected in that record) 6) Convert *I* into *R* if :option:`EnergyMethod` = **I2R** or **I2RFITTED** (:cpp:func:`convertI2R`) diff --git a/doc/_build/html/SIRENAfunctions.html b/doc/_build/html/SIRENAfunctions.html index 4fd4d1f..0f7462c 100644 --- a/doc/_build/html/SIRENAfunctions.html +++ b/doc/_build/html/SIRENAfunctions.html @@ -3457,14 +3457,15 @@

SIRENA functions -
-int procRecord(ReconstructInitSIRENA **reconstruct_init, double tstartRecord, double samprate, fitsfile *dtcObject, gsl_vector *record, gsl_vector *recordWithoutConvert2R, PulsesCollection *foundPulses)
+
+int procRecord(ReconstructInitSIRENA **reconstruct_init, double tstartRecord, double samprate, fitsfile *dtcObject, gsl_vector *record, gsl_vector *recordWithoutConvert2R, PulsesCollection *foundPulses, int oscillations)

Located in file: tasksSIRENA.cpp

This function processes the input record (detecting the pulses):

  1. Declare and initialize variables

  2. Allocate GSL vectors

  3. (Low-pass filtering and) differentiation

    4. +

  4. If there are weird oscillations in the record, it is not processed => numPulses = 0

  5. Find the events (pulses) in the record

    • If production mode (opmode = 1):

      @@ -3534,6 +3535,10 @@

      SIRENA functions @@ -3664,7 +3669,7 @@

      SIRENA functionsrunDetect()).

      • +

    • Detect pulses in input record (runDetect()).

    • If reconstruction (opmode = 1) and not PCA:
      diff --git a/libsixt/tasksSIRENA.cpp b/libsixt/tasksSIRENA.cpp index 8485ac5..a439f98 100755 --- a/libsixt/tasksSIRENA.cpp +++ b/libsixt/tasksSIRENA.cpp @@ -249,11 +249,14 @@ void runDetect(TesRecord* record, int trig_reclength, int lastRecord, int nrecor cout<<"sgTEST*100/meanTEST: "< sgTES>1% is a weird oscillation // Noise records are not important if no pulses are detected in them + int oscillations = 0; if ((gsl_vector_max(invector)-meanTEST) < (meanTEST-gsl_vector_min(invector)) && ((gsl_vector_max(invector)-meanTEST>sgTEST) || (meanTEST-gsl_vector_min(invector)>sgTEST)) && (meanTEST-gsl_vector_min(invector)>2*sgTEST) && (sgTEST*100/meanTEST>1)) { + oscillations = 1; char str_nrecord[125]; snprintf(str_nrecord,125,"%d",nrecord); message = "Weird oscillations in record " + string(str_nrecord); - EP_EXIT_ERROR(message,EPFAIL); + //EP_EXIT_ERROR(message,EPFAIL); + EP_PRINT_ERROR(message,-999); // Only a warning } // Convert I into R if 'EnergyMethod' = I2R or I2RFITTED @@ -273,7 +276,7 @@ void runDetect(TesRecord* record, int trig_reclength, int lastRecord, int nrecor log_trace("Detecting..."); // Process each record - if (procRecord(reconstruct_init, tstartRecord, 1/record->delta_t, dtcObject, invector, invectorOriginal,*pulsesInRecord, pulsesAll->ndetpulses, record->pixid,record->phid_list->phid_array[0])) + if (procRecord(reconstruct_init, tstartRecord, 1/record->delta_t, dtcObject, invector, invectorOriginal,*pulsesInRecord, pulsesAll->ndetpulses, record->pixid,record->phid_list->phid_array[0], oscillations)) { message = "Cannot run routine procRecord for record processing"; EP_EXIT_ERROR(message,EPFAIL); @@ -872,11 +875,14 @@ void th_runDetect(TesRecord* record, int trig_reclength, int lastRecord, int nre } // sgTES~0.1% of meanTEST if the record is noise only => sgTES>1% is a weird oscillation // Noise records are not important if no pulses are detected in them + int oscillations = 0; if ((gsl_vector_max(invector)-meanTEST) < (meanTEST-gsl_vector_min(invector)) && ((gsl_vector_max(invector)-meanTEST>sgTEST) || (meanTEST-gsl_vector_min(invector)>sgTEST)) && (meanTEST-gsl_vector_min(invector)>2*sgTEST) && (sgTEST*100/meanTEST>1)) { + oscillations = 1; char str_nrecord[125]; snprintf(str_nrecord,125,"%d",nrecord); message = "Weird oscillations in record " + string(str_nrecord); - EP_EXIT_ERROR(message,EPFAIL); + //EP_EXIT_ERROR(message,EPFAIL); + EP_PRINT_ERROR(message,-999); // Only a warning } if ((strcmp((*reconstruct_init)->EnergyMethod,"I2R") == 0) || (strcmp((*reconstruct_init)->EnergyMethod,"I2RFITTED") == 0)) @@ -934,7 +940,7 @@ void th_runDetect(TesRecord* record, int trig_reclength, int lastRecord, int nre // Process each record // thread safe if (procRecord(reconstruct_init, tstartRecord, 1/record->delta_t, dtcObject, - invector, invectorOriginal, *pulsesInRecord, pulsesAll->ndetpulses,record->pixid,record->phid_list->phid_array[0])) + invector, invectorOriginal, *pulsesInRecord, pulsesAll->ndetpulses,record->pixid,record->phid_list->phid_array[0],oscillations)) { message = "Cannot run routine procRecord for record processing"; EP_EXIT_ERROR(message,EPFAIL); @@ -2153,6 +2159,7 @@ int loadRecord(TesRecord* record, double *time_record, gsl_vector **adc_double) * - Declare and initialize variables * - Allocate GSL vectors * - (Low-pass filtering and) differentiation + * - If there are weird oscillations in the record, it is not processed => numPulses = 0 * - Find the events (pulses) in the record * - Production mode: * - No detect: 'noDetect' if tStartPulse1!=0 @@ -2182,8 +2189,9 @@ int loadRecord(TesRecord* record, double *time_record, gsl_vector **adc_double) * - num_previousDetectedPulses: Number of previous detected pulses (to know the index to get the proper element from tstartPulse1_i in case tstartPulse1=nameFile) * - pixid: Pixel ID (from the input file) to be propagated * - phid: Photon ID (from the input file) to be propagated + * - oscillations: 1 (there are weird oscillations in the record) or 0 (record without weird oscillations) ****************************************************************************/ -int procRecord(ReconstructInitSIRENA** reconstruct_init, double tstartRecord, double samprate, fitsfile *dtcObject, gsl_vector *record, gsl_vector *recordWithoutConvert2R, PulsesCollection *foundPulses, long num_previousDetectedPulses, int pixid, int phid) +int procRecord(ReconstructInitSIRENA** reconstruct_init, double tstartRecord, double samprate, fitsfile *dtcObject, gsl_vector *record, gsl_vector *recordWithoutConvert2R, PulsesCollection *foundPulses, long num_previousDetectedPulses, int pixid, int phid, int oscillations) { int status = EPOK; string message = ""; @@ -2308,72 +2316,81 @@ gsl_vector_memcpy(recordDERIVATIVE,record);*/ gsl_vector *recordDERIVATIVEOriginal = gsl_vector_alloc(recordDERIVATIVE->size); // To be used in 'writeTestInfo' gsl_vector_memcpy(recordDERIVATIVEOriginal,recordDERIVATIVE); - // Find the events (pulses) in the record - if ((*reconstruct_init)->opmode == 1) // In PRODUCTION mode - { - int tstartFirstEvent = 0; - bool triggerCondition; - int flagTruncated; - int tstartProvided; - - if ((!isNumber((*reconstruct_init)->tstartPulse1)) || ((isNumber((*reconstruct_init)->tstartPulse1)) && (atoi((*reconstruct_init)->tstartPulse1) != 0))) - { - if (noDetect(recordDERIVATIVE, (*reconstruct_init), &numPulses, &tstartgsl, &qualitygsl, &maxDERgsl, &samp1DERgsl, num_previousDetectedPulses, samprate, tstartRecord)) - { - message = "Cannot run routine noDetect"; - EP_PRINT_ERROR(message,EPFAIL);return(EPFAIL); - } - } - else + // If there are weird oscillations in the record, it is not processed => numPulses = 0 + if (oscillations == 0) + { + // Find the events (pulses) in the record + if ((*reconstruct_init)->opmode == 1) // In PRODUCTION mode { - if (InitialTriggering (recordDERIVATIVE, nSgms, - scaleFactor, samprate, stopCriteriaMKC, kappaMKC, - &triggerCondition, &tstartFirstEvent, &flagTruncated, - &threshold)) - { - message = "Cannot run routine InitialTriggering"; - EP_PRINT_ERROR(message,EPFAIL);return(EPFAIL); - } + int tstartFirstEvent = 0; + bool triggerCondition; + int flagTruncated; + int tstartProvided; - if (strcmp((*reconstruct_init)->detectionMode,"STC") == 0) + if ((!isNumber((*reconstruct_init)->tstartPulse1)) || ((isNumber((*reconstruct_init)->tstartPulse1)) && (atoi((*reconstruct_init)->tstartPulse1) != 0))) { - if (FindSecondariesSTC ((*reconstruct_init)->maxPulsesPerRecord, recordDERIVATIVE, threshold, (*reconstruct_init), tstartFirstEvent, &numPulses, &tstartgsl, &qualitygsl, &maxDERgsl, &samp1DERgsl)) + if (noDetect(recordDERIVATIVE, (*reconstruct_init), &numPulses, &tstartgsl, &qualitygsl, &maxDERgsl, &samp1DERgsl, num_previousDetectedPulses, samprate, tstartRecord)) { - message = "Cannot run FindSecondariesSTC"; + message = "Cannot run routine noDetect"; EP_PRINT_ERROR(message,EPFAIL);return(EPFAIL); } } - else if (strcmp((*reconstruct_init)->detectionMode,"AD") == 0) - { - if (FindSecondaries ((*reconstruct_init)->maxPulsesPerRecord, - //recordDERIVATIVE, threshold,sigmaout, - recordDERIVATIVE, threshold, samprate, - (*reconstruct_init), - tstartFirstEvent, - &numPulses,&tstartgsl,&qualitygsl, &maxDERgsl,&samp1DERgsl,&lagsgsl)) + else + { + if (InitialTriggering (recordDERIVATIVE, nSgms, + scaleFactor, samprate, stopCriteriaMKC, kappaMKC, + &triggerCondition, &tstartFirstEvent, &flagTruncated, + &threshold)) { - message = "Cannot run routine FindSecondaries"; + message = "Cannot run routine InitialTriggering"; EP_PRINT_ERROR(message,EPFAIL);return(EPFAIL); } + + if (strcmp((*reconstruct_init)->detectionMode,"STC") == 0) + { + if (FindSecondariesSTC ((*reconstruct_init)->maxPulsesPerRecord, recordDERIVATIVE, threshold, (*reconstruct_init), tstartFirstEvent, &numPulses, &tstartgsl, &qualitygsl, &maxDERgsl, &samp1DERgsl)) + { + message = "Cannot run FindSecondariesSTC"; + EP_PRINT_ERROR(message,EPFAIL);return(EPFAIL); + } + } + else if (strcmp((*reconstruct_init)->detectionMode,"AD") == 0) + { + if (FindSecondaries ((*reconstruct_init)->maxPulsesPerRecord, + //recordDERIVATIVE, threshold,sigmaout, + recordDERIVATIVE, threshold, samprate, + (*reconstruct_init), + tstartFirstEvent, + &numPulses,&tstartgsl,&qualitygsl, &maxDERgsl,&samp1DERgsl,&lagsgsl)) + { + message = "Cannot run routine FindSecondaries"; + EP_PRINT_ERROR(message,EPFAIL);return(EPFAIL); + } + } } } - } - else if ((*reconstruct_init)->opmode == 0) // In CALIBRATION mode - { - if (findPulsesCAL (recordNOTFILTERED, recordDERIVATIVE, &tstartgsl, &qualitygsl, &pulseHeightsgsl, &maxDERgsl, - &numPulses, &threshold, - scaleFactor, samprate, - samplesUp, nSgms, - Lb, Lrs, - (*reconstruct_init), - stopCriteriaMKC, - kappaMKC)) + else if ((*reconstruct_init)->opmode == 0) // In CALIBRATION mode { - message = "Cannot run routine findPulsesCAL"; - EP_PRINT_ERROR(message,EPFAIL);return(EPFAIL); + if (findPulsesCAL (recordNOTFILTERED, recordDERIVATIVE, &tstartgsl, &qualitygsl, &pulseHeightsgsl, &maxDERgsl, + &numPulses, &threshold, + scaleFactor, samprate, + samplesUp, nSgms, + Lb, Lrs, + (*reconstruct_init), + stopCriteriaMKC, + kappaMKC)) + { + message = "Cannot run routine findPulsesCAL"; + EP_PRINT_ERROR(message,EPFAIL);return(EPFAIL); + } } + (*reconstruct_init)->threshold = threshold; + } + else // There are weird oscillations in the record => No processed + { + numPulses = 0; + (*reconstruct_init)->threshold = -999.0; } - (*reconstruct_init)->threshold = threshold; // Write test info if ((*reconstruct_init)->intermediate == 1) diff --git a/libsixt/tasksSIRENA.h b/libsixt/tasksSIRENA.h index 9d1d8aa..3c4c6db 100755 --- a/libsixt/tasksSIRENA.h +++ b/libsixt/tasksSIRENA.h @@ -63,7 +63,7 @@ int createLibrary(ReconstructInitSIRENA* reconstruct_init, bool *appendToLibrary int createDetectFile(ReconstructInitSIRENA* reconstruct_init, double samprate, fitsfile **dtcObject, int inputPulseLength); int filderLibrary(ReconstructInitSIRENA** reconstruct_init, double samprate); int loadRecord(TesRecord* record, double *time_record, gsl_vector **adc_double); -int procRecord(ReconstructInitSIRENA** reconstruct_init, double tstartRecord, double samprate, fitsfile *dtcObject, gsl_vector *record, gsl_vector *recordWithoutConvert2R, PulsesCollection *foundPulses,long num_previousDetectedPulses, int pixid, int phid); +int procRecord(ReconstructInitSIRENA** reconstruct_init, double tstartRecord, double samprate, fitsfile *dtcObject, gsl_vector *record, gsl_vector *recordWithoutConvert2R, PulsesCollection *foundPulses,long num_previousDetectedPulses, int pixid, int phid, int oscillations); int writePulses(ReconstructInitSIRENA** reconstruct_init, double samprate, double initialtime, gsl_vector *invectorNOTFIL, int numPulsesRecord, gsl_vector *tstart, gsl_vector *tend, gsl_vector *quality, gsl_vector *taurise, gsl_vector *taufall, fitsfile *dtcObject); int writeTestInfo(ReconstructInitSIRENA* reconstruct_init, gsl_vector *recordDERIVATIVE, double threshold, fitsfile *dtcObject); int calculateTemplate (ReconstructInitSIRENA *reconstruct_init, PulsesCollection *pulsesAll, PulsesCollection *pulsesInRecord, double samprate, gsl_vector **pulseaverage, double *pulseaverageHeight, gsl_matrix **covariance, gsl_matrix **weight, gsl_vector **pulseaverageMaxLengthFixedFilter); diff --git a/libsixt/versionSIRENA.h b/libsixt/versionSIRENA.h index b8b7d65..d0746e4 100644 --- a/libsixt/versionSIRENA.h +++ b/libsixt/versionSIRENA.h @@ -22,11 +22,11 @@ // CANTABRIA (CSIC-UC) with funding from the Spanish Ministry of Science and // Innovation (MICINN) // -// DATE: 2021/02/11, 10:38:37 +// DATE: 2021/02/18, 08:57:15 #ifndef SIRENA_VERSION_H #define SIRENA_VERSION_H -#define SIRENA_VERSION "3.8.5" +#define SIRENA_VERSION "3.8.6" #endif