GETSystemAnalyser.hh

Go to the documentation of this file.

//======================================================================
/*! \file GETSystemAnalyser.hh
 *
 *  Header file for the GETSystemAnalyser class.
 */
//======================================================================

#ifndef GET_SYSTEM_ANALYSER_HH
#define GET_SYSTEM_ANALYSER_HH

// includes from current library
#include "GETSystem.hh"
#include "GETSampleAnalysis.hh"

#include "RExtremaFinder.hh"

#include "TFitResultPtr.h"
#include "TFitResult.h"


//======================================================================
/*! \class GETSystemAnalyser
 *
 *  This is a base class for various analysis of experimental data
 *  from the GETSystem objects.
 *  When a GETSystemAnalyser object is associated to a GETSystem object,
 *  the analysis functions from this class are called by the GETSystem
 *  specific functions:
 *  - GETSystemAnalyser::AnalyseRunStart at the beginning of a run processing
 *  - GETSystemAnalyser::AnalyseRunStop at the end of a run processing
 *  - GETSystemAnalyser::AnalyseRawEvent when a new event has been read from
 *    a run file, \b before corrections are applied to the output samples
 *  - GETSystemAnalyser::AnalyseCorEvent when a new event has been read from
 *    a run file, \b after corrections are applied to the output samples
 *
 *  In addition, the function GETSystemAnalyser::ClearAnalysis may
 *  be called when all analysis data should be resetted.
 *  Note that this function is never called by the GETSystem class.
 *
 *  The GETSystemAnalyser::ClearEvent function is automatically called
 *  by the GETSystem object when its clears its own event data.
 *
 *  This class may be derived in order to perform specific analysis,
 *  eventually from a derived class of GETSystem: see for example
 *  the GETActarDemAnalyser for the GETActarDem system class.
 *
 *  \note
 *  When a GETSystemAnalyser (or derived) class object is used to perform
 *  events analysis while the events are not read using the GETSystem
 *  (or derived) object event processing functions, then the analyser
 *  functions (GETSystemAnalyser::AnalyseRunStart,
 *  GETSystemAnalyser::AnalyseRunStop, GETSystemAnalyser::AnalyseRawEvent and
 *  GETSystemAnalyser::AnalyseCorEvent) are not called automatically.\n
 *  This for example the case when events are read from simulation and
 *  imported in a GETSystem object.
 *  In this case, the functions must be explicitely called for the
 *  analysis to be performed and the histograms to be filled properly.
 *
 *
 *  The GETSystemAnalyser class can perform various types of analysis,
 *  depending on what has been defined in terms of analysis parameters and
 *  histograms.
 *
 *  The analysis available for the base class are:
 *
 *  - <B>automatic histograms</B> filling, if the histograms have been
 *    defined with the functions:
 *    - GETSystemAnalyser::CreateFullEventRawHisto and
 *      GETSystemAnalyser::CreateFullEventCorHisto for full
 *      event 2D histograms with all samples of all channels
 *      (raw or corrected output)
 *    - GETSystemAnalyser::CreateChannelsMaxRawHisto and
 *      GETSystemAnalyser::CreateChannelsMaxCorHisto for summary
 *      histograms with maximum amplitude of all channels
 *      (raw or corrected output)
 *
 *  - <B>predefined analysis</B> (see \ref get_analyser):
 *    - check of the <B>number of data (time buckets) per channel</B>
 *      (performed on raw outputs): \ref get_analyser_datanum
 *    - check of the <B>data continuity</B>
 *      (performed on raw outputs): \ref get_analyser_datacont
 *    - <B>noise analysis</B> (performed on corrected outputs):
 *      \ref get_analyser_noise
 *    - <B>baseline fluctuations</B> analysis: the baseline of current
 *      event is compared to the average baseline of previous events
 *      (performed on raw outputs): \ref get_analyser_baseline
 *    - <B>resolution</B> analysis: the maximum signal amplitude variation
 *      is computed: \ref get_analyser_resolution (meaningful only for
 *      constant signal such as pulser calibration).
 *
 *  In addition, it is possible to read the analysis parameters from a
 *  text file: see functions GETSystemAnalyser::ReadConfigFile and
 *  GETSystemAnalyser::ConfigCommand.
 *
 */
