SdFootprintAnalyzer.cc

Go to the documentation of this file.

#include "SdFootprintAnalyzer.h"

#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CentralConfig.h>

#include <sdet/SDetector.h>
#include <sdet/Station.h>

#include <evt/ShowerRecData.h>
#include <evt/ShowerSRecData.h>

#include <sevt/SdFootprintData.h>
#include <sevt/PMT.h>

#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>

using namespace fwk;
using namespace utl;
using namespace evt;
using namespace sevt;

using namespace SdFootprintAnalyzerNS;


ByIncreasingDistanceToAxis::ByIncreasingDistanceToAxis(const CoordinateSystemPtr& cs):
  fShowerCS(cs)
{
}


bool
ByIncreasingDistanceToAxis::operator()(const sevt::Station* const s1, const sevt::Station* const s2)
  const
{
  static const sdet::SDetector& sdetector = det::Detector::GetInstance().GetSDetector();
  // put the stations with no rec data at the end
  if (!s1->HasRecData())
    return false;
  else if (!s2->HasRecData())
    return true;
  else {
    const double d1 = sdetector.GetStation(*s1).GetPosition().GetRho(fShowerCS);
    const double d2 = sdetector.GetStation(*s2).GetPosition().GetRho(fShowerCS);
    return d1 < d2;
  }
}

bool ByIncreasingTime::operator()(SdFootprintAnalyzerNS::Station i, SdFootprintAnalyzerNS::Station j) const
{
  return (i.fTime<j.fTime);
}


SdFootprintAnalyzer::SdFootprintAnalyzer():
  fMinProjDistForSpeed(1000),
  fRecomputeBarycenter(false)
{

}


SdFootprintAnalyzer::~SdFootprintAnalyzer()
{

}


VModule::ResultFlag
SdFootprintAnalyzer::Init()
{
  INFO("SdFootprintAnalyzer::Init()");
  
  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("SdFootprintAnalyzerOG");
  
  topB.GetChild("Verbose").GetData(fVerbose);
  fVerbose=fVerbose<0 ? 0 : fVerbose;
  
  topB.GetChild("MinProjDistForSpeed").GetData(fMinProjDistForSpeed);
  fMinProjDistForSpeed=fMinProjDistForSpeed<0 ? abs(fMinProjDistForSpeed) : fMinProjDistForSpeed;
  
  topB.GetChild("RecomputeBarycenter").GetData(fRecomputeBarycenter);
  
  topB.GetChild("MaxPhiDiffLine").GetData(fMaxPhiDiffLine);
  fMaxPhiDiffLine=fMaxPhiDiffLine<0 ? abs(fMaxPhiDiffLine) : fMaxPhiDiffLine;
  
  if(fVerbose)
    cout << " Config parameters: " << endl
         << "\t                    Verbose: " << fVerbose << endl
         << "\t  Min Proj. Dist. for Speed: " << fMinProjDistForSpeed << endl
         << "\t       Recompute Barycenter: " << fRecomputeBarycenter << endl
         << "\t    Max. Phi Diff. for Line: " << fMaxPhiDiffLine << endl
         << endl << endl;

  return eSuccess;
}


