DataCheck_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////
/// \brief   Commissioning files to look at the quality of our data
/// \author  lein@physics.umn.edu, zirnstein@umn.edu
///////////////////////////////////////////////////////////////////////
#include <cmath>
#include <string> 
#include <time.h>
#include <vector>
#include <iostream>
#include <ostream>
#include <fstream>
#include <sstream>

// ROOT includes
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TF1.h"

// Framework includes
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Core/EDAnalyzer.h"   
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Principal/Handle.h"

// NOvA includes
#include "ChannelInfo/BadChanList.h"
#include "CMap/service/CMapService.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "NovaDatabase/Table.h"
#include "RawData/RawDigit.h"
#include "RawData/RawTrigger.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "TrackFit/RLFit.h"
#include "RunHistory/RunHistory.h" //....tx

namespace comi {
  
  class DataCheck : public art::EDAnalyzer {  
  public:
    explicit DataCheck(fhicl::ParameterSet const& pset);
    virtual ~DataCheck();
    
    void beginJob();
    virtual void beginSubRun(art::SubRun const& sr);
    void analyze(art::Event const& evt);  
    void reconfigure(const fhicl::ParameterSet& p);
    void endJob();      
    virtual void endSubRun(art::SubRun const& sr);
    
  private:
    
    TH1F*        fNSlices;
    TH1F*        fNHitsPerSlice;
    TH1F*        fNHitsPerSliceNoNoise;
    TH1F*        fSliceLength;
    TH1F*        fNHitsPerEvent;
    TH2F*        fPlaneVCell;
    TH2F*        fPlaneVCellRaw;
    TH2F*        fChannelMap;
    TH3F*        fBlocDcmFeb;
    TH1F*        fOccupancy;
    TH1F*        fEmptyEvents;
    TH1F*        fNoiseEvents;
    TH1F*        fNEvents;
    
    TH2F*        fDCM_temp;
    std::vector<TH2F*>   vDCMs;
    
    int          fRun;
    int          fSubRun;
    unsigned int fNevents; 
    unsigned int fNempty;
    unsigned int fNoise;
    unsigned int fTriggerType;
    unsigned int fNGoodChannels;
    unsigned int fNRawChannels;
    
    int          fDCM_FEB[4][3][64];
    int          fDCM_nFEB[4][3];
 
    time_t       fFirstEventTimeStamp;
    time_t       fFirstNonEmptyEventTimeStamp;
    time_t       fLastNonEmptyEventTimeStamp;
    time_t       fLastEventTimeStamp;
    
    int          fWriteToDB;
    int          fIsNearline;
    std::vector<float> fNHits;
    unsigned int fNDCM;
    
  private:
    std::string  fSlicerModule;     ///< Label of the slicer module in the fcl file used 
    std::string  fCalHitModule;     ///< Label of the calhit module in the fcl file used 
    std::string  fRawDigitModule;   ///< Label of the raw digit module; set to default
    std::string  fRawDataLabel;
    std::string  fSoftwareRelease;
    std::string  fReportDir;
  }; // end of class DataCheck

  // NOvA uses a 64 Mhz clock
  //static const double kUSEC_PER_TDC = (1.0E6/64.0E6);
  //static const int    kTDC_PER_USEC = 64;
  
  //......................................................................
  
  DataCheck::DataCheck(fhicl::ParameterSet const& p):
    EDAnalyzer(p),
    fRun(-999),
    fSubRun(-999),
    fNevents(0),
    fNempty(0),
    fNoise(0),
    fTriggerType(-999),
    fNGoodChannels(0),
    fNRawChannels(0),
    fFirstEventTimeStamp(0),
    fFirstNonEmptyEventTimeStamp(0),
    fLastNonEmptyEventTimeStamp(0),
    fLastEventTimeStamp(0),
    fNHits(0),
    fNDCM(0)

  {
    this->reconfigure(p);
  }
  
  //......................................................................
  
