sbndcode/sbndcode/CRT/CRTTzeroProducer_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       CRTTzeroProducer
// Module Type: producer
// File:        CRTTzeroProducer_module.cc
// Description: Module for constructiong over-simplified CRT tracks.
// Copied from CRTTrackProducer by David Lorca Galindo 
//  Edited by Michelle Stancari April 3, 2018
////////////////////////////////////////////////////////////////////////

#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/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/Assns.h"
#include "canvas/Persistency/Provenance/ProductID.h"
#include "art/Persistency/Common/PtrMaker.h"

//#include "lardata/Utilities/PtrMaker.h"
#include "lardata/Utilities/AssociationUtil.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "art_root_io/TFileService.h"

#include "sbnobj/Common/CRT/CRTHit.hh"
#include "sbnobj/Common/CRT/CRTTrack.hh"
#include "sbnobj/Common/CRT/CRTTzero.hh"

#include "TTree.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3S.h"
#include "TProfile.h"
#include "TF1.h"
#include "TDatime.h"

#include <iostream>
#include <stdio.h>
#include <sstream>
#include <vector>
#include <map>
#include <utility>
#include <cmath> 
#include <memory>

namespace sbnd{

class CRTTzeroProducer : public art::EDProducer {

public:

  explicit CRTTzeroProducer(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.
  CRTTzeroProducer(CRTTzeroProducer const &) = delete;
  CRTTzeroProducer(CRTTzeroProducer &&) = delete;
  CRTTzeroProducer & operator = (CRTTzeroProducer const &) = delete;
  CRTTzeroProducer & operator = (CRTTzeroProducer &&) = delete;

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

  // Selected optional functions.
  void beginJob() override;
  void endJob() override;

private:

  // Declare member data here.
  //  art::ServiceHandle<art::TFileService> tfs;

  std::string  data_label_;
  double max_time_difference_ ;//max time for coincidence 
  int store_tzero_;
  int verbose_ = 0;

};

void vmanip(std::vector<double> v, double* ave, double* rms);

void set_def(sbn::crt::CRTTzero tz);

CRTTzeroProducer::CRTTzeroProducer(fhicl::ParameterSet const & p)
  : EDProducer(p),

