CrossCorrelator.h

/* -*- C++ -*-.*********************************************************************************************
 Author: Ben Strutt
 Email: b.strutt.12@ucl.ac.uk

 Description: 
             A Cross Correlator to interact with the ROOTified ANITA-3 data and do some interferometry. 
***********************************************************************************************************/

#ifndef CROSSCORRELATOR_H
#define CROSSCORRELATOR_H

// Ryan things
#include "UsefulAnitaEvent.h"
#include "AnitaEventCalibrator.h"
#include "AnitaGeomTool.h"
#include "UsefulAdu5Pat.h"

// #include "FFTtools.h"
#include "FancyFFTs.h"

// My things
#include "RootTools.h"

// ROOT things
#include "TGraph.h"
#include "TH2D.h"
#include "TThread.h"

// standard c++ things
#include <iostream>

// Number of phi-sectors to cross correlate between
#define DELTA_PHI_SECT 2

// Offline reconstruction definitions
#define NUM_COMBOS 336
#define NUM_THREADS 4

// Typical number of samples in waveform
#define NUM_SAMPLES 256
#define UPSAMPLE_FACTOR 40
#define NOMINAL_SAMPLING_DELTAT (1./2.6f)
#define PAD_FACTOR 2
#define GETNUMFREQS(n)((n)/2+1)

// Image definitions
#define NUM_BINS_THETA 100
#define NUM_BINS_PHI 15
#define THETA_RANGE 150
#define PHI_RANGE 22.5

#define NUM_BINS_THETA_ZOOM 100
#define NUM_BINS_PHI_ZOOM 100
// #define NUM_BINS_THETA_ZOOM 200
// #define NUM_BINS_PHI_ZOOM 200
// #define NUM_BINS_THETA_ZOOM 60
// #define NUM_BINS_PHI_ZOOM 60
#define ZOOM_BIN_SIZE_PHI 0.2
#define ZOOM_BIN_SIZE_THETA 0.2
// #define ZOOM_BIN_SIZE_PHI 0.05
// #define ZOOM_BIN_SIZE_THETA 0.05
#define THETA_RANGE_ZOOM (NUM_BINS_THETA_ZOOM*ZOOM_BIN_SIZE_THETA)
#define PHI_RANGE_ZOOM (NUM_BINS_PHI_ZOOM*ZOOM_BIN_SIZE_PHI)

// #define NUM_BINS_PHI_ZOOM_TOTAL 7200
// #define NUM_BINS_THETA_ZOOM_TOTAL 3000
#define NUM_BINS_PHI_ZOOM_TOTAL 7269
#define NUM_BINS_THETA_ZOOM_TOTAL 3069
// #define NUM_BINS_PHI_ZOOM_TOTAL 7229
// #define NUM_BINS_THETA_ZOOM_TOTAL 3029

// Anita Geometry definitions, shouldn't really be here
#define NUM_POL AnitaPol::kNotAPol
#define NUM_RING AnitaRing::kNotARing
#define DEGREES_IN_CIRCLE 360

#define MAX_NUM_PEAKS 2
#define PEAK_PHI_DEG_RANGE 10
#define PEAK_THETA_DEG_RANGE 10

#define SPEED_OF_LIGHT 2.99792458e8
#define SPEED_OF_LIGHT_NS 0.299792458

#define ALL_PHI_TRIGS 0xffff



class CrossCorrelator{

public:

  enum mapMode_t{
    kGlobal,
    kTriggered,
    kNumMapModes
  };

  enum zoomMode_t{
    kZoomedOut,
    kZoomedIn,
    kNumZoomModes
  };

  struct threadArgs{
    Long_t threadInd; 
    CrossCorrelator* ptr; 
  };

  CrossCorrelator();
  ~CrossCorrelator();

  void initializeVariables();
  void printInfo();

  void getNormalizedInterpolatedTGraphs(UsefulAnitaEvent* realEvent, AnitaPol::AnitaPol_t pol);
  void simple260MHzSatelliteNotch(AnitaPol::AnitaPol_t pol, Int_t ant);
  void simple370MHzSatelliteNotch(AnitaPol::AnitaPol_t pol, Int_t ant);
  void simpleNotch(AnitaPol::AnitaPol_t pol, Int_t ant, Double_t notchLowEdgeMHz, Double_t notchHighEdgeMHz);
  
  void renormalizeFourierDomain(AnitaPol::AnitaPol_t pol, Int_t ant);
  
