AnitaEventSummary.cc

#include "AnitaEventSummary.h"

#include "RawAnitaHeader.h" 
#include "UsefulAdu5Pat.h"
#include "TruthAnitaEvent.h"
#include "TBuffer.h"
#include "TClass.h"
//are these even necessary anymore? Who knows! 
//ClassImp(AnitaEventSummary)
//ClassImp(AnitaEventSummary::PointingHypothesis)
//ClassImp(AnitaEventSummary::SourceHypothesis)


const double C_IN_M_NS = 0.299792;
static double impulsivityFractionThreshold = -1; //modifies behavior of mostImpulsive() set of utils


//---------------------------------------------------------------------------------------------------------
  AnitaEventSummary::AnitaEventSummary()
: anitaLocation()
{
  zeroInternals();
}




//---------------------------------------------------------------------------------------------------------
  AnitaEventSummary::AnitaEventSummary(const RawAnitaHeader* header)
: anitaLocation()
{
  zeroInternals();

  setTriggerInfomation(header);
  eventNumber = header->eventNumber;
  run = header->run;
  realTime = header->realTime;

}





//---------------------------------------------------------------------------------------------------------
  AnitaEventSummary::AnitaEventSummary(const RawAnitaHeader* header, UsefulAdu5Pat* pat, const TruthAnitaEvent * truth)
: anitaLocation(dynamic_cast<Adu5Pat*>(pat))
{

  zeroInternals();
  setTriggerInfomation(header);
  setSourceInformation(pat,truth);
  eventNumber = header->eventNumber;
  run = header->run;
  realTime = header->realTime;
}








void AnitaEventSummary::zeroInternals(){

  fLastEventNumber = -1; 
  run = 0;
  eventNumber = 0;
  resetNonPersistent();


  for(Int_t polInd=0; polInd < AnitaPol::kNotAPol; polInd++){
    AnitaPol::AnitaPol_t pol = (AnitaPol::AnitaPol_t) polInd;
    nPeaks[polInd] = 0; 
    for(Int_t dir=0; dir < maxDirectionsPerPol; dir++)
    {
      memset(&peak[polInd][dir],0,sizeof(PointingHypothesis)); 
      memset(&coherent[polInd][dir],0,sizeof(WaveformInfo)); 
      memset(&deconvolved[polInd][dir],0,sizeof(WaveformInfo)); 
      memset(&coherent_filtered[polInd][dir],0,sizeof(WaveformInfo)); 
      memset(&deconvolved_filtered[polInd][dir],0,sizeof(WaveformInfo));

    }
    for(Int_t ant=0; ant < 48; ant++){
      memset(&channels[polInd][ant],0,sizeof(ChannelInfo));
      channels[polInd][ant].pol = pol;
      channels[polInd][ant].ant = ant;
    }
  }
  memset(&flags,0, sizeof(EventFlags)); 
  flags.pulser = EventFlags::NONE; 

  sun.reset();
  wais.reset(); 
  ldb.reset(); 
  mc.reset();

}



AnitaPol::AnitaPol_t AnitaEventSummary::highestPol() const{
  findHighestPeak();
  return fHighestPol;
}


int AnitaEventSummary::highestPolAsInt() const{
  findHighestPeak();
  return int(fHighestPol);
}


int AnitaEventSummary::highestPeakInd() const{
  findHighestPeak();
  return fHighestPeakIndex;
}



const AnitaEventSummary::PointingHypothesis& AnitaEventSummary::highestPeak() const{
  findHighestPeak();
  return peak[fHighestPol][fHighestPeakIndex];
}


const AnitaEventSummary::WaveformInfo& AnitaEventSummary::highestCoherent() const{
  findHighestPeak();
  return coherent[fHighestPol][fHighestPeakIndex];
}


const AnitaEventSummary::WaveformInfo& AnitaEventSummary::highestDeconvolved() const{
  findHighestPeak();
  return deconvolved[fHighestPol][fHighestPeakIndex];
}




const AnitaEventSummary::WaveformInfo& AnitaEventSummary::highestCoherentFiltered() const{
  findHighestPeak();
  return coherent_filtered[fHighestPol][fHighestPeakIndex];
}


const AnitaEventSummary::WaveformInfo& AnitaEventSummary::highestDeconvolvedFiltered() const{
  findHighestPeak();
  return deconvolved_filtered[fHighestPol][fHighestPeakIndex];
}









