EvtAddRawDigitBlockModule.cxx

Go to the documentation of this file.

// $Id: EvtAddRawDigitBlockModule.cxx,v 1.1 2014/06/24 20:25:53 rhatcher Exp $
//
// A JobControl Module for adding a RawDigitDataBlock from a separate
// stream (served by EvtOverlaySingleton) into the current record
//
// rhatcher@fnal.gov

#include <cassert>
#include <csignal>
#include <cmath>
#include <sstream>
#include <algorithm>

#include <Conventions/Munits.h>

#include <EvtDataOverlay/EvtAddRawDigitBlockModule.h>
#include <EvtDataOverlay/EvtOverlaySingleton.h>
#include <Util/WhiteBoardRegistry.h>

#include <MinosObjectMap/MomNavigator.h>
#include <MessageService/MsgService.h>
#include <JobControl/JobCModuleRegistry.h>

#include <Record/SimSnarlHeader.h>
#include <Record/SimSnarlRecord.h>
#include <RawData/RawRecord.h>
#include <RawData/RawDaqSnarlHeader.h>
#include <RawData/RawDigitDataBlock.h>
#include <RawData/RawDigit.h>
#include <RawData/RawVtmTimeInfoBlock.h>

#include <HistMan/HistMan.h>
#include <Calibrator/Calibrator.h>
#include <Plex/PlexHandle.h>

#include <BeamDataUtil/BeamMonSpill.h>

#include "TMath.h"
#include "TString.h"

ClassImp(EvtAddRawDigitBlockModule)

//......................................................................

CVSID("$Id: EvtAddRawDigitBlockModule.cxx,v 1.1 2014/06/24 20:25:53 rhatcher Exp $");
JOBMODULE(EvtAddRawDigitBlockModule, "EvtAddRawDigitBlock",
         "Adds RawDigitDataBlock to DaqSnarl");

//......................................................................

// QIE electronics, tick = 18.832 ns (RF clock)
const double tdc_convert =  1000./ 53.1;  // clock is 53.1MHz


//......................................................................

EvtAddRawDigitBlockModule::EvtAddRawDigitBlockModule() 
  : fMaxDeltaT(1.3)
  , fSmallT0(20)   // less than 20n crateT0ns is "small"
{
   // construct a new "EvtAddRawDigitBlockModule" JobControl module

}

//......................................................................

EvtAddRawDigitBlockModule::~EvtAddRawDigitBlockModule() 
{
  MSG("Hepevt", Msg::kVerbose) << "EvtAddRawDigitBlockModule::Destructor\n";

}

//......................................................................

JobCResult EvtAddRawDigitBlockModule::Reco(MomNavigator *mom)
{

  // get what we can
  RawRecord* rawrec = 
    dynamic_cast<RawRecord*>(mom->GetFragment("RawRecord",0,"DaqSnarl"));
  if ( ! rawrec ) return JobCResult::kFailed;

  const RawDaqSnarlHeader* snarlHdr =
    dynamic_cast<const RawDaqSnarlHeader*>(rawrec->GetRawHeader());
  if ( ! snarlHdr ) return JobCResult::kFailed;

  VldTimeStamp vts = snarlHdr->GetVldContext().GetTimeStamp();

  EvtOverlaySingleton& eos = EvtOverlaySingleton::Instance();
  eos.FindClosestDaqSnarl(vts);  // find it
  if ( ! eos.IsValidSnarl() ) {
    MSG("EvtDataOverlay",Msg::kError) << "invalid DaqSnarl match for "
                                      << vts.AsString("c") << std::endl;
    return JobCResult::kFailed;
  }
  double deltat = eos.GetSnarlTimeDifference();
  if ( std::abs(deltat) > fMaxDeltaT ) {
    MSG("EvtDataOverlay",Msg::kWarning) << "found closest DaqSnarl record was " 
                                        << deltat << " away from " 
                                        << vts.AsString("c") << std::endl;
  }

  // get an unowned pointer to the raw block in the .mdaq.root file
  const RawDigitDataBlock* origRDDB = eos.GetClosestRawDigitDataBlock();

  // now get our very own copy ... to love and keep forever
  RawDigitDataBlock* copyRDDB = new RawDigitDataBlock(*origRDDB);

  // oh, wait, I can't keep this poor little thing
  // I have neither the means nor wherewithal to care for it
  // guess I'll look around for a good home
  rawrec->AdoptRawBlock(copyRDDB);

  // Minerva also expects NtpSRTimeStatus to be filled
  // ... that is filled from the DaqMonitor stream RawVtmTimeInfoBlock
  // we're not faking up a DaqMonitor stream, so just find the right
  // one of those blocks and copy it into the DaqSnarl stream
  eos.FindLowerBoundDaqMonitor(vts);  // find it
  if ( ! eos.IsValidMonitor() ) {
    MSG("EvtDataOverlay",Msg::kError) << "invalid DaqMonitor match for "
                                      << vts.AsString("c") << std::endl;
    return JobCResult::kFailed;
  }
  double deltatm = eos.GetMonitorTimeDifference();
  // record time should preceed requested time (req-found>0)
  if ( deltatm > fMaxDeltaT || deltatm < 0 ) {
    MSG("EvtDataOverlay",Msg::kWarning) << "found lowerbound DaqMonitor record was " 
                                        << deltatm << " away from " 
                                        << vts.AsString("c") << std::endl;
  }

  // get an unowned pointer to the raw block in the .mdaq.root file
  const RawVtmTimeInfoBlock* origRVTIB = 
    eos.GetLowerBoundRawVtmTimeInfoBlock();
  RawVtmTimeInfoBlock* copyRVTIB = new RawVtmTimeInfoBlock(*origRVTIB);
  rawrec->AdoptRawBlock(copyRVTIB);

  return JobCResult::kAOK; // All Ok
}