VModule::ResultFlag
SdFootprintAnalyzer::Run(Event & event)
{
  INFO("SdFootprintAnalyzer::Run()");

  if(!event.HasSEvent()){
    INFO("No sd event. Skipping.");
    return eSuccess;
  }
  SEvent & theSEvent = event.GetSEvent();
  
  if(!event.HasRecShower())
    event.MakeRecShower();
  
  if(!event.GetRecShower().HasSRecShower()){
    event.GetRecShower().MakeSRecShower();
    fRecomputeBarycenter=true; 
    if(fVerbose) cout << " Barycenter info not available. Will be recomputed." << endl;
  }
  
  ShowerSRecData & theSRecShower = event.GetRecShower().GetSRecShower();
  
  if(!theSRecShower.HasFootprintData())
    theSRecShower.MakeFootprintData();
  
  SdFootprintData & theFootprint = theSRecShower.GetFootprintData();
  
  // If not exist, compute the Barycenter
  if(fRecomputeBarycenter) theSRecShower.SetBarycenter(ComputeBarycenter(theSEvent));

  // Get collection of valid ground stations
  StationCollection AllStations=GetFootprintStations(theSEvent);
  
  // Compute Ground Variables
  ComputeGroundVariables(event, AllStations, theFootprint);
  if(fVerbose) cout << "   Footprint Variables: " << endl
                    << "               Length : " << theFootprint.GetLength() << endl
                    << "                Width : " << theFootprint.GetWidth() << endl
                    << "         Ground Speed : " << theFootprint.GetSpeed() << endl
                    << "   Speed Standard Dev.: " << theFootprint.GetSpeedStandardDeviation() << endl
                    << "            Alignment : " << theFootprint.GetAlignment() << endl;
  
  
  // Compute AoP Asymmetry
  theFootprint.SetAreaOverPeakAsymmetry(ComputeAoPAsymmetry(AllStations));
  if(fVerbose) cout << "  Area over Peak Asymm: " << theFootprint.GetAreaOverPeakAsymmetry() << endl << endl;
                 
  return eSuccess;
}


VModule::ResultFlag
SdFootprintAnalyzer::Finish()
{
  
  return eSuccess;
}


Point
SdFootprintAnalyzer::ComputeBarycenter(const SEvent& theSEvent)
  const
{
  static const sdet::SDetector& theSDet =
    det::Detector::GetInstance().GetSDetector();
  static CoordinateSystemPtr siteCs =
    det::Detector::GetInstance().GetSiteCoordinateSystem();
  
  Point origin = Point(0,0,0,siteCs);
  Vector barycenter = Vector(0,0,0,siteCs);
  double wSum = 0;
  
  for (SEvent::ConstCandidateStationIterator sIt = theSEvent.CandidateStationsBegin();
       sIt != theSEvent.CandidateStationsEnd();++sIt){
    const sdet::Station& detStation = theSDet.GetStation(sIt->GetId());
    if (!sIt->IsCandidate() || !detStation.IsInGrid()) // Should we think on Infill studies?
      continue;
    barycenter += sIt->GetRecData().GetTotalSignal()*(detStation.GetPosition()-origin);
    wSum += sIt->GetRecData().GetTotalSignal();
  }

  if(fVerbose) cout << "  Computed Barycenter (X, Y, Z): (" 
                    << (origin + barycenter/wSum).GetX(siteCs) << " , "
                    << (origin + barycenter/wSum).GetY(siteCs) << " , "
                    << (origin + barycenter/wSum).GetZ(siteCs) << ")"
                    << endl;
  
  return origin + barycenter/wSum;
}


void
SdFootprintAnalyzer::SortStationsAroundBarycenter(Event& event, Point barycenter)
  const
{
  CoordinateSystemPtr baryCS = LocalCoordinateSystem::Create(barycenter);
  ByIncreasingDistanceToAxis criterion(baryCS);
  event.GetSEvent().SortStations(criterion);
}