  TGraph* interpolateWithStartTimeAndZeroMean(TGraph* grIn, Double_t startTime, Double_t dt, Int_t nSamp);
  void doFFTs(AnitaPol::AnitaPol_t pol);
  void correlateEvent(UsefulAnitaEvent* realEvent);
  void correlateEvent(UsefulAnitaEvent* realEvent, AnitaPol::AnitaPol_t pol);
  void reconstructEvent(UsefulAnitaEvent* usefulEvent, Int_t numFinePeaks=MAX_NUM_PEAKS, Int_t numCoarsePeaks=MAX_NUM_PEAKS);
  void findPeakValues(AnitaPol::AnitaPol_t pol, Int_t numPeaks, Double_t* peakValues,
              Double_t* phiDegs, Double_t* thetaDegs);


  void getMaxCorrelationTimeValue(AnitaPol::AnitaPol_t pol, Int_t combo,Double_t& time, Double_t& value);
  void getMaxCorrelationTimeValue(AnitaPol::AnitaPol_t pol, Int_t ant1, Int_t ant2, Double_t& time, Double_t& value);
  void getMaxUpsampledCorrelationTimeValue(AnitaPol::AnitaPol_t pol, Int_t combo, Double_t& time, Double_t& value);
  void getMaxUpsampledCorrelationTimeValue(AnitaPol::AnitaPol_t pol, Int_t ant1, Int_t ant2, Double_t& time, Double_t& value);


  
  void doAllCrossCorrelationsThreaded(AnitaPol::AnitaPol_t pol);
  void doUpsampledCrossCorrelationsThreaded(AnitaPol::AnitaPol_t pol, Int_t phiSector);

  
  Double_t getDeltaTExpected(AnitaPol::AnitaPol_t pol, Int_t ant1, Int_t ant2,
                 Double_t phiWave, Double_t thetaWave);
  // Double_t getDeltaTExpectedFast(AnitaPol::AnitaPol_t pol, Int_t ant1, Int_t ant2,
  //                    Int_t phiIndex, Int_t thetaIndex);

  Double_t relativeOffAxisDelay(AnitaPol::AnitaPol_t pol, Int_t ant1, Int_t ant2, Double_t phiDeg);
  Double_t singleAntennaOffAxisDelay(Double_t deltaPhiDeg);


  Bool_t useCombo(Int_t ant1, Int_t ant2, Int_t phiSector, Int_t deltaPhiSect);
  void fillCombosToUse();
  void do5PhiSectorCombinatorics();
  TH2D* initializeNewCombinatorics();


  void fillDeltaTLookup();
  Double_t getBin0PhiDeg();








  void getCoarsePeakInfo(AnitaPol::AnitaPol_t pol, Int_t peakIndex, Double_t& value,
             Double_t& phiDeg, Double_t& thetaDeg);
  void getFinePeakInfo(AnitaPol::AnitaPol_t pol, Int_t peakIndex, Double_t& value,
               Double_t& phiDeg, Double_t& thetaDeg);  


  
  Double_t getInterpolatedUpsampledCorrelationValue(AnitaPol::AnitaPol_t pol, Int_t combo, Double_t deltaT);

  TH2D* getMap(AnitaPol::AnitaPol_t pol, Double_t& peakValue,
           Double_t& peakPhiDeg, Double_t& peakThetaDeg,
           UShort_t l3TrigPattern=ALL_PHI_TRIGS);
    

  TH2D* getZoomMap(AnitaPol::AnitaPol_t pol);
  
  void reconstruct(AnitaPol::AnitaPol_t pol, Double_t& imagePeak,
           Double_t& peakPhiDeg, Double_t& peakThetaDeg);
  void reconstructZoom(AnitaPol::AnitaPol_t pol,
               Double_t& imagePeak, Double_t& peakPhiDeg,
               Double_t& peakThetaDeg, Double_t zoomCenterPhiDeg=0, Double_t zoomCenterThetaDeg=0);
  
  TH2D* makeGlobalImage(AnitaPol::AnitaPol_t pol, Double_t& imagePeak,
            Double_t& peakPhiDeg, Double_t& peakThetaDeg);
  TH2D* makeGlobalImage(AnitaPol::AnitaPol_t pol);

  TH2D* makeTriggeredImage(AnitaPol::AnitaPol_t pol, Double_t& imagePeak, Double_t& peakPhiDeg,
               Double_t& peakThetaDeg, UShort_t l3TrigPattern);

