AnitaEventSummary.h

#ifndef _ANITA_EVENT_SUMMARY_H_
#define _ANITA_EVENT_SUMMARY_H_

#include "TObject.h"
#include "AnitaConventions.h"
#include "Adu5Pat.h"
#include <iostream>
#include "TVector3.h" 

class UsefulAdu5Pat;
class RawAnitaHeader;
class TruthAnitaEvent;

class AnitaEventSummary : public TObject
{
 public:


  //------------------------------------------------------------------------------------
  /*************************************************************************************
   * Static members (for array sizes inside class)
   *************************************************************************************/
  static const Int_t maxDirectionsPerPol = 5; 
  static const Int_t peaksPerSpectrum = 3; 
  static const Int_t numFracPowerWindows = 5;
  static const Int_t numBlastPowerBands = 3;

  //------------------------------------------------------------------------------------
  /*************************************************************************************
   * Sub-classes
   *************************************************************************************/


  class SourceHypothesis; // Forward class declarations for PointingHypothesis utils
  class PayloadLocation; // Forward class declarations for PointingHypothesis utils

  class PointingHypothesis
  {
   public:
    PointingHypothesis() : fContainer(NULL) { ; }
    Double_t phi;  
    Double_t theta; 
    Double_t value; 
    Double_t snr; 
    Double_t mapRMS; 
    Double_t mapHistoryVal; 
    Double_t hwAngle; 
    Double_t hwAngleXPol; 
    Double_t latitude;
    Double_t longitude;
    Double_t altitude;
    Double_t distanceToSource; 
    Double_t sigma_theta;  
    Double_t sigma_phi;  
    Double_t rho;  
    Double_t chisq; 
    Double_t theta_adjustment_needed; 
    Double_t phi_separation; 
    Double_t dphi_rough;  
    Double_t dtheta_rough; 
    Bool_t triggered; 
    Bool_t triggered_xpol; 
    Bool_t masked; 
    Bool_t masked_xpol; 
    Double_t antennaPeakAverage; 

    // Most basic resolution utility functions in payload coordinates relative to ADU5-aft-fore line
    Double_t dPhi(Double_t phi) const;
    Double_t dTheta(Double_t theta, bool different_sign_conventions = false) const;

    // Peak direction relative to north (N->E->S->W)
    Double_t bearing() const;
    Double_t dPhiNorth() const;

    // Find angle from stored source hypotheses
    Double_t dPhiWais() const;
    Double_t dThetaWais(bool different_sign_conventions = true) const;
    Double_t dPhiSun() const;
    Double_t dThetaSun(bool different_sign_conventions = true) const;
    Double_t dPhiLDB() const;
    Double_t dThetaLDB(bool different_sign_conventions = true) const;
    Double_t dPhiMC() const;
    Double_t dThetaMC(bool different_sign_conventions = true) const;
    Double_t dPhiTagged() const;    
    Double_t dThetaTagged(bool different_sign_conventions = true) const;  

    // Are you within this theta/phi of a stored hypothesis?
    Bool_t closeToMC(Double_t deltaPhiDeg, Double_t deltaThetaDeg) const;
    Bool_t closeToWais(Double_t deltaPhiDeg, Double_t deltaThetaDeg) const;
    Bool_t closeToLDB(Double_t deltaPhiDeg, Double_t deltaThetaDeg) const;
    Bool_t closeToSun(Double_t deltaPhiDeg, Double_t deltaThetaDeg) const;
    Bool_t closeToTagged(Double_t deltaPhiDeg, Double_t deltaThetaDeg) const;

    // peak direction relative to trigger information
    Double_t minAbsHwAngle() const;
    Bool_t absHwAngleLessThanAbsHwAngleXPol() const;

