CosmicMetrics_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////
/// \brief   Module that plots metrics from reconstructed cosmic ray data
/// \author  gsdavies@iastate.edu
////////////////////////////////////////////////////////////////////////////

#include <string>

// ROOT includes
#include "TH1D.h"
#include "TVector3.h"

// 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/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// NOvA includes
#include "GeometryObjects/CellGeo.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/PlaneGeo.h"
#include "NovaDatabase/Table.h"
#include "RawData/RawTrigger.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Track.h"
#include "Utilities/func/MathUtil.h"

static const double ns = 1.0E-9;

namespace comi {
  // A module to plot cosmic metrics from reconstructed cosmic ray muons
  class CosmicMetrics : public art::EDAnalyzer {
  public:
    explicit CosmicMetrics(fhicl::ParameterSet const& pset);
    virtual ~CosmicMetrics();

    void analyze(art::Event const& evt);
    void reconfigure(const fhicl::ParameterSet& pset);

    void beginJob();
    void endJob();
        
  private:
  
    bool IsEntering(double x, double y, double z); ///< yes/no is "track" entering within a buffer window
    std::string   fRawDataLabel;      ///< Where to find RawData
    std::string   fCellHitsListLabel; ///< Where to find CellHits
    std::string   fTrackLabel;        ///< Where to find Tracks
    std::string   fSoftwareRelease;   ///< Software release
    
    double        fTrkLengthCut;      ///< Minimum length of track to accept [cm]
    double        fMinTDC;            ///< Minimum TDC value to accept track
    double        fMaxTDC;            ///< Maximum TDC value to accept track
    unsigned int  fMuonNhit;          ///< Required # hits on muon track
    unsigned int  fMuonNhitX;         ///< # required in x view
    unsigned int  fMuonNhitY;         ///< # required in y view
    double        fDetEdgeBuffer;     ///< Detector edge buffer zone for entering cosmic ray track decision [cm]

    unsigned int  fNevents;           ///< total number of events counter
    int           fNTracks;           ///< total number of reconstructed tracks counter, using time cuts
    int           fRun;               ///< Run number
    int           fSubRun;            ///< Subrun number
    unsigned int  fTriggerType;
    int           fWriteToDB;   

    //......................................................................
    
    TH1D *fTrackLength;               ///< length of track [cm]
    TH1D *fCosCosmic;                 ///< cosine theta of cosmic track with respect to zenith
    TH1D *fCosmicRate;                ///< histogram holding single value for cosmic rate at end of job [Hz]
    TH1D *fNEvents;                   ///< total number of events
    TH1D *fTNS[2];                    ///< distribution of XZ/YZ CellHit times in nano seconds
    
  };
}

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

namespace comi
{  
  //...................................................................... 
  CosmicMetrics::CosmicMetrics(fhicl::ParameterSet const& pset):
    EDAnalyzer(pset),
    fRawDataLabel     (pset.get< std::string >("RawDataLabel")                  ),
    fCellHitsListLabel(pset.get< std::string >("CellHitsListLabel")             ),
    fTrackLabel       (pset.get< std::string >("TrackLabel")                    ),
    fTrkLengthCut     (pset.get< double      >("TrkLengthCut")                  ),
    fMinTDC           (pset.get< double      >("MinTDC")                        ),
    fMaxTDC           (pset.get< double      >("MaxTDC")                        ),
    fMuonNhit         (pset.get< unsigned int>("MuonNhit")                      ),
    fMuonNhitX        (pset.get< unsigned int>("MuonNhitX")                     ),
    fMuonNhitY        (pset.get< unsigned int>("MuonNhitY")                     ),
    fDetEdgeBuffer    (pset.get< double      >("DetEdgeBuffer")                 ),
    fNevents(0), fNTracks(0), fRun(-999), fSubRun(-999), fTriggerType(-999)
  {
    reconfigure(pset);    
  }
  
  //..................................................................................
  
  CosmicMetrics::~CosmicMetrics()
  {
  }
  
  //......................................................................
  
  void CosmicMetrics::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    
    fNEvents        = tfs->make<TH1D>("fNEvents",";;Number of events",2,0,2);
    fTrackLength    = tfs->make<TH1D>("trackLength",   ";Track Length (cm);"  , 100, 0., 1000.);
    fCosCosmic      = tfs->make<TH1D>("coscosmic",";Cos #theta wrt zenith",1000,0.,1.);
    fCosmicRate     = tfs->make<TH1D>("cosmicrate",";Cosmic Rate (Hz)",1,0.,100000);
    
