PEResponse_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////                         
// Class:       PEResponse                                                                                          
// Module Type: Analyzer                                                                                              
// File:        PEResponse_module.cc
// \brief       A module to study channel response at a very basic level.
//
////////////////////////////////////////////////////////////////////////                         

// art and FHiCL includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindManyP.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/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/ParameterSet.h"

// NOvASoft includes
#include "RecoBase/Track.h"
#include "Geometry/Geometry.h"
#include "CMap/service/CMapService.h"
#include "RunHistory/RunHistory.h"
#include "CalibrationDataProducts/PCHit.h"
#include "SummaryData/CosmicExposure.h"
#include "DAQChannelMap/DAQChannelMap.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "ChannelInfo/BadChanList.h"

// C++                                                                                                               
#include <memory>
#include <string>
#include <iostream>
#include <vector>

// ROOT
#include "TFile.h"
#include "TTree.h"
#include "TVector3.h"
#include "TH1.h"
#include "TH2.h"

////////////////////////////////////////////////////////////////////////                         
////////////////////////////////////////////////////////////////////////                         


namespace calib 
{
  class PEResponse : public art::EDAnalyzer 
  {
  public:
    explicit PEResponse ( fhicl::ParameterSet const& p );
    virtual ~PEResponse ();
    
    void analyze     ( art::Event const& e );
    void reconfigure ( fhicl::ParameterSet const& p );
    void beginJob    ( );
    void beginRun    ( art::Run const& r );
    void endSubRun   ( art::SubRun const&  sr);
  private:

    void SetDefault();
    std::string fCellHitLabel;
    std::string fTrkLabel;     ///< Track module label
    std::string fPCHitLabel;   ///< PCHit module label
    std::string fQualXYLabel;  ///< Instance label for hits with xz quality                                   
    std::string fQualZLabel;   ///< Instance label for hits with z quality                                      
    std::string fQualAvgLabel; ///< Instance label for hits with avg quality    
    std::string fQualBadLabel; ///< Instance label for hits with bad track quality    
    std::string fExpLabel;     ///< Cosmic exposure module label name

    TTree* fHitTree;           ///< Tree to save hit information
    std::vector<TH2F*> fThreshHists;
    TH1F* hx[15];
    TH1F* hy[15];

    TH2F* hthresh;

    TH1F* hall;
    TH1F* hexp;
    //TH1F* ncellperplane;
    float fPe;
    float fW;
    float fPath;
    float fTns;
    float fReadDist;
    float fFlightDist;
    bool  fIsBad;
    bool  fIsXYAdj;
    bool  fIsZAdj;
    bool  fIsAvg;
    bool  fIsBadTrk;
    //bool  fHasApd;
    //bool  fIsEnabled;
    int   fEvent;
    int   fRun;
    int   fSubRun;
    int   fPlane;
    int   fCell;
    int   fDcm;
    int   fDiblock;
    int   fApd;
    int   fPixel;
    int   fView;
    int   fTrkPlanes;
    bool  fFillHitTree;

    art::ServiceHandle< geo::Geometry > fGeom;
    nova::dbi::RunHistory fRunHist;  //Run history object        
    

  }; // end of PEResponse class definition

  ////////////////////////////////////////////////////////////////////////  
  ////////////////////////////////////////////////////////////////////////                         

  ///< Constructor
  PEResponse::PEResponse( fhicl::ParameterSet const& pset )
    :EDAnalyzer( pset )
  {
    this->reconfigure(pset);
  }

  ////////////////////////////////////////////////////////////////////////  
  ////////////////////////////////////////////////////////////////////////                         


  ///< Destructor
  PEResponse::~PEResponse()
  {
    
  }

  ////////////////////////////////////////////////////////////////////////  
  ////////////////////////////////////////////////////////////////////////                         