StationCollection
SdFootprintAnalyzer::GetFootprintStations(SEvent& theSEvent)
  const
{
  static const sdet::SDetector& theSDet =
    det::Detector::GetInstance().GetSDetector();
  
  StationCollection AllStations;
  
  if(fVerbose) cout << "  Selected Station Id (Id, GPSNs, Signal, AoP, TOT) " << endl;
  
  // Fill candidate Stations
  for (SEvent::CandidateStationIterator sIt = theSEvent.CandidateStationsBegin();
       sIt != theSEvent.CandidateStationsEnd();
       ++sIt) {
    
    if(!sIt->HasRecData()){
      if(fVerbose) cerr << "  Warning! Candidate Station without Rec Data. Skiping..." << endl;
      continue;
    }  
    
    // Compute AoP
    if(!sIt->GetRecData().GetAreaOverPeak()){
      double AoP=0;
      double Peak=0;
      utl::TraceD vemPeakTrace = sIt->GetVEMTrace();
      utl::TraceD vemChargeTrace = GetVEMChargeTrace(*sIt);
      unsigned int Start = sIt->GetRecData().GetSignalStartSlot();
      unsigned int Stop = sIt->GetRecData().GetSignalEndSlot();
      for(unsigned int slot=Start; slot<Stop; slot++){
        AoP += vemChargeTrace[slot];
        Peak = Peak > vemPeakTrace[slot] ? Peak : vemPeakTrace[slot];
      }
      AoP = Peak != 0 ? AoP/Peak : 0;
      sIt->GetRecData().SetAreaOverPeak(AoP);
    }

    if(!sIt->HasTriggerData()){
      if(fVerbose) cerr << "  Warning! Candidate Station without Trigger Data. Skiping..." << endl;
      continue;
    }  
    
    int Id=sIt->GetId();
    SdFootprintAnalyzerNS::Station thisSt(Id,
                                          sIt->GetRecData().GetSignalStartTime(),
                                          theSDet.GetStation(Id).GetPosition(),
                                          sIt->GetRecData().GetTotalSignal(),
                                          sIt->GetRecData().GetAreaOverPeak(),
                                          sIt->GetTriggerData().IsTimeOverThreshold(),
                                          0);

    // cout << "           PMT AoP: (PMT1, PMT2, PMT3, Mean) ="
    //   << " ( " << sIt->GetPMT(1).GetRecData().GetAreaOverPeak()
    //   << " , " << sIt->GetPMT(2).GetRecData().GetAreaOverPeak()
    //   << " , " << sIt->GetPMT(3).GetRecData().GetAreaOverPeak()
    //   << " , " << sIt->GetPMT(1).GetRecData().GetAreaOverPeak()/3.
    //   + sIt->GetPMT(2).GetRecData().GetAreaOverPeak()/3. + sIt->GetPMT(3).GetRecData().GetAreaOverPeak()/3. << " )" << endl;
    
    if(fVerbose) cout << "                   (" << thisSt.fId
                      << " , " << thisSt.fTime.GetGPSNanoSecond()
                      << " , " << thisSt.fSignal
                      << " , " << thisSt.fAoP
                      << " , " << thisSt.fTOT
                      << ")" << endl;
    
    AllStations.push_back(thisSt);
  }
  return AllStations;
}