AnitaPol::AnitaPol_t AnitaEventSummary::mostImpulsivePol(int whichMetric) const{
  findMostImpulsive(whichMetric);
  return fMostImpulsivePol;
}


int AnitaEventSummary::mostImpulsivePolAsInt(int whichMetric) const{
  findMostImpulsive(whichMetric);
  return int(fMostImpulsivePol);
}


int AnitaEventSummary::mostImpulsiveInd(int whichMetric) const{
  findMostImpulsive(whichMetric);
  return fMostImpulsiveIndex;
}



const AnitaEventSummary::PointingHypothesis& AnitaEventSummary::mostImpulsivePeak(int whichMetric) const{
  findMostImpulsive(whichMetric);
  return peak[fMostImpulsivePol][fMostImpulsiveIndex];
}


const AnitaEventSummary::WaveformInfo& AnitaEventSummary::mostImpulsiveCoherent(int whichMetric) const{
  findMostImpulsive(whichMetric);
  return coherent[fMostImpulsivePol][fMostImpulsiveIndex];
}


const AnitaEventSummary::WaveformInfo& AnitaEventSummary::mostImpulsiveDeconvolved(int whichMetric) const{
  findMostImpulsive(whichMetric);
  return deconvolved[fMostImpulsivePol][fMostImpulsiveIndex];
}




const AnitaEventSummary::WaveformInfo& AnitaEventSummary::mostImpulsiveCoherentFiltered(int whichMetric) const{
  findMostImpulsive(whichMetric);
  return coherent_filtered[fMostImpulsivePol][fMostImpulsiveIndex];
}


const AnitaEventSummary::WaveformInfo& AnitaEventSummary::mostImpulsiveDeconvolvedFiltered(int whichMetric) const{
  findMostImpulsive(whichMetric);
  return deconvolved_filtered[fMostImpulsivePol][fMostImpulsiveIndex];
}


int AnitaEventSummary::countChannelAboveThreshold(int threshold) const{
  int count = 0;
  for(Int_t polInd=0; polInd < AnitaPol::kNotAPol; polInd++){
    for(Int_t ant=0; ant < 48; ant++){
      if (channels[polInd][ant].rms>threshold){
        count++;
      };
    }
  }
  return count;
}




void AnitaEventSummary::setTriggerInfomation(const RawAnitaHeader* header){  

  // setting to zero, for testing this function.
  flags.isVPolTrigger = 0;
  flags.isHPolTrigger = 0;

  // need two adjacent L3 phi-sectors to trigger ANITA-3
  for(Int_t phi=0; phi<NUM_PHI; phi++){
    Int_t phi2 = (phi+1)%NUM_PHI; // adjacent phi-sector

    Int_t trigBit = (header->l3TrigPattern >> phi) & 1;
    Int_t trigBit2 = (header->l3TrigPattern >> phi2) & 1;

    // std::cout << header->l3TrigPattern << "\t"
    //        << phi << "\t" << phi2 << "\t"
    //        << trigBit << "\t" << trigBit2 << std::endl;

    if(trigBit > 0 && trigBit2 > 0){
      flags.isVPolTrigger = 1;
    }

    Int_t trigBitH = (header->l3TrigPatternH >> phi) & 1;
    Int_t trigBitH2 = (header->l3TrigPatternH >> phi2) & 1;

    if(trigBitH > 0 && trigBitH2 > 0){
      flags.isHPolTrigger = 1;      
    }
  }

  flags.isRF = header->getTriggerBitRF();
  flags.isSoftwareTrigger = header->getTriggerBitADU5() | header->getTriggerBitG12() | header->getTriggerBitSoftExt();
  flags.isAdu5Trigger = header->getTriggerBitADU5();
  flags.isG12Trigger = header->getTriggerBitG12();
  flags.isSoftwareTrigger = header->getTriggerBitSoftExt();
  flags.isMinBiasTrigger = flags.isAdu5Trigger || flags.isG12Trigger || header->getTriggerBitG12() || flags.isSoftwareTrigger;


}




