DDTEvd_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       CosmicUniformity
// Module Type: analyzer
// File:        CosmicUniformity_module.cc
//
// Generated at Aug 12 2015 by Matthew Tamsett
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.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 "canvas/Persistency/Common/FindOneP.h"
#include "canvas/Utilities/InputTag.h"

#include "fhiclcpp/ParameterSet.h"


#include "DDTBaseDataProducts/BaseProducts.h"
#include "DDTBaseDataProducts/CompareDAQHit.h"
#include "DDTBaseDataProducts/DAQHit.h"
#include "DDTBaseDataProducts/HitList.h"
#include "DDTBaseDataProducts/TriggerDecision.h"
#include "DDTBaseDataProducts/Track3D.h"
#include "DDTBaseDataProducts/EventHeader.h"

#include "art/Framework/Services/Optional/TFileService.h"
#include <iomanip>
#include <stdlib.h>
#include "TH1.h"
#include "TFile.h"
#include <fstream>
#include <iostream>

// To convert TDC to NOvA time
#include <NovaTimingUtilities/TimingUtilities.h>
// To convert vectors to lenght
#include "SlowMonopoleTrigger/Constants.h"

//---------------------------------------------------------------------
namespace novaddt {
    class DDTEvd;
}
//---------------------------------------------------------------------
class novaddt::DDTEvd : public art::EDAnalyzer {
    public:
        explicit DDTEvd(fhicl::ParameterSet const & p);
        virtual ~DDTEvd();
        void analyze(art::Event const & e) override;
        void endJob() override;
        
    private:
        // Get 3D track
        std::string fTrack3DLabel;
        std::string fInputHitLabel;
        std::string fInputHitInstance;
        
        // Handle file
        ofstream fOutputEvd;
        
        // Funcitons
        int GetX(DAQHit& h, TVector3 start, TVector3 end);
        int GetY(DAQHit& h, TVector3 start, TVector3 end);
        void WriteEventInfo(const art::Event& event);
        void WriteTrackInfo(const Track3D& track, long timeStart, const HitList& hitlist);
        void WriteHitListInfo(const HitList& hitlist, long timeStart, TVector3 start, TVector3 end);
};
//---------------------------------------------------------------------
novaddt::DDTEvd::DDTEvd(fhicl::ParameterSet const & p)
  : EDAnalyzer(p),
  fTrack3DLabel(p.get<std::string>("Track3DLabel")),
  fInputHitLabel(p.get<std::string>("InputHitLabel")),
  fInputHitInstance(p.get<std::string>("InputHitInstance")) 
{
    std::cout << "--- novaddt::DDTEvd instantiate" << std::endl;
    fOutputEvd.open("DDTEvd_Data.txt");
}
//---------------------------------------------------------------------
novaddt::DDTEvd::~DDTEvd()
{
    // Clean up dynamic memory and other resources here.
}


void novaddt::DDTEvd::WriteEventInfo(const art::Event& event)
{
        // Event info
        art::Handle<novaddt::EventHeader> header;
        event.getByLabel(fInputHitLabel, fInputHitInstance, header);
        fOutputEvd << "E " << header->eventNumber << " ";
        fOutputEvd << header->timeStart/novadaq::timeutils::NOVA_TIME_FACTOR + novadaq::timeutils::NOVA_EPOCH << std::endl;
        
        // Tracks
        art::Handle<std::vector<Track3D>> tracks;
        event.getByLabel(fTrack3DLabel, tracks);

        // We need the hitlist just so we can get the times for the trigger
        art::FindOneP<HitList> hls(tracks, event, fTrack3DLabel);
        
        int i = 0;
    
        // Write info for this event to file
    
    
        for(const Track3D& track: *tracks)
        {
                const HitList& hlist = *hls.at(i);
                i++;
                
                WriteTrackInfo(track, header->timeStart, hlist);
        }
}

void novaddt::DDTEvd::WriteTrackInfo(const Track3D& track, long timeStart, const HitList& hitlist)
{
        TVector3 start = track.Start();
        TVector3 end = track.End();
        if(start.Y() < end.Y()) std::swap(start, end); // all tracks are down-going 
      
        //Check whether 3D track vector makes sense
        if(!track.Is3D() ||
        (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;
        
        WriteHitListInfo(hitlist, timeStart, start, end);

}
void novaddt::DDTEvd::WriteHitListInfo(const HitList& hitlist, long timeStart, TVector3 start, TVector3 end)
{
        for(uint j = 0; j < hitlist.size(); j++)
        {
                if(j != 0)
                        fOutputEvd << ",";
                        
                DAQHit h = hitlist.at(j);
                int view = (int)h.View().val;
                fOutputEvd << ((long)h.TDC().val - timeStart)/64 << "/";
                fOutputEvd << (view == 1 ? (int)h.Cell().val : GetX(h, start, end)) << "/";
                fOutputEvd << (view == 2 ? (int)h.Cell().val : GetY(h, start, end)) << "/";
                fOutputEvd << (int)h.Plane().val << "/";
                fOutputEvd << view;
        }
        
        fOutputEvd << std::endl;
}
int novaddt::DDTEvd::GetX(DAQHit& h, TVector3 start, TVector3 end)
{
        int z = h.Plane().val;
        
        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 novaddt::DDTEvd::GetY(DAQHit& h, TVector3 start, TVector3 end)
{
        int z = h.Plane().val;
        
        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);
}

//---------------------------------------------------------------------
void novaddt::DDTEvd::analyze(art::Event const & event)
{                          
        WriteEventInfo(event);
}
//---------------------------------------------------------------------
void novaddt::DDTEvd::endJob()
{
    std::cout << "=== novaddt::DDTEvd endJob  " << std::endl;
    fOutputEvd.close();
}

DEFINE_ART_MODULE(novaddt::DDTEvd)