DLECorrectionGG/DLECorrection.cc

Go to the documentation of this file.

#include "DLECorrection.h"

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

#include <sevt/SEvent.h>
#include <sevt/Station.h>
#include <sevt/PMTCalibData.h>
#include <sevt/PMTRecData.h>
#include <sevt/StationRecData.h>
#include <sevt/StationConstants.h>

#include <det/Detector.h>

#include <sdet/SDetector.h>

#include <utl/ErrorLogger.h>
#include <utl/MathConstants.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/AxialVector.h>
#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/AugerUnits.h>
#include <utl/PhysicalConstants.h>
#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>

#define max(a,b) ((a) > (b) ? (a) : (b))
#define max3(a,b,c) (max((a),max((b),(c))))

using namespace sdet;
using namespace det;
using namespace sevt;
using namespace evt;
using namespace std;
using namespace utl;
using namespace fwk;


namespace DLECorrection{
  
  const double VEM_Threshold_DLE = 0.2; //Threshold for direct light correction
  const StationConstants::SignalComponent source1 = StationConstants::eDirectLightSubtracted;
  
  //Parametrized function to remove average azimuthal dependence
  double DLECorrection::b(double theta){
    
    //Fit Result:
    const double a = 1.60484e-03;
    const double a0 = 2.00000e-05;
    const double a1 = 1.70001e-03;
    const double a2 = -1.74851e-05;
    const double a3 = -1.35945e-01;
    const double a4 = 5.79209e+01;
    const double B = a*theta + (a0 + a1*theta + a2*theta*theta)
      / (1 + exp( a3 * (theta -a4)));
    
    return B;
  }
  