  TH2D* makeZoomedImage(AnitaPol::AnitaPol_t pol, UShort_t l3TrigPattern,
            Double_t zoomCenterPhiDeg, Double_t zoomCenterThetaDeg);
  TH2D* makeZoomedImage(AnitaPol::AnitaPol_t pol, Double_t& imagePeak, Double_t& peakPhiDeg,
            Double_t& peakThetaDeg, UShort_t l3TrigPattern, Double_t zoomCenterPhiDeg,
            Double_t zoomCenterThetaDeg);
  TH2D* makeZoomedImage(AnitaPol::AnitaPol_t pol, Double_t& imagePeak, Double_t& peakPhiDeg,
            Double_t& peakThetaDeg, Double_t zoomCenterPhiDeg, Double_t zoomCenterThetaDeg);

  Int_t getPhiSectorOfAntennaClosestToPhiDeg(AnitaPol::AnitaPol_t pol, Double_t phiDeg);
  TGraph* makeCoherentlySummedWaveform(AnitaPol::AnitaPol_t pol, Double_t phiDeg,
                       Double_t thetaDeg, Int_t maxDeltaPhiSect, Double_t& snr);
  TGraph* makeUpsampledCoherentlySummedWaveform(AnitaPol::AnitaPol_t pol, Double_t phiDeg,
                        Double_t thetaDeg, Int_t maxDeltaPhiSect, Double_t& snr);  
  TGraph* makeCoherentWorker(AnitaPol::AnitaPol_t pol, Double_t phiDeg, Double_t thetaDeg,
                 Int_t maxDeltaPhiSect, Double_t& snr, Int_t nSamp);
  static void* doSomeCrossCorrelationsThreaded(void* voidPtrArgs);
  static void* doSomeUpsampledCrossCorrelationsThreaded(void* voidPtrArgs);
  static void* makeSomeOfImageThreaded(void* voidPtrArgs);
  static void* makeSomeOfZoomImageThreaded(void* voidPtrArgs);
  
  void deleteAllWaveforms(AnitaPol::AnitaPol_t pol);

  TGraph* getCrossCorrelationGraph(AnitaPol::AnitaPol_t pol, Int_t ant1, Int_t ant2);
  TGraph* getUpsampledCrossCorrelationGraph(AnitaPol::AnitaPol_t pol, Int_t ant1, Int_t ant2);
  TGraph* getCrossCorrelationGraphWorker(Int_t numSamps, AnitaPol::AnitaPol_t pol,
                     Int_t ant1, Int_t ant2);
  static Int_t directlyInsertGeometry(TString pathToLindasFile, AnitaPol::AnitaPol_t pol);
  void insertPhotogrammetryGeometry();
  
  TString mapModeNames[kNumMapModes];
  TString zoomModeNames[kNumZoomModes];
  
  // Double_t interpPreFactors[NUM_POL][NUM_COMBOS][NUM_BINS_THETA][NUM_PHI*NUM_BINS_PHI]; //!< The interpolation factor for neighbouring samples
  Double_t interpPreFactors[NUM_POL][NUM_COMBOS][NUM_PHI*NUM_BINS_PHI][NUM_BINS_THETA]; 
  // Int_t offsetLows[NUM_POL][NUM_COMBOS][NUM_BINS_THETA][NUM_PHI*NUM_BINS_PHI]; //!< The interpolation factor for neighbouring samples
  Int_t offsetLows[NUM_POL][NUM_COMBOS][NUM_PHI*NUM_BINS_PHI][NUM_BINS_THETA]; 

  Double_t crossCorrelations[NUM_POL][NUM_COMBOS][NUM_SAMPLES*PAD_FACTOR]; 
  Double_t coarseMap[NUM_POL][NUM_BINS_PHI*NUM_PHI][NUM_BINS_THETA]; 

  std::set<std::vector<int> >::iterator comboSetIterator[NUM_BINS_THETA][NUM_BINS_PHI*NUM_PHI]; 
  
  Double_t partBAsZoom[NUM_POL][NUM_COMBOS][NUM_BINS_THETA_ZOOM_TOTAL]; 
  Double_t part21sZoom[NUM_POL][NUM_COMBOS][NUM_BINS_PHI_ZOOM_TOTAL]; 

