TimingCalibration_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       TimingCalibration
// Module Type: producer
// File:        TimingCalibration_module.cc
//
// Generated at Mon Feb 11 09:01:32 2013 by Evan Niner using artmod
// from art v1_02_06.
////////////////////////////////////////////////////////////////////////

//NOvASoft includes
#include "Geometry/Geometry.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "Utilities/func/MathUtil.h"
#include "CalibrationDataProducts/PCHit.h"
#include "CalibrationDataProducts/TCTrack.h"
#include "CalibrationDataProducts/DCMStat.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "DAQChannelMap/DAQChannelMap.h"
#include "Calibrator/Calibrator.h"
#include "MCCheater/BackTracker.h"

//framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "Utilities/AssociationUtil.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"


//ROOT includes
#include "TTree.h"
#include "TH2F.h"
#include "TH1F.h"

namespace calib {

  class TimingCalibration : public art::EDProducer {
  public:
    explicit TimingCalibration(fhicl::ParameterSet const & p);
    virtual ~TimingCalibration();

    void produce(art::Event & e) override;

    void reconfigure(const fhicl::ParameterSet& pset);
    void beginJob() override;
    void clearNTuple();
    void beginRun(art::Run& run) override;


  private:

    std::string   fTCTrackLabel;        ///< Where to find calibration tracks
    std::string   fTrackLabel;        ///< Where to find reco Tracks
            
    //variables to select track quality
    int    fMinDCMHits;      ///<minimum number of hits on a dcm to keep results
    int    fMinDCM;         ///<minimum number of good dcms on a track
    bool   fDebug;          ///<Add exta debugging info to ntuple

    float fFiberSpeed;      ///<fixed fiberspeed
    bool fRemoveOutlier;    ///<remove outlier hits in DCMs based on timing resolution
    int  fRemoveThreshold;  ///<number of sigma away from mean to remove hits
    bool fCorrectTimewalk;  ///<adjust hit times based on timewalk effect for distance to readout
    bool fAdjustSpeed; ///<adjust fiberspeed based on distance to readout instead of assuming flat speed


    //map relating unique dcm to collection of pchits
    std::map<int, std::vector<art::Ptr<caldp::PCHit> > > fDCMHitMap;

    //ntuple variables
    //event statistics
    int run;
    int subrun;
    int event;
    int cosmictr;
    //track start coordinate
    double sX;
    double sY;
    double sZ;
    //track end coordinate
    double eX;
    double eY;
    double eZ;

    double plength;

    //min, max plane
    int minPl;
    int maxPl;
    //min, max cell
    int minPlX;
    int minPlY;
    int maxPlX;
    int maxPlY;
    //number of cells used in final fiber speed fit
    int ncellsfitT;
    int ncellsfitX;
    int ncellsfitY;

    uint32_t evtTime;

    //vectors filled for each track, one entry per dcm
    std::vector<float> dcmnum; ///<dcm index
    std::vector<float> dcmhits; ///<number of hits in dcm
    std::vector<float> dcmsimtime; ///<simultaneous hit time in dcm

    //vectors filled for each track, one entry per hit, sorted by dcm
    std::vector<float> cellnum; ///<cell number
    std::vector<float> planenum; ///<plane number
    std::vector<float> dcmsimresid; ///<difference between simultaneous hit time and average
    std::vector<float> dcmtime; ///<uncorrected times of hits
    std::vector<float> dcmstime; ///<times with distance to readout, time-of-flight, timewalk corrections
    std::vector<float> dcmrd;  ///<distance to readout
    std::vector<float> dcmpe;  ///<number of photoelectrons
    std::vector<float> dcmpl; ///<path length along track
    

    //tree to store ntuple statistics
    TTree* fOutTree;

    fhicl::ParameterSet fPSetND, fPSetFD, fPSetTB;



   
    

  };
}

namespace calib {

  //define speed of light in cm/ns
  static const double sCmPerNsec = 29.9792458;

  TimingCalibration::TimingCalibration(fhicl::ParameterSet const & pset)
  // :
  // Initialize member data here.
  {
    this->reconfigure(pset); 

    produces< std::vector<caldp::DCMStat> >();
    produces< art::Assns<caldp::DCMStat, rb::Track> >();
  }

//---------------------------------------------------------------------------
  void TimingCalibration::reconfigure(const fhicl::ParameterSet& pset)
  {
    fTCTrackLabel     = pset.get< std::string         >("TCTrackLabel");
    fTrackLabel       = pset.get< std::string         >("TrackLabel");
    fMinDCMHits       = 0;
    fMinDCM           = 0;
    fDebug            = pset.get< bool                >("Debug");
    fFiberSpeed       = pset.get< float               >("FiberSpeed");
    fRemoveOutlier    = pset.get< bool                >("RemoveOutlier");
    fRemoveThreshold  = pset.get< int                 >("RemoveThreshold");
    fCorrectTimewalk  = pset.get< bool                >("CorrectTimewalk");
    fAdjustSpeed = pset.get< bool                >("AdjustSpeed");
    fPSetND           = pset.get< fhicl::ParameterSet >("nd");
    fPSetFD           = pset.get< fhicl::ParameterSet >("fd");
    fPSetTB           = pset.get< fhicl::ParameterSet >("tb");    

  }  

