AcclaimClustering.h

/* -*- C++ -*-.*********************************************************************************************
 Author: Ben Strutt
 Email: strutt@physics.ucla.edu

 Description:
             A class to cluster. What were you expecting?
***********************************************************************************************************/

#ifndef ACCLAIM_CLUSTERING_H
#define ACCLAIM_CLUSTERING_H

#include "UsefulAdu5Pat.h"
#include "AnitaConventions.h"
#include "BaseList.h"
#include "AnitaEventSummary.h"
#include "TKDTree.h"
#include "TCut.h"
#include "TEntryList.h"
#include "Math/Minimizer.h"
#include "Math/Functor.h"
#include "AcclaimOpenMP.h"
#include "TRandom3.h"
#include "TChain.h"

class TTree;
class TGraphAntarctica;
class TFile;
class TH2DAntarctica;
class TCanvas;
class TArrowAntarctica;


namespace Acclaim{

  namespace Clustering {

    const Double_t default_sigma_theta = 0.25;
    const Double_t default_sigma_phi = 0.5;
    const Double_t default_range_easting_northing = 750e3;
    const Double_t default_horizon_distance = 750e3;

    void getAngularResolution(const AnitaEventSummary* sum, AnitaPol::AnitaPol_t pol, Int_t peakInd, double& sigma_theta, double& sigma_phi);
    void getAngularResolution(double snr, double& sigma_theta, double& sigma_phi);
    TCanvas* drawAngularResolutionModel(double maxSnr);

    template <class T>
    Bool_t isVaguelyNearMcMurdo(const T& t){
      return t.longitude >= 150 && t.longitude < 180 && t.latitude >= -85 && t.latitude < -75;
    }
    template <class T>
    Bool_t inSandbox(const T& t){
      return t.longitude >= -90 && t.longitude < 90;
    }

    class Event{
    public:

      //--------------------------------------------------------------------------------
      // determined by reconstruction
      //--------------------------------------------------------------------------------
      UInt_t eventNumber;               
      Int_t run;                    
      AnitaPol::AnitaPol_t pol;             
      Int_t peakIndex;                  

      Double_t centre[3];
      Double_t latitude;                
      Double_t longitude;               
      Double_t altitude;                
      Double_t easting;                 
      Double_t northing;                
      AnitaEventSummary::PayloadLocation anita;     

      Double_t theta;                   
      Double_t phi;                 
      Double_t thetaAdjustmentRequired;         

      Double_t sigmaTheta;              
      Double_t sigmaPhi;                

      //--------------------------------------------------------------------------------
      // determined by clustering
      //--------------------------------------------------------------------------------
      Int_t nThresholds;                        
      Int_t* cluster;//[nThresholds]            /// which cluster am I associated with?
      Double_t* dThetaCluster;//[nThresholds]       /// theta distance to cluster
      Double_t* dPhiCluster;//[nThresholds]     /// phi distance to cluster

      Bool_t eventEventClustering;          
      Double_t nearestKnownBaseLogLikelihood;       
      Double_t nearestKnownBaseSurfaceSeparationKm; 
      Int_t nearestKnownBaseCluster;            
      Int_t nearestKnownBaseClusterSurface;     
      Double_t selfLogLikelihood;           
      Double_t nearestEventSurfaceDistanceKm;       
      UInt_t nearestEventSurfaceEventNumber;        
      Double_t nearestEventSurfaceLogLikelihood;    
      UInt_t nearestEventSurfaceLLEventNumber;      

      Int_t antarcticaHistBin;              
      UsefulAdu5Pat usefulPat; 
      mutable Bool_t fDebug; 

      Event(Int_t nT=0);
      Event(const AnitaEventSummary* sum, AnitaPol::AnitaPol_t pol, Int_t peakInd, Int_t nT=0);
      Event(const Event& event);
      Event& operator=(const Event& event);
      Event(int pol, int peakInd, double peak_phi, double peak_theta, int nT, UInt_t eventNumber, Int_t run,
        double anita_longitude, double anita_latitude, double anita_altitude, double anita_heading, double coherent_filtered_snr);