//......................................................................
JobCResult EvtAddRawDigitBlockModule::Ana(const MomNavigator *mom)
{
   // Check out current RecSimSnarl Record

  RawRecord* rawrec = 
    dynamic_cast<RawRecord*>(mom->GetFragment("RawRecord",0,"DaqSnarl"));
  if ( ! rawrec ) {
    MSG("Overlay",Msg::kWarning)
      << "No RawRecord found" << std::endl;
    return JobCResult::kFailed;
  }

  const RawDaqSnarlHeader* snarlHdr =
    dynamic_cast<const RawDaqSnarlHeader*>(rawrec->GetRawHeader());
  if ( ! snarlHdr ) {
    MSG("Overlay",Msg::kWarning)
      << "No RawDaqSnarlHeader found" << std::endl;
    return JobCResult::kFailed;
  }

  // book histograms after we know what run/subrun
  // can't necessarily use BeginRun() or BeginSubrun() because
  // there might not be marker records if MC
  static bool first = true;
  if ( first ) {
    int run = snarlHdr->GetRun();
    int subrun = snarlHdr->GetSubRun();
    BookHist(run,subrun);
    first = false;
  }

  VldContext vldc = snarlHdr->GetVldContext();
  Calibrator& cal = Calibrator::Instance();
  cal.Reset(vldc);
  PlexHandle ph(vldc);

  EvtOverlaySingleton& eos = EvtOverlaySingleton::Instance();

  // fill some HistMan histograms
  HistMan hm("overlay");

  std::vector<int> tdc0v, tdc1v;
  VldTimeStamp t0_0(VldTimeStamp::GetBOT()), t0_1(VldTimeStamp::GetBOT());

  int nrddb = -1;  // which RawDigitDataBlock do we have 
  TIter rbitr = rawrec->GetRawBlockIter();
  RawDataBlock* rdb = 0;
  // need to do cast from TObject
  // also rbitr() advances iterator as side effect
  while ( (rdb = dynamic_cast<RawDataBlock*>(rbitr())) ) {
    RawDigitDataBlock* rddb = dynamic_cast<RawDigitDataBlock*>(rdb);
    if ( ! rddb ) continue; // no what we want?

    ++nrddb;  // got one

    // figure out what histograms we might fill
    const char* name_small = ((nrddb==0) ? "smallT0 block0":"smallT0 block1");
    const char* name_big   = ((nrddb==0) ? "bigT0 block0"  :"bigT0 block1");
    TH1D* hsmall = hm.Get<TH1D>(name_small);
    TH1D* hbig   = hm.Get<TH1D>(name_big);

    std::vector<int>& tdcv = ((nrddb==0)? tdc0v : tdc1v);

    //RWH test
    //if ( nrddb == 0 ) rddb->FormatToOStream(std::cout,"");


    std::set<VldTimeStamp> t0s = eos.GetCrateT0s(rddb,vldc);
    if ( ! t0s.empty() ) {
      if (nrddb==0) t0_0 = *(t0s.begin());
      else          t0_1 = *(t0s.begin());
      if ( t0s.size() > 1 ) {
        MSG("Overlay",Msg::kWarning) 
          << "!!!!!!! more than one crateT0 nrddb=" << nrddb << endl;
      }
    } else {
        MSG("Overlay",Msg::kWarning) 
          << "!!!!!!! no crateT0 nrddb=" << nrddb << endl;      
    }

    // loop over the digits
    TIter rditer = rddb->GetDatumIter();
    RawDigit* digit = 0;
    int lastCrate = -1;
    VldTimeStamp crateT0;
    bool isSmallT0 = true;
    while ( ( digit = (RawDigit*)rditer() ) ) {
      RawChannelId rcid = digit->GetChannel();
      ReadoutType::Readout_t digitType = ph.GetReadoutType(rcid);
      if ( ReadoutType::kScintStrip != digitType ) continue;

      // new crate ? check T0 
      int crate = rcid.GetCrate();
      if ( crate != lastCrate ) {
        lastCrate = crate;
        crateT0 = digit->GetCrateT0();
        isSmallT0 = ( abs(crateT0.GetNanoSec()) < fSmallT0 );
      }

      tdcv.push_back(digit->GetTDC());

      // calculate time like AlgDigit would ... (hopefully)
      // calibrator returns ticks to seconds (not ns like it says)
      double digitTimeNS = cal.GetTimeFromTDC(digit->GetTDC(),rcid)*1.0e9;
      //  double calns = digitTimeNS;  // testing
      // adjust for skew
      digitTimeNS += ( digit->GetCrateT0().GetSec() -
                       vldc.GetTimeStamp().GetSec()    ) * 1.0e9;

      double tadj = digitTimeNS - 
        ( vldc.GetTimeStamp().GetNanoSec() - crateT0.GetNanoSec() );
      tadj *= 1.0e-9;  // convert to seconds

      MSG("Overlay",Msg::kDebug)
        << "nrddb " << nrddb << "|" << ((isSmallT0)?"s":"B")
        << " tadj = " << setw(8) << tadj/Munits::microsecond << " usec " 
        << " TDC " << digit->GetTDC()
        //  << " ditiTimeNS " << calns << " / " << digitTimeNS 
        //  << " vldc.ns " << vldc.GetTimeStamp().GetNanoSec()
        //  << " crateT0.ns " << crateT0.GetNanoSec()
        << "  " << rcid.AsString("c")
        << endl;

      if ( isSmallT0 ) hsmall->Fill(tadj);
      else             hbig->Fill(tadj);

    } // loop over digits

    if ( lastCrate != -1 ) {
      // did we see anything?  if not then crateT0 & isSmallT0 are meaningless
      TH1I* hseen = hm.Get<TH1I>("rddb seen");
      if ( isSmallT0 ) hseen->Fill(name_small,1.0);
      else             hseen->Fill(name_big,1.0);

      /*
      int ibtype = ((isSmallT0)?0:2) + ((nrddb==0)?0:1);
      //hseen->Fill(ibtype);  // how many blocks of each type have we seen
      if (   isSmallT0 && nrddb==0 ) hseen->Fill("smallT0 block0",1);
      if (   isSmallT0 && nrddb==1 ) hseen->Fill("smallT0 block1",1);
      if ( ! isSmallT0 && nrddb==0 ) hseen->Fill("bigT0 block0",1);
      if ( ! isSmallT0 && nrddb==1 ) hseen->Fill("bigT0 block1",1);

      MSG("Overlay",Msg::kWarning)
        << "#RWH ibtype " << ibtype << " nrddb " << nrddb 
        << " isSmallT0 " << isSmallT0 << " crateT0 " << crateT0.AsString("c")
        << endl;
      */

    }

  } // loop over blocks

  // make some snarl specific histograms

  if ( t0_0 != t0_1 ) {
    MSG("Overlay",Msg::kWarning)
      << "!!!!!!! RDDB blocks had different crateT0 "
      << t0_0.AsString("c") << " " << t0_1.AsString("c") << endl;
  }

  std::sort(tdc0v.begin(),tdc0v.end());
  std::sort(tdc1v.begin(),tdc1v.end());

  int maxtdc = *(tdc1v.rbegin());
  int mintdc = *(tdc1v.begin());

  int snarl = snarlHdr->GetSnarl();
  std::ostringstream s;
  s << "snarl" << snarl << "_";
  std::string name0 = s.str() + "0";
  std::string name1 = s.str() + "1";

  MSG("Overlay",Msg::kSynopsis)
    << "Snarl " << snarl << " TDC min/max " << mintdc << " " << maxtdc 
    << " range " << (maxtdc-mintdc) << endl;

#undef HIST_INDIVIDUAL
#ifdef HIST_INDIVIDUAL
  HistMan hmMany("individual");
  TH1I* htdc0 = hmMany.Book<TH1I>(name0.c_str(),name0.c_str(),(maxtdc-mintdc)+1,mintdc,maxtdc+1);
  TH1I* htdc1 = hmMany.Book<TH1I>(name1.c_str(),name1.c_str(),(maxtdc-mintdc)+1,mintdc,maxtdc+1);

  for (size_t i=0; i<tdc0v.size(); ++i) htdc0->Fill(tdc0v[i]);
  for (size_t i=0; i<tdc1v.size(); ++i) htdc1->Fill(tdc1v[i]);
#endif

  // some info about POTs and beam spill structure 
  int index = 0;
  minos::GetFromWhiteBoard("EvtInfo","TimeIndex",index);
  double pots             = eos.GetClosestSpillPOTs(true);
  const BeamMonSpill* bms = eos.GetClosestBeamMonSpill();

  HistMan hmPOTs(fHM_pots_dir.c_str());
  hmPOTs.Get<TH1D>("POTs")->Fill(index,pots);
  TH2D* hbi =  hmPOTs.Get<TH2D>("BatchIntensity");
  if ( tdc0v.size() > 0 ) {
    double x = 100. + log((double)(tdc0v.size()));
    hbi->Fill(index,7,x);  // scale up so colors work
  }
  if ( bms ) {
    for (int ib = 0; ib < 6; ++ib) {
      double batchi = bms->fBpmInt[ib];
     hbi->Fill(index,ib,batchi);
    }
  }

  return JobCResult::kPassed; // All Ok
}
//......................................................................

