MergeWireCell_module.cc

Go to the documentation of this file.

// Class:       MergeWireCell
// Module Type: producer
// File:        MergeWireCell_module.cc
//
// Generated at Sat Sep  5 12:50:05 2015 by Tingjun Yang using artmod
// from cetpkgsupport v1_08_06.

#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 "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/Shower.h"
#include "larcore/Geometry/Geometry.h"
#include "lardata/Utilities/AssociationUtil.h"

#include "TTree.h"
#include "TBranch.h"
#include "TLeaf.h"
#include "TFile.h"

#include <memory>
#include <string>
#include <dirent.h>
#include <iostream>

namespace wc {
  class MergeWireCell;
}

class wc::MergeWireCell : public art::EDProducer {
public:
  explicit MergeWireCell(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.
  MergeWireCell(MergeWireCell const &) = delete;
  MergeWireCell(MergeWireCell &&) = delete;
  MergeWireCell & operator = (MergeWireCell const &) = delete;
  MergeWireCell & operator = (MergeWireCell &&) = delete;

  // Required functions.
  void produce(art::Event & e) override;

  // Selected optional functions.
  void beginJob() override;
  void reconfigure(fhicl::ParameterSet const & p) ;

  // Make recob::Tracks
  void MakeTracks(art::Event &evt,
                  std::unique_ptr<std::vector<recob::Track>> &trk_coll,
                  std::unique_ptr<art::Assns<recob::Track, recob::Hit>> &trkhassn,
                  std::unique_ptr<art::Assns<recob::Track, recob::SpacePoint>> &trksassn,
                  std::unique_ptr<std::vector<recob::Hit>> &hit_coll,
                  std::unique_ptr<std::vector<recob::SpacePoint>> &spt_coll,              
                  std::vector<int> &index,
                  TTree *tree);

  // Make recob::Showers
  void MakeShowers(art::Event &evt,
                  std::unique_ptr<std::vector<recob::Shower>> &shw_coll,
                  std::unique_ptr<art::Assns<recob::Shower, recob::Hit>> &shwhassn,
                  std::unique_ptr<art::Assns<recob::Shower, recob::SpacePoint>> &shwsassn,
                  std::unique_ptr<std::vector<recob::Hit>> &hit_coll,
                  std::unique_ptr<std::vector<recob::SpacePoint>> &spt_coll,              
                  std::vector<int> &index,
                  TTree *tree);

private:

  // Declare member data here.
  std::string fInput;

  static constexpr float Sqrt2Pi = 2.5066;
  static constexpr float SqrtPi  = 1.7725;

};


wc::MergeWireCell::MergeWireCell(fhicl::ParameterSet const & p) :
  fInput(p.get<std::string>("WireCellInput"))
{
  // Call appropriate produces<>() functions here.
  produces<std::vector<recob::Hit>                      >();
  produces<std::vector<recob::SpacePoint>               >();
  produces<art::Assns<recob::Hit, recob::SpacePoint>    >();
  produces<std::vector<recob::Track>                    >();
  produces<std::vector<recob::Shower>                   >();
  produces<art::Assns<recob::Track, recob::Hit>         >();
  produces<art::Assns<recob::Track, recob::SpacePoint>  >();
  produces<art::Assns<recob::Shower, recob::Hit>        >();
  produces<art::Assns<recob::Shower, recob::SpacePoint> >();
}

void wc::MergeWireCell::produce(art::Event & evt){

  art::ServiceHandle<geo::Geometry> geom;

  int run = evt.run();
  int subrun = evt.subRun();
  int event = evt.id().event();

  std::unique_ptr<std::vector<recob::Hit>> hit_coll(new std::vector<recob::Hit>);
  std::unique_ptr<std::vector<recob::SpacePoint>> spt_coll(new std::vector<recob::SpacePoint>);
  std::unique_ptr<art::Assns<recob::SpacePoint, recob::Hit> >   shassn(new art::Assns<recob::SpacePoint, recob::Hit>);
  std::unique_ptr<std::vector<recob::Track>> trk_coll(new std::vector<recob::Track>);
  std::unique_ptr<art::Assns<recob::Track, recob::Hit>> trkhassn(new art::Assns<recob::Track, recob::Hit>);
  std::unique_ptr<art::Assns<recob::Track, recob::SpacePoint>> trksassn(new art::Assns<recob::Track, recob::SpacePoint>);
  std::unique_ptr<std::vector<recob::Shower>> shw_coll(new std::vector<recob::Shower>);
  std::unique_ptr<art::Assns<recob::Shower, recob::Hit>> shwhassn(new art::Assns<recob::Shower, recob::Hit>);
  std::unique_ptr<art::Assns<recob::Shower, recob::SpacePoint>> shwsassn(new art::Assns<recob::Shower, recob::SpacePoint>);

  std::string path(fInput);
  path=(path+"/");

  DIR *pDIR;
  struct dirent *entry;
  if( (pDIR=opendir(path.c_str())) != NULL ){
     
    while((entry = readdir(pDIR)) != NULL){

      if( strcmp(entry->d_name, ".")==0 || strcmp(entry->d_name, "..")==00) continue;

      std::string filename(entry->d_name);    
      if((int)filename.find(".root")==-1)
      continue;

      std::string file = (path + entry->d_name);
      TFile *f = new TFile(file.c_str());
      TTree *Trun = (TTree*)f->Get("Trun");
      int eventNo, runNo, subRunNo;
      Trun->SetBranchAddress("eventNo",&eventNo);
      Trun->SetBranchAddress("runNo",&runNo);
      Trun->SetBranchAddress("subRunNo",&subRunNo);
      if (!Trun->GetEntries()) continue;
      Trun->GetEntry(0);
      //std::cout<<run<<" "<<subrun<<" "<<event<<" "<<runNo<<" "<<subRunNo<<" "<<eventNo<<std::endl;
      if (runNo!=run||eventNo!=event||subRunNo!=subrun) {
        f->Close();
        continue;
      }
      TTree *TC = (TTree*)f->Get("TC");
      Int_t           time_slice;
      Double_t        charge;
      Double_t        xx;
      Double_t        yy;
      Double_t        zz;
      Int_t           u_index;
      Int_t           v_index;
      Int_t           w_index;
      Double_t        u_charge;
      Double_t        v_charge;
      Double_t        w_charge;
      Double_t        u_charge_err;
      Double_t        v_charge_err;
      Double_t        w_charge_err;
      Int_t           tpc_no;
      Int_t           cryostat_no;
      Int_t           mcell_id;

      TC->SetBranchAddress("time_slice", &time_slice);
      TC->SetBranchAddress("charge", &charge);
      TC->SetBranchAddress("xx", &xx);
      TC->SetBranchAddress("yy", &yy);
      TC->SetBranchAddress("zz", &zz);
      TC->SetBranchAddress("u_index", &u_index);
      TC->SetBranchAddress("v_index", &v_index);
      TC->SetBranchAddress("w_index", &w_index);
      TC->SetBranchAddress("u_charge", &u_charge);
      TC->SetBranchAddress("v_charge", &v_charge);
      TC->SetBranchAddress("w_charge", &w_charge);
      TC->SetBranchAddress("u_charge_err", &u_charge_err);
      TC->SetBranchAddress("v_charge_err", &v_charge_err);
      TC->SetBranchAddress("w_charge_err", &w_charge_err);
      TC->SetBranchAddress("cryostat_no", &cryostat_no);
      TC->SetBranchAddress("mcell_id",&mcell_id);
        
        TC->GetDirectory()->Print();
        
//      TC->SetBranchAddress("cell.tpc_no", &tpc_no);
        
        TObjArray* leafList = TC->GetListOfLeaves();
        
        for(int idx = 0; idx < leafList->GetEntries(); idx++)
        {
            leafList->At(idx)->Print();
        }
        
        TObject* candObj = leafList->FindObject("tpc_no");
        
        std::cout << "candObj: " << candObj << std::endl;

        TLeaf* candLeaf = (TLeaf*)candObj;
        candLeaf->Print();
        candLeaf->SetAddress(&tpc_no);

      std::vector<int> vcell_id;

      for (int i = 0; i<TC->GetEntries(); ++i){
        TC->GetEntry(i);
    
    tpc_no = candLeaf->GetTypedValue<Int_t>(i) - 1;
          
        vcell_id.push_back(mcell_id);
        double xyz[3] = {xx,yy,zz};
        double err[3] = {0,0,0};
        spt_coll->push_back(recob::SpacePoint(xyz,err,charge));
        geo::WireID wireu(cryostat_no,tpc_no,0,u_index);
        raw::ChannelID_t chanu = geom->PlaneWireToChannel(wireu);
        geo::WireID wirev(cryostat_no,tpc_no,1,v_index);
        raw::ChannelID_t chanv = geom->PlaneWireToChannel(wirev);
        geo::WireID wirew(cryostat_no,tpc_no,2,w_index);
        raw::ChannelID_t chanw = geom->PlaneWireToChannel(wirew);
        size_t hitStart = hit_coll->size();
//      raw::TDCtick_t start_tick = time_slice*4;
//      raw::TDCtick_t end_tick = (time_slice+1)*4;
//      float peak_time = (time_slice+0.5)*4;
//      float sigma_peak = 2.0;
//      float rms = 2.0;
//      float peak_ampu = u_charge/(Sqrt2Pi*2);
//      float sigma_peak_a
        hit_coll->push_back(recob::Hit(chanu,
                                       time_slice*4,
                                       (time_slice+1)*4,
                                       (time_slice+0.5)*4,
                                       2.,
                                       2.,
                                       u_charge/(Sqrt2Pi*2),
                                       u_charge_err/(Sqrt2Pi*2),
                                       u_charge,
                                       u_charge,
                                       u_charge_err,
                                       1,
                                       0,
                                       0,
                                       0,
                                       geom->View(chanu),
                                       geom->SignalType(chanu),
                                       wireu));
        hit_coll->push_back(recob::Hit(chanv,
                                       time_slice*4,
                                       (time_slice+1)*4,
                                       (time_slice+0.5)*4,
                                       2.,
                                       2.,
                                       v_charge/(Sqrt2Pi*2),
                                       v_charge_err/(Sqrt2Pi*2),
                                       v_charge,
                                       v_charge,
                                       v_charge_err,
                                       1,
                                       0,
                                       0,
                                       0,
                                       geom->View(chanv),
                                       geom->SignalType(chanv),
                                       wirev));
        hit_coll->push_back(recob::Hit(chanw,
                                       time_slice*4,
                                       (time_slice+1)*4,
                                       (time_slice+0.5)*4,
                                       2.,
                                       2.,
                                       w_charge/(Sqrt2Pi*2),
                                       w_charge_err/(Sqrt2Pi*2),
                                       w_charge,
                                       w_charge,
                                       w_charge_err,
                                       1,
                                       0,
                                       0,
                                       0,
                                       geom->View(chanw),
                                       geom->SignalType(chanw),
                                       wirew));
        size_t hitEnd = hit_coll->size();
        util::CreateAssn(*this, evt, *spt_coll, *hit_coll, *shassn, hitStart, hitEnd);
      }//Loop over TC
/*
      TTree *T_goodtrack = (TTree*)f->Get("T_goodtrack");
      MakeTracks(evt, trk_coll, trkhassn, trksassn, hit_coll, spt_coll, vcell_id, T_goodtrack);
      TTree *T_shorttrack = (TTree*)f->Get("T_shorttrack");
      MakeTracks(evt, trk_coll, trkhassn, trksassn, hit_coll, spt_coll, vcell_id, T_shorttrack);
      TTree *T_paratrack = (TTree*)f->Get("T_paratrack");
      MakeTracks(evt, trk_coll, trkhassn, trksassn, hit_coll, spt_coll, vcell_id, T_paratrack);

      TTree *T_shower = (TTree*)f->Get("T_shower");
      MakeShowers(evt, shw_coll, shwhassn, shwsassn, hit_coll, spt_coll, vcell_id, T_shower);
*/
      f->Close();
      break;
    }
    closedir(pDIR);

  }
    
    std::cout << "space point size: " << spt_coll->size() << ", hit coll size: " << hit_coll->size() << std::endl;
    
  evt.put(std::move(spt_coll));
  evt.put(std::move(hit_coll));
  evt.put(std::move(shassn));
  evt.put(std::move(trk_coll));
  evt.put(std::move(trkhassn));
  evt.put(std::move(trksassn));
  evt.put(std::move(shw_coll));
  evt.put(std::move(shwhassn));
  evt.put(std::move(shwsassn));
}

void wc::MergeWireCell::MakeTracks(art::Event &evt,
                                   std::unique_ptr<std::vector<recob::Track>> &trk_coll,
                                   std::unique_ptr<art::Assns<recob::Track, recob::Hit>> &trkhassn,
                                   std::unique_ptr<art::Assns<recob::Track, recob::SpacePoint>> &trksassn,
                                   std::unique_ptr<std::vector<recob::Hit>> &hit_coll,
                                   std::unique_ptr<std::vector<recob::SpacePoint>> &spt_coll,             
                                   std::vector<int> &index,
                                   TTree *tree){
  Int_t           npoints;
  Int_t           trackid;
  Int_t           msc_id[1000];   //[npoints]
  Double_t        x[1000];   //[npoints]
  Double_t        y[1000];   //[npoints]
  Double_t        z[1000];   //[npoints]
  Double_t        theta[1000];   //[npoints]
  Double_t        phi[1000];   //[npoints]

  tree->SetBranchAddress("npoints", &npoints);
  tree->SetBranchAddress("trackid", &trackid);
  tree->SetBranchAddress("msc_id", msc_id);
  tree->SetBranchAddress("x", x);
  tree->SetBranchAddress("y", y);
  tree->SetBranchAddress("z", z);
  tree->SetBranchAddress("theta", theta);
  tree->SetBranchAddress("phi", phi);

  Long64_t nentries = tree->GetEntriesFast();
  
  for (Long64_t i=0; i<nentries; ++i) {
    tree->GetEntry(i);
    std::vector<TVector3> xyz;
    std::vector<TVector3> dircos;
    std::vector< std::vector<double> > dQdx;
    std::vector<double> mom(2, util::kBogusD);
    std::vector<size_t> hits;
    std::vector<size_t> spts;
    for (int i = 0; i<npoints; ++i){
      xyz.push_back(TVector3(x[i],y[i],z[i]));
      dircos.push_back(TVector3(sin(theta[i])*cos(phi[i]),sin(theta[i])*sin(phi[i]),cos(theta[i])));
      for (size_t j = 0; j<index.size(); ++j){
        if (index[j]==msc_id[i]){
          spts.push_back(j);
          hits.push_back(3*j);
          hits.push_back(3*j+1);
          hits.push_back(3*j+2);
        }
      }
    }
    trk_coll->push_back(recob::Track(xyz, dircos, dQdx, mom, trackid));
    // make associations between the track and space points
    util::CreateAssn(*this, evt, *trk_coll, *spt_coll, *trksassn, spts);
    // make associations between the track and hits
    util::CreateAssn(*this, evt, *trk_coll, *hit_coll, *trkhassn, hits);
  }
}    

void wc::MergeWireCell::MakeShowers(art::Event &evt,
                                   std::unique_ptr<std::vector<recob::Shower>> &shw_coll,
                                   std::unique_ptr<art::Assns<recob::Shower, recob::Hit>> &shwhassn,
                                   std::unique_ptr<art::Assns<recob::Shower, recob::SpacePoint>> &shwsassn,
                                   std::unique_ptr<std::vector<recob::Hit>> &hit_coll,
                                   std::unique_ptr<std::vector<recob::SpacePoint>> &spt_coll,             
                                   std::vector<int> &index,
                                   TTree *tree){
  Int_t           npoints;
  Int_t           showerid;
  Int_t           msc_id[1000];   //[npoints]

  tree->SetBranchAddress("npoints", &npoints);
  tree->SetBranchAddress("showerid", &showerid);
  tree->SetBranchAddress("msc_id", msc_id);

  Long64_t nentries = tree->GetEntriesFast();
  
  for (Long64_t i=0; i<nentries; ++i) {
    tree->GetEntry(i);
    TVector3 v1(0,0,0); 
    std::vector< double > tmp(2, util::kBogusD);
    int bestplane = -1;
    std::vector<size_t> hits;
    std::vector<size_t> spts;
    for (int i = 0; i<npoints; ++i){
      for (size_t j = 0; j<index.size(); ++j){
        if (index[j]==msc_id[i]){
          spts.push_back(j);
          hits.push_back(3*j);
          hits.push_back(3*j+1);
          hits.push_back(3*j+2);
        }
      }
    }
    shw_coll->push_back(recob::Shower(v1, v1, v1, v1, tmp, tmp, tmp, tmp, bestplane, showerid));
    // make associations between the track and space points
    util::CreateAssn(*this, evt, *shw_coll, *spt_coll, *shwsassn, spts);
    // make associations between the track and hits
    util::CreateAssn(*this, evt, *shw_coll, *hit_coll, *shwhassn, hits);
  }
}    

void wc::MergeWireCell::beginJob()
{
  // Implementation of optional member function here.
}

void wc::MergeWireCell::reconfigure(fhicl::ParameterSet const & p)
{
  // Implementation of optional member function here.
}

DEFINE_ART_MODULE(wc::MergeWireCell)