void AnitaEventSummary::setSourceInformation(UsefulAdu5Pat* pat, const TruthAnitaEvent * truth){


  pat->getSunPosition(sun.phi, sun.theta);
  sun.distance = 150e9;  // I guess in theory we could compute this! 

  pat->getThetaAndPhiWaveWaisDivide(wais.theta, wais.phi);
  wais.theta *= 180/ TMath::Pi(); 
  wais.phi *= 180/ TMath::Pi(); 
  wais.distance = pat->getWaisDivideTriggerTimeNs() * C_IN_M_NS; 

  pat->getThetaAndPhiWaveLDB(ldb.theta, ldb.phi);
  ldb.distance = pat->getLDBTriggerTimeNs() * C_IN_M_NS; 
  ldb.theta *= 180/ TMath::Pi(); 
  ldb.phi *= 180/ TMath::Pi();


  if (truth) 
  {
    pat->getThetaAndPhiWave(truth->sourceLon, truth->sourceLat, truth->sourceAlt, mc.theta,mc.phi);
    mc.theta*=TMath::RadToDeg();
    mc.phi*=TMath::RadToDeg();

    // Calculate information about the associated positions from mc
    // Store lat, long, alt and angles.
    // Note: lat, long etc are already calculated for some positions in icemc using position.cc
    //             This method converts the cartesian coords to lon/lat using usefulAdu5Pat
    //             The results are very similar
    AnitaGeomTool* geom  = AnitaGeomTool::Instance();
    // nu interaction pos
    Double_t nuPos[3];
    nuPos[0] = truth->nuPos[0];
    nuPos[1] = truth->nuPos[1];
    nuPos[2] = truth->nuPos[2];
    geom->getLatLonAltFromCartesian(nuPos,mc.interactionLat,mc.interactionLon,mc.interactionAlt);
    // Note: This gets the angles from interaction pos to balloon pos, NOT interaction pos to rfExit pos.
    pat->getThetaAndPhiWave(mc.interactionLon, mc.interactionLat, mc.interactionAlt, mc.interactionTheta,mc.interactionPhi);
    mc.interactionTheta*=TMath::RadToDeg();
    mc.interactionPhi*=TMath::RadToDeg();
    // rfExit pos (in summary tree) calculated from raw cartesian coords from truth tree
    // truth tree version = sourceLat/sourceLon/sourceAlt
    Double_t rfExitPos[3];
    rfExitPos[0] = truth->rfExitPos[2][0];
    rfExitPos[1] = truth->rfExitPos[2][1];
    rfExitPos[2] = truth->rfExitPos[2][2];
    geom->getLatLonAltFromCartesian(rfExitPos,mc.rfExitLat,mc.rfExitLon,mc.rfExitAlt);
    pat->getThetaAndPhiWave(mc.rfExitLon, mc.rfExitLat, mc.rfExitAlt, mc.rfExitTheta,mc.rfExitPhi);
    mc.rfExitTheta*=TMath::RadToDeg();
    mc.rfExitPhi*=TMath::RadToDeg();
    // balloon pos (in summary tree) calculated from raw cartesian coords from truth tree
    // truth tree version = latitude/longitude/altitude
    Double_t balloonPos[3];
    balloonPos[0] = truth->balloonPos[0];
    balloonPos[1] = truth->balloonPos[1];
    balloonPos[2] = truth->balloonPos[2];
    geom->getLatLonAltFromCartesian(balloonPos,mc.balloonLat,mc.balloonLon,mc.balloonAlt);


    if (truth->payloadPhi > 0 ) 
    {
      mc.theta = truth->payloadTheta; 
      mc.phi = truth->payloadPhi; 
    }
    else
    {
      pat->getThetaAndPhiWave(truth->sourceLon, truth->sourceLat, truth->sourceAlt, mc.theta,mc.phi);
      mc.theta*=TMath::RadToDeg();
      mc.phi*=TMath::RadToDeg();
    }

    mc.weight = truth->weight; 
    mc.distance = pat->getTriggerTimeNsFromSource(truth->sourceLat, truth->sourceLon, truth->sourceAlt);
    mc.energy = truth->nuMom; // I guess this won't be true if icemc ever simulates non-relativistic neutrinos :P
    mc.triggerSNR[0] = truth->maxSNRAtTriggerH; 
    mc.triggerSNR[1] = truth->maxSNRAtTriggerV; 
    TVector3 v(truth->nuDir[0], truth->nuDir[1], truth->nuDir[2]); 
    TVector3 p(truth->nuPos[0], truth->nuPos[1], truth->nuPos[2]); 
    pat->getThetaAndPhiWaveOfRayAtInfinity(p,v,mc.nuTheta,mc.nuPhi, true, 1e-5*TMath::DegToRad(), 1e7,&mc.nuDirection ); 
    mc.nuTheta *= TMath::RadToDeg(); 
    mc.nuPhi *= TMath::RadToDeg(); 
  }
}



