MoonShadowAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
//// Class:       MoonShadowAna
//// Module Type: analyzer
//// File:        MoonShadowAna_module.cc
////
//////////////////////////////////////////////////////////////////////////


#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/FindOneP.h"

#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"

#include <art/Framework/Services/Registry/ServiceHandle.h>
#include "art_root_io/TFileService.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"

#include <iostream>
#include <set>
#include <cmath>
#include <algorithm>
#include <boost/math/special_functions.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include "TNtuple.h"

// To convert TDC to NOvA time  
#include <NovaTimingUtilities/TimingUtilities.h>

#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "nusimdata/SimulationBase/MCNeutrino.h"
#include "MEFinder/MEClusters.h"
#include "MCCheater/BackTracker.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Track.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Vertex.h"
#include "Calibrator/Calibrator.h"
#include "Utilities/func/MathUtil.h"
#include "Simulation/FLSHit.h"
#include "Simulation/FLSHitList.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleNavigator.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include <messagefacility/MessageLogger/MessageLogger.h>

#include "CelestialLocator/CelestialLocator.h"

#include <TH2D.h>
#include <iostream>
#include <stdlib.h>
#include <fstream>

namespace moonshadowana {
  class MoonShadowAna;
  
  const double PLANE_WIDTH               =  6.6522; // cm
      const double CELL_WIDTH                =  3.9887; // cm
      const double TDC_TICK                  = 15.625;  // ns
      const double SPEED_OF_LIGHT            = 30.0;    // cm / ns

      const unsigned N_CELLS                 = 384; // cells
      const unsigned N_PLANES                = 896; // planes
}

class moonshadowana::MoonShadowAna : public art::EDAnalyzer {
public:
  explicit MoonShadowAna(fhicl::ParameterSet const & p);
  virtual ~MoonShadowAna();

  void analyze(art::Event const & e) override;

  void beginJob() override;

  TVector3 vec(double zen, double azi);
  double GetDot(double zen1, double azi1, double zen2, double azi2) const;
  TVector3 AnglesToVector(double zen, double azi) const;
  void InitLookupTable();
  double GetAngularPrescale(double zen1, double azi1, double zen2, double azi2);
  bool GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2);
  double GetRate(double zen, double azi);
  double GetSmartPrescale(double zen1, double azi1, double zen2, double azi2);
  double GetSmartPrescale(double zen, double azi);
  bool GetSmartPrescaleDecision(double prescale);
  void GetHisto(std::string filename, int xmax, int ymax, TH2D* hist);
  std::string FindSmartPrescaleLookup() const;
  
      void WriteEventInfo(const art::Event& event);
    void WriteTrackInfo(const rb::Track& track, long timeStart);
    void WriteHitListInfo(const art::PtrVector<rb::CellHit> hlist, long timeStart, TVector3 start, TVector3 end);
    int GetX(art::Ptr<rb::CellHit> h, TVector3 start, TVector3 end);
    int GetY(art::Ptr<rb::CellHit> h, TVector3 start, TVector3 end);
    TVector3 GetDetVec(TVector3 vec) const;

private:
  TNtuple* fNtuple;
  std::string fTrackModuleLabel;
  art::ServiceHandle<locator::CelestialLocator> loc;
  // Lookup table for introducing "smart prescale"
  TH2D* rate_collected_angular;
  int bin_size;
  
  // cut on angle
  const double fDotCut = cos(5*M_PI/180);
  // Limit on what rate is safe for detector hardware
  double max_rate;
  // We will cut more tracks on the border of 10 degree window than in the middle.
  double angular_stdev;
  // file name of lookup table
 std::string fLookupFile;
 
  std::ofstream fOutputEvd;
 bool init;
  
};

moonshadowana::MoonShadowAna::MoonShadowAna(fhicl::ParameterSet const & p)
  : EDAnalyzer           (p)
  , fTrackModuleLabel    (p.get<std::string>("TrackModuleLabel")),
    max_rate(p.get<double>("AvgRate")),
    angular_stdev(p.get<double>("AngularSuppression")),
    fLookupFile(FindSmartPrescaleLookup())
{
  InitLookupTable();
   bin_size = 1;  
   
   init = false;
   
   //fOutputEvd.open("OFFEvd_Data.txt");
}