      TArrowAntarctica* makeArrowFromAnitaToEvent();
      void setupUsefulPat(bool calculateNow = true);
      void resetClusteringNumbers();
      void deleteArrays();
      void setNThresholds(int n);

      double logLikelihoodFromPoint(double sourceLon, double sourceLat, double sourceAlt, bool addOverHorizonPenalty=false) const;
      template<class T>
      double logLikelihoodFromPoint(const T& point, bool addOverHorizonPenalty=false) const {
    return logLikelihoodFromPoint(point.longitude, point.latitude, point.altitude, addOverHorizonPenalty);
      }
      template<class T>
      inline double cartesianSeparation(const T& event2){
    double d0 = centre[0] - event2.centre[0];
    double d1 = centre[1] - event2.centre[1];
    double d2 = centre[2] - event2.centre[2];
    return TMath::Sqrt(d0*d0 + d1*d1 + d2*d2);
      }

      inline bool operator < (const Event& other) const {
        return (eventNumber < other.eventNumber);
      }

      virtual ~Event();
      ClassDef(Event, 14)
    };



    class McEvent : public Event{
    public:
      Double_t weight;                      
      Double_t energy;                      

      McEvent(Int_t nT = 0);
      McEvent(const AnitaEventSummary* sum, AnitaPol::AnitaPol_t pol, Int_t peakInd, Int_t nT=0);
      McEvent(double weight, double energy, int pol, int peakInd, double peak_phi, double peak_theta, int nT, UInt_t eventNumber, Int_t run,
          double anita_longitude, double anita_latitude, double anita_altitude, double anita_heading, double coherent_filtered_snr);

      virtual ~McEvent(){;}
    
      ClassDef(McEvent, 1)
    };    





    //--------------------------------------------------------------------------------------------------------
    class Cluster{
    public:
      Cluster(Int_t i=-1);
      Cluster(const Event& seedEvent, Int_t i=-1);
      Cluster(const BaseList::base& base, Int_t i=-1);

      virtual ~Cluster(){ ;}    

      Double_t centre[3]; 
    
      Double_t latitude;                    
      Double_t longitude;                   
      Double_t altitude;                    
    
      Int_t numDataEvents;                  
      Double_t sumMcWeights;                    
      Int_t knownBase;                      
      Int_t index;                                              
      Int_t llEventCutInd; 
      Double_t llEventCut; 

      Int_t antarcticaHistBin; 
      Int_t seedEvent; 

      void resetClusteringNumbers();

      ClassDef(Cluster, 5)
    };

    
    class LogLikelihoodMethod {
    public:

      static const Int_t SmallClusterSizeThreshold = 100;

      LogLikelihoodMethod();
      virtual ~LogLikelihoodMethod();

      void doClustering(const char* dataGlob, const char* mcGlob, const char* outFileName, bool useAcclaimFiles=true);
      void testSmallClusters(const char* dataGlob, const char* outFileName, int clusterSizeMin=5, int clusterSizeMax=16, int nAttempts=100);
      
      bool getUseBaseList(){return fUseBaseList;}
      void setUseBaseList(bool useBaseList){ // *TOGGLE *GETTER=GetUseBaseList
    fUseBaseList = useBaseList;
      }

      void setCut(TCut cut){ fCut = cut; }
      void setCutHical(bool hc){ fCutHical = hc; } // *TOGGLE *GETTER=GetCutHical
      void setSelfLLMax(double llmax){ fSelfLLMax = llmax; } 
      void setSurfaceDistThresholdKm(Double_t dist){ surfaceDistThresholdKm = dist; }
      void setPercentOfMC(Int_t percent){ fPercentOfMC = percent; }

      bool getDebug(){return fDebug;}
      void setDebug(bool db){fDebug = db;} // *TOGGLE *GETTER=GetDebug

      size_t addCut(const TCut c){
    fThermalChainCuts.push_back(c);
    return fThermalChainCuts.size();
      }

      bool fStoreUnclusteredHistograms;
      
    private:
      Double_t getDistSqEventCluster(const Event& event, const Cluster& cluster);
      Double_t getAngDistSqEventCluster(const Event& event, const Cluster& cluster);
      void getDeltaThetaDegDeltaPhiDegEventCluster(const Event& event, const Cluster& cluster, Double_t& deltaThetaDeg, Double_t& deltaPhiDeg);