void AnitaEventSummary::MCTruth::reset()
{
  SourceHypothesis::reset();
  memset(&wf[0],0,sizeof(WaveformInfo));
  memset(&wf[1],0,sizeof(WaveformInfo));
  weight = 0;
  energy = 0;
  triggerSNR[0] = 0;
  triggerSNR[1] = 0; 
  nuDirection.SetXYZ(0,0,0); 
}



void AnitaEventSummary::WaveformInfo::cacheQuantitiesDerivedFromNarrowestWidths() const{
  if(!(fContainer && fContainer->eventNumber == fLastEventNumberCache)){

    double fracPowerWindow[numFracPowerWindows];
    for(int i=0; i < numFracPowerWindows; i++){
      fracPowerWindow[i] = fracPowerWindowEnds[i] - fracPowerWindowBegins[i];
    }

    // first do the mean
    nwMeanCache = TMath::Mean(numFracPowerWindows, fracPowerWindow);

    // then the gradient
    // The power fractions x[numFracPowerWindows] = {0.1, 0.2, 0.3, 0.4, 0.5} (sum = 1.5, mean = 0.3)
    // the widths are y[numFracPowerWindows] = fracPowerWindow
    // want least squres gradient = sum over (x - xbar)(y - ybar) / ((x - xbar)^{2})
    const double gradDenom = 0.1; //0.04 + 0.01 + 0 + 0.01 + 0.04 = sum over (x[i] - mean_x)^{2}

  // gradNumerator = sum over (x[i] - mean_x)*(y[i] - mean_y), skip third term since x[i] - mean_x  = 0
  double gradNumerator = - 0.2*(fracPowerWindow[0] - nwMeanCache) + -0.1*(fracPowerWindow[1] - nwMeanCache)
    + 0.1*(fracPowerWindow[3] - nwMeanCache) + 0.2*(fracPowerWindow[4] - nwMeanCache);
  nwGradCache = gradNumerator/gradDenom;

  // then the intercept
  nwInterceptCache = nwMeanCache - nwGradCache*0.3; // 0.3 = mean power fraction

  // finally the chisquare
  nwChisquareCache = 0;
  for(int i=0; i < numFracPowerWindows; i++){
    double f = 0.1*(i+1);
    double y = nwGradCache*f + nwInterceptCache;
    double dy = (y - fracPowerWindow[i]);
    nwChisquareCache += dy*dy;
  }
  }
}


double AnitaEventSummary::WaveformInfo::fracPowerWindowMean() const{
  cacheQuantitiesDerivedFromNarrowestWidths();
  return nwMeanCache;
}

double AnitaEventSummary::WaveformInfo::fracPowerWindowGradient() const {
  cacheQuantitiesDerivedFromNarrowestWidths();
  return nwGradCache;
}

double AnitaEventSummary::WaveformInfo::fracPowerWindowIntercept() const {
  cacheQuantitiesDerivedFromNarrowestWidths();
  return nwInterceptCache;
}


double AnitaEventSummary::WaveformInfo::fracPowerWindowChisquare() const {
  cacheQuantitiesDerivedFromNarrowestWidths();
  return nwChisquareCache;
}



double AnitaEventSummary::WaveformInfo::linearPolFrac() const {

  double value = TMath::Sqrt( pow(Q,2) + pow(U,2) ) / I;

  return value;
}


double AnitaEventSummary::WaveformInfo::linearPolAngle() const {

  double value = (TMath::ATan(U/Q)/2)*TMath::RadToDeg();

  return value;

}



double AnitaEventSummary::WaveformInfo::circPolFrac() const {

  double value = TMath::Abs(V)/I;

  return value;

}

double AnitaEventSummary::WaveformInfo::totalPolFrac() const {

  double value = TMath::Sqrt(pow(Q,2) + pow(U,2) + pow(V,2))/I;

  return value;

}

double AnitaEventSummary::WaveformInfo::instantaneousLinearPolFrac() const {

  double value = TMath::Sqrt( pow(max_dQ,2) + pow(max_dU,2) ) / max_dI;

  return value;
}


double AnitaEventSummary::WaveformInfo::instantaneousLinearPolAngle() const {

  double value = (TMath::ATan(max_dU/max_dQ)/2)*TMath::RadToDeg();

  return value;

}



double AnitaEventSummary::WaveformInfo::instantaneousCircPolFrac() const {

  double value = TMath::Abs(max_dV)/max_dI;

  return value;

}

double AnitaEventSummary::WaveformInfo::instantaneousTotalPolFrac() const {

  double value = TMath::Sqrt(pow(max_dQ,2) + pow(max_dU,2) + pow(max_dV,2))/max_dI;

  return value;

}