moonshadowana::MoonShadowAna::~MoonShadowAna()
{
fOutputEvd.close();
}

void moonshadowana::MoonShadowAna::beginJob()
{
  art::ServiceHandle<art::TFileService> tfs;
  fNtuple = tfs->make<TNtuple>("moon_ntuple", "Moon Shadow Ntuple", "zen_moon:azi_moon:zen_trk:azi_trk:nhit:len:dot:smart_pre:ang_pre");
}

  // write event info
  void moonshadowana::MoonShadowAna::WriteEventInfo(const art::Event& event)
{
        if(!init)
        {
                fOutputEvd.open(std::string("OFFEvd_Data_") + std::string(std::to_string(event.run())) + std::string("_") + std::string(std::to_string(event.subRun())) + std::string(".txt"));
                init = true;
        }

        // Event info
        fOutputEvd << "E " << event.event() << " ";
        fOutputEvd << event.time().timeHigh() << std::endl;
        
        // Tracks       
        art::Handle<std::vector<rb::Track> > tracks;
        event.getByLabel(fTrackModuleLabel, "", tracks);
        
        for(unsigned int trkIdx = 0; trkIdx < tracks->size(); ++trkIdx)
        {
                const rb::Track& trk = (*tracks)[trkIdx];
                if(trk.Is2D())continue;
                
                
                
                WriteTrackInfo(trk, event.time().timeHigh());
        }
}

void moonshadowana::MoonShadowAna::WriteTrackInfo(const rb::Track& track, long timeStart)
{
        TVector3 start = GetDetVec(track.Start());
        TVector3 end = GetDetVec(track.Stop());
        if(start.Y() < end.Y()) std::swap(start, end); // all tracks are down-going 
      
        //Check whether 3D track vector makes sense
        if(
        (start.X() == 0 && end.X() == 0) || 
        (start.Y() == 0 && end.Y() == 0) || 
        (start.Z() == 0 && end.Z() == 0))
        {
                return;
        }

        // Track info
        fOutputEvd << "T ";
        fOutputEvd << "(" << start.X() << "," << start.Y() << "," << start.Z() << ") ";
        fOutputEvd << "(" << end.X() << "," << end.Y() << "," << end.Z() << ")";
        fOutputEvd << std::endl;
        
        art::PtrVector<rb::CellHit> hlist = track.AllCells();
        
        WriteHitListInfo(hlist, timeStart, start, end);

}
void moonshadowana::MoonShadowAna::WriteHitListInfo(const art::PtrVector<rb::CellHit> hlist, long timeStart, TVector3 start, TVector3 end)
{
        for (size_t j=0; j < hlist.size(); ++j) {
                art::Ptr<rb::CellHit> h = hlist.at(j);

                if(j != 0)
                        fOutputEvd << ",";
                        
                int view = (int)h->View();
                fOutputEvd << ((long)h->TDC())/64 << "/";
                fOutputEvd << (view == 0 ? (int)h->Cell() : GetX(h, start, end)) << "/";
                fOutputEvd << (view == 1 ? (int)h->Cell() : GetY(h, start, end)) << "/";
                fOutputEvd << (int)h->Plane() << "/";
                fOutputEvd << view + 1;
        }
        
        fOutputEvd << std::endl;
}
int moonshadowana::MoonShadowAna::GetX(art::Ptr<rb::CellHit> h, TVector3 start, TVector3 end)
{
        int z = h->Plane();
        
        double xSpread = end.X() - start.X();
        double zSpread = end.Z() - start.Z();
        
        double zSpreadLoc = z - start.Z();
        
        double x = start.X() + xSpread * zSpreadLoc/zSpread;
        
        return round(x);
}
int moonshadowana::MoonShadowAna::GetY(art::Ptr<rb::CellHit> h, TVector3 start, TVector3 end)
{
        int z = h->Plane();
        
        double ySpread = end.Y() - start.Y();
        double zSpread = end.Z() - start.Z();
        
        double zSpreadLoc = z - start.Z();
        
        double y = start.Y() + ySpread * zSpreadLoc/zSpread;
        
        return round(y);
}