  ///< Reconfigure parameter set
  void PEResponse::reconfigure( const fhicl::ParameterSet& pset)
  {
    fTrkLabel     = pset.get< std::string >("TrkLabel");
    fExpLabel     = pset.get< std::string >("ExpLabel");
    fCellHitLabel = pset.get< std::string >("CellHitLabel");
    fPCHitLabel   = pset.get< std::string >("PCHitLabel");
    fQualXYLabel  = pset.get< std::string >("QualXYLabel");
    fQualZLabel   = pset.get< std::string >("QualZLabel");
    fQualAvgLabel = pset.get< std::string >("QualAvgLabel");
    fQualBadLabel = pset.get< std::string >("QualBadLabel");
    fFillHitTree  = pset.get< bool        >("FillHitTree");
  }


  ////////////////////////////////////////////////////////////////////////  
  ////////////////////////////////////////////////////////////////////////                         

  ///< beginJob function
  void PEResponse::beginJob()
  {
    art::ServiceHandle< art::TFileService > tfs;
    
    // declare all the histogram necessary 

    hall = tfs->make<TH1F>("hall", "All W; PE; Number of Hits", 1000, 0, 1000);
    hexp = tfs->make<TH1F>("hexp", "Total Exposure;;",1,0,1);

    
    for( int i = 0; i < 15; i++){

      char *hxname = new char[10];
      sprintf(hxname, "hx%d", i);
      char *hyname = new char[10];
      sprintf(hyname, "hy%d", i);

      char *htitle = new char[100];
      int offset = -750;
      int lower  = (offset+(i*100));
      int upper  = offset+(100+(i*100));
      sprintf(htitle, "%d cm < W < %d cm ; PE; Number of Hits", lower, upper );
      
      hx[i]    = tfs->make<TH1F>(hxname, htitle, 1000, 0, 1000);
      hy[i]    = tfs->make<TH1F>(hyname, htitle, 1000, 0, 1000);
    }

    
    fHitTree = tfs->make<TTree>("HitTree","Hit Tree");
    // define branches of the tree
    fHitTree->Branch("run",        &fRun);
    fHitTree->Branch("subrun",     &fSubRun);
    fHitTree->Branch("event",      &fEvent);
    fHitTree->Branch("plane",      &fPlane);
    fHitTree->Branch("cell",       &fCell);
    fHitTree->Branch("diblock",    &fDiblock);
    fHitTree->Branch("dcm",        &fDcm);
    fHitTree->Branch("apd",        &fApd);
    fHitTree->Branch("pixel",      &fPixel);
    fHitTree->Branch("view",       &fView);
    fHitTree->Branch("pe",         &fPe);
    fHitTree->Branch("w",          &fW);
    fHitTree->Branch("path",       &fPath);
    fHitTree->Branch("tns",        &fTns);
    fHitTree->Branch("readdist",   &fReadDist);
    fHitTree->Branch("flightdist", &fFlightDist);
    fHitTree->Branch("isxyadj",    &fIsXYAdj);
    fHitTree->Branch("iszadj",     &fIsZAdj);
    fHitTree->Branch("isavg",      &fIsAvg);
    fHitTree->Branch("isbadtrk",   &fIsBadTrk);
    fHitTree->Branch("isbad",      &fIsBad);

  }// end of beginJob 