//======================================================================
/*! \page get_analyser GET system analyser predefined analysis
 *
 *  This page presents the procedures defined for the various analysis
 *  available in the GET system analyser class.
 *  The base class for analysis is GETSystemAnalyser.
 *
 *  These analyses should be defined \b after the GET system has been created
 *  with its number of modules and channels.
 *
 *  - \ref get_analyser_noise
 *  - \ref get_analyser_baseline
 *  - \ref get_analyser_resolution
 *  - \ref get_analyser_datanum
 *  - \ref get_analyser_datacont
 *
 *
 *
 *  \section get_analyser_noise Noise monitoring analysis
 *
 *  The noise analysis evaluates, on an event by event basis, the signal
 *  RMS of each channel.
 *
 *  The analysis is created when defining the samples interval \f$[i_0;i_1]\f$
 *  on which the noise is estimated (function GETSystemAnalyser::SetNoiseRange).
 *  The interval must be selected so that it has very little chance that
 *  a real signal appears in this range.
 *  This is in principle achieved at the beginning of the sample, if the
 *  trigger delay from GET electronics is set properly.
 *
 *  The analysis is performed by the analyser on corrected outputs
 *  (FPN, baseline...), namely in the GETSystemAnalyser::AnalyseCorEvent
 *  function (automatically called by the GETSystem::ReadEvent function).
 *
 *  For a given channel \e j, the average of the channel signal \f$S_{j}\f$
 *  is computed in the \f$[i_0;i_1]\f$ interval:
 *  \f[\left< S_{j} \right> = \frac{1}{i_1 - i_0 + 1}
                              \cdot \sum_{i=i_0}^{i_1} S_j[i]\f]
 *
 *  Then the corresponding RMS:
 *  \f[\Delta S_{j} = \sqrt{ \frac{1}{i_1 - i_0 + 1}
                                \cdot \sum_{i=i_0}^{i_1}
                                  \left(S_j[i] - \left< S_{j} \right>\right)^2}\f]
 *
 *  This RMS is the estimated noise for the channel.
 *  The analysis is performed on corrected outputs (FPN, baseline...)
 *  in the GETSystemAnalyser::AnalyseCorEvent function.
 *
 *  If the noise analysis is defined, when an event is analysed, the RMS
 *  computed for each channel is stored in an histogram
 *  (GETSystemAnalyser::GetChannelsNoiseHisto).
 *
 *  This analysis does not work when the zero suppress mode is active.
 *
 *
 *  \image html "AnalyserNoise_Event.png"
 *  This picture shows the noise histogram for one event, indicating
 *  an abnormal noise level for channels around number 2000, compared
 *  to the other system channels.
 *
 *  \image html "AnalyserNoise_Pads.png"
 *  The noise analysis functions are overloaded in the GETActarDemAnalyser
 *  class to define also a noise histogram for the ACTAR TPC (demonstrator)
 *  pad layout (following the GETSystem lookup table).
 *  This picture shows the same event as the previous picture, for pads
 *  instead of electronics channels.
 *
 *  \image html "AnalyserNoise_Cumul.png"
 *  It is possible to cumulate the results for all events in a 2D histogram
 *  using the GETSystemAnalyser::FillCumulNoiseHisto.
 *
 *
 *
 *  \section get_analyser_baseline Baseline fluctuations monitoring analysis
 *
 *  While the noise analysis is defined for a single event, the
 *  <B>baseline fluctuation</B> analysis compares the current event
 *  to the previous ones, to check that the system is in stable
 *  conditions.
 *  It appears to be a good monitor for the AsAd-AGET calibration and
 *  clock alignment.
 *
 *  The analysis is created by the GETSystemAnalyser::InitBaseFluctuations
 *  function.
 *  Like for the noise analysis, the time buckets interval  \f$[i_0;i_1]\f$
 *  must be chosen so that real signal should not appear in the analysis
 *  range.
 *  Note that the interval may be different than the one used for the
 *  noise analysis.
 *  In addition, a warnig threshold is also defined: if the baseline
 *  fluctuation overpasses this threshold for one or more channels, a
 *  warning is issued on the terminal output.
 *
 *  Since the analysis compares current event to the previous ones,
 *  it requires several events before being effective.
 *  The analysis works as follows:
 *  - for the first few events, only the average (over events) baseline
 *    is built;
 *  - for event number \e n, the average baseline (from previous
 *    events) of channel \e j is:
 *  \f[\left< B_{j}^{(n)} \right> = \frac{1}{n - 1}
                              \cdot \sum_{k=0}^{n-1}
                                      \left< S_j^{(k)} \right>_{[i_0;i_1]}\f]
 *    where \f$\left< S_{j}^{(k)} \right>\f$ is computed, for event \e k in
 *    the same way than for noise analysis;
 *  - for a given channel, the signal is taken into account only if the
 *    channel contains data;
 *  - the baseline variation for channel \e j of event \e n is:
 *  \f[\Delta B_{j}^{(n)} = \left< B_{j}^{(n)} \right> - \left< S_{j}^{(n)} \right>\f]
 *
 *  This baseline fluctuations analysis is performed on the raw data, in
 *  GETSystemAnalyser::AnalyseRawEvent function.
 *
 *  If the analysis is defined, when an event is analysed, the baseline
 *  fluctuation for each channel is stored in an histogram
 *  (GETSystemAnalyser::GetBaseFluctuationsHisto).
 *
 *  This analysis does not work when the zero suppress mode is active.
 *
 *
 *  \image html "AnalyserBaseline_Event.png"
 *  This picture shows the baseline fluctuation histogram for one event,
 *  indicating a problem for channels around 1000-1100.
 *
 *  \image html "AnalyserBaseline_Cumul.png"
 *  It is possible to cumulate the results for all events in a 2D histogram
 *  using the GETSystemAnalyser::FillCumulFluctuationsHisto.
 *
 *
 *
 *  \section get_analyser_resolution Maximum amplitude resolution analysis
 *
 *  The resolution analysis proposes a tool to estimate the channel
 *  resolution for run files with constant signal on each channel,
 *  such as pulser calibration measurements.
 *
 *  The analysis stores, on an event by event basis, the maximum amplitude
 *  of all channel in a 2D histogram (channel number versus amplitude).
 *  Then the RMS is computed for each channel, and stored in the
 *  \b resolution histogram, that is recomputed for each event.
 *
 *  Note: it has been designed to optimize the GET clocks offsets.
 *
 *
 *  \section get_analyser_datanum Channels data count monitoring analysis
 *
 *  This analysis checks the number of data (from the circular buffer
 *  memory) have been read from the data frame, for each channel.
 *  In principle, for a sample depth of 512 time buckets, this
 *  number should be exactly 512 or 0 (if partial readout mode is on).
 *  This is true if the zero suppress mode is off.
 *
 *  The analysis is defined by the GETSystemAnalyser::InitDataNumberCheck
 *  function, that sets the interval \f$[n_0;n_1]\f$ in which the number
 *  of read data is considered bad.
 *  This interval should be for example \f$[0;511]\f$ in full readout mode
 *  and \f$[1;511]\f$ in partial readout mode.
 *
 *  When analysing an event, if some channels have a number of data in
 *  this interval, a warning is issued.
 *
 *
 *
 *  \section get_analyser_datacont Data continuity (or accidents) monitoring analysis
 *
 *  A possible result of problems in the GET system hardware configuration
 *  process is to get corrupted data.
 *  This may happen in several ways, like the example below or with
 *  absolutly inconsistant channels data.
 *
 *  The data continuity proposes a test based on a maximum accepted variation
 *  \e dmax between values of contiguous time buckets.
 *  This limit is set by the GETSystemAnalyser::InitDataContinuityCheck
 *  function.
 *
 *  In order to accept fast signal variations that may happen (for short
 *  peaking times), the analysis test checks an abnormal increase and
 *  decrease of the signal (or the inverse).
 *
 *  \image html "AnalyserContinuity_Event.png"
 *  This picture illustrate the type of problem that can be detected
 *  by the continuity analysis.
 *
 */

