SuperDDTEva_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// \file    SuperDDTEva.cxx
//
// \brief   First official DDT analyzer in offline!
//
// \author  ricken@fnal.gov
////////////////////////////////////////////////////////////////////////
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "MCCheater/BackTracker.h"
#include "nusimdata/SimulationBase/MCTruth.h"

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

#include "DAQChannelMap/DAQChannelMap.h"
#include "RecoBase/Track.h"
#include "RecoBase/Cluster.h"
#include "Simulation/FLSHit.h"
#include "Simulation/FLSHitList.h"


#include "TVector3.h"
#include "TLorentzVector.h"
#include "TNtuple.h"
#include "TH1I.h"

#include <string>
#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <fstream>

class TTree;

namespace ddt {
  class SuperDDTEva;
  struct Chit;
}

struct ddt::Chit 
{
  Chit(unsigned short p, unsigned short c, uint32_t t, int16_t a, int ssid, int tsid, int rnsid, int trid) :
    plane(p), cell(c), tdc(t), adc(a), ssliceid(ssid), tsliceid(tsid), rnsliceid(rnsid),  trackid(trid) {}

  unsigned short plane, cell;
  uint32_t tdc;
  int16_t adc;
  int ssliceid, tsliceid, rnsliceid, trackid;
};

struct trueHit
{
  trueHit(unsigned short p, unsigned short c, double t, bool f, double e, int _pdg) :
    plane(p), cell(c), texit(t), found(f), energy(e), P_D_G(_pdg) {}
  unsigned short plane, cell;
  double texit;
  bool found; 
  double energy;
  int P_D_G;
};

class ddt::SuperDDTEva : public art::EDAnalyzer
{
public:
  explicit SuperDDTEva(fhicl::ParameterSet const &pset);
  virtual ~SuperDDTEva();
                        
  void analyze(const art::Event &evt);
  void beginJob();
  void reconfigure(fhicl::ParameterSet const &pset);
      
private:
  TTree*      MCTruth_Tree;
  TTree*      HitList_Tree;
  TTree*      Track_Tree;
  TTree*      TrueHits_Tree;
  /*
  TH1I*    nhits_hist;
  TH1I*    ntracks_hist;
  */
  TH1I*       n_RNSlice;
  TH1I*       n_SpaceSlice;

  //Truth info    
  int pdg;
  int nFLShits;
  double Px_truth;
  double Py_truth;
  double Pz_truth;
  double energydump_truth;
  double pathlength_truth;
  double x0_truth;
  double y0_truth;
  double z0_truth;
  double t0_truth;
  double x1_truth;
  double y1_truth;
  double z1_truth;
  double t1_truth;

  //Track info
  int view;
  double v0_ddt;
  double z0_ddt;
  double t0_ddt;
  double v1_ddt;
  double z1_ddt;
  double t1_ddt;

  //HitList Info:
  int Plane, Cell,ADC;
  uint32_t TDC;
  int trackID, TsliceID, RNsliceID, SsliceID,eventID;
  bool isNoise, isMM;

  std::string TimeSlice_ModuleLabel;
  std::string TimeSlice_InstanceLabel;

  std::string RNSlice_ModuleLabel;
  std::string RNSlice_InstanceLabel;

  std::string SpaceSlice_ModuleLabel;
  std::string SpaceSlice_InstanceLabel;

  std::string OnlineTracker_ModuleLabel;
  std::string OnlineTracker_InstanceLabel;
  std::string fFLSHitListLabel;

  //True Hits Info
  int Plane_true, Cell_true;
  double Time_true, Energy_true;
  bool Found_true; 

};



ddt::SuperDDTEva::SuperDDTEva(fhicl::ParameterSet const& pset) 
  : EDAnalyzer(pset)
{
  reconfigure(pset);
}


//......................................................................
void ddt::SuperDDTEva::reconfigure(fhicl::ParameterSet const &pset)
{
  TimeSlice_ModuleLabel          = pset.get< std::string >("TimeSlice_ModuleLabel" );
  TimeSlice_InstanceLabel        = pset.get< std::string >("TimeSlice_InstanceLabel" );
  RNSlice_ModuleLabel            = pset.get< std::string >("RNSlice_ModuleLabel" );
  RNSlice_InstanceLabel          = pset.get< std::string >("RNSlice_InstanceLabel" );
  SpaceSlice_ModuleLabel         = pset.get< std::string >("SpaceSlice_ModuleLabel" );
  SpaceSlice_InstanceLabel       = pset.get< std::string >("SpaceSlice_InstanceLabel" );
  OnlineTracker_ModuleLabel      = pset.get< std::string >("OnlineTracker_ModuleLabel" );
  OnlineTracker_InstanceLabel    = pset.get< std::string >("OnlineTracker_InstanceLabel" );
  fFLSHitListLabel               = pset.get< std::string >("FLSHitListLabel") ;
}