  ////////////////////////////////////////////////////////////////////////  
  ////////////////////////////////////////////////////////////////////////                         
  ///< beginRun function
  void PEResponse::beginRun( art::Run const&r )
  {
    art::ServiceHandle< art::TFileService > tfs;
    int run = r.run();
    daqchannelmap::DAQChannelMap* cmap = daqchannelmap::DAQChannelMap::getInstance(fGeom->DetId());
    novadaq::cnv::DetId           det = fGeom->DetId();
    nova::dbi::RunHistory         runhist(run,det);

    fRunHist = runhist;
    fRunHist.LoadPixelInfo();

    char *htitle = new char[100];
    sprintf(htitle, "Threshold Map For Run %d;Plane;Cell", run);

    char *hname = new char[20];
    sprintf(hname, "thresh-%d", run);
    hthresh = tfs->make<TH2F>(hname, htitle, 1792, -0.5, 1791.5, 384, -0.5, 383.5);     

    int ndb = fRunHist.NDiBlocks();

    for( int idb = 0; idb < ndb; ++idb){
      nova::dbi::RunHistory::DiBlock thisdb = fRunHist.GetDiBlock( idb );
      std::vector< nova::dbi::RunHistory::DCM > dcms = thisdb.dcm;
      int ndcm = dcms.size();

      for( int idcm = 0; idcm < ndcm; ++idcm){
        nova::dbi::RunHistory::DCM thisdcm = dcms[ idcm];
        std::vector< nova::dbi::RunHistory::FEB > febs = thisdcm.feb;
        int nfeb = febs.size();
        
        for(int ifeb =0; ifeb < nfeb; ++ifeb){
          nova::dbi::RunHistory::FEB thisfeb = febs[ ifeb];
          std::vector< int16_t > thresh = thisfeb.pixelThresh;
          int npix = thresh.size();
          
          for(int ipix = 0; ipix < npix; ++ipix){
            
            uint32_t daqchan = cmap->encodeDChan( det, thisdb.num, idcm+1, ifeb, ipix);
            uint32_t logchan = cmap->encodeLChan(daqchan);
            int cell  = cmap->getCell( logchan );
            int plane = cmap->getPlane( logchan );
            // If threshold is 4095, the pixel was masked out
            // Don't fill the histogram for it, otherwise
            // it dominates the colz scale.
            if (thresh[ipix] != 4095)
              hthresh->Fill( plane, cell, thresh[ipix]);                                                        
            
          }// end loop over pixels                                                                                                   
        }// end loop over febs                                                                                                       
      }// end loop over dcms                                                                                                        
    }// end loop over diblocks                                

    fThreshHists.push_back( hthresh );
  }// end of beginRun

  ////////////////////////////////////////////////////////////////////////  
  ////////////////////////////////////////////////////////////////////////                         
  