  void DataCheck::beginJob() 
  {
    art::ServiceHandle<art::TFileService> f;
    art::ServiceHandle<geo::Geometry> geom;
   //    const unsigned int maxCell = std::max(geom->Plane(0)->Ncells(), geom->Plane(1)->Ncells());
    Int_t NPlanes = 910;// geom->NPlanes();
    Int_t NCells  =  394;//std::max(geom->Plane(0)->Ncells(), geom->Plane(1)->Ncells());


    // Some event-level histograms
    
    fNSlices               = f->make<TH1F>("fNSlices",";NSlices;Events",100,0.0,100.0);
    fNHitsPerSlice         = f->make<TH1F>("fNHitsPerSlice",";NHits;Slices",100,0.0,100.0);
    fNHitsPerSliceNoNoise  = f->make<TH1F>("fNHitsPerSliceNoNoise",";NHits;Slices",100,0.0,100.0);
    fSliceLength           = f->make<TH1F>("fSliceLength",";Slice Duration; NSlices", 550,0.0, 5500);
    
    fNHitsPerEvent         = f->make<TH1F>("fNHitsPerEvent",";NHits/NGoodChannels;Events",100,0.0,0.5);
    fPlaneVCell            = f->make<TH2F>("fPlaneVCell",";Plane Number;Cell Number", NPlanes,-0.5,NPlanes-0.5, NCells,-0.5,NCells  -0.5);
    fPlaneVCellRaw         = f->make<TH2F>("fPlaneVCellRaw",";Plane Number;Cell Number",NPlanes,-0.5,NPlanes-0.5, NCells,-0.5,NCells  -0.5);
    fChannelMap            = f->make<TH2F>("fChannelMap","Channel Map;Plane Number;Cell Number",NPlanes,-0.5,NPlanes-0.5, NCells,-0.5,NCells  -0.5);

    fBlocDcmFeb            = f->make<TH3F>("fBlockDcmFeb","Channel Map;Diblock; DCM;FEB",14,-0.5, 13.5, 12,-0.5,11.5 ,64,-0.5,63.5);
    
    fOccupancy             = f->make<TH1F>("fOccupancy",";log of Occupancy;Channels",200,-3.0,0.0);
    fEmptyEvents           = f->make<TH1F>("fEmptyEvents",";;Number of empty events",2,0,2);
    fNoiseEvents           = f->make<TH1F>("fNoiseEvents",";;Number of events with noise only slice",2,0,2);
    fNEvents               = f->make<TH1F>("fNEvents",";;Number of events",2,0,2);
    
    
    for(int ndiblock=0;ndiblock<4;ndiblock++){
      for(int ndcm=0;ndcm<3;ndcm++){
        for(int nfeb=0;nfeb<64;nfeb++){
          fDCM_FEB[ndiblock][ndcm][nfeb]=0;
        }
        fDCM_nFEB[ndiblock][ndcm]=0;
        
      }
    }
  }  //End of beginJob 
  
  //......................................................................  
  DataCheck::~DataCheck() 
  { 
    delete fDCM_temp;
    vDCMs.clear();
  }
  
  //............................begin SubRun..............................
  void DataCheck::beginSubRun(art::SubRun const& sr)
  {
    art::ServiceHandle<art::TFileService> f;
    art::ServiceHandle<geo::Geometry> geom;

    novadaq::cnv::DetId detectorid = geom->DetId();
    if(fRun < 0) fRun = sr.run();
    nova::dbi::RunHistory rh(fRun, detectorid); 
    rh.LoadDAQRunHistory();

    unsigned int nDB = rh.NDiBlocks();
    const nova::dbi::RunHistory::DiBlock tempDiBlock = rh.GetDiBlock(0);
    float nFEB = (float)tempDiBlock.dcm[0].feb.size();
    
    int ic=0;
    for (unsigned int i=0; i<nDB; ++i) {
      const nova::dbi::RunHistory::DiBlock db = rh.GetDiBlock(i);
      for( unsigned int j=0;j<db.dcm.size(); j++, ++ic ){
        std::string dcmName = db.dcm[j].appname; 
        char histname[30];
        strcpy(histname, dcmName.c_str());
        vDCMs.push_back(fDCM_temp); 
        vDCMs.at(ic) = f->make<TH2F>(histname,";Event;FEB",100,0.0,38000.0,nFEB,0,nFEB);
      } 
    }
  }