void
SdFootprintAnalyzer::ComputeGroundVariables(Event& event, StationCollection station, SdFootprintData& theFootprint)
{
  
  utl::CoordinateSystemPtr fBaryCS = LocalCoordinateSystem::Create(event.GetRecShower().GetSRecShower().GetBarycenter());
  
  
  
  double wSum = 0;
  double Ixx = 0;
  double Iyy = 0;
  double Ixy = 0;
  double ToTFrac=0;
  int NStations=station.size();
  
  if(fVerbose>3) cout << "    Computing Length and Width for " << NStations << endl;
  
  for (int i = 0; i<NStations; i++) {
    
    double wi = station[i].fSignal;
    wSum+=wi;
    
    Point StPos = station[i].fPos;
    double xDiff = StPos.GetX(fBaryCS);
    Ixx += wi*xDiff*xDiff;
    double yDiff = StPos.GetY(fBaryCS);
    Iyy += wi*yDiff*yDiff;
    
    Ixy += wi*xDiff*yDiff;
    
    ToTFrac+=station[i].fTOT/(double)NStations;
  }
  
  if(fVerbose>4) cout << "     Computed TOT Fraction " << ToTFrac << endl;
  
  theFootprint.SetTOTFraction(ToTFrac);
  
  Ixx /= wSum;
  Iyy /= wSum;
  Ixy /= wSum;
  
  double lambda1 = (Ixx + Iyy + sqrt(pow(Ixx-Iyy,2.)+4.*Ixy*Ixy)) /2.; //First eigenvalue of the covariance matrix
  double lambda2 = (Ixx + Iyy - sqrt(pow(Ixx-Iyy,2.)+4.*Ixy*Ixy)) /2.; //Second eigenvalue of the covariance matrix
  double l = sqrt(lambda1/wSum);// As defined in tau limit long paper
  double w = sqrt(lambda2/wSum);// As defined in tau limit long paper
  
  // First eigenvector e1 (1,(lambda1-Ixx)/Ixy) This vector has the direction of length
  double q  = (lambda1-Ixx)/Ixy; //Second coordinate of first eigenvector, first coordinate is 1
  double q2 = q*q;
  // First normalized eigenvector v
  double v1 = 1./sqrt(1.+ q2); // Normalized first coordinate of first eigenvector
  double v2 = q/sqrt(1.+ q2); // Normalized second coordinate of first eigenvector
  
  // Second eigenvector e2 (1,(lambda2-Ixx)/Ixy)
  // double r  = (lambda2-Ixx)/Ixy; //Second coordinate of eigenvector, first coordinate is 1
  // double r2 = q*q;
  // Second normalized eigenvector w
  // double w1 = 1./sqrt(1.+ q2); // Normalized first coordinate of second eigenvector
  // double w2 = q/sqrt(1.+ q2); // Normalized second coordinate of second eigenvector
  
  if(fVerbose>3) cout << "     Momentum matrix of footprint Stations: " << endl
                      << "       wSum : " << wSum  << endl
                      << "        Ixx : " << Ixx << endl
                      << "        Iyy : " << Iyy << endl
                      << "        Ixy : " << Ixy << endl
                      << "    lambda1 : " << lambda1 << endl
                      << "    lambda2 : " << lambda2 << endl
                      << "          q : " << q << endl
                      << "         q2 : " << q2 << endl
                      << "         v1 : " << v1 << endl
                      << "         v2 : " << v2 << endl
                      << "     Length : " << l << endl
                      << "      Width : " << w << endl
                      << endl;
  
  theFootprint.SetWidth(w);
  theFootprint.SetLength(l);
  // theFootprint.SetElliPhi=acos(v1);
  
  
  for (int i = 0; i<NStations; i++) {
    Point StPosi = station[i].fPos;
    station[i].fLp=v1*StPosi.GetX(fBaryCS) + v2*StPosi.GetY(fBaryCS);
  }
  
  double speSum=0;
  double spe2Sum=0;
  double speN=0;
  wSum=0;
  
  for (int i=0;i<NStations;i++){
    TimeStamp ti = station[i].fTime;
    for (int j=0;j<i;j++) {
      double dLPij = station[i].fLp-station[j].fLp;
      if (abs(dLPij) < fMinProjDistForSpeed)
        continue;
      
      TimeStamp tj = station[j].fTime;
      double dTij  = (tj - ti).GetInterval();
      double speedij=dLPij/dTij;
      double w=station[i].fSignal*station[j].fSignal;
      speSum+=speedij*w;
      spe2Sum+=speedij*speedij*w;
      wSum+=w;
      speN++;
    }
  }
  
  double speMean=speSum/wSum;
  double spe2Mean=spe2Sum/wSum;
  double Sspe=spe2Mean-speMean*speMean;

  if(fVerbose>3) cout << "  Variables for speed computing: " << endl
                      << "               Weight : " << wSum  << endl
                      << "            Speed Sum : " << speSum  << endl
                      << "     Square Speed Sum : " << spe2Sum << endl
                      << "                 Sspe : " << Sspe << endl
                      << endl;
  
  
  theFootprint.SetSpeed(abs(speMean), Sspe);
  

  double phiDeg=acos(v1)/deg;

  theFootprint.SetAlignment(SdFootprintEnumerations::eAligned);
  
  for (int i=0;i<NStations;i++){
    double xi = station[i].fPos.GetX(fBaryCS);
    double yi = station[i].fPos.GetY(fBaryCS);
    TimeStamp ti = station[i].fTime;
    
    for (int j=0;j<NStations;j++){
      if(i==j) continue;
      double xj = station[j].fPos.GetX(fBaryCS);
      double yj = station[j].fPos.GetY(fBaryCS);
      TimeStamp tj = station[j].fTime;
      
      double dXij = xj - xi;
      double dYij = yj - yi;
      double dTij = (tj - ti).GetInterval();
      
      double vXij = dXij/dTij;
      double vYij = dYij/dTij;
      
      double Phi = 180. + atan2(vYij,vXij)/deg;
      if(Phi>180) Phi = Phi-360;
      
      double alfa = fabs(phiDeg - Phi);
      if(alfa > 360 - fMaxPhiDiffLine) alfa -= 360;
      if ( alfa > fMaxPhiDiffLine){
        if(fVerbose>3) cout << "    Station " << station[i].fId << " and " << station[j].fId  << " disaligned with angle " << alfa << endl;
        if (theFootprint.GetAlignment() == SdFootprintEnumerations::eIll) {
          theFootprint.SetAlignment(SdFootprintEnumerations::eNotAligned);
          return;
        }
        else 
          theFootprint.SetAlignment(SdFootprintEnumerations::eIll);
      }
    }
  }
  
  return;
}


