HitEfficiency_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       HitEfficiency
// Module Type: Analyzer
// File:        HitEfficiency_module.cc
//
////////////////////////////////////////////////////////////////////////

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

// novasoft stuff
#include "RecoBase/Track.h"
#include "RecoBase/CellHit.h"
#include "Calibrator/Calibrator.h"
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "GeometryObjects/CellGeo.h"
#include "Utilities/func/MathUtil.h"
#include "CMap/service/CMapService.h"
#include "RunHistory/RunHistory.h"
#include "ChannelInfo/BadChanList.h"


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

// ROOT
#include "TVector3.h"
#include "TH2F.h"


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

  private:

    std::string fTrackLabel;
    double fNTrkHits;
    double fNPlanes;
    double fMaxPlane;
    double fDEdge;

    unsigned int fplanemax; // maximum plane to calculate efficiencies on
    double zmaxinstr; // maximum instrumented z position
    double zmininstr; // minimum instrumented z position

    //int fDet; // Detector ID
    nova::dbi::RunHistory frunhist;
    
    TH2F* hitsVwVlenx;
    TH2F* cellsVwVlenx;
    
    TH2F* hitsVwVleny;
    TH2F* cellsVwVleny;

    TH2F* hitsVapddistVlenx;
    TH2F* cellsVapddistVlenx;
    
    TH2F* hitsVapddistVleny;
    TH2F* cellsVapddistVleny;
    
    // efficiency plots
    TH2F* hitsVwcell;
    TH2F* cellsVwcell;

    // plots after applying track quality cuts plotted in relation to distance from apd, not w
    TH2F* hitsVapddistcell;
    TH2F* cellsVapddistcell;

    // plots after applying path length cut
    TH2F* hitsVwcelllencut;
    TH2F* cellsVwcelllencut;

    TH2F* hitsVapddistcelllencut;
    TH2F* cellsVapddistcelllencut;


    // Histograms to keep track of the number of events used
    TH1F* NEvents;
    int fnevents;

    // define length cuts
    double minpath = 6.0;
    double maxpath = 7.0;

  };
  
  
  HitEfficiency::HitEfficiency(fhicl::ParameterSet const & p)
    : EDAnalyzer(p)
  {
    this->reconfigure(p);
  }
  
  void HitEfficiency::reconfigure(const fhicl::ParameterSet& pset)
  {
    fTrackLabel = pset.get< std::string >("TrackLabel");
    fNTrkHits = pset.get< double >("NTrkHits");
    fNPlanes = pset.get< double >("NPlanes");
    fMaxPlane = pset.get< double >("MaxPlane");
    fDEdge = pset.get< double >("DEdge");
  }


  HitEfficiency::~HitEfficiency()
  {
    // Clean up dynamic memory and other resources here.
  }

  void HitEfficiency::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    
    fplanemax = 896;

    fnevents = 0;

    zmaxinstr = 0.0; // maximum instrumented z position
    zmininstr = 0.0;
    
    // Define my histograms binning here
    double nwbins = 32;
    double wmin = -800;
    double wmax = 800;

    double napdbins = nwbins;
    double apdmin = 0.0;
    double apdmax = 1600;
    
    double nlbins = 25;
    double lmin = 0.0;
    double lmax = 10.0;
    
    int ncellbins = 384*fplanemax+1;
    double cellmin = -0.5;
    double cellmax = 384*(double)fplanemax+0.5;

    // efficiency plots 
    hitsVwVlenx = tfs->make<TH2F>("hitsVwVlenx",";W (cm);Length (cm)",
                                  nwbins,wmin,wmax,nlbins,lmin,lmax);
    cellsVwVlenx = tfs->make<TH2F>("cellsVwVlenx",";W (cm);Length (cm)",
                                   nwbins,wmin,wmax,nlbins,lmin,lmax);

    
    hitsVwVleny = tfs->make<TH2F>("hitsVwVleny",";W (cm);Length (cm)",
                                  nwbins,wmin,wmax,nlbins,lmin,lmax);
    cellsVwVleny = tfs->make<TH2F>("cellsVwVleny",";W (cm);Length (cm)",
                                   nwbins,wmin,wmax,nlbins,lmin,lmax);


    hitsVapddistVlenx = tfs->make<TH2F>("hitsVapddistVlenx",";Distance from APD (cm);Length (cm)",
                                        napdbins,apdmin,apdmax,nlbins,lmin,lmax);
    cellsVapddistVlenx = tfs->make<TH2F>("cellsVapddistVlenx",";Distance from APD (cm);Length (cm)",
                                         napdbins,apdmin,apdmax,nlbins,lmin,lmax);

    
    hitsVapddistVleny = tfs->make<TH2F>("hitsVapddisVleny",";Distance from APD (cm);Length (cm)",
                                        napdbins,apdmin,apdmax,nlbins,lmin,lmax);
    cellsVapddistVleny = tfs->make<TH2F>("cellsVapddisVleny",";Distance from APD (cm);Length (cm)",
                                         napdbins,apdmin,apdmax,nlbins,lmin,lmax);

    
    // efficiency plots as a function of w
    hitsVwcell = tfs->make<TH2F>("hitsVwcell",";Cell;W (cm)",
                                 ncellbins,cellmin,cellmax,nwbins,wmin,wmax);
    cellsVwcell = tfs->make<TH2F>("cellsVwcell",";Cell;W (cm)",
                                  ncellbins,cellmin,cellmax,nwbins,wmin,wmax);


    // plots as a function of apd distance
    hitsVapddistcell = tfs->make<TH2F>("hitsVapddistcell",";Cell;Distance from APD (cm)",
                                       ncellbins,cellmin,cellmax,napdbins,apdmin,apdmax);
    cellsVapddistcell = tfs->make<TH2F>("cellsVapddistcell",";Cell;Distance from APD (cm)",
                                        ncellbins,cellmin,cellmax,napdbins,apdmin,apdmax);

    
    // plots after applying path length cut
    hitsVwcelllencut = tfs->make<TH2F>("hitsVwcelllencut",";Cell;W (cm)",
                                       ncellbins,cellmin,cellmax,nwbins,wmin,wmax);
    cellsVwcelllencut = tfs->make<TH2F>("cellsVwcelllencut",";Cell;W (cm)",
                                        ncellbins,cellmin,cellmax,nwbins,wmin,wmax);

        
    hitsVapddistcelllencut = tfs->make<TH2F>("hitsVapddistcelllencut",";Cell;Distance from APD (cm)",
                                       ncellbins,cellmin,cellmax,napdbins,apdmin,apdmax);
    cellsVapddistcelllencut = tfs->make<TH2F>("cellsVapddistcelllencut",";Cell;Distance from APD (cm)",
                                        ncellbins,cellmin,cellmax,napdbins,apdmin,apdmax);

    NEvents = tfs->make<TH1F>("NEvents",";NEvents;",1,-0.5,0.5);


  }
  
  void HitEfficiency::beginRun(art::Run const& r)
  {
    // make a geometry service handle
    art::ServiceHandle<geo::Geometry> geom;
    int planemin = 0;
    int planemax = geom->NPlanes()-1;
    if(fMaxPlane >= planemin && fMaxPlane <= planemax){
      fplanemax = fMaxPlane;
    }
    else{ fplanemax = planemax; }

    art::ServiceHandle<geo::LiveGeometry> lgeom;
    double zfront1,zback1,zfront2,zback2;
    int ndets = lgeom->InstrumentedDetEnds(zfront1,zback1,zfront2,zback2);
    
    // take the most downstream zfront and zback
    if(ndets > 0){
      zmaxinstr = zback1;
      zmininstr = zfront1; 
    }

    // leaving this disabled for now
//     fDet = geom->DetId();
//     // get the run history
//     run = r.run();
//     nova::dbi::RunHistory runhist(run,fDet);
//     frunhist = runhist;
//     frunhist.LoadDAQRunHistory();
//     frunhist.LoadPixelInfo();

  }
  
  void HitEfficiency::analyze(art::Event const& e)
  {

    ++fnevents;

    // make a geometry service handle
    art::ServiceHandle<geo::Geometry> geom;
    // make a cmap service handle
    art::ServiceHandle<cmap::CMap> cmap;


    //get the cosmic tracks associated with each event
    art::Handle<std::vector<rb::Track> > trackcol;
    e.getByLabel(fTrackLabel,trackcol);
    art::PtrVector<rb::Track> tracks;
    for(unsigned int i = 0; i<trackcol->size(); ++i){
      art::Ptr<rb::Track>track(trackcol,i);
      tracks.push_back(track);
    }

    // loop over the tracks
    for(unsigned int itr = 0; itr < tracks.size(); ++itr){
      art::Ptr<rb::Track> track = tracks[itr];

      if(track->View() != geo::kXorY){ continue; }

      // get basic track information here to store in tree
      int nxhit = track->NXCell();
      int nyhit = track->NYCell();    
      if(nxhit < fNTrkHits || nyhit < fNTrkHits){ continue; }

      TVector3 dir = track->Dir();
      double dx = dir.X();
      double dy = dir.Y();
      double dz = dir.Z();
      if( dx > 1.0 || dx < -1.0 ||
          dy > 1.0 || dy < -1.0 ||
          dz > 1.0 || dz < -1.0){ continue; }
            
      // find the minimum, maximum planes
      unsigned int minplane = track->MinPlane();
      unsigned int maxplane = track->MaxPlane();

      if((maxplane - minplane) < fNPlanes){ continue; }
      if(maxplane > fplanemax || minplane > fplanemax){ continue; }

      // get the track start/stop point info
      TVector3 start = track->Start();
      TVector3 stop = track->Stop();
      double ix = start.X();
      double iy = start.Y();
      double iz = start.Z();
      double fx = stop.X();
      double fy = stop.Y();
      double fz = stop.Z();

      double zmin = zmininstr+fDEdge;
      double zmax = zmaxinstr-fDEdge;
      double xmax = geom->DetHalfWidth() - fDEdge;
      double ymax = geom->DetHalfHeight() - fDEdge;

      // through going track cut
      if(fabs(ix) < xmax && fabs(iy) < ymax && iz > zmin && iz < zmax){ continue; } //start
      if(fabs(fx) < xmax && fabs(fy) < ymax && fz > zmin && fz < zmax){ continue; } // stop

      // get the hits that belong to this track
      art::PtrVector<rb::CellHit> trkhits = track->AllCells();

      // sort the hits by plane
      rb::SortByPlane(trkhits);

      // get the list of cells that the track passes through 
      std::vector<geo::OfflineChan> allcellsontrack, ycellsontrack;
      geom->CountCellsOnLine(start,stop,allcellsontrack,ycellsontrack);

      for(unsigned int i = 0; i < ycellsontrack.size(); ++i){
        allcellsontrack.push_back(ycellsontrack[i]);
      }

      // make a list of tricells
      std::vector<geo::OfflineChan> tricells;
      for(unsigned int icell = 0; icell < allcellsontrack.size(); ++icell){
        // check if the z position of the cell falls in the instrumented region
        double zcent[3];
        geom->Plane(allcellsontrack[icell].Plane())->Cell(allcellsontrack[icell].Cell())->GetCenter(zcent);
        if(zcent[2] < zmininstr || zcent[2] > zmaxinstr){ continue; }
        // find the maximum and minimum cell hit on this plane
        int plane = allcellsontrack[icell].Plane();
        int cellmin = allcellsontrack[icell].Cell()-1;
        int cellmax = allcellsontrack[icell].Cell()+1;
        bool topok = false;
        bool botok = false;
        for(unsigned int ihit = 0; ihit < trkhits.size(); ++ihit){
          if(trkhits[ihit]->Plane() < plane){ continue; }
          if(trkhits[ihit]->Plane() > plane){ break; }
          // now have a hit on the same plane
          if(trkhits[ihit]->Cell() == cellmax){ topok = true; }
          else if(trkhits[ihit]->Cell() == cellmin){ botok = true; }
        }
        if(topok && botok){ tricells.push_back(allcellsontrack[icell]); }
      }

      // loop over all the cells this track passes through
      for(unsigned int icell = 0; icell < tricells.size(); ++icell){

        unsigned int cell = tricells[icell].Cell();
        unsigned int plane = tricells[icell].Plane();

        // if plane is not less than plane max skip
        if(plane > fplanemax-1){ continue; }

        // check if this channel is enabled
        bool enable = true;
        // need to double check this is working, leave disabled for now
        // get dcm
        //      const daqchannelmap::lchan LogicalChannel = cmap->encodeLChan(fDet,plane,cell);
        // //   get the dcm id
        //      const uint32_t dcmid = cmap->getDCMOffline(LogicalChannel);
        // //   get the feb id
        //      const int febid = cmap->getDCMLinkOffline(LogicalChannel);
        // //   get the run history dcm object for this dcm id
        //      nova::dbi::RunHistory::DCM dcm = frunhist.GetDCM(dcmid);
        //      // get enable information for this febid
        //      enable = dcm.feb[febid].isEnabled;
        // skip if not enabled
        if(!enable){ continue; }

        // get the center of this cell
        double xyz[3];
        geom->Plane(plane)->Cell(cell)->GetCenter(xyz);
        
        double deltax,deltay,deltaz,w;
        if(geom->Plane(plane)->View() == geo::kX){
          deltax =  2.0*geom->Plane(plane)->Cell(cell)->HalfW();
          deltay = (dy/dx)*deltax;
          deltaz = (dz/dx)*deltax;
          // find the w position (y) for this cell
          w = (xyz[2]-iz)*dy/dz + iy;
        }
        else{
          deltay = 2.0*geom->Plane(plane)->Cell(cell)->HalfW();
          deltax = (dx/dy)*deltay;
          deltaz = (dz/dy)*deltay;
          // find the w position (x) for this cell
          w = (xyz[2]-iz)*dx/dz + ix;
        }

        double tripath = sqrt(deltax*deltax+deltay*deltay+deltaz*deltaz);
        double apddist = geom->DistToElectronics(w,*(geom->Plane(plane)->Cell(cell)));

        // loop over hits and see if this is a hit or not
        bool hit = false;
        for(unsigned int ihit = 0; ihit < trkhits.size(); ++ihit){
          if(trkhits[ihit]->Cell() == cell && trkhits[ihit]->Plane() == plane){
            hit = true;
            break;
          }
        }
    
        // Fill some histograms
        int cellbin = 384*plane+cell;

        if(geom->Plane(plane)->View() == geo::kX){
          cellsVwVlenx->Fill(w,tripath);
          cellsVapddistVlenx->Fill(apddist,tripath);    
        }
        else{
          cellsVwVleny->Fill(w,tripath);
          cellsVapddistVleny->Fill(apddist,tripath);
        }
        cellsVwcell->Fill(cellbin,w);
        cellsVapddistcell->Fill(cellbin,apddist);       

        if(hit){
          hitsVwcell->Fill(cellbin,w);
          hitsVapddistcell->Fill(cellbin,apddist);
          if(geom->Plane(plane)->View() == geo::kX){
            hitsVwVlenx->Fill(w,tripath);
            hitsVapddistVlenx->Fill(apddist,tripath);
          }
          else{
            hitsVwVleny->Fill(w,tripath);
            hitsVapddistVleny->Fill(apddist,tripath);
          }
        }
        
        // make plots with path length cut
        if(tripath > minpath && tripath < maxpath){
          cellsVwcelllencut->Fill(cellbin,w);
          cellsVapddistcelllencut->Fill(cellbin,apddist);
          if(hit){
            hitsVwcelllencut->Fill(cellbin,w);
            hitsVapddistcelllencut->Fill(cellbin,apddist);
          }
        }
        
      } // end of loop over tricells
            
    } // end of loop over tracks
  
  } // end of analyze
  
  void HitEfficiency::endJob()
  {
    NEvents->SetBinContent(1,fnevents);
  }
  
  DEFINE_ART_MODULE(HitEfficiency)
  
} // end of namespace