const Registry& EvtAddRawDigitBlockModule::DefaultConfig() const
{
//======================================================================
// Create a registry which holds the default configuration and return it
//======================================================================
  static Registry r;
 
  // Set name of config
  std::string name = this->GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  // Set values of config
  r.UnLockValues();

  r.Set("MaxDeltaT",1.3);  // warn if found snarl data is more than this away

  r.LockValues();

  return r;
}

//......................................................................

void EvtAddRawDigitBlockModule::Config(const Registry& r)
{
//======================================================================
// Configure the module given the registry r
//======================================================================
  //const char*  tmps;
  double tmpd;
  //int    tmpi;

  if (r.Get("MaxDeltaT", tmpd)) { fMaxDeltaT = tmpd; }

}

//......................................................................

void EvtAddRawDigitBlockModule::BeginJob()
{

}

void EvtAddRawDigitBlockModule::BookHist(int run, int subrun)
{
  // protect against being called more than once
  static bool first = true;
  if ( ! first ) return;
  first = false;

  // book some HistMan histograms
  HistMan hm("overlay");

  double upper = 20000. * 1.0e-9;  // look over 20usec, in seconds
  int ngang = 1;          // 1 or 4 tdc ticks per bin
  int nbins = TMath::Nint((double)upper/ (tdc_convert*1.0e-9)) / ngang;

  TH1I* hseen = hm.Book<TH1I>("rddb seen","rddb seen",4,0,4);
  hseen->SetBit(TH1::kCanRebin);

  const char* hnames[4] = { "smallT0 block0", "smallT0 block1",
                            "bigT0 block0",   "bigT0 block1" };

  for (int ihn=0; ihn<4; ++ihn) {
    const char* hname = hnames[ihn];
    std::string htitle(hname);
    htitle += "; (sec); # digits";
    hm.Book<TH1D>(hname,htitle.c_str(),nbins,0,upper);
    // can only alphabetize labels after some filling
    hseen->Fill(hname,0.0);
  }
  hseen->LabelsOption("a");

  MSG("Overlay",Msg::kInfo) 
    <<"EvtAddRawDigitBlockModule BOOK block[0|1] [small|big]T0 "
    << "(" << nbins << ",0," << upper <<")" << endl;

  // set the directory for the POTs info
  // separate out by run/subrun so "hadd" doesn't smack them together
  fHM_pots_dir = Form("r%ds%d POTs",run,subrun);
  HistMan hmPOTs(fHM_pots_dir.c_str());

  // keep track of POTs and 5-6 batch intensities for each snarl
  // use "7" for (scaled) # digits in RawDigitDataBlock 0 (Minerva)
  EvtOverlaySingleton& eos = EvtOverlaySingleton::Instance();
  size_t ntimes = eos.NumUserTimes();

  hmPOTs.Book<TH1D>("POTs","POTs",ntimes,0,ntimes);
  hmPOTs.Book<TH2D>("BatchIntensity","BatchIntensity",ntimes,0,ntimes,8,0,8);

}

//......................................................................

void EvtAddRawDigitBlockModule::EndJob()
{

  MSG("Hepevt",Msg::kInfo)
    << "EvtAddRawDigitBlockModule Added " << endl;

}

//......................................................................

---