double AnitaEventSummary::WaveformInfo::standardizedPeakMoment(int k) const {

  // the first moment k=1, is stored in peakMoments[0]...  
  if(k <= 0){
    return -1;
  }
  else if (k <= 5){
    return peakMoments[k-1]/pow(peakMoments[1], 0.5*double(k));
  }
  else{
    return -1;
  }
}



double AnitaEventSummary::ChannelInfo::getPhi() const {
  return AnitaGeomTool::Instance()->getPhiFromAnt(ant);
}





void AnitaEventSummary::SourceHypothesis::reset() {

  theta = -999;
  phi = -999;
  distance = -999;

  memset(mapHistoryVal,0,NUM_POLS*sizeof(Double_t));
  memset(mapValue,0,NUM_POLS*sizeof(Double_t));
}




AnitaEventSummary::PayloadLocation::PayloadLocation(const Adu5Pat* pat){
  update(pat);
}


void AnitaEventSummary::PayloadLocation::update(const Adu5Pat* pat){
  if(pat){
    latitude = pat->latitude;
    longitude = pat->longitude;
    altitude = pat->altitude;

    prevHeading = heading;
    heading = pat->heading;

  }
  else{
    reset();
  }
}


double AnitaEventSummary::PointingHypothesis::dPhiSource(const SourceHypothesis& source) const{
  return dPhi(source.phi);
}


double AnitaEventSummary::PointingHypothesis::dPhi(double phi2) const{
  double dPhi = phi - phi2;
  if(dPhi < -180){
    dPhi += 360;
  }
  else if(dPhi >= 180){
    dPhi -= 360;
  }

  if(dPhi < -180 || dPhi >= 180){
    std::cerr << "Warning in " << __PRETTY_FUNCTION__ << " dPhi = " << dPhi
      << " is outside expected range, peak phi = " << phi << ", source phi = "
      << phi2 << std::endl;
  }

  return dPhi;
}




double AnitaEventSummary::PointingHypothesis::dTheta(double theta2, bool different_sign_conventions) const{
  int factor = different_sign_conventions ? -1 : 1;
  double dTheta = theta - factor*theta2; // + instead of - due to sign convention difference
  return dTheta;
}


double AnitaEventSummary::PointingHypothesis::dThetaSource(const SourceHypothesis& source, bool different_sign_conventions) const{
  return dTheta(source.theta, different_sign_conventions);
}



double AnitaEventSummary::PointingHypothesis::bearing() const{
  // heading increases clockwise, payload phi increases anti-clockwise so we subtract it from heading.
  const AnitaEventSummary* sum = getContainer(__PRETTY_FUNCTION__);
  double bearing = -9999;
  if(sum){
    bearing = double(sum->anitaLocation.heading) - phi;

    bearing = bearing < 0 ? bearing + 360 : bearing;
    bearing = bearing >= 360 ? bearing - 360 : bearing;

    if(bearing < 0 || bearing >= 360){
      std::cerr << "Warning in " << __PRETTY_FUNCTION__ << ", peak bearing = "
        << bearing << ", phi = " << phi << ", heading = " << sum->anitaLocation.heading << std::endl;
    }
  }
  return bearing;
}

double AnitaEventSummary::PointingHypothesis::dPhiNorth() const{
  // heading increases clockwise, payload phi increases anti-clockwise so we subtract it from heading.
  Double_t b = bearing();
  return  b > 180 ? b - 360 : b;
}



AnitaPol::AnitaPol_t AnitaEventSummary::trainingPol() const {
  findTrainingPeak();
  return fTrainingPol;
}

int AnitaEventSummary::trainingPolAsInt() const {
  findTrainingPeak();
  return int(fTrainingPol);
}

int AnitaEventSummary::trainingPeakInd() const {
  findTrainingPeak();
  return fTrainingPeakIndex;
}

const AnitaEventSummary::PointingHypothesis& AnitaEventSummary::trainingPeak() const {
  findTrainingPeak();
  return peak[fTrainingPol][fTrainingPeakIndex];
}

const AnitaEventSummary::WaveformInfo& AnitaEventSummary::trainingCoherent() const {
  findTrainingPeak();
  return coherent[fTrainingPol][fTrainingPeakIndex];
}

const AnitaEventSummary::WaveformInfo& AnitaEventSummary::trainingDeconvolved() const {
  findTrainingPeak();
  return deconvolved[fTrainingPol][fTrainingPeakIndex];
}