class GETSystemAnalyser
{
  //------------------------------------------------------------
  /*! \object_doc */
  GObject(GETSystemAnalyser);
  //------------------------------------------------------------

  friend class GETSystem;

  protected:

    GETSystem   * get_ptr;            ///< Pointer to the GET system object
    GETHitMask  * hit_mask;           ///< Pointer to the hit mask (same as base class hit mask, but known as GETHitMaskXY)

    u_short     selected_data;        ///< Information on the channel data to analyse: GET::signalOut (default), GET::signalTst or GET::signalRec

    bool      * channel_bad;          ///< Whether the channel contains bad data
    bool        channel_bad_flag;     ///< Whether there are bad channels

    u_int       events_count;         ///< Number of events since analysis cleard
    u_int       events_skip;          ///< Number of first events to skip for analysis

    // Full events histograms
    TH2D      * h_full_event_raw;     ///< Pointer to the histogram with full raw data (before processing)
    TH2D      * h_full_event_cor;     ///< Pointer to the histogram with full data after corrections (if any)

    // Channels maximum analysis
    TH1D      * h_channels_max_raw;   ///< Pointer to the channels summary histogram with maximum raw output values
    TH1D      * h_channels_max_cor;   ///< Pointer to the channels summary histogram with maximum corrected output values
    TH1I      * h_channels_data_num;  ///< Pointer to the channels summary histogram with number of raw data per channel