  void DataCheck::endSubRun(art::SubRun const& sr)
  {
  }
  //........................................................................
  void DataCheck::analyze(art::Event const& evt)
  {

    art::ServiceHandle<cmap::CMap> cmap;

    art::ServiceHandle<geo::Geometry> geom;
    double detid = geom->DetId(); 

    art::Timestamp ts = evt.time();
    
    time_t uevt_sec = ts.timeHigh(); //utval >> 32 & 0xFFFFFFFFF;
    
    art::Handle< std::vector<rawdata::RawTrigger> > trigv;
    evt.getByLabel(fRawDataLabel, trigv);
    const rawdata::RawTrigger& trig = (*trigv)[0];
    fTriggerType = trig.fTriggerMask_TriggerType;  
    
    
    if(fRun < 0) fRun = evt.run();
    if(fSubRun < 0) fSubRun = evt.subRun();    
    
    if( fFirstEventTimeStamp == 0 ) {
      fFirstEventTimeStamp=uevt_sec;
    }
    if (fNevents > 0 ) {
      fLastEventTimeStamp=uevt_sec;
    }
    
    
    ++fNevents; //Updating counter of number of events 
    
    fNEvents->Fill(1);
    
    //-------------------------------------------------------------------------------------------------
    
    //Looking at Slices
    
    //Load the slicer list from the event
    art::Handle<std::vector<rb::Cluster> > slicecol;
    evt.getByLabel(fSlicerModule,slicecol);
    
    // put in a mf::LogWarning why this statement returns if true
    if(slicecol->empty()){
      mf::LogWarning ("No Slices")<<"No Slices in the input file";
      return;
    }
    
    //Fill histogram with number of slices
    fNSlices->Fill(slicecol->size());
   
    //Loop over slices  
    for(unsigned int iSlice = 0; iSlice < slicecol->size(); ++iSlice){ 
      const rb::Cluster &slice = slicecol->at(iSlice);  
      
      fNHitsPerSlice->Fill(slice.NCell()); //Hits per slice (with Noise Slice)      
      
      //branching for noise slice
      if(slice.IsNoise()){
        continue;
      }
      
      // Initialize the beginning and ending times to first hit times.
      if(slice.NCell() != 0){
        fSliceLength->Fill(slice.ExtentTNS());   //Fill histogram with slice length
      }
      
      fNHitsPerSliceNoNoise->Fill(slice.NCell());       //Hits per slice (without Noise Slice) 
      
    }
    
    //-------------------------------------------------------------------------------------------------
    
    //Looking at CellHits

    
    // Load hit list (CellHits)
    art::Handle< std::vector<rb::CellHit> > hitcol;
    evt.getByLabel(fCalHitModule, hitcol);
    
    //Make it an art::PtrVector
    art::PtrVector<rb::CellHit> hitlist;
    for(unsigned int i = 0; i < hitcol->size(); ++i){
      art::Ptr<rb::CellHit> hit(hitcol,i);
      hitlist.push_back(hit);
    }
    
    //Loop over CellHits
    int cell, plane;
    for(unsigned int iHit = 0; iHit < hitlist.size(); iHit++){
      
      cell =    hitlist[iHit]->Cell();
      plane =   hitlist[iHit]->Plane(); 
      
      fPlaneVCell->Fill(plane,cell);
      
    } //End of loop over CellHits
    
    //-------------------------------------------------------------------------------------------------
    
    //Looking at RawDigits
    
    // Load hit list (RawDigits)
    art::Handle< std::vector<rawdata::RawDigit> > rawhitcol;
    evt.getByLabel(fRawDigitModule, rawhitcol);
    
    //Make it an art::PtrVector
    art::PtrVector<rawdata::RawDigit> rawhitlist;
    for(unsigned int i = 0; i < rawhitcol->size(); ++i){
      art::Ptr<rawdata::RawDigit> rawhit(rawhitcol,i);
      rawhitlist.push_back(rawhit);
    }
    
    
    //Loop over RawDigits 
    
    uint32_t daq_channel = 0;
    uint32_t feb         = 0;
    uint32_t dcm         = 0;
    uint32_t diblock     = 0;
    uint32_t rawplane    = 0;
    uint32_t rawcell     = 0;
    
    for(unsigned int iRawHit = 0; iRawHit < rawhitlist.size(); iRawHit++){ 
      
      daq_channel = rawhitlist[iRawHit]->DaqChannel();
      
      feb         = cmap->Map()->getFEB(daq_channel);
      dcm         = cmap->Map()->getDCM(daq_channel);
      diblock     = cmap->Map()->getDiBlock(daq_channel);
      rawplane    = cmap->Map()->getPlane(daq_channel);
      rawcell     = cmap->Map()->getCell(daq_channel);

      fPlaneVCellRaw->Fill(rawplane,rawcell);
      fBlocDcmFeb->Fill(diblock,dcm, feb);
      if( detid==1 ){
        vDCMs.at(diblock*12+dcm-13)->Fill(fNevents,feb);
      }
      else if( detid==2 ){
        vDCMs.at(diblock*12+dcm-13)->Fill(fNevents,feb);
      }
      else if( detid==3 ){
        vDCMs.at(diblock*3+dcm-4)->Fill(fNevents,feb);
      }
      else{
        std::cout << "no such detector!!!!!!!!!!!!!" << std::endl;
      }  
        

        if(diblock<2 && dcm <12) fDCM_FEB[diblock-1][dcm-1][feb]++;
      
    } //End of loop over RawDigits
    
    //Timestamps of non-empty events; Counters for empty events and noise only events
    
    if( (fFirstNonEmptyEventTimeStamp == 0) && (rawhitlist.size() !=0) ) {
      fFirstNonEmptyEventTimeStamp=uevt_sec;
        }
    
    if ( (fNevents > 0) && (rawhitlist.size() !=0) ) {
      fLastNonEmptyEventTimeStamp=uevt_sec;
    }
    
    if(hitlist.size()==0) {
      ++fNempty;
      fEmptyEvents->Fill(1);
    }
    else if (slicecol->size()==1) {
      ++fNoise;
      fNoiseEvents->Fill(1);
    }
    
    fNHits.push_back(hitlist.size());
    
  } //End of DataCheck::analyze
  