const AnitaEventSummary::WaveformInfo& AnitaEventSummary::trainingCoherentFiltered() const {
  findTrainingPeak();
  return coherent_filtered[fTrainingPol][fTrainingPeakIndex];
}

const AnitaEventSummary::WaveformInfo& AnitaEventSummary::trainingDeconvolvedFiltered() const {
  findTrainingPeak();
  return deconvolved_filtered[fTrainingPol][fTrainingPeakIndex];
}


const AnitaEventSummary* AnitaEventSummary::PointingHypothesis::getContainer(const char* funcName) const{
  if(!fContainer){
    std::cerr << "Error in " << funcName << ", don't have access to AnitaEventSummary that contains me!\n" 
      << "To fix this error, try calling AnitaEventSummary()::update() as the first argument of the cuts in TTree::Draw(),\n"
      << "Note: update() always returns true, so TTree::Draw(\"sum.peak[0][0].dPhiWais()\", "
      << "\"sum.update() && TMath::Abs(sum.peak[0][0].dPhiWais()) < 5\") will select all events within 5 degrees of phi of WAIS.\n"
      << std::endl;
  }
  return fContainer;
}


double AnitaEventSummary::PointingHypothesis::minAbsHwAngle() const {
  if(TMath::Abs(hwAngle) < TMath::Abs(hwAngleXPol)){
    return hwAngle;
  }
  else{
    return hwAngleXPol;
  }
}

Bool_t AnitaEventSummary::PointingHypothesis::absHwAngleLessThanAbsHwAngleXPol() const {
  return (TMath::Abs(hwAngle) < TMath::Abs(hwAngleXPol));
}




double AnitaEventSummary::PointingHypothesis::dPhiWais() const {
  return getContainer(__PRETTY_FUNCTION__) ? dPhiSource(fContainer->wais) : dPhi(-9999); // should trigger warning message
}
double AnitaEventSummary::PointingHypothesis::dThetaWais(bool different_sign_conventions) const {
  return getContainer(__PRETTY_FUNCTION__) ? dThetaSource(fContainer->wais, different_sign_conventions) : dPhi(-9999); // should trigger warning message
}
double AnitaEventSummary::PointingHypothesis::dPhiSun() const {
  return getContainer(__PRETTY_FUNCTION__) ? dPhiSource(fContainer->sun) : dPhi(-9999); // should trigger warning message
}
double AnitaEventSummary::PointingHypothesis::dThetaSun(bool different_sign_conventions) const {
  return getContainer(__PRETTY_FUNCTION__) ? dThetaSource(fContainer->sun, different_sign_conventions) : dPhi(-9999); // should trigger warning message
}
double AnitaEventSummary::PointingHypothesis::dPhiLDB() const {
  return getContainer(__PRETTY_FUNCTION__) ? dPhiSource(fContainer->ldb) : dPhi(-9999); // should trigger warning message
}
double AnitaEventSummary::PointingHypothesis::dThetaLDB(bool different_sign_conventions) const {
  return getContainer(__PRETTY_FUNCTION__) ? dThetaSource(fContainer->ldb, different_sign_conventions) : dPhi(-9999); // should trigger warning message
}
double AnitaEventSummary::PointingHypothesis::dPhiMC() const {
  return getContainer(__PRETTY_FUNCTION__) && fContainer->mc.weight > 0 ? dPhiSource(fContainer->mc) : -9999;
}
double AnitaEventSummary::PointingHypothesis::dThetaMC(bool different_sign_conventions) const {
  return getContainer(__PRETTY_FUNCTION__) && fContainer->mc.weight > 0 ? dThetaSource(fContainer->mc, different_sign_conventions) : -9999;
}
double AnitaEventSummary::PointingHypothesis::dPhiTagged() const {
  return getContainer(__PRETTY_FUNCTION__) && fContainer->sourceFromTag() ? dPhiSource(*fContainer->sourceFromTag()) : -9999;
}
double AnitaEventSummary::PointingHypothesis::dThetaTagged(bool different_sign_conventions) const {
  return getContainer(__PRETTY_FUNCTION__) && fContainer->sourceFromTag() ? dThetaSource(*fContainer->sourceFromTag(), different_sign_conventions) : -9999;
}