//......................................................................

ddt::SuperDDTEva::~SuperDDTEva(){}
  

//......................................................................

void ddt::SuperDDTEva::beginJob()
{
  art::ServiceHandle<art::TFileService> tfs;
  n_RNSlice    = tfs->make<TH1I>("n_RNSlice","number of remove noise slice",1000,-0.5,999.5);
  n_SpaceSlice = tfs->make<TH1I>("n_SpaceSlice","number of space slice",100,-0.5,99.5);

  MCTruth_Tree = tfs->make<TTree>("MCTruth","MC_info");   
  //  MCTruth_Tree -> Branch("PDG",&pdg,"PDG/I");
  MCTruth_Tree -> Branch("nhits",&nFLShits,"nhits/I");
  MCTruth_Tree -> Branch("Px",&Px_truth,"Px/D");
  MCTruth_Tree -> Branch("Py",&Py_truth,"Py/D");
  MCTruth_Tree -> Branch("Pz",&Pz_truth,"Pz/D");
  MCTruth_Tree -> Branch("EnergyDump",&energydump_truth,"energydump/D");
  MCTruth_Tree -> Branch("PathLength",&pathlength_truth,"pathlength/D");
  MCTruth_Tree -> Branch("x0",&x0_truth,"x0/D");
  MCTruth_Tree -> Branch("y0",&y0_truth,"y0/D");
  MCTruth_Tree -> Branch("z0",&z0_truth,"z0/D");
  MCTruth_Tree -> Branch("t0",&t0_truth,"t0/D");
  MCTruth_Tree -> Branch("x1",&x1_truth,"x1/D");
  MCTruth_Tree -> Branch("y1",&y1_truth,"y1/D");
  MCTruth_Tree -> Branch("z1",&z1_truth,"z1/D");
  MCTruth_Tree -> Branch("t1",&t1_truth,"t1/D");
  MCTruth_Tree -> Branch("evtID",&eventID,"evtID/I");

  HitList_Tree = tfs->make<TTree>("HitList","HitList_info");   
  HitList_Tree -> Branch("Plane",&Plane,"Plane/I");
  HitList_Tree -> Branch("Cell",&Cell,"Cell/I");
  HitList_Tree -> Branch("TDC",&TDC,"TDC/i");
  HitList_Tree -> Branch("ADC",&ADC,"ADC/I");
  HitList_Tree -> Branch("trackID",&trackID,"trackID/I");
  HitList_Tree -> Branch("TimesliceID",&TsliceID,"TimesliceID/I");
  HitList_Tree -> Branch("RNsliceID",&RNsliceID,"RNsliceID/I");
  HitList_Tree -> Branch("SsliceID",&SsliceID,"SsliceID/I");
  HitList_Tree -> Branch("eventID",&eventID,"eventID/I");
  HitList_Tree -> Branch("IsNoise",&isNoise,"IsNoise/O");
  HitList_Tree -> Branch("IsMM",&isMM,"IsMM/O");

  Track_Tree = tfs->make<TTree>("Track_Tree","Track_info");   
  Track_Tree -> Branch("View",&view,"View/I");
  Track_Tree -> Branch("v0",&v0_ddt,"v0/D");
  Track_Tree -> Branch("z0",&z0_ddt,"z0/D");
  Track_Tree -> Branch("t0",&t0_ddt,"t0/D");
  Track_Tree -> Branch("v1",&v1_ddt,"v1/D");
  Track_Tree -> Branch("z1",&z1_ddt,"z1/D");
  Track_Tree -> Branch("t1",&t1_ddt,"t1/D");
  Track_Tree -> Branch("eventID",&eventID,"eventID/I");

  TrueHits_Tree = tfs->make<TTree>("TrueHits_Tree","TrueHits_info");   
  TrueHits_Tree -> Branch("Plane",&Plane_true,"Plane/I");
  TrueHits_Tree -> Branch("Cell",&Cell_true,"Cell/I");
  TrueHits_Tree -> Branch("Energy",&Energy_true,"Energy/D");
  TrueHits_Tree -> Branch("Time",&Time_true,"Time/D");
  TrueHits_Tree -> Branch("Found",&Found_true,"Found/O");
  TrueHits_Tree -> Branch("eventID",&eventID,"eventID/I");


  /*
  nhits_hist   = tfs->make<TH1I>("nhitsHist","Number of Hits per Slice",1000,0,1000);
  ntracks_hist = tfs->make<TH1I>("ntracks","Number of Tracks per Slice",50,0,50);
  */
}