double
SdFootprintAnalyzer::ComputeAoPAsymmetry(StationCollection stations)
{
  double AoPAsym=0; 
  sort(stations.begin(), stations.end(), SdFootprintAnalyzerNS::ByIncreasingTime());  
  
  double AoPEarly=0;
  double AoPLate=0;
  unsigned int st=0;
  
  int NStations = stations.size();
  if(NStations==4 || NStations==5)
    st=2;
  else if(NStations==6 || NStations==7)
    st=3;
  else if(NStations>7)
    st=4;
  
  for(unsigned int i=0; i<st; i++){
    AoPEarly+=stations[i].fAoP;
    AoPLate+=stations[NStations-1-st].fAoP;
  }
  
  AoPAsym=AoPEarly-AoPLate;
  
  return AoPAsym;
}


TraceD 
SdFootprintAnalyzer::GetVEMChargeTrace(const sevt::Station& station)
  const
{
  vector<const PMT*> validPMTs;

  const sdet::Station& dStation =
    det::Detector::GetInstance().GetSDetector().GetStation(station);
  
  TraceD vemTrace(dStation.GetFADCTraceLength(),dStation.GetFADCBinSize(), 
                  sevt::StationConstants::eMySource0);
  
  const int traceLength = vemTrace.GetSize();
  for (int i = 0; i < traceLength; ++i)
    vemTrace[i]=0;
  
  for (sevt::Station::ConstPMTIterator pmtIt = station.PMTsBegin();
       pmtIt != station.PMTsEnd(); ++pmtIt)
    if (pmtIt->HasCalibData() && pmtIt->GetCalibData().IsTubeOk() && 
        pmtIt->HasRecData() && pmtIt->GetRecData().HasVEMTrace() )
      validPMTs.push_back(&(*pmtIt));
  
  
  if (!validPMTs.empty()) {
    const int nPMTs = validPMTs.size();
    for (vector<const PMT*>::const_iterator pmtIt = validPMTs.begin();
         pmtIt != validPMTs.end(); ++pmtIt) {
      const PMTRecData & pmtRec=(*pmtIt)->GetRecData();
      const TraceD& pmtTrace = pmtRec.GetVEMTrace();
      const double PeakChargeFactor = pmtRec.GetVEMPeak()/pmtRec.GetVEMCharge(); 
      for (int i = 0; i < traceLength; ++i)
        vemTrace[i] += (pmtTrace[i] * PeakChargeFactor) / nPMTs;
    }
  }
  
  return vemTrace;
}