SimTestPulseWire_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       SimTestPulseWire
// Module Type: producer
// File:        SimTestPulseWire_module.cc
//
// Imported to ICARUS on October 21, 2018
//
// Generated at Sat Feb 11 04:35:44 2017 by Kazuhiro Terao using artmod
// from cetpkgsupport v1_11_00.
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDProducer.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 "art_root_io/TFileService.h"

#include "lardata/DetectorInfoServices/DetectorClocksService.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardataobj/Simulation/SimEnergyDeposit.h"
#include "fhiclcpp/ParameterSet.h"
#include "larcoreobj/SummaryData/RunData.h"
#include <memory>

#include "lardataobj/Simulation/SimChannel.h"
#include "lardataobj/RawData/TriggerData.h"
#include "larcore/Geometry/Geometry.h"
#include "larcorealg/Geometry/GeometryCore.h"

#include "ParamHolder.h"

#include <TTree.h>
#include <TFile.h>

class SimTestPulseWire;

class SimTestPulseWire : public art::EDProducer
{
public:
    explicit SimTestPulseWire(fhicl::ParameterSet const & p);
    // The destructor generated by the compiler is fine for classes
    // without bare pointers or other resource use.
    
    // Plugins should not be copied or assigned.
    SimTestPulseWire(SimTestPulseWire const &) = delete;
    SimTestPulseWire(SimTestPulseWire &&) = delete;
    
    SimTestPulseWire & operator = (SimTestPulseWire const &) = delete;
    SimTestPulseWire & operator = (SimTestPulseWire &&) = delete;
    
    // Required functions.
    void produce(art::Event & e) override;
    void beginRun(art::Run& run) override;
    
    void beginJob() override;
    void endJob()   override;

private:

    bool fVerbose; ///< verbosity

    using WireVec        = std::vector<int>;
    using WireVecVec     = std::vector<WireVec>;

    using Position       = std::vector<double>;
    using PositionVec    = std::vector<Position>;
    using PositionVecVec = std::vector<PositionVec>;
    
    //
    // Parameters to be configured
    //
    double fTriggerTime;                      ///< Trigger timing in electronics clock [us]
    unsigned int        fCryostat;            ///< Cryostat to consider
    unsigned int        fTPC;                 ///< TPC in that cryostat (0-3 for split wire)
    std::vector<double> fSimTime_v;           ///< Charge timing in electronics clock [us]
    WireVecVec          fWiresByPlaneVec;     ///< For each plane the wires to deposit charge on
    std::vector<double> fNumElectrons_v;      ///< Charge amount in electron count
    
    //
    // Parameters to be calculated
    //
    WireVec             fTick_v;              ///< Corresponding tick
    WireVecVec          fPlaneChannelVec;     ///< Keep track of wires for each plane

    PositionVecVec      fPlanePositionVecVec; ///< Keep track of positions for each wire on each plane
    
    int _run, _subrun, _event;
    
    TTree*              fTupleTree;           ///< output analysis tree
};

SimTestPulseWire::SimTestPulseWire(fhicl::ParameterSet const & p)
  : EDProducer{p}
// Initialize member data here.
{
    produces< std::vector<sim::SimChannel> >();
    produces< std::vector<sim::SimEnergyDeposit> >();
    produces< std::vector<raw::Trigger> >();
    produces< sumdata::RunData, art::InRun >();
    
    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataForJob();
    fTriggerTime = clockData.TriggerTime();

    // Assume 3 planes
    fWiresByPlaneVec.resize(3);
    fPlaneChannelVec.resize(3);
    fPlanePositionVecVec.resize(3);

//    fTriggerTime    = p.get< double>              ("TriggerTime_us");
    fSimTime_v          = p.get< std::vector<double> >("SimTimeArray_us");
    fCryostat           = p.get< unsigned int        >("Cryostat");
    fTPC                = p.get< unsigned int        >("TPC");
    fWiresByPlaneVec[0] = p.get< WireVec             >("Plane0Wires");
    fWiresByPlaneVec[1] = p.get< WireVec             >("Plane1Wires");
    fWiresByPlaneVec[2] = p.get< WireVec             >("Plane2Wires");
    fNumElectrons_v     = p.get< std::vector<double> >("NumElectronsArray");
    fVerbose            = p.get< bool>                ("Verbose",false);
    
    assert( fSimTime_v.size() == fWiresByPlaneVec[0].size() &&
            fSimTime_v.size() == fWiresByPlaneVec[1].size() &&
            fSimTime_v.size() == fWiresByPlaneVec[2].size() &&
            fSimTime_v.size() == fNumElectrons_v.size() );
}