void ddt::SuperDDTEva::analyze(const art::Event& evt)
{
  eventID = evt.id().event();
  std::vector<Chit> chits;
  std::vector<trueHit> thits;


  //Get Truth Info
  //Make a backtracker
  art::ServiceHandle<cheat::BackTracker> bt;
  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();

  for (sim::ParticleNavigator::const_iterator i = pnav.begin(); i != pnav.end(); ++i){
    const sim::Particle *p = (*i).second;
    pdg = p->PdgCode();
    Px_truth = p->Px();
    Py_truth = p->Py();
    Pz_truth = p->Pz();
 
    if (bt->ParticleToFLSHit(p->TrackId()).size()==0) continue;

    //Later info needs to be obtained from the FLSHit:
    const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(p->TrackId());

    pathlength_truth = flshits[0].GetTotalPathLength();
    //pathlength_truth = 0;
    energydump_truth = flshits[0].GetTotalEnergyLoss();
    x0_truth = flshits[0].GetEntryX();
    y0_truth = flshits[0].GetEntryY();
    z0_truth = flshits[0].GetEntryZ();
    t0_truth = flshits[0].GetEntryT();
    x1_truth = flshits[0].GetExitX();
    y1_truth = flshits[0].GetExitY();
    z1_truth = flshits[0].GetExitZ();
    t1_truth = flshits[0].GetExitT();

    thits.emplace_back(flshits[0].GetPlaneID(),flshits[0].GetCellID(),t1_truth,false,energydump_truth, pdg); 

    nFLShits = flshits.size();
    for (int h =1; h!= nFLShits; ++h) {

      //     std::cout<<"plane: "<<flshits[h].GetPlaneID()<<" cell: "<<flshits[h].GetCellID()<<" time: "<<t1_truth<<" pdg: "<<pdg<<std::endl;
      thits.emplace_back(flshits[h].GetPlaneID(),flshits[h].GetCellID(),t1_truth,false,flshits[h].GetTotalEnergyLoss(), pdg); 
      //      std::cout<<"t_true:"<<t1_truth<<std::endl;
      //      std::cout<<"energy: "<<flshits[h].GetTotalEnergyLoss()<<std::endl;

      pathlength_truth+=flshits[h].GetTotalPathLength();
      energydump_truth+=flshits[h].GetTotalEnergyLoss();
      if (t0_truth > flshits[h].GetEntryT() ) {
        x0_truth = flshits[h].GetEntryX();
        y0_truth = flshits[h].GetEntryY();
        z0_truth = flshits[h].GetEntryZ();
        t0_truth = flshits[h].GetEntryT();
      }
      else if (t1_truth < flshits[h].GetExitT() ) {
        x1_truth = flshits[h].GetExitX();
        y1_truth = flshits[h].GetExitY();
        z1_truth = flshits[h].GetExitZ();
        t1_truth = flshits[h].GetExitT();
      }
    }
    MCTruth_Tree->Fill();
  }


  //Get TimeSlice info:
  art::Handle<novaddt::GroupedHitList<novaddt::TDC> > Tghits;
  evt.getByLabel(TimeSlice_ModuleLabel, TimeSlice_InstanceLabel, Tghits);
  novaddt::GroupedHitList<novaddt::TDC> ghl(*Tghits);
  if (ghl.n_groups()==0) return;
  //  std::cout<<"number of original time slices: "<< ghl.n_groups()<<std::endl;
  for (unsigned i =0; i!= ghl.n_groups();++i) {
    novaddt::HitList hits = ghl.get(i);
 
    for (auto hit_ptr = hits.begin(); hit_ptr!=hits.end();++hit_ptr){
      chits.emplace_back(hit_ptr->Plane().val, hit_ptr->Cell().val,
                         hit_ptr->TDC().val, hit_ptr->ADC().val,
                         0,i+1, 0, 0                     //spacesliceID, timesliceID, rnsliceID, trackID
                         );
    }

  }

  //Flip the RNSliceID:
  art::Handle<novaddt::GroupedHitList<novaddt::TDC> > RNghits;
  evt.getByLabel(RNSlice_ModuleLabel, RNSlice_InstanceLabel, RNghits);
  novaddt::GroupedHitList<novaddt::TDC> RNghl(*RNghits);
  n_RNSlice -> Fill(RNghl.n_groups());
  //  std::cout<<"number of removenoise slice: "<< RNghl.n_groups()<<std::endl;
  for (unsigned i =0; i!= RNghl.n_groups(); ++i) {
    novaddt::HitList hits = RNghl.get(i);
    for (auto hit_ptr = hits.begin(); hit_ptr!=hits.end();++hit_ptr){
      for (auto chit_ptr= chits.begin(); chit_ptr!=chits.end();++chit_ptr ) {
        if (chit_ptr->rnsliceid) continue;
        if (chit_ptr->plane == hit_ptr->Plane().val && 
            chit_ptr->cell == hit_ptr->Cell().val &&
            chit_ptr->tdc == hit_ptr->TDC().val) {
          chit_ptr->rnsliceid =i+1;
          break;
        }
      }
    }
  }


  //Flip the SpaceSliceID:
  art::Handle<std::vector<novaddt::HitList> > hit_lists; // from HitList
  evt.getByLabel(SpaceSlice_ModuleLabel, SpaceSlice_InstanceLabel, hit_lists);
  n_SpaceSlice->Fill(hit_lists->size());
  for (unsigned i=0; i!=hit_lists->size(); ++i) {
    novaddt::HitList const &hits = hit_lists->at(i);
    for (auto hit_ptr = hits.begin(); hit_ptr!=hits.end();++hit_ptr){
      for (auto chit_ptr= chits.begin(); chit_ptr!=chits.end();++chit_ptr ) {
        if (chit_ptr->ssliceid > 0) continue;
        if (chit_ptr->plane == hit_ptr->Plane().val && 
            chit_ptr->cell == hit_ptr->Cell().val &&
            chit_ptr->tdc == hit_ptr->TDC().val) {
          chit_ptr->ssliceid =i+1;
          break;
        }
      }
    }
  }
  /*
  //Get Track Info
  art::Handle<std::vector<novaddt::Track> > tracks;     
  evt.getByLabel(OnlineTracker_ModuleLabel,OnlineTracker_InstanceLabel, tracks);

  art::FindOneP<novaddt::HitList> fohl(tracks, evt, OnlineTracker_ModuleLabel);
  for (unsigned i =0; i != tracks->size(); ++i) {
    view   = tracks->at(i).View();
    v0_ddt = tracks->at(i).StartV();
    z0_ddt = tracks->at(i).StartZ();
    v1_ddt = tracks->at(i).EndV();
    z1_ddt = tracks->at(i).EndZ();
    t0_ddt = tracks->at(i).StartT();
    t1_ddt = tracks->at(i).EndT();

    Track_Tree->Fill();
    //Get associated hit list:
    art::Ptr<novaddt::HitList> hl_ptr = fohl.at(i);
    for (auto const& hit : (*hl_ptr) ) {
      //loop over the vector of chits:
      for (auto chit_ptr= chits.begin(); chit_ptr!=chits.end();++chit_ptr ) {
        if (chit_ptr->trackid) continue;
        if (chit_ptr->plane == hit.Plane().val && 
            chit_ptr->cell == hit.Cell().val &&
            chit_ptr->tdc == hit.TDC().val) {
          chit_ptr->trackid =i+1;
          break;
        }
      }
    }

  }
  */
  //  std::cout<<"Here is the list of the found true hits:"<<std::endl;

  //Fill the hitlist tree:
  for (auto const& chit : chits) {
    Plane = chit.plane;
    Cell = chit.cell;
    TDC = chit.tdc;
    ADC = chit.adc;
    TsliceID = chit.tsliceid;
    RNsliceID = chit.rnsliceid;
    SsliceID = chit.ssliceid;
    trackID = chit.trackid;
    isNoise = true;
    isMM = false; 
    for (auto true_hit_ptr= thits.begin(); true_hit_ptr != thits.end(); ++true_hit_ptr) {
      if (true_hit_ptr->plane == Plane && true_hit_ptr->cell== Cell) {

        if ( std::max(TDC,uint32_t(true_hit_ptr->texit / 15.625)) 
             -std::min(TDC,uint32_t(true_hit_ptr->texit / 15.625)) <70  ) {
          isNoise = false;

          //      std::cout<<"plane: "<<Plane<<"   cell:"<<Cell<<std::endl;

          //      if (true_hit_ptr->found) std::cout<<"This hit has already been assigned to a DAQHit!"<<std::endl;
          true_hit_ptr->found = true;
          if (true_hit_ptr->P_D_G == 42) {
            isMM = true;
            break;
          }         
        }
      } 
    }
    
    HitList_Tree -> Fill();
  }

  for (auto const& true_hit: thits) {
    Plane_true = true_hit.plane;
    Cell_true = true_hit.cell;
    Energy_true = true_hit.energy;
    Found_true = true_hit.found;
    Time_true = true_hit.texit;
    //    std::cout<<"energy dump: "<< true_hit.energy<<std::endl;
    TrueHits_Tree-> Fill();
  } 

  return;
}


DEFINE_ART_MODULE(ddt::SuperDDTEva)