Bool_t AnitaEventSummary::PointingHypothesis::closeToMC(double deltaPhiDeg, double deltaThetaDeg) const {
  return TMath::Abs(dThetaMC()) < deltaThetaDeg && TMath::Abs(dPhiMC()) < deltaPhiDeg;
}
Bool_t AnitaEventSummary::PointingHypothesis::closeToWais(double deltaPhiDeg, double deltaThetaDeg) const {
  return TMath::Abs(dThetaWais()) < deltaThetaDeg && TMath::Abs(dPhiWais()) < deltaPhiDeg;
}
Bool_t AnitaEventSummary::PointingHypothesis::closeToLDB(double deltaPhiDeg, double deltaThetaDeg) const {
  return TMath::Abs(dThetaLDB()) < deltaThetaDeg && TMath::Abs(dPhiLDB()) < deltaPhiDeg;
}
Bool_t AnitaEventSummary::PointingHypothesis::closeToSun(double deltaPhiDeg, double deltaThetaDeg) const {
  return TMath::Abs(dThetaSun()) < deltaThetaDeg && TMath::Abs(dPhiSun()) < deltaPhiDeg;
}
Bool_t AnitaEventSummary::PointingHypothesis::closeToTagged(double deltaPhiDeg, double deltaThetaDeg) const {
  const AnitaEventSummary* s = getContainer(__PRETTY_FUNCTION__);
  if(s){
    const AnitaEventSummary::SourceHypothesis* taggedSource = s->sourceFromTag();
    if(taggedSource){
      return TMath::Abs(dThetaSource(*taggedSource)) < deltaThetaDeg && TMath::Abs(dPhiSource(*taggedSource)) < deltaPhiDeg;
    }
  }
  return false;
}



void AnitaEventSummary::findHighestPeak() const {
  resetNonPersistent();
  if(fHighestPeakIndex < 0){ // then we've not done this before
    double highestVal = -1e99;
    for(int polInd=0; polInd < AnitaPol::kNotAPol; polInd++){
      AnitaPol::AnitaPol_t pol = (AnitaPol::AnitaPol_t) polInd;
      for(int peakInd=0; peakInd < nPeaks[polInd]; peakInd++){
        if(peak[polInd][peakInd].value > highestVal){
          highestVal = peak[polInd][peakInd].value;
          fHighestPeakIndex = peakInd;
          fHighestPol = pol;
        }
      }
    }
  }
}


void AnitaEventSummary::findMostImpulsive(int whichMetric) const {
  resetNonPersistent();
  if(fMostImpulsiveIndex < 0){ // then we've not done this before
    double highestVal = -1e99;

    for(int polInd=0; polInd < AnitaPol::kNotAPol; polInd++){
      AnitaPol::AnitaPol_t pol = (AnitaPol::AnitaPol_t) polInd;
      for(int peakInd=0; peakInd < nPeaks[polInd]; peakInd++){

        const WaveformInfo& wave = deconvolved_filtered[polInd][peakInd];

        // select which impulsivity measure to use...
        // 0 -> Cosmin's impulsivityMeasure
        // 1 -> Ben's fracPowerWindowGradient()

        // VERY IMPORTANT FACTOR OF -1 HERE
        // as lower == better for this metric
        // 2 -> Cosmins's impulsivityMeasure * peak.value 
        double val;
        switch (whichMetric)
        {
          case 1:
            val = -1*wave.fracPowerWindowGradient();
            break;
          case 2:
            val = wave.impulsivityMeasure * peak[polInd][peakInd].value;
            break;
          case 0:
          default:
            val = wave.impulsivityMeasure;
            break;
        };
        if(val > highestVal){
          highestVal = val;
          fMostImpulsiveIndex = peakInd;
          fMostImpulsivePol = pol;
        }
      }
    }

    if(whichMetric != 2 && impulsivityFractionThreshold > 0. && impulsivityFractionThreshold < 1.)
    {
      int prevPol = int(fMostImpulsivePol);
      int prevInd = fMostImpulsiveIndex;
      const WaveformInfo& wave = deconvolved_filtered[prevPol][prevInd];
      double highestVal = whichMetric <= 0 ? wave.impulsivityMeasure : -1*wave.fracPowerWindowGradient();
      double bright = peak[prevPol][prevInd].value;
      for(int polInd=0; polInd < AnitaPol::kNotAPol; polInd++){
        AnitaPol::AnitaPol_t pol = (AnitaPol::AnitaPol_t) polInd;
        for(int peakInd=0; peakInd < nPeaks[polInd]; peakInd++)
        {
          if((peakInd == prevInd) && (prevPol == polInd)) continue;
          const WaveformInfo& wave2 = deconvolved_filtered[polInd][peakInd];

          double val = whichMetric <= 0 ? wave.impulsivityMeasure : -1*wave.fracPowerWindowGradient();

          if(val/highestVal < impulsivityFractionThreshold) continue;
          if(bright < peak[polInd][peakInd].value)
          {
            bright = peak[polInd][peakInd].value;
            fMostImpulsiveIndex = peakInd;
            fMostImpulsivePol = pol;
          } 
        }
      }
    }
  }
}