              // the followings are defined to avoid calculating several times
              // the same data
    u_short     amplitude_mode;         ///< Signal amplitude analysis mode
    double      amplitude_window;       ///< Signal amplitude analysis window size (for fit or average)

    u_short     timing_mode;            ///< Signal timing analysis mode
    double      timing_cfd_fraction;    ///< CFD fraction for timing analysis
    double      timing_cfd_delay;       ///< CFD offset for timing analysis

    bool        amax_computed;          ///< Whether the maximum amplitude / time is computed
    double    * amax_amplitude;         ///< Array of maximum amplitude values
    double    * amax_time;              ///< Array of maximum amplitude times

    double      ch_max_signal;          ///< Maximum channel amplitude value
    u_int       ch_max_index;           ///< Channel with maximum amplitude

    // multiple maxima analysis (only for pads, not for FPN)
    RExtremaFinder  mmax_finder;        ///< Multiple maxima analyser
    bool        mmax_analysis;          ///< Whether multiple maxima analysis is requested
    bool        mmax_computed;          ///< Whether multiple maxima analysis is performed
    u_int       mmax_max;               ///< Maximum number of maxima per channel
    u_int     * mmax_num;               ///< Number of maxima for each channel
    double   ** mmax_amplitude;         ///< Maxima values arrays for each channel
    double   ** mmax_time;              ///< Maxima times arrays for each channel
    double      mmax_scale;             ///< Scaling factor for amplitude

    // Number of data per channel
    bool        data_num_on;            ///< Whether the analysis of number of data per channel is ON
    int         data_num_min;           ///< Lower value for bad range of data / channel
    int         data_num_max;           ///< Upper value for bad range of data / channel
    int         data_num_bad;           ///< Number of bad channels in analysis of the data per channel

    // Data continuity check
    bool        data_cont_on;           ///< Whether the analysis of data discontinuity is ON
    double      data_cont_max;          ///< Maximum variation for continuity
    int         data_cont_bad;          ///< Number of bad channels in analysis of the data continuity

    // Amplitude resolution analysis
    bool        resol_on;               ///< Whether amplitude resolution analysis is ON
    int         resol_dim;              ///< Dimension of amplitude axis
    double      resol_min;              ///< Lower limit of amplitude axis
    double      resol_max;              ///< Upper limit of amplitude axis
    double      resol_step;             ///< Lower limit of amplitude axis
    TH2I      * h_resol_ampl;           ///< Pointer to the channels versus amplitude maximum histogram
    TH1D      * h_resol_rms;            ///< Pointer to the channels amplitude variations histogram

    // Noise analysis
    bool        rms_on;                 ///< Whether noise analysis is ON
    u_int       rms_imin;               ///< Lower limit for RMS analysis
    u_int       rms_imax;               ///< Upper limit for RMS analysis
    TH1D      * h_channels_rms;         ///< Pointer to the channels RMS analysis summary histogram