      Long64_t readInSummaries(const char* summaryGlob);
      Long64_t readInTMVATreeSummaries(const char* summaryGlob, bool isMC);
      void readInSummariesForTesting(const char* summaryGlob);
      void pickEventsFromList(int n_in_cluster);
      size_t addEvent(const AnitaEventSummary* sum, AnitaPol::AnitaPol_t pol, Int_t peakInd);
      size_t addMcEvent(const AnitaEventSummary* sum, AnitaPol::AnitaPol_t pol, Int_t peakInd);      
      void assignSingleEventToCloserCluster(Int_t eventInd, Int_t isMC, Cluster& cluster, Int_t z, double llEventCut = -1);
      void readInBaseList();

      TRandom3* tr3;

      void assignMcEventsToClusters();
      void setInitialBaseClusters();
      void makeSummaryTrees();
      void addToHistograms(TH2D* h, TH2D* h2);
      void resetClusters();
      Double_t getSumOfMcWeights();
      void initKDTree();
      Acclaim::Clustering::Event* nextEvent();

      Double_t dMin(const Event* event1, const Event* event2);
      Double_t dSum(const Event* event1, const Event* event2);
      Double_t dAsym(const Event* event1, const Event* event2);
      Double_t dFit(const Event* event1, const Event* event2);

      void testTriangleInequality();
      void testSmallClustersFromPointSource();
      Int_t removeLargeBasesNearMcMurdo();

      void doBaseEventClustering();
      void doEventEventClustering();

      void doMcEventClustering();
      void doMcBaseClustering();
      bool considerBin(const Event& event, Int_t bx, Int_t by, double& easting, double& northing);
      void nearbyEvents2(UInt_t eventInd, std::vector<UInt_t>& nearbyEvents);

      void makeAndWriteNSquaredEventEventHistograms();
      Double_t evalPairLogLikelihoodAtLonLat(const Double_t* params);




      std::vector<const Acclaim::Clustering::Event*> fFitEvent1s; 
      std::vector<const Acclaim::Clustering::Event*> fFitEvent2s; 
      Int_t fMaxFitterAttempts; 
      TCut fCut; 
      bool fCutHical; 
      bool fSelfLLMax; 
      TEntryList* fEntryList; 

      UInt_t fTestEvent1; 
      UInt_t fTestEvent2; 
      std::vector<Double_t> fFitEastings;  
      std::vector<Double_t> fFitNorthings; 
      const Int_t numMcDivisions; 
      Int_t mcDivision; // Which of the MC divisions should I read in? (runs from 0 to numMcDivisions-1)

      std::vector<Double_t> llEventCuts;                        
      Double_t surfaceDistThresholdKm;
      Bool_t fEventsAlreadyClustered;

      std::vector<std::vector<Acclaim::Clustering::Cluster> >clusters;  
      // std::vector<Acclaim::Clustering::Cluster> clusters;    /// Vector of clusters,      
      std::vector<Acclaim::Clustering::Event> events;       
      std::vector<Acclaim::Clustering::McEvent> mcEvents;   
      std::vector<std::vector<std::vector<UInt_t> > > fLookupEN;
      AntarcticaBackground fMyBackground;


      TKDTreeID* fKDTree;                     
      std::vector<Double_t> fEventEastings;
      std::vector<Double_t> fEventNorthings;

      std::vector<TH2DAntarctica*> hUnclusteredEvents;  

      bool fDebug;
      bool fUseBaseList;
      TChain* fChain;
      Int_t fPercentOfMC;
      TGraph* grTestMinimizerWalk;
      TGraph* grTestMinimizerValue;
      TGraph* grTestMinimumPosition;
      std::vector<ROOT::Math::Minimizer*> fMinimizers;
      std::vector<ROOT::Math::Functor> fFunctors;
      Int_t fROOTgErrorIgnoreLevel;
      bool fDrawNewNearbyEventsHistograms;
      bool fReadInBaseList;

      std::vector<TCut> fThermalChainCuts;
    };
  }
}

#endif