  ///< analyze function
  void PEResponse::analyze( art::Event const& evt )
  {
    SetDefault();
    fSubRun = evt.subRun();
    fEvent  = evt.event();
    fRun    = evt.run();


    art::ServiceHandle< cmap::CMap >            cmap;
    art::ServiceHandle< chaninfo::BadChanList > badc;

    art::Handle< std::vector< rb::Track > > trkCol;
    evt.getByLabel( fTrkLabel, trkCol );


    art::FindManyP<caldp::PCHit> fmpcxy(trkCol, evt, art::InputTag(fPCHitLabel, fQualXYLabel) );
    art::FindManyP<caldp::PCHit> fmpcz (trkCol, evt, art::InputTag(fPCHitLabel, fQualZLabel) );
    art::FindManyP<caldp::PCHit> fmpcav(trkCol, evt, art::InputTag(fPCHitLabel, fQualAvgLabel) );
    art::FindManyP<caldp::PCHit> fmpcbt(trkCol, evt, art::InputTag(fPCHitLabel, fQualBadLabel) );

    int nTrks = trkCol->size();

    
    for( int iTrk = 0; iTrk < nTrks; ++iTrk){
      
      art::Ptr< rb::Track > track(trkCol, iTrk);
      // grab  the pc hits associated with this track
      if ((track->NXCell()==0) || (track->NYCell()==0)) continue;

      // get the pc hits associated with the cosmic track                                                              
      std::vector<art::Ptr<caldp::PCHit> > pcxy = fmpcxy.at(iTrk);
      std::vector<art::Ptr<caldp::PCHit> > pcz  = fmpcz.at(iTrk);
      std::vector<art::Ptr<caldp::PCHit> > pcav = fmpcav.at(iTrk);
      std::vector<art::Ptr<caldp::PCHit> > pcbt = fmpcbt.at(iTrk);

      fTrkPlanes = track->ExtentPlane();
      // Now we will loop over all the pchits of different types
      fIsXYAdj   = true;
      fIsZAdj    = false;
      fIsAvg     = false;
      fIsBadTrk  = false;
      
      
      for( auto const& iPc: pcxy ){
        fW         = iPc->W();
        fPath      = iPc->Path();
        fTns       = iPc->TNS();
        fReadDist  = iPc->ReadDist();
        fFlightDist= iPc->FlightLen();
        fPlane     = iPc->Plane();
        fCell      = iPc->Cell();
        fDcm       = iPc->DCM();
        fDiblock   = iPc->Diblock();
        fApd       = iPc->APD();
        fPixel     = iPc->Pixel();
        fView      = iPc->View();
        fIsBad     = badc->IsBad( fPlane, fCell);
        fPe        = iPc->PE();
        
        int offset = -750;
        if( fPlane < 128 && fIsXYAdj && fTrkPlanes > 20){

          for( int iH = 0; iH < 15; ++iH){
            hall->Fill(fPe);
            if( fW > offset+(iH*100) && fW < offset+(100+(iH*100)) ){
              if( fView== geo::kX )
                hx[iH]->Fill( fPe);
              if( fView== geo::kY )
                hy[iH]->Fill(fPe);
            }
          }// end loop over iH
        }

        fHitTree->Fill();
      }// end loop over pcxy
      

      fIsXYAdj   = false;
      fIsZAdj    = true;
      fIsAvg     = false;
      fIsBadTrk  = false;
      for( auto const& iPc: pcz ){
        fW         = iPc->W();
        fPath      = iPc->Path();
        fTns       = iPc->TNS();
        fReadDist  = iPc->ReadDist();
        fFlightDist= iPc->FlightLen();
        fPlane     = iPc->Plane();
        fCell      = iPc->Cell();
        fDcm       = iPc->DCM();
        fDiblock   = iPc->Diblock();
        fApd       = iPc->APD();
        fPixel     = iPc->Pixel();
        fView      = iPc->View();
        fIsBad     = badc->IsBad( fPlane, fCell);
        fPe        = iPc->PE();
        fHitTree->Fill();
      }// end loop over pcz
      

      fIsXYAdj   = false;
      fIsZAdj    = false;
      fIsAvg     = true;
      fIsBadTrk  = false;
      for( auto const& iPc: pcav ){
        fW         = iPc->W();
        fPath      = iPc->Path();
        fTns       = iPc->TNS();
        fReadDist  = iPc->ReadDist();
        fFlightDist= iPc->FlightLen();
        fPlane     = iPc->Plane();
        fCell      = iPc->Cell();
        fDcm       = iPc->DCM();
        fDiblock   = iPc->Diblock();
        fApd       = iPc->APD();
        fPixel     = iPc->Pixel();
        fView      = iPc->View();
        fIsBad     = badc->IsBad( fPlane, fCell);
        fPe        = iPc->PE();
        fHitTree->Fill();
      }// end loop over pcav
      

      fIsXYAdj   = false;
      fIsZAdj    = false;
      fIsAvg     = false;
      fIsBadTrk  = true;
      for( auto const& iPc: pcbt ){
        fW         = iPc->W();
        fPath      = iPc->Path();
        fTns       = iPc->TNS();
        fReadDist  = iPc->ReadDist();
        fFlightDist= iPc->FlightLen();
        fPlane     = iPc->Plane();
        fCell      = iPc->Cell();
        fDcm       = iPc->DCM();
        fDiblock   = iPc->Diblock();
        fApd       = iPc->APD();
        fPixel     = iPc->Pixel();
        fView      = iPc->View();
        fIsBad     = badc->IsBad( fPlane, fCell);
        fPe        = iPc->PE();
        fHitTree->Fill();
      }// end loop over pcbt
      
    }// end loop over tracks      

  }// end of analyze

  void PEResponse::SetDefault(){

    fIsBad   = false;
    fIsXYAdj = false;
    fIsZAdj  = false;
    fIsAvg   = false;
    fIsBadTrk= false;
    fPe      = -5;
    fW       = -5;
    fPath    = -5;
    fTns     = -5;
    fReadDist= -5;
    fEvent   = -5;
    fRun     = -5;
    fSubRun  = -5;
    fPlane   = -5;
    fCell    = -5;
    fDcm     = -5;
    fDiblock = -5;
    fApd     = -5;
    fPixel   = -5;
    fView    = -5;
    fFlightDist= -5;
    fTrkPlanes = -5;
  }// end of SetDefault


  void PEResponse::endSubRun(const art::SubRun& sr)
  {
    art::Handle< sumdata::CosmicExposure  > exp;
    sr.getByLabel(fExpLabel, exp);
    if( exp.failedToGet())
      return;
    hexp->Fill( 0.5, exp->totlivetime );
  }
  
  DEFINE_ART_MODULE(PEResponse)

}// end of namespace calib