    for(int view = geo::kX; view <= geo::kY; ++view){   
      const TString viewStr = (view == geo::kX) ? "X" : "Y";
      fTNS[view]            = tfs->make<TH1D>("HitTNS"+viewStr, ";Time (ns);hits", 70, -100000.0, 600000.);
    }       
  }  
  
  //......................................................................
  
  ///
  /// \brief The analysis of CosmicTrack reconstructed tracks in roder to
  ///        calculate cosmic ray muon rates, track lengths and more 
  ///        cosmic ray metrics to be monitored for commissioning
  ///
  /// \param evt - Read-only access to the event data
  ///
  
  void CosmicMetrics::analyze(art::Event const& evt)
  { 
  
    art::ServiceHandle<geo::Geometry> geom;   

    ++fNevents; //Updating counter of number of events in sub-run
    if(fRun < 0)    fRun = evt.run();
    if(fSubRun < 0) fSubRun = evt.subRun();
    fNEvents->Fill(1);
    
    art::Handle< std::vector<rawdata::RawTrigger> > trigv;
    evt.getByLabel(fRawDataLabel, trigv);
    const rawdata::RawTrigger& trig = (*trigv)[0];      
    fTriggerType = trig.fTriggerMask_TriggerType;  
    
    // now get the tracks that were made by CosmicTrack
    // and check that there are tracks in the event
    art::Handle< std::vector<rb::Track> > trks;
    evt.getByLabel(fTrackLabel,trks);
    if (trks->size()==0) {
      mf::LogWarning("CosmicMetricsWarn") << "There are " << trks->size() 
                                          << " Tracks in this event, skip";
      return;
    }
    LOG_DEBUG("CosmicMetrics") << "There are " << trks->size() << " Tracks in this event";
    
    const unsigned int trkMax = trks->size(); // maximum number of tracks in event
    double cosCosmic   = 0.;
    int numtracks      = 0 ;
    
    // Loop over tracks
    for(unsigned int trkIdx = 0; trkIdx < trkMax; ++trkIdx){
      
      bool isgoodXZtrack  = false;
      bool isgoodYZtrack  = false;
      
      art::Ptr<rb::Track> track(trks, trkIdx); // ART-pointer to each track in event
      
      // Get vectors for the end points of the track
      TVector3 start(track->Start());
      TVector3 end(track->Stop());
      TVector3 dir(track->Dir());
      TVector3 length = end - start;
      
      if(end.Y() > start.Y()) std::swap(start, end); // force start point in Y to be greater than end point in Y for cosmics

      // Does this look like an entering track?         
      if (this->IsEntering(start.X(),start.Y(),start.Z()) == false) continue;
  
      const unsigned int hitMax = track->NCell();  // maximum number of CellHits on track
      const unsigned int hitX   = track->NXCell(); // maximum number of CellHits on track in XZ
      const unsigned int hitY   = track->NYCell(); // maximum number of CellHits on track in YZ
      
      // Do NOT select tracks if max number of hits or XZ- YZ-hits less than configurable values
      if ( (hitMax < fMuonNhit) || (hitX < fMuonNhitX) || (hitY < fMuonNhitY) ) continue;
      assert((track->NXCell() > 0) && (track->NYCell() > 0));
      
      double TNS[2]        = {0.}; // TNS() value in each view
      double TNStotal[2]   = {0.}; // total TNS() in each view
      double TNSaverage[2] = {0.}; // average TNS() in each view
      
      // Only select tracks that cross more than fTrkLengthCut [cm] in Z 
      if(fabs(end.Z() - start.Z()) > fTrkLengthCut){
      
        fTrackLength->Fill(length.Mag());
        
        // Loop through both views and grab every CellHit in this track
        for(int view = geo::kX; view <= geo::kY; ++view){       
          // C++ is a bit strict about the conversions. Save typing.    
          const geo::View_t geoview = geo::View_t(view);
          const int vhitMax =  track->NCell(geoview);
        
          for(int hitIdx = 0; hitIdx < vhitMax; ++hitIdx){
          
            art::Ptr<rb::CellHit> chit = track->Cell(geoview, hitIdx);
            fTNS[view]->Fill(chit->TNS());
            TNS[view]      = chit->TNS();
            TNStotal[view]+=TNS[view];
          } // hitIdx
          
          TNSaverage[view] = TNStotal[view]/track->NCell(geoview); // average TNS() CellHits in each view   
        } // view
      } // track length cut
      
      if(TNSaverage[geo::kX] >= fMinTDC && TNSaverage[geo::kX] <= fMaxTDC){
        isgoodXZtrack = true;
      }
      if(TNSaverage[geo::kY] >= fMinTDC && TNSaverage[geo::kY] <= fMaxTDC){
        isgoodYZtrack = true;
      }  
      
      // Only count 3D tracks within configurable window size and with matching XZ/YZ 2D tracks
      if(isgoodXZtrack == true && isgoodYZtrack == true){
        numtracks++;
      }  
      
      // Cosmic track angular distribution with respect to zenith
      cosCosmic = TMath::Abs(-dir.Y()/util::pythag(dir.X(),dir.Y(),dir.Z()));
      fCosCosmic->Fill(cosCosmic);
      
    } // loop over tracks  
    fNTracks += numtracks;
    
    return;
  } 
  
  //......................................................................
    
  void CosmicMetrics::reconfigure(const fhicl::ParameterSet& pset)
  {
    fTrkLengthCut = pset.get< double >("TrkLengthCut"  );
    fMinTDC       = pset.get< double >("MinTDC"        );
    fMaxTDC       = pset.get< double >("MaxTDC"        );
    fDetEdgeBuffer= pset.get< double >("DetEdgeBuffer");
    fWriteToDB    = pset.get< int    >("WriteToDB");
  }  
           
  //......................................................................          
  
  void CosmicMetrics::endJob()
  {


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

    float fLivetime   = fNevents*((fMaxTDC-fMinTDC)*ns);
    float fCosmicrate = fNTracks/(fNevents*((fMaxTDC-fMinTDC)*ns));
    std::cout<< "There are: " << fNevents << " events in this job.\n";
    std::cout<< "There are: " << fNTracks << " tracks selected in this job.\n";
    std::cout<< "The livetime is " << fLivetime << "!\n";
    std::cout<< "The cosmic rate is therefore " << fCosmicrate << " Hz \n";
    std::cout<< "tracklengthmean="<<fTrackLength->GetMean(1)<<std::endl;
    std::cout<< "tracklengthrms=" <<fTrackLength->GetRMS(1) <<std::endl;
    fCosmicRate->Fill(fCosmicrate);
    
    //Get software release
    fSoftwareRelease = getenv("SRT_BASE_RELEASE");
    
    if (fWriteToDB) { 

      nova::database::Table* t = new nova::database::Table(std::string("/Commissioning/tables/CosmicMetrics.xml"));     

      t->SetDetector(geom->DetId());
      t->SetValidityRange("detector", novadaq::cnv::DetInfo::GetName(geom->DetId()));
      t->SetValidityRange("run",fRun);
      t->SetValidityRange("subrun",fSubRun);
      t->SetValidityRange("triggertype",fTriggerType);
      t->SetValidityRange("softwarerelease", fSoftwareRelease);
      t->LoadFromDB();
      if ( t->NRow() ){
        std::cerr<< " Subrun metrics are already in the database "<<std::endl;
      }
      else{
        std::cout<< " Writting subrun metrics to the database "<<std::endl;
        
        nova::database::Row* r = t->NewRow();

        r->Set("detector",                   novadaq::cnv::DetInfo::GetName(geom->DetId()));
        r->Set("run",                        fRun);
        r->Set("subrun",                     fSubRun);
        r->Set("triggertype",                fTriggerType);
        r->Set("nevents",                    fNevents);
        r->Set("ntracks",                    fNTracks);
        r->Set("livetime",                   fLivetime);
        r->Set("cosmicrate",                 fCosmicrate);
        r->Set("tracklengthmean",            fTrackLength->GetMean(1));
        r->Set("tracklengthrms",             fTrackLength->GetRMS(1));
        r->Set("softwarerelease",            fSoftwareRelease);

        t->AddRow(r);   
        delete r;
        t->WriteToDB();
      }
    }        
  }

  //......................................................................          
              
  bool CosmicMetrics::IsEntering(double x, double y, double z)
  {     
    art::ServiceHandle<geo::Geometry> geom;     
    if (y > geom->DetHalfHeight()-fDetEdgeBuffer)      return true;
    if (fabs(x) > geom->DetHalfWidth()-fDetEdgeBuffer) return true;
    if (z < fDetEdgeBuffer)                            return true;     
    if (z > geom->DetLength()-fDetEdgeBuffer)          return true;
    return false;       
  } 
} // end namespace comi

//////////////////////////////////////////////////////////////////////////
namespace comi 
{
  DEFINE_ART_MODULE(CosmicMetrics)
}
////////////////////////////////////////////////////////////////////////