TVector3 moonshadowana::MoonShadowAna::GetDetVec(TVector3 vec) const
{
        int x = int((vec.X() + CELL_WIDTH * N_CELLS /2) / CELL_WIDTH);
        int y = int((vec.Y() + CELL_WIDTH * N_CELLS /2) / CELL_WIDTH);
        int z = int(vec.Z() / PLANE_WIDTH);
    
        return TVector3(x,y,z);
}

void moonshadowana::MoonShadowAna::analyze(art::Event const & e){
  
  if(e.event() == 486130)
        WriteEventInfo(e);
  
  //fOutputEvd << "E " << e.event() << std::endl;
  art::Timestamp ts = e.time();
  time_t timeSec    = ts.timeHigh();
  // std::cout<<"timeSec "<<timeSec<<std::endl;

  double zen_moon, azi_moon,  zen_sun, azi_sun;
  loc->GetMoonPosition_FD(timeSec,zen_moon, azi_moon);
  loc->GetSunPosition_FD(timeSec,zen_sun, azi_sun);

  art::Handle<std::vector<rb::Track> > tracks;
  e.getByLabel(fTrackModuleLabel, "", tracks);

  //  float vars[10];
  float vars[9];
  vars[0] =  zen_moon; // zen_moon
  vars[1] =  azi_moon; // azi_moon

  double smartPrescale = GetSmartPrescale(zen_moon, azi_moon, zen_sun, azi_sun);
  
  double zen_trk, azi_trk, cos_trk_moon, cos_trk_sun;
  for(unsigned int trkIdx = 0; trkIdx < tracks->size(); ++trkIdx){
    const rb::Track& trk = (*tracks)[trkIdx];
    if(trk.Is2D())continue;
    
    //const TVector3 dir = (trk.Start().Y() > trk.Stop().Y()) ? -trk.Dir() : trk.StopDir();
    
    TVector3 start = trk.Start();
    TVector3 end = trk.Stop();
    if(start.Y() < end.Y()) std::swap(start, end);
    const TVector3 dir = (start-end).Unit();
    
      // test print out params of the track
      //std::cout << "Track #" << trkIdx << std::endl;
      //std::cout << start.X() << "," << start.Y() << "," << start.Z() << " || ";
      //std::cout << end.X() << "," << end.Y() << "," << end.Z() << std::endl;
    
    zen_trk = acos(dir.Y()) * 180 / M_PI;      
    azi_trk = atan2(-dir.X(), dir.Z()) * 180 / M_PI + 332.066111;
    if(azi_trk < 0)   azi_trk += 360;
    if(azi_trk > 360) azi_trk -= 360;
    
    cos_trk_moon = vec(zen_moon, azi_moon).Dot(vec(zen_trk,azi_trk));
    cos_trk_sun  = vec(zen_sun,  azi_sun).Dot(vec(zen_trk,azi_trk));
    (void)cos_trk_sun; // Suppress unused variable warning

    vars[2] = zen_trk;           // zen_trk
    vars[3] = azi_trk;           // azi_trk
    vars[4] = trk.NCell();       // nhit
    vars[5] = trk.TotalLength(); // len
    vars[6] = cos_trk_moon;      // cosine angle trk and moon
    
    double angularPrescale = GetAngularPrescale(zen_moon, azi_moon, zen_trk, azi_trk);
    
    vars[7] = smartPrescale;
    vars[8] = angularPrescale;
    
    /*if(true)//cos_trk_moon > cos(5*3.14159/180))
    {
        fOutputEvd << "T " << "(" << start.X() << "," << start.Y() << "," << start.Z() << ") ";
        fOutputEvd << "(" << end.X() << "," << end.Y() << "," << end.Z() << ") ";
        fOutputEvd << std::endl;
    }*/
    
    fNtuple->Fill(vars);
  }
    
}