  Double_t crossCorrelationsUpsampled[NUM_POL][NUM_COMBOS][NUM_SAMPLES*PAD_FACTOR*UPSAMPLE_FACTOR*PAD_FACTOR]; 
  Double_t fineMap[NUM_POL][NUM_BINS_THETA_ZOOM][NUM_BINS_PHI_ZOOM]; 
  
  std::complex<Double_t> fftsPadded[NUM_POL][NUM_SEAVEYS][GETNUMFREQS(NUM_SAMPLES*PAD_FACTOR*UPSAMPLE_FACTOR)]; 

  std::complex<Double_t> ffts[NUM_POL][NUM_SEAVEYS][GETNUMFREQS(NUM_SAMPLES*PAD_FACTOR)]; 
  TGraph* grs[NUM_POL][NUM_SEAVEYS]; 
  TGraph* grsResampled[NUM_POL][NUM_SEAVEYS]; 
  Double_t interpRMS[NUM_POL][NUM_SEAVEYS]; 
  Double_t interpRMS2[NUM_POL][NUM_SEAVEYS]; 
  Int_t comboIndices[NUM_SEAVEYS][NUM_SEAVEYS]; 

  UInt_t eventNumber[NUM_POL]; 
  UInt_t lastEventNormalized[NUM_POL]; 
  Double_t nominalSamplingDeltaT; 
  Double_t correlationDeltaT; 
  Int_t numSamples; 
  Int_t numSamplesUpsampled; 
  Int_t numCombos; 
  
  std::vector<Int_t> comboToAnt1s; 
  std::vector<Int_t> comboToAnt2s; 
  std::vector<Int_t> combosToUseGlobal[NUM_PHI]; 

  std::vector<Double_t> rArray[NUM_POL]; 
  std::vector<Double_t> phiArrayDeg[NUM_POL]; 
  std::vector<Double_t> zArray[NUM_POL]; 

  Double_t coarseMapPeakValues[NUM_POL][MAX_NUM_PEAKS]; 
  Double_t coarseMapPeakPhiDegs[NUM_POL][MAX_NUM_PEAKS]; 
  Double_t coarseMapPeakThetaDegs[NUM_POL][MAX_NUM_PEAKS]; 

  Double_t fineMapPeakValues[NUM_POL][MAX_NUM_PEAKS]; 
  Double_t fineMapPeakPhiDegs[NUM_POL][MAX_NUM_PEAKS]; 
  Double_t fineMapPeakThetaDegs[NUM_POL][MAX_NUM_PEAKS]; 

  
  Double_t zoomedThetaWaves[NUM_BINS_THETA_ZOOM_TOTAL]; 
  Double_t zoomedTanThetaWaves[NUM_BINS_THETA_ZOOM_TOTAL]; 
  Double_t zoomedCosThetaWaves[NUM_BINS_THETA_ZOOM_TOTAL]; 
  Double_t zoomedPhiWaveLookup[NUM_BINS_PHI_ZOOM_TOTAL]; 
  Double_t zoomedCosPartLookup[NUM_POL][NUM_SEAVEYS][NUM_BINS_PHI_ZOOM_TOTAL]; 
  Double_t offAxisDelays[NUM_POL][NUM_COMBOS][NUM_BINS_PHI_ZOOM_TOTAL]; 
  
  Double_t thetaWaves[NUM_BINS_THETA]; 
  Double_t phiWaveLookup[NUM_BINS_PHI*NUM_PHI]; 
  
  AnitaPol::AnitaPol_t threadPol; 
  UInt_t threadL3TrigPattern; 
  Int_t threadPhiSector; 
  Double_t threadImagePeak[NUM_THREADS]; 
  Double_t threadPeakPhiDeg[NUM_THREADS]; 
  Double_t threadPeakThetaDeg[NUM_THREADS]; 
  
  std::vector<TThread*> mapThreads; 
  std::vector<TThread*> corrThreads; 
  std::vector<TThread*> upsampledCorrThreads; 

  Int_t kOnlyThisCombo; 
  Int_t kDeltaPhiSect; 
  Int_t kUseOffAxisDelay; 
  Double_t maxDPhiDeg; 
  Int_t kDoSimpleSatelliteFiltering; 

private:
  
  std::vector<threadArgs> threadArgsVec; 

  Double_t aftForeOffset; 
  Double_t minThetaDegZoom; 
  Double_t minPhiDegZoom; 
  Double_t zoomPhiMin; 
  Double_t zoomThetaMin; 

  
};
#endif