  //correcting for direct light...
  //zenithal & azimuthal dependence:
  void DLECorrection::Direct_Light_Azim_Depen_PMT(int pmtId, double ang, 
                                                  double phi, 
                                                  utl::TraceD &times)
  {
    if (pmtId == 1)
      f = b(ang)*cos((phi-30.+180.)*3.14159/180.);
    if (pmtId == 2)
      f = b(ang)*cos((phi-270.+180.)*3.14159/180.);
    if (pmtId == 3)
      f = b(ang)*cos((phi-150.+180.)*3.14159/180.);
    
    for (int j = 0; j < fTraceSize; j++){
      if (times[j] > VEM_Threshold_DLE)
        times[j] -= f;
    }
  }
  //correcting individual instances of Direct Light
  void DLECorrection:: PMT2_Correction(utl::TraceD &t1,utl::TraceD &t2)
  {
    for (int j = 0; j < fTraceSize; j++){
      const double s1=t1[j];
      const double s2=t2[j];
      const double smean=(s1+s2)/2.;
      const double sigma=sqrt(((s1-smean)*(s1-smean)+(s2-smean)*(s2-smean))/2.);
      if(sigma>=1.)
        (s1>s2) ? t1[j]=t2[j] : t2[j]=t1[j] ;
    }
  }
  void DLECorrection::Individual_Direct_Light_Corr(int NPMTs,sevt::Station& station, const StationConstants::SignalComponent source)
  {
    
    if(NPMTs==2){
      
      vector<unsigned int> workingPMTsIds;
      for (sevt::Station::PMTIterator pIt = station.PMTsBegin();
           pIt != station.PMTsEnd(); ++pIt) {
        if (!pIt->HasRecData())
          continue;
        if (!pIt->GetRecData().HasVEMTrace(source))
          continue;
        workingPMTsIds.push_back(pIt->GetId());
      }
      
      if (workingPMTsIds.size() !=2 ){
        INFO("DLECorrection::Individual_Direct_Light_Corr" 
             " => Something is wrong, less than 2 PMTS ");
        return;
      }
      TraceD &times1 = station.GetPMT(workingPMTsIds[0]).GetRecData().GetVEMTrace(source);
      TraceD &times2 = station.GetPMT(workingPMTsIds[1]).GetRecData().GetVEMTrace(source);
      PMT2_Correction(times1,times2);
      
    }
    else //Correction for the case the 3 PMTs work normally
      {
        TraceD &times1 = station.GetPMT(1).GetRecData().GetVEMTrace(source);
        TraceD &times2 = station.GetPMT(2).GetRecData().GetVEMTrace(source);
        TraceD &times3 = station.GetPMT(3).GetRecData().GetVEMTrace(source);
        
        for ( int j = 0; j < fTraceSize; j++){
          double s1 = times1[j];
          double s2 = times2[j];
          double s3 = times3[j];
          double smean = (s1+s2+s3)/3;
          double sigma = sqrt(((s1-smean)*(s1-smean)+(s2-smean)*(s2-smean)+(s3-smean)*(s3-smean))/3);
          if(sigma>=1.){
            double peak = max3(s1,s2,s3);
            
            if(peak==s1)
              times1[j] = (s2+s3)/2.;
            else if(peak==s2)
              times2[j] = (s1+s3)/2.;
            else
              times3[j] = (s1+s2)/2.;
          }
        }
      }
  }
  
  
  void DLECorrection::CorrectDLE(Event & event, const StationConstants::SignalComponent component, const StationConstants::SignalComponent source){
    INFO("DLECorrection::CorrectDLE()");
    sevt::SEvent&  sevent = event.GetSEvent(); 
    const det::Detector& detector = det::Detector::GetInstance();
    const CoordinateSystemPtr siteCS = detector.GetSiteCoordinateSystem();        
    evt::ShowerSRecData& showersrec  = event.GetRecShower().GetSRecShower();
    const Vector& axis        = showersrec.GetAxis();           
    double theta=axis.GetTheta(siteCS)/deg;
    double phi=axis.GetPhi(siteCS)/deg;
    
    //Loop over stations
    for (SEvent::StationIterator sIt = sevent.StationsBegin();
         sIt != sevent.StationsEnd(); ++sIt) {
      
      sevt::Station& station = *sIt;
      
      if (!station.IsCandidate()) // do not work with bad stations
        continue;
      if (station.IsLowGainSaturation()) // do not work with saturated stations
        continue;

      const sdet::Station& dStation =
        det::Detector::GetInstance().GetSDetector().GetStation(station);
      
      fTraceSize = dStation.GetFADCTraceLength();
      
      int NumberOfPMTs=0;
      //Loop over 3 PMTs
      for (sevt::Station::PMTIterator pIt = sIt->PMTsBegin();
           pIt != sIt->PMTsEnd(); ++pIt) {
        if (!pIt->HasRecData())
          continue;
        if (pIt->GetRecData().HasVEMTrace(component))
          ++NumberOfPMTs;
      }
      
      for (sevt::Station::PMTIterator pIt = sIt->PMTsBegin();
           pIt != sIt->PMTsEnd(); ++pIt) {
        
        if (!pIt->HasRecData())
          continue;
        
        PMTRecData& pmtRec = pIt->GetRecData();
        
        if (!pmtRec.HasVEMTrace(component))
          continue;
        
        const int pmtId = pIt->GetId();
        
        const TraceD& vemtrace = pmtRec.GetVEMTrace(component);
        
        if (!pmtRec.HasVEMTrace(source))
          pmtRec.MakeVEMTrace(source);
        
        for ( int i =0 ; i < fTraceSize; ++i)
          pmtRec.GetVEMTrace(source)[i] = vemtrace[i];
        
        //Correct for azimuthal dependence of direct light
        Direct_Light_Azim_Depen_PMT(pmtId,theta,phi,pmtRec.GetVEMTrace(source));
        
        
      }//loop over the pmts
      //Correct individual instances of direct light    
      if(NumberOfPMTs>1)
        Individual_Direct_Light_Corr(NumberOfPMTs,station,source);
      
    } //Loop over stations 
    
    
  }
  
  DLECorrection::DLECorrection():f(0){
  }
  
  DLECorrection::~DLECorrection(){
  }
  
  
  VModule::ResultFlag DLECorrection::Init()
  {
    
    INFO("DLECorrection::Init()");
    
    return eSuccess;
  }
  
  
  VModule::ResultFlag DLECorrection::Run(evt::Event& event)
  {
    
    INFO("DLECorrection::Run()");
    if (!(event.HasSEvent() && event.HasRecShower() &&
          event.GetRecShower().HasSRecShower()))
      return eContinueLoop;
    
    CorrectDLE(event,StationConstants::eTotal,source1); 
    
    
    return eSuccess;
  }
  
  
  VModule::ResultFlag DLECorrection::Finish()
  {
    
    INFO("DLECorrection::Finish()");
    
    return eSuccess;
  }
  
}