  //---------------------------------------------------------------------------
  TimingCalibration::~TimingCalibration()
  {
  }

  //---------------------------------------------------------------------------
  void TimingCalibration::beginJob()
  {

    art::ServiceHandle<art::TFileService> tfs;

    // define the tree
    fOutTree = tfs->make<TTree>("TimingCalib","Timing Calibaration");

    //define ntuple branches
    fOutTree->Branch("run",    &run,    "run/I");
    fOutTree->Branch("subrun", &subrun, "subrun/I");
    fOutTree->Branch("event",  &event,  "event/I");
    fOutTree->Branch("cosmictr",  &cosmictr,  "cosmictr/I");
    fOutTree->Branch("sX", &sX, "sX/D");
    fOutTree->Branch("sY", &sY, "sY/D");
    fOutTree->Branch("sZ", &sZ, "sZ/D");
    fOutTree->Branch("eX", &eX, "eX/D");
    fOutTree->Branch("eY", &eY, "eY/D");
    fOutTree->Branch("eZ", &eZ, "eZ/D");
    fOutTree->Branch("plength", &plength, "plength/D");
    fOutTree->Branch("minPl", &minPl, "minPl/I");
    fOutTree->Branch("maxPl", &maxPl, "maxPl/I");
    fOutTree->Branch("minPlX", &minPlX, "minPlX/I");
    fOutTree->Branch("minPlY", &minPlY, "minPlY/I");
    fOutTree->Branch("maxPlX", &maxPlX, "maxPlX/I");
    fOutTree->Branch("maxPlY", &maxPlY, "maxPlY/I");
    fOutTree->Branch("ncellsfitT", &ncellsfitT, "ncellsfitT/I");
    fOutTree->Branch("ncellsfitX", &ncellsfitX, "ncellsfitX/I");
    fOutTree->Branch("ncellsfitY", &ncellsfitY, "ncellsfitY/I");
    fOutTree->Branch("evtTime", &evtTime, "evtTime/i");

    fOutTree->Branch("dcmnum", &dcmnum);
    fOutTree->Branch("dcmhits", &dcmhits);
    fOutTree->Branch("dcmsimtime", &dcmsimtime);
    if (fDebug){
      fOutTree->Branch("planenum", &planenum);
      fOutTree->Branch("cellnum", &cellnum);
      fOutTree->Branch("dcmsimresid", &dcmsimresid);
      fOutTree->Branch("dcmtime", &dcmtime);
      fOutTree->Branch("dcmstime", &dcmstime);
      fOutTree->Branch("dcmrd", &dcmrd);
      fOutTree->Branch("dcmpe", &dcmpe);
      fOutTree->Branch("dcmpl",&dcmpl);
    }


  }

  //----------------------------------------------------------------------
  void TimingCalibration::beginRun(art::Run&)
  {
    art::ServiceHandle<geo::Geometry> geom;

    fhicl::ParameterSet pset;
    switch(geom->DetId()){
    case novadaq::cnv::kNEARDET:
      pset = fPSetND;
      break;
    case novadaq::cnv::kFARDET:
      pset = fPSetFD;
      break;
    case novadaq::cnv::kTESTBEAM:
      pset = fPSetTB;
      break;
    default:
      assert(0 && "Unknown detector");
    }

    fMinDCMHits     = pset.get< int         >("MinDCMHits");
    fMinDCM         = pset.get< int         >("MinDCM");

  }