TVector3 moonshadowana::MoonShadowAna::vec(double zen, double azi){
  zen *= M_PI / 180;
  azi *= M_PI / 180;
  return TVector3(sin(zen)*cos(azi), sin(zen)*sin(azi), cos(zen));      
}

 double moonshadowana::MoonShadowAna::GetRate(double zen, double azi)
  {
        if(zen > 90)
                return 0;
        // Smart prescale is ratio of maximal rate to expected rate divided by angular prescale factor for given position
        int bin_azi = round(azi/bin_size);
        int bin_zen = round(zen/bin_size);
        
        //return 0;
        return rate_collected_angular->GetBinContent(bin_azi, bin_zen);
  }

 double moonshadowana::MoonShadowAna::GetAngularPrescale(double zen1, double azi1, double zen2, double azi2)
  {
    double dTheta = acos(GetDot(zen1, azi1, zen2, azi2)) * 180 / M_PI;
    return exp(-dTheta*dTheta/(angular_stdev*angular_stdev));
  }
  bool moonshadowana::MoonShadowAna::GetAngularPrescaleDecision(double zen1, double azi1, double zen2, double azi2)
  {
      double prescale = GetAngularPrescale(zen1, azi1, zen2, azi2);
      
      double r = (double)rand() / (double)RAND_MAX;
      
      return r < prescale;
  }

void moonshadowana::MoonShadowAna::InitLookupTable()
  {
     // Init lookup table and make 2d hist with expected rate
     // Get histogram with full rates in 2degx2deg bins
     int azi_max = 360;
     int zen_max = 90;
     
     rate_collected_angular = new TH2D("LookupPrescale", "", azi_max, 0, 360, zen_max, 0, 90);
     
     GetHisto(fLookupFile, azi_max, zen_max, rate_collected_angular);
  }

void moonshadowana::MoonShadowAna::GetHisto(std::string filename, int xmax, int ymax, TH2D* hist)
{
  std::ifstream data;
  data.open(filename.c_str());
  for(int j=0; j<ymax;j++)
  {
    for(int i=0; i<xmax; i++)
    {
      double val;
      data >> val;
      hist->SetBinContent(i,j, val);
    }
  }
 art::ServiceHandle<art::TFileService> tfs;
  TH2D* htmp = tfs->make<TH2D>("LookupPrescale_stored", "", 360, 0, 360, 90, 0, 90);
  htmp = (TH2D*)hist->Clone();
  htmp->Write();
  data.close();
}

double moonshadowana::MoonShadowAna::GetSmartPrescale(double zen1, double azi1, double zen2, double azi2)
  {
        double rate1 = GetRate(zen1, azi1);
        double rate2 = GetRate(zen2, azi2);
                
        if(rate1 == 0 && rate2 == 0)
                return 0; 
        
        return max_rate / (rate1 + rate2);
  }
  double moonshadowana::MoonShadowAna::GetSmartPrescale(double zen, double azi)
  {
        double rate = GetRate(zen, azi);
        
        if(rate == 0)
                return 0;
        
        return max_rate / rate;
  }

bool moonshadowana::MoonShadowAna::GetSmartPrescaleDecision(double prescale)
  {
        // Is it good enough for small prescale factors ~1/100?
        double r = (double)rand() / (double)RAND_MAX;
        return r < prescale;
  }

double moonshadowana::MoonShadowAna::GetDot(double zen1, double azi1, double zen2, double azi2) const
  {
    if(zen1 == zen2 && azi1 == azi2)
      return 1;

    zen1 *= M_PI / 180;
    zen2 *= M_PI / 180;
    azi1 *= M_PI / 180;
    azi2 *= M_PI / 180;
    double dot = sin(zen1)*cos(azi1)*sin(zen2)*cos(azi2) + sin(zen1)*sin(azi1)*sin(zen2)*sin(azi2) + cos(zen1)*cos(zen2);
    
    return dot;
  }

// find lookup table file
std::string moonshadowana::MoonShadowAna::FindSmartPrescaleLookup() const
{
  cet::search_path sp("FW_SEARCH_PATH");
  std::string lookupFname;
  sp.find_file(("Eval/PrescaleLookup.txt"), lookupFname);
  assert(!lookupFname.empty());
  return lookupFname;
}
      
DEFINE_ART_MODULE(moonshadowana::MoonShadowAna)