   private:
    //----------------------------------------------------------------------------------------------------
    // WARNING! This private nonsense is fragile, and a bit hacky.
    // The //! comment after the fContainer member means it does not persist in ROOT.
    // Please do not edit that comment.
    //----------------------------------------------------------------------------------------------------
    mutable AnitaEventSummary* fContainer; 
    const AnitaEventSummary* getContainer(const char* funcName) const; 
    Double_t dPhiSource(const SourceHypothesis& source) const;   // Won't work inside TTree::Draw due to limitations in TTreeFormula, so are private, use e.g. dPhiWais() instead.
    Double_t dThetaSource(const SourceHypothesis& source, bool different_sign_conventions = true) const; // Won't work inside TTree::Draw due to limitations in TTreeFormula, so are private, use e.g. dThetaWais() instead.
    void printEvent() const;
    friend class AnitaEventSummary;

    ClassDefNV(PointingHypothesis,15);
  };

  class WaveformInfo
  {

   public:
    WaveformInfo() : fContainer(NULL), fLastEventNumberCache(0), nwMeanCache(-1),
                     nwGradCache(-1), nwInterceptCache(-1), nwChisquareCache(-1) {; }
    Double_t snr; //[0,100,16]  /// Signal to Noise of waveform
    Double_t peakHilbert;//[0,4096,21]  /// peak of hilbert envelope
    Double_t peakVal;  //[0,4096,21] /// peak value
    Double_t xPolPeakVal; //[0,4096,21]  /// Peak of xpol trace
    Double_t xPolPeakHilbert; //[0,4096,21]  /// Peak of xpol hilbert Envelope

    Double_t I,Q,U,V;  
    Double_t max_dI,max_dQ,max_dU,max_dV,polErr; 
    Int_t NPointsMaxStokes; 

    Double_t totalPower;  
    Double_t totalPowerXpol;  

    //spectrum info
    Double_t bandwidth[peaksPerSpectrum]; //[0,2,16] /// bandwidth of each peak (implementation defined, may not be comparable between analyses)
    Double_t peakFrequency[peaksPerSpectrum]; //[0,2,16] /// peak frequency of power spectrum
    Double_t peakPower[peaksPerSpectrum]; //power within +/- bandwidth of each peak
    Double_t spectrumSlope; // [-100,100,16] ///  Slope of line fit to spectrum (in log-space, so this is spectral-index)
    Double_t spectrumIntercept; // [-200,200,16] /// Intercept of line fit to spectrum (in log-space)

    //Shape parameters, computed using hilbert envelope
    // This should probably taken out into its own class
    Double_t riseTime_10_90; //[0,128,16]  /// Rise time of hilbert env from 10% to 90% of peak
    Double_t riseTime_10_50; //[0,128,16] /// Rise time of hilbert env from 10% to 50% of peak
    Double_t fallTime_90_10; //[0,128,16]/// Fall time of hilbert env from 90% to 10% of peak
    Double_t fallTime_50_10; //[0,128,16] /// Fall time of hilbert env from 50% to 10% of peak
    Double_t width_50_50;  //[0,128,16] /// Width from first envelope crossing of 50 percent of peak to last
    Double_t width_10_10;  //[0,128,16]/// Width from first envelope crossing of 10 percent of peak to last
    Double_t power_10_10;  
    Double_t power_50_50;  
    Double_t peakTime;  //[-128,384,18] /// Time that peak hilbert env occurs
    Double_t peakMoments[5];  

    Double_t impulsivityMeasure; //[-1,1, 16]  /// A number that has something to do with how impulsive it is
    Double_t bandwidthMeasure; 
    Double_t fracPowerWindowBegins[numFracPowerWindows]; //[0,128,16] /// Narrowest width containing {10%, 20%, 30%, 40%, 50%} of the total power
    Double_t fracPowerWindowEnds[numFracPowerWindows]; //[0,128,16] /// Narrowest width containing {10%, 20%, 30%, 40%, 50%} of the total power

    Int_t numAntennasInCoherent; 