  // Initialize member data here.
  data_label_(p.get<std::string>("data_label")),
  max_time_difference_(p.get<double>("max_time_difference")),
  store_tzero_(p.get<int>("store_tzero")),
  verbose_(p.get<int>("verbose"))

{

  // Call appropriate produces<>() functions here.
  if(store_tzero_ == 1)  {
    produces< std::vector<sbn::crt::CRTTzero>   >();
    produces<art::Assns<sbn::crt::CRTTzero, sbn::crt::CRTHit> >();
  }

}

void CRTTzeroProducer::produce(art::Event & evt)
{
  // Implementation of required member function here.
  
  art::Handle< std::vector<sbn::crt::CRTHit> > rawHandle;
  evt.getByLabel(data_label_, rawHandle);   
  //check to make sure the data we asked for is valid                                                                       
  if(!rawHandle.isValid()){
    std::cout << "Run " << evt.run() << ", subrun " << evt.subRun()
              << ", event " << evt.event() << " has zero"
              << " CRTHits " << " in module " << data_label_ << std::endl;
    std::cout << std::endl;
    return;
  }
  
  //get better access to the data               
  std::vector<sbn::crt::CRTHit> const& CRTHitCollection(*rawHandle);

  //CRTTzero collection on this event                                                              
  std::unique_ptr<std::vector<sbn::crt::CRTTzero> > CRTTzeroCol(new std::vector<sbn::crt::CRTTzero>);
  
  // Output collections  
  std::unique_ptr<art::Assns<sbn::crt::CRTTzero, sbn::crt::CRTHit>> outputHits(new art::Assns<sbn::crt::CRTTzero, sbn::crt::CRTHit>);
  //  auto outputHits    = std::make_unique<art::Assns<sbn::crt::CRTTzero, sbn::crt::CRTHit>>();
  art::PtrMaker<sbn::crt::CRTHit> hitPtrMaker(evt, rawHandle.id());
  art::PtrMaker<sbn::crt::CRTTzero> tzeroPtrMaker(evt);
 
  int N_CRTHits = CRTHitCollection.size();

  //int iflag[1000] = {};
  std::vector<int> iflag;
  for(int i = 0; i < N_CRTHits; i++) iflag.push_back(0);
  int nTzero = 0;
  uint planeA, planeB;
  for(int  i = 0; i < N_CRTHits; i++) {//A 
    if (iflag[i]==0) {  // new tzero
      //temporary hit collection for each tzero
      std::vector<art::Ptr<sbn::crt::CRTHit>> CRTHitCol;
      sbn::crt::CRTHit CRTHiteventA = CRTHitCollection[i];
      art::Ptr<sbn::crt::CRTHit> hptr = hitPtrMaker(i);
      CRTHitCol.push_back(hptr);

      double time_s_A = 0; //CRTHiteventA.ts0_s;
      double time_ns_A = CRTHiteventA.ts1_ns;

      iflag[i]=1;

      // create and initialize, ugly code :(
      sbn::crt::CRTTzero CRTcanTzero;
      CRTcanTzero.ts0_ns=0;
      CRTcanTzero.ts1_ns=0;
      for (int j=0; j<7 ;++j) { // NUMBER OF PLANES, CHANGE TO 7
        CRTcanTzero.nhits[j]=0;
        CRTcanTzero.pes[j]=0;
      }

      // 
      CRTcanTzero.ts0_s=CRTHiteventA.ts0_s;
      CRTcanTzero.ts0_s_err=0;
      planeA = CRTHiteventA.plane;
      CRTcanTzero.nhits[planeA]=1;      

      int icount=1;

      CRTcanTzero.pes[planeA]=CRTHiteventA.peshit;
      for(int j = i+1; j < N_CRTHits; j++) {//B
        if (iflag[j]==0) {
          sbn::crt::CRTHit CRTHiteventB = CRTHitCollection[j];
          //look for coincidences
          double time_s_B = 0; //CRTHiteventB.ts0_s;
          double time_ns_B = CRTHiteventB.ts1_ns;
          double time_diff = time_ns_B - time_ns_A;

          if( (time_s_A == time_s_B) && (abs(time_diff)<max_time_difference_)  ){//D
            art::Ptr<sbn::crt::CRTHit> hptr = hitPtrMaker(j);
            CRTHitCol.push_back(hptr);
            planeB = CRTHiteventB.plane; 
            CRTcanTzero.nhits[planeB]+=1;
            CRTcanTzero.pes[planeB]+=CRTHiteventB.peshit;
            iflag[j]=1;
            CRTcanTzero.ts1_ns+=(int)(time_diff);
            CRTcanTzero.ts0_ns+=(CRTHiteventB.ts0_ns-CRTHiteventA.ts0_ns);
            icount++;
          }
        }
      }      // done with this tzero

      // Make a tzero data product
      CRTcanTzero.ts1_ns/=icount; 
      CRTcanTzero.ts1_ns+=(int)time_ns_A;
      CRTcanTzero.ts0_ns/=icount;
      CRTcanTzero.ts0_ns+=(int)CRTHiteventA.ts0_ns;
      CRTcanTzero.ts1_ns_err=0.;
      CRTcanTzero.ts0_ns_err=0.;
      CRTTzeroCol->push_back(CRTcanTzero);
      nTzero++;

      //associate hits to this Tzero
      art::Ptr<sbn::crt::CRTTzero> aptz = tzeroPtrMaker(CRTTzeroCol->size()-1);
      util::CreateAssn(*this,evt,aptz,CRTHitCol,*outputHits);
     
    }//B
    
  }//A

  //store tzero collection into event
  if(store_tzero_ == 1) {
    evt.put(std::move(CRTTzeroCol));
    evt.put(std::move(outputHits));
  }
}

void CRTTzeroProducer::beginJob()
{
  
}

void CRTTzeroProducer::endJob()
{
  // Implementation of optional member function here.
}

DEFINE_ART_MODULE(CRTTzeroProducer)

} //namespace sbnd