void SimTestPulseWire::beginJob()
{
    art::ServiceHandle<art::TFileService> tfs;

    fTupleTree = tfs->make<TTree>("SimTestPulseWire", "Tree by SimTestPulseWire_module");
    fTupleTree->Branch("run",&_run,"run/I");
    fTupleTree->Branch("subrun",&_subrun,"subrun/I");
    fTupleTree->Branch("event",&_event,"event/I");
    fTupleTree->Branch("trigger_time",&fTriggerTime,"trigger_time/D");
    fTupleTree->Branch("charge_time_v","std::vector<double>",&fSimTime_v);
    fTupleTree->Branch("tick_v","std::vector<int>",&fTick_v);
//    fTupleTree->Branch("y_v","std::vector<double>",&fY_v);
//    fTupleTree->Branch("z_v","std::vector<double>",&fZ_v);
    fTupleTree->Branch("e_v","std::vector<double>",&fNumElectrons_v);
    
    fTupleTree->Branch("ch_plane0","std::vector<int>",&fPlaneChannelVec[0]);
    fTupleTree->Branch("ch_plane1","std::vector<int>",&fPlaneChannelVec[1]);
    fTupleTree->Branch("ch_plane2","std::vector<int>",&fPlaneChannelVec[2]);
    
    fTupleTree->Branch("wire_plane0","std::vector<int>",&fWiresByPlaneVec[0]);
    fTupleTree->Branch("wire_plane1","std::vector<int>",&fWiresByPlaneVec[1]);
    fTupleTree->Branch("wire_plane2","std::vector<int>",&fWiresByPlaneVec[2]);
}

void SimTestPulseWire::endJob()
{
    alternative::ParamHolder::destroy();
}

void SimTestPulseWire::beginRun(art::Run& run)
{
  // grab the geometry object to see what geometry we are using                                                                   
    art::ServiceHandle<geo::Geometry> geo;
    
    std::unique_ptr<sumdata::RunData> runData(new sumdata::RunData(geo->DetectorName()));
    
    run.put(std::move(runData));
    
    return;
}