    Double_t localMaxToMin; //[0,4096,21] /// Largest value of local max to neighbouring local min (see Acclaim::RootTools::getLocalMaxToMin)
    Double_t localMaxToMinTime; //[0,100,16] /// Time between local maxima and minima +ve means max is before min, -ve means min is before max
    Double_t globalMaxToMin; //[0,4096,21] /// Difference between maximum and minimum voltage
    Double_t globalMaxToMinTime; //[0,128,16] /// Time between maximum and minimum volts, +ve means max is before min, -ve means min is before max


    //some utilities for polarization info
    Double_t linearPolFrac() const;
    Double_t linearPolAngle() const;
    Double_t circPolFrac() const;
    Double_t totalPolFrac() const;
    
    //utilities for instantaneous polarization info
    Double_t instantaneousLinearPolFrac() const;
    Double_t instantaneousLinearPolAngle() const;
    Double_t instantaneousCircPolFrac() const;
    Double_t instantaneousTotalPolFrac() const;

    Double_t standardizedPeakMoment(Int_t i) const;
    inline Double_t skewness(){return standardizedPeakMoment(3);}
    inline Double_t kurtosis(){return standardizedPeakMoment(4);}

    Double_t fracPowerWindowMean() const;
    Double_t fracPowerWindowGradient() const;
    Double_t fracPowerWindowIntercept() const;
    Double_t fracPowerWindowChisquare() const;

    ClassDefNV(WaveformInfo, 18);

   private:
    friend class AnitaEventSummary;
    void cacheQuantitiesDerivedFromNarrowestWidths() const;
    mutable AnitaEventSummary* fContainer; 
    mutable UInt_t fLastEventNumberCache; 
    mutable Double_t nwMeanCache;           
    mutable Double_t nwGradCache;           
    mutable Double_t nwInterceptCache;      
    mutable Double_t nwChisquareCache;      
  };

  class ChannelInfo
  {

   public:
    ChannelInfo() : pol(AnitaPol::kNotAPol), ant(-1) {; }

    Double_t rms; //[0,1024,20]
    Double_t avgPower;
    Double_t snr; //[0,100,16]/// Signal to Noise of waveform
    Double_t peakHilbert;//[0,4096,21]  /// peak of hilbert envelope

    Double_t getPhi() const;
    inline Int_t getAnt() const {return ant;}                // could add some errors on -1 here...
    inline AnitaPol::AnitaPol_t getPol() const {return pol;} // could add some errors on 2 here...

   private:
    friend class AnitaEventSummary;
    // WARNING! THE COMMENTS TRAILING THESE AFFECT WHAT HAPPENS IN ROOT!
    // The //! comment initializer means ROOT does not store these variables when writing to files.
    // These variables are filled in the AnitaEventSummary constructor, so each ChannelInfo knows
    // its location in the ChannelInfo channels[2][48] array in AnitaEventSummary
    AnitaPol::AnitaPol_t     pol; 
    Int_t                    ant; 

    ClassDefNV(ChannelInfo, 5);
  };


  class EventFlags
  {
   public:
    EventFlags() {; }
    enum CalPulser
    {
      NONE,
      WAIS,  // is actually Hpol wais in both A3 and A4
      LDB,
      SIPLE,
      TD,
      WAIS_V, // the Vpol wais in A4
      HICAL // HICAL2 flag
    };

    Int_t isGood;
    Int_t isRF;
    Int_t isAdu5Trigger;
    Int_t isG12Trigger;
    Int_t isSoftwareTrigger;
    Int_t isMinBiasTrigger;
    Int_t isPayloadBlast;
    Int_t nadirFlag;
    Int_t strongCWFlag;
    Int_t isHPolTrigger;
    Int_t isVPolTrigger;
    Int_t isStepFunction;
    Int_t hasGlitch;

    CalPulser pulser;
    Bool_t isVarner;
    Bool_t isVarner2;