    // Baseline fluctuations analysis
    bool        fluct_on;               ///< Whether baseline fluctuations analysis is ON
    u_int       fluct_imin;             ///< Lower limit for baseline fluctuations analysis
    u_int       fluct_imax;             ///< Upper limit for baseline fluctuations analysis
    double      fluct_warn;             ///< Baseline fluctuation value to issue a warning
    int         fluct_bad;              ///< Number of bad channels in fluctuations analysis
    double    * fluct_sum;              ///< Cumulated values for baseline check
    int       * fluct_cnt;              ///< Number of events for baseline check
    TH1D      * h_channels_fluct;       ///< Pointer to the channels baseline fluctuations analysis summary average histogram


 public:
  //----------------------------------------------------------
  /** @name Constructors, affectation, destructor */
  //@{
    GETSystemAnalyser ( GETSystem * get = NULL );
    virtual ~GETSystemAnalyser ( );
  //@}

  //----------------------------------------------------------
  /** @name General functions */
  //@{
    virtual void  Reset     ( );
    virtual int   SetSystem ( GETSystem * gptr );
    u_int         GetChannelNumber ( ) const;   // inline

    void          SetEventsSkip ( u_int n );      // inline
    u_int         GetEventsSkip ( );              // inline
    u_int         GetAnalysisEventsCounter ( );   // inline

    void          SetSelectedData ( u_short id ); // inline
    int           GetSelectedData ( ) const;      // inline
  //@}

  //----------------------------------------------------------
  /** @name Main analysis functions, called by GETSystem */
  //@{
    virtual void  ClearAnalysis ( );
    virtual int   AnalyseRunStart ( );
    virtual int   AnalyseRunStop  ( );

    virtual void  ClearEvent ( );
    virtual int   AnalyseRawEvent ( );
    virtual int   AnalyseCorEvent ( );

    virtual void  BadChannelsWarning ( );

    virtual double  GetMaximumSignal  ( ) const;  // inline
    virtual u_int   GetMaximumChannel ( ) const;  // inline
  //@}

  //----------------------------------------------------------
  /** @name Predefined analysis configuration files */
  //@{
    virtual int   ReadConfigFile ( const string & cfg );
    virtual int   ReadConfigFile ( FILE * fp );
    virtual bool  ConfigCommand  ( const GString & code, const GString & args = "", FILE * fptr = NULL );
  //@}

  //----------------------------------------------------------
  /** @name Full event data */
  //@{
    TH2D  * CreateFullEventRawHisto ( const string & hname );
    TH2D  * CreateFullEventCorHisto ( const string & hname );

    TH2D  * GetFullEventRawHisto ( ) const;     // inline
    TH2D  * GetFullEventCorHisto ( ) const;     // inline
  //@}

  //----------------------------------------------------------
  /** @name Channels amplitude analysis */
  //@{
    void            SetAmplitudeMode   ( u_short mode, double win = 0. );
    void            SetAmplitudeWindow ( double win );

    void            SetTimingMode    ( u_short mode );
    void            SetTimingModeCFD ( double frac, double delay );

    virtual int     CalcChannelOutput ( u_int ic, double & vmax, double & tmax );
    virtual int     CalcChannelMaxima ( u_int ic );
    void            SetMaximaScaling  ( double fact );  // inline

    void            InitAmplitudeData  ( );
    void            InitMaximaData     ( );
    void            ResetMaximaData    ( );
    void            ClearAmplitudeData ( );
    bool            IsAmplitudeDataComputed ( ) const;  // inline
    bool            IsMaximaDataComputed    ( ) const;  // inline
    virtual void    ComputeAmplitudeData ( );

    TH1D  * CreateChannelsMaxRawHisto ( const string & hname );
    TH1D  * CreateChannelsMaxCorHisto ( const string & hname );

    TH1D  * GetChannelsMaxRawHisto ( ) const;   // inline
    TH1D  * GetChannelsMaxCorHisto ( ) const;   // inline
  //@}

  //----------------------------------------------------------
  /** @name Analysis of number of data per channel */
  //@{
    virtual void    InitDataNumberCheck ( u_int n0, u_int n1 );
    bool            IsDataNumberCheckON ( ) const;    // inline