void SimTestPulseWire::produce(art::Event & e)
{
    std::unique_ptr<std::vector<raw::Trigger> > trigger_v(new std::vector<raw::Trigger> );
    trigger_v->push_back(raw::Trigger(0,fTriggerTime,fTriggerTime,1));
    
    std::unique_ptr<std::vector<sim::SimChannel> > simch_v(new std::vector<sim::SimChannel> );

    std::unique_ptr<std::vector<sim::SimEnergyDeposit>> simDep_v(new std::vector<sim::SimEnergyDeposit>);
    
    geo::Point_t chargeDepCoords = {0., 0., 0.};

    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService>()->DataFor(e);
    art::ServiceHandle<geo::Geometry> geo;

    auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob(clockData);

    for(size_t idx = 0; idx < fPlaneChannelVec.size(); idx++)
    {
        fPlaneChannelVec[idx].clear();
        fPlanePositionVecVec[idx].clear();
    }

    fTick_v.clear();
    
    auto& pholder = alternative::ParamHolder::get();
    pholder.Clear();

    for(size_t index=0; index < fSimTime_v.size(); ++index) {
        
        int tdc = fSimTime_v[index] / clockData.TPCClock().TickPeriod();
        fTick_v.push_back(clockData.TPCTDC2Tick(tdc));
        
        if(fVerbose)
            std::cout << "[BUFFOON!] Charge injection id " << index << " @ TDC=" << tdc
//            << " @ Y=" << fY_v[index] << " @ Z=" << fZ_v[index]
            << " with " << fNumElectrons_v[index] << " electrons " << std::endl;
        
        if(tdc<0) {
            std::cerr << "\033[93m[WARNING BUFFOON!]\033[00m ignoring fSimTime " << fSimTime_v[index]
            << " as it results in negative TDC (invalid)" << std::endl;
            continue;
        }
        
        const geo::PlaneGeo& planeGeo = geo->Plane(geo::PlaneID(fCryostat,fTPC,0)); // Get the coordinates of the first wire plane

        TVector3 planeCoords = planeGeo.GetCenter();
        TVector3 planeNormal = planeGeo.GetNormalDirection();

        // Assume C=0, T=1
        chargeDepCoords = geo::Point_t(planeCoords[0] + 1. * planeNormal[0],planeCoords[1],planeCoords[2]);

        alternative::TruthHit pulse_record;
        pulse_record.tdc = tdc;
        pulse_record.num_electrons = fNumElectrons_v[index];
        pulse_record.tick = clockData.TPCTDC2Tick(tdc);

        double nElecADC = detProp.ElectronsToADC() * fNumElectrons_v[index];

        std::cout << "==> x position of plane 0: " << planeCoords[0] << ", normal: " << planeNormal[0] << ", nElec: " << fNumElectrons_v[index] << ", nElecADC: " << nElecADC << std::endl;

        simDep_v->emplace_back(0,fNumElectrons_v[index],0.,nElecADC,chargeDepCoords,chargeDepCoords);

        for(size_t plane=0; plane<3; ++plane) 
        {
            int wire = fWiresByPlaneVec[plane][index];

            geo::WireID wireID(fCryostat,fTPC,plane,wire);

            raw::ChannelID_t channel = geo->PlaneWireToChannel(wireID);

            fPlaneChannelVec[plane].emplace_back(channel);

            // Let's go backwards...
            std::vector<geo::WireID> wireIDVec = geo->ChannelToWire(channel);

            std::cout << ">>Channel " << channel << ", C/T/P/W: " << wireID.Cryostat << "/" << wireID.TPC << "/" << wireID.Plane << "/" << wireID.Wire 
                      << " returns " << wireIDVec.size() << " IDs, C/T/P/W: " << wireIDVec.front().Cryostat << "/" << wireIDVec.front().TPC << "/" << wireIDVec.front().Plane << "/" << wireIDVec.front().Wire << std::endl;

            // Recover the positions 
            const geo::WireGeo& wireGeo = geo->Wire(wireID);

            double xyz[3];

            wireGeo.GetCenter(xyz);

            xyz[0] = chargeDepCoords.X();

            Position position = {xyz[0],xyz[1],xyz[2]};

            fPlanePositionVecVec[plane].emplace_back(position);

            if(fVerbose) std::cout << "[BUFFOON!]    plane " << plane << " channel " << channel << " ... wire " << wire << std::endl;

            pulse_record.channel_list[plane] = channel;

            sim::SimChannel sch(channel);

            sch.AddIonizationElectrons(1,   /// track id, keep 0 = invalid
                                       (unsigned int)tdc,
                                       fNumElectrons_v[index],
                                       xyz,
                                       100.); /// energy, keep -1 = invalid
            simch_v->emplace_back(std::move(sch));
        }

        pholder.Register(std::move(pulse_record));
    }

    _run = e.id().run();
    _subrun = e.id().subRun();
    _event = e.id().event();
    
    fTupleTree->Fill();
    e.put(std::move(simch_v));
    e.put(std::move(simDep_v));
    e.put(std::move(trigger_v));
    
    return;
}

DEFINE_ART_MODULE(SimTestPulseWire)