    Double_t meanPower[1+AnitaRing::kNotARing];
    Double_t medianPower[1+AnitaRing::kNotARing];
    Double_t meanPowerFiltered[1+AnitaRing::kNotARing];
    Double_t medianPowerFiltered[1+AnitaRing::kNotARing];

    Double_t maxBottomToTopRatio[AnitaPol::kNotAPol];
    Int_t maxBottomToTopRatioSector[AnitaPol::kNotAPol];
    Double_t minBottomToTopRatio[AnitaPol::kNotAPol];
    Int_t minBottomToTopRatioSector[AnitaPol::kNotAPol];

    Int_t nSectorsWhereBottomExceedsTop;

    Double_t blastFraction;

    Double_t middleOrBottomPower[numBlastPowerBands];
    Double_t topPower[numBlastPowerBands];
    Int_t middleOrBottomAnt[numBlastPowerBands];
    Int_t middleOrBottomPol[numBlastPowerBands];

    ClassDefNV(EventFlags,15);
  };





  class SourceHypothesis
  {
   public:
    SourceHypothesis() { reset(); }
    Double_t theta;
    Double_t phi;
    Double_t distance;

    Double_t mapValue[NUM_POLS];  
    Double_t mapHistoryVal[NUM_POLS]; 

    void reset(); 
    virtual ~SourceHypothesis() { ; } 

    ClassDef(SourceHypothesis,4);
  };




  class MCTruth : public SourceHypothesis
  {
   public:
    MCTruth() { reset(); }
    WaveformInfo wf[AnitaPol::kNotAPol];
    Double_t weight;
    Double_t energy;
    Double_t triggerSNR[2]; 
    TVector3 nuDirection; //in earth centered fixed coordinates (i.e. Cartesian coordinates of neutrino far away)
    double nuTheta, nuPhi; //neutrino (not RF!) direction in payload coordinates -> ray at infinity
    double interactionTheta, interactionPhi; //neutrino interaction angle to ANITA in payload coordinates
    double interactionLon, interactionLat, interactionAlt;
    double rfExitTheta, rfExitPhi; //rf exit
    double rfExitLon, rfExitLat, rfExitAlt; 
    double balloonLon, balloonLat, balloonAlt; //balloon co-ords
    void reset();
    virtual ~MCTruth() { ; } 

    ClassDefNV(MCTruth,9);
  };




  class PayloadLocation
  {
   public:
    PayloadLocation() { reset(); }
    PayloadLocation(const Adu5Pat* pat); 

    Float_t latitude;
    Float_t longitude;
    Float_t altitude;
    Float_t heading;

    Float_t prevHeading; //useful for determining rotation rate

    void reset() { latitude = -999; longitude = -999; altitude = -999; heading = -999; prevHeading = -999;};
    void update(const Adu5Pat* pat); 

    Adu5Pat pat () const {
      Adu5Pat pat;
      pat.longitude = longitude;
      pat.latitude = latitude;
      pat.altitude = altitude;
      pat.heading = heading;
      return pat;
    }

    ClassDefNV(PayloadLocation,2);
  };






  //------------------------------------------------------------------------------------
  /*************************************************************************************
   * Public member variables
   *************************************************************************************/
  Int_t run; 
  UInt_t eventNumber; 
  UInt_t realTime; 
  Int_t nPeaks[AnitaPol::kNotAPol]; 
  PointingHypothesis peak[AnitaPol::kNotAPol][maxDirectionsPerPol]; 
  WaveformInfo coherent[AnitaPol::kNotAPol][maxDirectionsPerPol]; 
  WaveformInfo deconvolved[AnitaPol::kNotAPol][maxDirectionsPerPol]; 
  WaveformInfo coherent_filtered[AnitaPol::kNotAPol][maxDirectionsPerPol]; 
  WaveformInfo deconvolved_filtered[AnitaPol::kNotAPol][maxDirectionsPerPol]; 
  ChannelInfo channels[AnitaPol::kNotAPol][NUM_SEAVEYS]; 
  EventFlags flags; 
  SourceHypothesis sun; 
  SourceHypothesis wais; 
  SourceHypothesis ldb; 
  MCTruth mc; 
  PayloadLocation anitaLocation; 