    virtual int     ProcessDataNumberCheck ( );
    virtual int     GetDataNumberBadCount  ( ) const; // inline
  //@}

  //----------------------------------------------------------
  /** @name Analysis of data continuity */
  //@{
    virtual void    InitDataContinuityCheck ( double dmax );
    bool            IsDataContinuityCheckON ( ) const;    // inline

    virtual bool    ContinuityCheck ( GETSample & sample );
    virtual int     ProcessDataContinuityCheck ( );
    virtual int     GetDataContinuityBadCount  ( ) const; // inline
  //@}

  //----------------------------------------------------------
  /** @name Noise analysis */
  //@{
    virtual void    SetNoiseRange ( u_int i0 = 2, u_int i1 = 509 );
    u_int           GetNoiseIndexMin ( ) const;       // inline
    u_int           GetNoiseIndexMax ( ) const;       // inline
    TH1D          * GetChannelsNoiseHisto ( ) const;  // inline
    bool            IsNoiseAnalysisON ( ) const;      // inline

    virtual int     ProcessNoiseAnalysis ( );
    virtual int     FillCumulNoiseHisto  ( TH2 * hptr );
  //@}

  //----------------------------------------------------------
  /** @name Amplitude resolution analysis */
  //@{
    virtual void    InitResolutionAnalysis  ( int n, double amin, double amax );
    virtual void    ClearResolutionAnalysis ( );
    bool            IsResolutionAnalysisON  ( ) const;   // inline
    TH2I          * GetResolutionMaxHisto   ( ) const;   // inline
    TH1D          * GetResolutionHisto      ( ) const;   // inline

    virtual int     ProcessResolutionAnalysis ( );

  //@}

  //----------------------------------------------------------
  /** @name Baseline fluctuations analysis */
  //@{
    virtual void    InitBaseFluctuations ( u_int i0 = 2, u_int i1 = 32, double warn = 0. );
    virtual void    ClearBaseFluctuations ( );
    TH1D          * GetBaseFluctuationsHisto ( ) const;   // inline
    bool            IsFluctuationsAnalysisON ( ) const;   // inline

    virtual int     ProcessFluctuationsAnalysis ( );
    virtual int     GetFluctuationsBadCount     ( ) const; // inline
    virtual int     FillCumulFluctuationsHisto  ( TH2 * hptr );
  //@}


  //----------------------------------------------------------
  /** @name ROOT related functions */
  //@{
  public:

    TH1D  * CreateChannelsHisto1D ( const string & name, const string & ytit = "", bool fpn = true );
    TH2D  * CreateChannelsHisto2D ( const string & name, int ny, double ymin, double ymax, const string & ytit = "", bool fpn = true );
    TH1I  * CreateChannelsHisto1I ( const string & name, const string & ytit = "", bool fpn = true );
    TH2I  * CreateChannelsHisto2I ( const string & name, int ny, double ymin, double ymax, const string & ytit = "", bool fpn = true );

    // the filling functions may be used for Tst, Out or Rec signals from
    // the GET system channels

    TH2D  * CreateFullEventHisto ( const string & name, bool fpn = true );
    TH2D  * CreateFullEventHisto ( const string & name,
                                   int ny, double ymin, double ymax, bool fpn = true );  // inline
    int     FillFullEventHisto   ( TH2D * hptr, bool reset = true, bool fpn = true, double scale = 1. );

    TH1D  * CreateChannelSummaryHisto ( const string & name, bool fpn = true );  // inline
    int     FillChannelsSummaryHisto  ( TH1D * hptr, u_short mode = GET::dataMaxSignal, bool fpn = true, bool reset = true );
    int     CumulChannelsSummaryHisto ( TH2  * hptr, u_short mode = GET::dataMaxSignal, bool fpn = true );

    TH1D ** CreateAllChannelsHisto ( const string & prefix = "GET_", const string & suffix = "" );
    int     FillAllChannelsHisto   ( TH1D ** htab, bool reset = true );

    /*! for use within ROOT.*/
    ClassDef(GETSystemAnalyser,0);
  //@}
};


//----------------------------------------------------------------------
//  Inline functions
#include "icc/GETSystemAnalyser.icc"


//======================================================================
#endif