  //---------------------------------------------------------------------------
  void TimingCalibration::produce(art::Event & evt)
  {
    art::ServiceHandle<calib::Calibrator> cal;
    art::ServiceHandle<cheat::BackTracker> bt;

    std::unique_ptr< std::vector<caldp::DCMStat>  > dcmstats(new std::vector<caldp::DCMStat>);
      
    std::unique_ptr<art::Assns<caldp::DCMStat, rb::Track> > dcmAssns(new art::Assns<caldp::DCMStat, rb::Track>);

    art::ServiceHandle<geo::Geometry>         geom;
    daqchannelmap::DAQChannelMap* cmap = daqchannelmap::DAQChannelMap::getInstance(geom->DetId());
    novadaq::cnv::DetId detId = geom->DetId();

    //get the cosmic tracks associated with each event
    art::Handle<std::vector<rb::Track> > trackcol;
    evt.getByLabel(fTrackLabel,trackcol);

    //get associationgs between slices and tracks
    art::FindManyP<caldp::TCTrack> fmtc(trackcol, evt, fTCTrackLabel);

    std::vector< art::Ptr<rb::Track> > tracks;
    art::fill_ptr_vector(tracks, trackcol);

    //store event info in ntuple
    run = evt.run();
    subrun = evt.subRun();
    event = evt.id().event();

    evtTime  = evt.time().timeHigh();

    //now loop over the tracks
    for (unsigned int i=0; i<tracks.size(); ++i){

      this->clearNTuple();

      const art::Ptr<rb::Track> tr = tracks.at(i);

      //if ((tr->NXCell()==0) || (tr->NYCell()==0)) continue;
      //get the tctrack associated with the cosmic track
      std::vector<art::Ptr<caldp::TCTrack> >  tct = fmtc.at(i);
      //if we failed to get tctracks skip s
      if (!(fmtc.isValid())) continue;
      //there should only be one tctrack associated with a reco track
      if (tct.size() != 1) continue;

      //make selection cuts on tracks
      fDCMHitMap.clear();


      cosmictr = i; //store track number
      //passed initial round of cuts, now get track information
      TVector3 start(tr->Start());
      TVector3 end(tr->Stop());
      sX = start.X();
      sY = start.Y();
      sZ = start.Z();
      eX = end.X();
      eY = end.Y();
      eZ = end.Z();
      plength = (end - start).Mag();
      minPl = tct[0]->MinPlane(geo::kXorY);
      maxPl = tct[0]->MaxPlane(geo::kXorY);

      minPlX = tct[0]->MinPlane(geo::kX);
      minPlY = tct[0]->MinPlane(geo::kY);
      maxPlX = tct[0]->MaxPlane(geo::kX);
      maxPlY = tct[0]->MaxPlane(geo::kY);

      //make other selection cuts if desired


      //put hits into map connected to DCM
      for(unsigned int j=0; j<tct[0]->NPCHit(); ++j){
        art::Ptr<caldp::PCHit> phit = tct[0]->PCHit(j);
        int key = cmap->getNumberOfDCM(daqchannelmap::AFULLAFULL_DIBLOCK)*(phit->Diblock()-1) + phit->DCM();
        fDCMHitMap[key].push_back(phit);
      }

      //for each dcm calculate timing offsets
      caldp::DCMStat dcms; //create object to hold dcm info for the track
      for(std::map<int, std::vector<art::Ptr<caldp::PCHit> > >::iterator itr = fDCMHitMap.begin(); itr!=fDCMHitMap.end(); ++itr){
        std::vector<art::Ptr<caldp::PCHit> > dcmcells = itr->second;
        if ((int)dcmcells.size() < fMinDCMHits) continue;

        caldp::DCMTime dcmTime;

        std::vector<float> cellnumtemp;
        std::vector<float> planenumtemp;
        std::vector<float> dcmtimetemp;
        std::vector<float> dcmstimetemp;
        std::vector<float> dcmrdtemp;
        std::vector<float> dcmpetemp;
        std::vector<float> dcmpltemp;
        std::vector<double> weights;
        std::vector<double> sigma;
        for(unsigned int k=0; k<dcmcells.size(); ++k){
          double tns = dcmcells[k]->TNS();
          double rd = dcmcells[k]->ReadDist();
          double pl = dcmcells[k]->FlightLen();
          if (fCorrectTimewalk) tns += tct[0]->ReadoutDistCorrection(dcmcells[k]->ReadDist(), detId);
          //reverse any timing calibration that has been baked in
          tns -= cal->GetTimingOffset(dcmcells[k]->Plane(), dcmcells[k]->Cell(), !bt->HaveTruthInfo());
          float tmpFiberSpeed = fFiberSpeed;
          if (fAdjustSpeed) tmpFiberSpeed = tct[0]->CalcFiberVelocity(rd, detId);
          dcmstimetemp.push_back(tns - pl/sCmPerNsec - rd/tmpFiberSpeed);
          dcmtimetemp.push_back(tns);
          dcmrdtemp.push_back(rd);
          dcmpetemp.push_back(dcmcells[k]->PE());
          dcmpltemp.push_back(pl);
          cellnumtemp.push_back(dcmcells[k]->Cell());
          planenumtemp.push_back(dcmcells[k]->Plane());
          double stdev = tct[0]->TNSUncertainty(dcmcells[k]->PE(), dcmcells[k]->GoodTime(), detId);
          sigma.push_back(stdev);
          double w = 1.0/(stdev*stdev);
          weights.push_back(w);
        }

        while (true){ 
          double st = 0.0;
          double weight = 0.0;
          for(unsigned int k=0; k<dcmstimetemp.size(); ++k){
            st += dcmstimetemp[k]*weights[k];
            weight += weights[k];
          }
          st /= weight;
          double maxresid = 0;
          unsigned int maxpos = 0;
          for(unsigned int k=0; k<dcmstimetemp.size(); ++k){
            float recresid = st - dcmstimetemp[k];
            if (fabs(recresid/sigma[k]) > fabs(maxresid)){
              maxresid = fabs(recresid/sigma[k]);
              maxpos = k;
            }
          }
          if (dcmcells.size() <= 2) break;
          if (maxresid > 2.0){
            if (fRemoveOutlier){
              dcmcells.erase(dcmcells.begin()+maxpos);
              cellnumtemp.erase(cellnumtemp.begin()+maxpos);
              planenumtemp.erase(planenumtemp.begin()+maxpos);
              dcmtimetemp.erase(dcmtimetemp.begin()+maxpos);
              dcmstimetemp.erase(dcmstimetemp.begin()+maxpos);
              dcmrdtemp.erase(dcmrdtemp.begin()+maxpos);
              dcmpetemp.erase(dcmpetemp.begin()+maxpos);
              dcmpltemp.erase(dcmpltemp.begin()+maxpos);
              weights.erase(weights.begin()+maxpos);
              sigma.erase(sigma.begin()+maxpos);
              continue;
            }
          }
          if ((int)dcmcells.size() < fMinDCMHits) break;
          dcmnum.push_back(itr->first);
          dcmhits.push_back(dcmstimetemp.size());
          dcmsimtime.push_back(st);
          dcmTime.dcmTime = st;
          dcmTime.dcmNum = itr->first;  
          dcmTime.dcmHits = dcmcells.size();
          double sum = 0.0;
          for(unsigned int k=0; k<dcmstimetemp.size(); ++k){
            double recresid = st - dcmstimetemp[k];
            dcmsimresid.push_back(recresid);
            sum += weights[k]*(dcmstimetemp[k] - st)*(dcmstimetemp[k]-st);
            cellnum.push_back(cellnumtemp[k]);
            planenum.push_back(planenumtemp[k]);
            dcmtime.push_back(dcmtimetemp[k]);
            dcmstime.push_back(dcmstimetemp[k]);
            dcmrd.push_back(dcmrdtemp[k]);
            dcmpe.push_back(dcmpetemp[k]);
            dcmpl.push_back(dcmpltemp[k]);
          }
          dcmTime.sigma = sqrt(sum/weight);
          dcms.dcmList.push_back(dcmTime);
          break;
        }

      }

      if ((int)dcms.dcmList.size() < fMinDCM) continue;

      dcmstats->push_back(dcms);

      util::CreateAssn(*this, evt, *dcmstats, tr, *dcmAssns);

      
        
      ncellsfitT = tct[0]->NPCHit();
      ncellsfitX = tct[0]->NXPCHit();
      ncellsfitY = tct[0]->NYPCHit();

      
      fOutTree->Fill();
      
    
    }//end loop over tracks

    evt.put(std::move(dcmstats));
    evt.put(std::move(dcmAssns));

  } 

  //---------------------------------------------------------------------------
  void TimingCalibration::clearNTuple(){

    //track start coordinate
    sX = -999.0;
    sY = -999.0;
    sZ = -999.0;
    //track end coordinate
    eX = -999.0;
    eY = -999.0;
    eZ = -999.0;
    plength = -5.0;
    //min, max plane
    minPl = -5;
    maxPl = -5;
    //min, max cell
    minPlX = -5;
    minPlY = -5;
    maxPlX = -5;
    maxPlY = -5;
    //number of cells used in final fiber speed fit
    ncellsfitT = 0;
    ncellsfitX = 0;
    ncellsfitY = 0;

  
    dcmnum.clear();
    dcmhits.clear();
    dcmsimtime.clear();
    planenum.clear();
    cellnum.clear();
    dcmsimresid.clear();
    dcmtime.clear();
    dcmstime.clear();
    dcmrd.clear();
    dcmpe.clear();
    dcmpl.clear();

  }

} //end namespace

//-----------------------------------------------------------------------------

DEFINE_ART_MODULE(calib::TimingCalibration)