  //------------------------------------------------------------------------------------
  /*************************************************************************************
   * Public member functions
   *************************************************************************************/
  // See source file for doxygen comments
  AnitaEventSummary();
  AnitaEventSummary(const RawAnitaHeader* header);
  AnitaEventSummary(const RawAnitaHeader* header, UsefulAdu5Pat* pat, const TruthAnitaEvent * truth = 0);
  void setTriggerInfomation(const RawAnitaHeader* header);
  void setSourceInformation(UsefulAdu5Pat* pat, const TruthAnitaEvent * truth = 0);
  void zeroInternals();
  Bool_t update() const {resetNonPersistent(); return true;}

  // Utilities to find interesting entries in AnitaEventSummary
  AnitaPol::AnitaPol_t highestPol() const;
  Int_t highestPolAsInt() const;
  Int_t highestPeakInd() const;
  const PointingHypothesis& highestPeak() const;
  const WaveformInfo& highestCoherent() const;
  const WaveformInfo& highestDeconvolved() const;
  const WaveformInfo& highestCoherentFiltered() const;
  const WaveformInfo& highestDeconvolvedFiltered() const;
    
  static void setThresholdForMostImpulsive(double threshold); 
  AnitaPol::AnitaPol_t mostImpulsivePol(int whichMetric=0) const;
  Int_t mostImpulsivePolAsInt(int whichMetric=0) const;
  Int_t mostImpulsiveInd(int whichMetric=0) const;
  const PointingHypothesis& mostImpulsivePeak(int whichMetric=0) const;
  const WaveformInfo& mostImpulsiveCoherent(int whichMetric=0) const;
  const WaveformInfo& mostImpulsiveDeconvolved(int whichMetric=0) const;
  const WaveformInfo& mostImpulsiveCoherentFiltered(int whichMetric=0) const;
  const WaveformInfo& mostImpulsiveDeconvolvedFiltered(int whichMetric=0) const;

  inline Double_t weight(){return mc.weight > 0 ? mc.weight : 1;} 
  AnitaPol::AnitaPol_t trainingPol() const;
  Int_t trainingPolAsInt() const;
  Int_t trainingPeakInd() const;
  const PointingHypothesis& trainingPeak() const;
  const WaveformInfo& trainingCoherent() const;
  const WaveformInfo& trainingDeconvolved() const;
  const WaveformInfo& trainingCoherentFiltered() const;
  const WaveformInfo& trainingDeconvolvedFiltered() const;
  Int_t countChannelAboveThreshold(int threshold=100) const;

  //------------------------------------------------------------------------------------
 private:

  /*************************************************************************************
   * Private member variables/functions
   *************************************************************************************/
  // WARNING! THE COMMENTS TRAILING THESE AFFECT WHAT HAPPENS IN ROOT!
  // The //! comment initializer means ROOT does not store these variables when writing to files.
  // We want to keep it this way, since they are only used to cache the results of the utility
  // funtions to find the highest peak and peak nearest the monte carlo truth info.
  mutable Int_t                fHighestPeakIndex;  
  mutable AnitaPol::AnitaPol_t fHighestPol;        
  mutable Int_t                fTrainingPeakIndex; 
  mutable AnitaPol::AnitaPol_t fTrainingPol;       
  mutable Int_t                fMostImpulsiveIndex;
  mutable AnitaPol::AnitaPol_t fMostImpulsivePol;  
  mutable UInt_t               fLastEventNumber;   


  void findHighestPeak() const;
  void findTrainingPeak() const;
  void findMostImpulsive(int whichMetric) const;
  void resetNonPersistent() const;
  const SourceHypothesis* sourceFromTag() const;

  ClassDefNV(AnitaEventSummary, 41);
};

#endif