const AnitaEventSummary::SourceHypothesis* AnitaEventSummary::sourceFromTag() const {
  switch(flags.pulser){
    case EventFlags::WAIS_V:
    case EventFlags::WAIS:
      // std::cerr << "wais" << std::endl;
      return &wais;
    case EventFlags::LDB:  
      // std::cerr << "ldb" << std::endl;
      return &ldb;      
    default:
      if(mc.weight > 0){
        // std::cerr << "mc" << std::endl;
        return &mc;
      }
      else{
        // std::cerr << "null" << std::endl;
        return NULL;
      }
  }
}


void AnitaEventSummary::findTrainingPeak() const {
  resetNonPersistent();

  if(fTrainingPeakIndex < 0){
    // then we've not done this before
    // and we need to figure out the peak of interest

    // Time to make this yet more complicated...
    // in the case we have a calPulser tagged event
    // return the peak closest to that source...
    const SourceHypothesis* peakOfInterest = sourceFromTag();
    if(peakOfInterest){
      // double lowestCloseFracPowWinGrad = DBL_MAX;
      double highestClosePeakVal = -1;
      for(int polInd=0; polInd < AnitaPol::kNotAPol; polInd++){
        AnitaPol::AnitaPol_t pol = (AnitaPol::AnitaPol_t) polInd;
        for(int peakInd=0; peakInd < nPeaks[polInd]; peakInd++){
          const double dPhiClose = 5.5;
          const double dThetaClose = 3.5;
          if(peak[polInd][peakInd].closeToTagged(dPhiClose, dThetaClose)){
            // double fpwg = deconvolved_filtered[polInd][peakInd].fracPowerWindowGradient();
            // if(fpwg  < lowestCloseFracPowWinGrad){
            double mp = peak[polInd][peakInd].value;
            if(mp > highestClosePeakVal){
              // lowestCloseFracPowWinGrad = fpwg;
              highestClosePeakVal = mp;
              fTrainingPeakIndex = peakInd;
              fTrainingPol = pol;
            }
          }
          }
        }
      }
      if(fTrainingPeakIndex < 0){
        // didn't find one or no tagged source
        // so, just do highest peak in map

        findHighestPeak();
        fTrainingPeakIndex = fHighestPeakIndex;
        fTrainingPol = fHighestPol;

        // const int metric = 1;
        // findMostImpulsive(metric);
        // fTrainingPeakIndex = fMostImpulsiveIndex;
        // fTrainingPol = fMostImpulsivePol;
      }
    }
  }



  void AnitaEventSummary::resetNonPersistent() const{
    if(fLastEventNumber!=eventNumber){
      fHighestPeakIndex = -1;
      fHighestPol = AnitaPol::kNotAPol;
      fMostImpulsiveIndex = -1;
      fMostImpulsivePol = AnitaPol::kNotAPol;
      fTrainingPeakIndex = -1;
      fTrainingPol = AnitaPol::kNotAPol;
      for(Int_t polInd=0; polInd < AnitaPol::kNotAPol; polInd++){
        for(Int_t dir=0; dir < maxDirectionsPerPol; dir++){
          peak[polInd][dir].fContainer = const_cast<AnitaEventSummary*>(this); // Set non-persistent pointer to container in hacky fashion.
          coherent[polInd][dir].fContainer = const_cast<AnitaEventSummary*>(this); // Set non-persistent pointer to container in hacky fashion.
          coherent_filtered[polInd][dir].fContainer = const_cast<AnitaEventSummary*>(this); // Set non-persistent pointer to container in hacky fashion.                    
          deconvolved[polInd][dir].fContainer = const_cast<AnitaEventSummary*>(this); // Set non-persistent pointer to container in hacky fashion.
          deconvolved_filtered[polInd][dir].fContainer = const_cast<AnitaEventSummary*>(this); // Set non-persistent pointer to container in hacky fashion.
        }
      }
      fLastEventNumber=eventNumber;
    }
  }

  void AnitaEventSummary::setThresholdForMostImpulsive(double threshold)
  {
    impulsivityFractionThreshold = threshold;
  }