  //......................................................................
  
  void DataCheck::endJob() 
  {
    //Rescaling fPlaneVCell by dividing by number of events
    art::ServiceHandle<geo::Geometry> geom;

    fPlaneVCell->Scale(1.0/(float)fNevents);
    fPlaneVCellRaw->Scale(1.0/(float)fNevents);
    
    //Occupancy Plots
    
    for (int xbin = 1; xbin <= fPlaneVCell->GetNbinsX(); xbin++) {
      for (int ybin = 1; ybin <= fPlaneVCell->GetNbinsY(); ybin++) {
        
        double cellcontent = fPlaneVCell->GetBinContent(xbin,ybin);
        double logcontent = -6;
        
        if (cellcontent > 0.0) {
          logcontent = log10(cellcontent);
          fNGoodChannels++;
          //fNGoodChannels->Fill(1);
          fChannelMap->SetBinContent(xbin,ybin,1);
        }// End of assigning log to cell content
        
        fOccupancy->Fill(logcontent);
        
      }// End of loop over y cells
    }// End of loop over x cells
    
    // Count number of active raw channels
    
    for (int xbin = 1; xbin <= fPlaneVCellRaw->GetNbinsX(); xbin++) {
      for (int ybin = 1; ybin <= fPlaneVCellRaw->GetNbinsY(); ybin++) {
        
        double cellcontent = fPlaneVCellRaw->GetBinContent(xbin,ybin);
        
        if (cellcontent > 0.0) {
          fNRawChannels++;
        }
        
      }// End of loop over y cells
    }// End of loop over x cells
    
    // occupancy fit
    TF1 *f1=new TF1("f1","landau");
    fOccupancy->Fit("f1");
    
    
    //rescale fNHits by the number of good channels
    for(unsigned int i = 0; i < fNHits.size(); i++){
      fNHitsPerEvent->Fill(fNHits[i]/(float)fNGoodChannels);
    }
    
    
   

   std::cout<< fBlocDcmFeb->GetXaxis()->GetNbins()<<std::endl;
   for(int nx=1;nx<=14;nx++){
       for(int ny=1;ny<=12;ny++){
          for(int nz=1;nz<=64;nz++){

      if( fBlocDcmFeb->GetBinContent(nx,ny,nz)) fNDCM++;

     }}}




    //Get software release
    fSoftwareRelease = getenv("SRT_BASE_RELEASE");
    
    //Figure out "Good Run Status" Standard cuts are below 100% empty events, less than 10% noise only events, hits per event between 0.00 and 0.06
    unsigned int stat = 0;
    float femptyp = fNempty/(1.0*fNevents);
    float fnoisep = 0.0;
    if (fNempty == fNevents) fnoisep = fNoise/(1.0*(fNevents-fNempty));
    float fhits   = fNHitsPerEvent->GetMean(1);
    
    if(fRun>=11570 && fRun<12978){
        if(femptyp < 1.0 && fnoisep < 0.1 && fhits < 0.06 && fhits > 0.0) stat=1;
        else stat=0;
    }
    if(fRun>12978){
        if(femptyp < 1.0 && fnoisep < 0.25 && fhits < 0.06 && fhits > 0.0) stat=1;
        else stat=0;
    }

    std::cout<<"detector="                   <<"NDOS"<<std::endl;
    std::cout<<"run="                        <<fRun<<std::endl;
    std::cout<<"subrun="                     <<fSubRun<<std::endl;
    std::cout<<"triggertype="                <<fTriggerType<<std::endl;            
    std::cout<<"nevents="                    <<fNevents<<std::endl;
    std::cout<<"nemptyevents="               <<fNempty<<std::endl;
    std::cout<<"nnoiseevents="               <<fNoise<<std::endl;
    std::cout<<"nactivechannels="            <<fNGoodChannels<<std::endl;
    std::cout<<"nrawactivechannels="         <<fNRawChannels<<std::endl;
    std::cout<<"nslices="                    <<fNSlices->GetMean(1)<<std::endl;
    std::cout<<"sliceduration="              <<fSliceLength->GetMean(1)<<std::endl;
    std::cout<<"hitsperslice="               <<fNHitsPerSlice->GetMean(1)<<std::endl;
    std::cout<<"hitsperslicenonoise="        <<fNHitsPerSliceNoNoise->GetMean(1)<<std::endl;
    std::cout<<"hitsperevent="               <<fNHitsPerEvent->GetMean(1)<<std::endl;
    std::cout<<"occupancympv="               <<f1->GetParameter(1)<<std::endl;
    std::cout<<"occupancysigma="             <<f1->GetParameter(2)<<std::endl;
    std::cout<<"software release="           <<fSoftwareRelease<<std::endl;
    std::cout<<"Good Run Status= "           <<stat<<std::endl;
    std::cout<<"Number of active DCMs= "     <<fNDCM<<std::endl;
    for(int ndiblock=0;ndiblock<4;ndiblock++){
      for(int ndcm=0;ndcm<3;ndcm++){
        std::cout<<"nfebfor DCM"<<ndiblock+1<<"-"<<ndcm+1<<"="<<fDCM_nFEB[ndiblock][ndcm]<<std::endl;
      }
    }
    
    char fFirstEventTimeStampStr[256];
    char fFirstNonEmptyEventTimeStampStr[256];
    char fLastNonEmptyEventTimeStampStr[256];
    char fLastEventTimeStampStr[256];
    
    struct tm * fFirstEventTimeStampinfo;
    struct tm * fFirstNonEmptyEventTimeStampinfo;
    struct tm * fLastNonEmptyEventTimeStampinfo;
    struct tm * fLastEventTimeStampinfo;
    
    fFirstEventTimeStampinfo = localtime(&fFirstEventTimeStamp);
    strftime(fFirstEventTimeStampStr,256,"%Y-%m-%d %H:%M:%S", fFirstEventTimeStampinfo);
    
    std::cout<<"Time of First Event=           "    <<fFirstEventTimeStampStr<<std::endl;
    
    fFirstNonEmptyEventTimeStampinfo = localtime(&fFirstNonEmptyEventTimeStamp);
    strftime(fFirstNonEmptyEventTimeStampStr,256,"%Y-%m-%d %H:%M:%S", fFirstNonEmptyEventTimeStampinfo);
    
    std::cout<<"Time of First Non Empty Event= "    <<fFirstNonEmptyEventTimeStampStr<<std::endl; 
    
    fLastNonEmptyEventTimeStampinfo = localtime(&fLastNonEmptyEventTimeStamp);
    strftime(fLastNonEmptyEventTimeStampStr,256,"%Y-%m-%d %H:%M:%S", fLastNonEmptyEventTimeStampinfo);
    
    std::cout<<"Time of Last Non Empty Event=  "    <<fLastNonEmptyEventTimeStampStr<<std::endl; 
    
    fLastEventTimeStampinfo  = localtime(&fLastEventTimeStamp);
    strftime(fLastEventTimeStampStr,256,"%Y-%m-%d %H:%M:%S", fLastEventTimeStampinfo);
    
    std::cout<<"Time of Last Event=            "    <<fLastEventTimeStampStr<<std::endl; 
    
    
    
    if (fWriteToDB && !fIsNearline){   
      
      nova::database::Table* t = new nova::database::Table(std::string("/Commissioning/tables/DataCheck.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 is 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("nemptyevents",                   fNempty);
        r->Set("nnoiseevents",                   fNoise); 
        r->Set("nactivechannels",                fNGoodChannels);
        r->Set("nrawactivechannels",             fNRawChannels);
        r->Set("firsteventtime",                 fFirstEventTimeStampStr); 
        r->Set("firstnonemptyeventtime",         fFirstNonEmptyEventTimeStampStr); 
        r->Set("lastnonemptyeventtime",          fLastNonEmptyEventTimeStampStr); 
        r->Set("lasteventtime",                  fLastEventTimeStampStr);
        r->Set("nslices",                        fNSlices->GetMean(1));
        r->Set("sliceduration",                  fSliceLength->GetMean(1));
        r->Set("hitsperslice",                   fNHitsPerSlice->GetMean(1)); 
        r->Set("hitsperslicenonoise",            fNHitsPerSliceNoNoise->GetMean(1));
        r->Set("hitsperevent",                   fNHitsPerEvent->GetMean(1));        
        r->Set("occupancympv"  ,                 f1->GetParameter(1));
        r->Set("occupancysigma" ,                f1->GetParameter(2));
        r->Set("softwarerelease",                fSoftwareRelease);
        r->Set("status",                         stat);
        r->Set("nDCMs",                          fNDCM);
        r->Set("DCM11",                      fDCM_nFEB[1-1][1-1]); 
        r->Set("DCM12",                      fDCM_nFEB[1-1][2-1]);
        r->Set("DCM13",                      fDCM_nFEB[1-1][3-1]);
        r->Set("DCM21",                      fDCM_nFEB[2-1][1-1]);
        r->Set("DCM22",                      fDCM_nFEB[2-1][2-1]);
        r->Set("DCM23",                      fDCM_nFEB[2-1][3-1]);
        r->Set("DCM31",                      fDCM_nFEB[3-1][1-1]);
        r->Set("DCM32",                      fDCM_nFEB[3-1][2-1]);
        r->Set("DCM33",                      fDCM_nFEB[3-1][3-1]);
        r->Set("DCM41",                      fDCM_nFEB[4-1][1-1]);
        r->Set("DCM42",                      fDCM_nFEB[4-1][2-1]);

        t->AddRow(r);
        delete r;
        t->WriteToDB();
      } 
      
    }
    
    if (fWriteToDB && fIsNearline){   
      
      nova::database::Table* t = new nova::database::Table(std::string("/Commissioning/tables/NearlineTable.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 is 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("nemptyevents",                   fNempty);
        r->Set("nnoiseevents",                   fNoise); 
        r->Set("nactivechannels",                fNGoodChannels);
        r->Set("nrawactivechannels",             fNRawChannels);
        r->Set("firsteventtime",                 fFirstEventTimeStampStr); 
        r->Set("firstnonemptyeventtime",         fFirstNonEmptyEventTimeStampStr); 
        r->Set("lastnonemptyeventtime",          fLastNonEmptyEventTimeStampStr); 
        r->Set("lasteventtime",                  fLastEventTimeStampStr);
        r->Set("nslices",                       fNSlices->GetMean(1));
        r->Set("sliceduration",              fSliceLength->GetMean(1));
        r->Set("hitsperslice",               fNHitsPerSlice->GetMean(1)); 
        r->Set("hitsperslicenonoise",        fNHitsPerSliceNoNoise->GetMean(1));
        r->Set("hitsperevent",               fNHitsPerEvent->GetMean(1));        
        r->Set("occupancympv"  ,             f1->GetParameter(1));
        r->Set("occupancysigma" ,            f1->GetParameter(2));
        r->Set("softwarerelease",            fSoftwareRelease);
        r->Set("status",                         stat);
        r->Set("nDCMs",                      fNDCM);
        r->Set("DCM11",                      fDCM_nFEB[1-1][1-1]); 
        r->Set("DCM12",                      fDCM_nFEB[1-1][2-1]);
        r->Set("DCM13",                      fDCM_nFEB[1-1][3-1]);
        r->Set("DCM21",                      fDCM_nFEB[2-1][1-1]);
        r->Set("DCM22",                      fDCM_nFEB[2-1][2-1]);
        r->Set("DCM23",                      fDCM_nFEB[2-1][3-1]);
        r->Set("DCM31",                      fDCM_nFEB[3-1][1-1]);
        r->Set("DCM32",                      fDCM_nFEB[3-1][2-1]);
        r->Set("DCM33",                      fDCM_nFEB[3-1][3-1]);
        r->Set("DCM41",                      fDCM_nFEB[4-1][1-1]);
        r->Set("DCM42",                      fDCM_nFEB[4-1][2-1]);
        
        t->AddRow(r);
        delete r;
        t->WriteToDB();
      }

    }   
    
  //.........................................................
    novadaq::cnv::DetId detectorid = geom->DetId();
    nova::dbi::RunHistory rh(fRun, detectorid); 
    rh.LoadDAQRunHistory();

    std::ostringstream oss;     
    oss << fReportDir << "feb_shutoff_r" << fRun << "_s" << fSubRun << ".txt";  
//    std::cout << "run: " << fRun << " subrun: " << fSubRun << std::endl;
    std::string fname(oss.str());
    std::ofstream fout(fname.c_str(),std::ios::app);    

    int nbinx = vDCMs.at(1)->GetNbinsX(); // number of events
    int nbiny = vDCMs.at(1)->GetNbinsY(); // number of febs
    double width = vDCMs.at(1)->GetBinWidth(1); //number of events per bin
    nbinx = (int)(fNevents/width);    
    // 2013-10-19 JMP: not sure if this is what is desired... orig. author(s)
    // should check this!
    int nDiBlocks = rh.NDiBlocks();
    int ic=0;
    for ( int i=0; i<nDiBlocks; ++i) {
      const nova::dbi::RunHistory::DiBlock db = rh.GetDiBlock(i);
      for (unsigned int idcm=0; idcm<db.dcm.size(); ++idcm, ++ic) {
        std::string dcmname = db.dcm[idcm].appname;
        for( int j=0;j<nbiny;j++ ){ 
          bool biscool = db.dcm[idcm].feb[j].isCooled; //feb name is from 0-63, the list starts from 0
          //        double APDv = tempDCM.feb[j].apdV;
          std::string strcool  = "hello";
          if( biscool ) strcool = "is cool";
          else strcool = "is not cool";
          int prev = 999;
          int curr = 999;
          for( int k=0;k<nbinx;k++ ){
            if( ( (vDCMs.at(ic))->GetBinContent(k+1,j+1))!=0.0 ) curr = 1;
            else curr = 0;
            if( (k!=0) & (prev>curr) ){
              if(fout){
                fout << "FEB" << j << " in " << dcmname << " shutoff in subrun" << fSubRun << ". This feb is " << strcool << std::endl;
              }
            }
            else if( (k!=0) & (prev<curr) ){
              if(fout){
                fout << "FEB" << j << " in " << dcmname << " came back in subrun" << fSubRun << ". This feb is " << strcool << std::endl;
              }
            }
            prev = curr;
          }
        } //loop over events 
      } //loop over febs
    } //loop over DCMs  

    fout.close();
    
  } //End of DataCheck::endJob
  
  //......................................................................
  
  void DataCheck::reconfigure(const fhicl::ParameterSet& p)
  {
    fSlicerModule   = p.get<std::string>("SlicerModule");
    fCalHitModule   = p.get<std::string>("CalHitModule");
    fRawDigitModule = p.get<std::string>("RawDigitModule", "daq");
    fRawDataLabel   = p.get<std::string>("RawDataLabel", "daq");
    fWriteToDB      = p.get<int>("WriteToDB");
    fIsNearline     = p.get<int>("Nearline");
    fReportDir      = p.get<std::string>("ReportDir");
  }

} // end namespace comi

DEFINE_ART_MODULE(comi::DataCheck)