CalibAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       CalibAna
// Plugin Type: analyzer (art v2_12_01)
// File:        CalibAna_module.cc
//
// Generated at Sat Apr 27 09:35:48 2019 by Tyler Alion using cetskelgen
// from cetlib version v3_06_00.
////////////////////////////////////////////////////////////////////////

#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 "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Persistency/Common/FindManyP.h"

#include "TTree.h"

#include "CalibrationDataProducts/PCHit.h"
#include "RecoBase/Track.h"
#include "RecoBase/CellHit.h"
//#include "MCCheater/BackTracker.h"
#include "Calibrator/Calibrator.h"

//#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/PlaneGeo.h"

namespace calib {
  struct CalibAnaParams
  {
    template<class T> using Atom = fhicl::Atom<T>;
    template<class T> using Table = fhicl::Table<T>;
    using Comment = fhicl::Comment;
    using Name = fhicl::Name;
    Atom<std::string> TrackLabel{
        Name("TrackLabel"),
        Comment("Label of track reconstruction module (all tracks or stopping seleciton)")
        };
    Atom<std::string> PCHitsLabel{
        Name("PCHitsLabel"),
        Comment("Label of PCHitsList module (pclist or pcliststop)")
        };
    Atom< std::string > WhichHitTrees{
        Name("WhichHitTrees"),
          Comment("Which Hit Trees to make. All: \"tricell, trajectory, cellhit\"")
        };
  };
}

namespace caldp {
  class PCHit;
}

namespace calib {
  class CalibAna;
}

class calib::CalibAna : public art::EDAnalyzer {
public:
  using Parameters = art::EDAnalyzer::Table<CalibAnaParams>;

  explicit CalibAna(const Parameters& params);
  // The compiler-generated destructor is fine for non-base
  // classes without bare pointers or other resource use.

  // Plugins should not be copied or assigned.
  CalibAna(CalibAna const &) = delete;
  CalibAna(CalibAna &&) = delete;
  CalibAna & operator = (CalibAna const &) = delete;
  CalibAna & operator = (CalibAna &&) = delete;

  // Required functions.
  void analyze(art::Event const & e) override;

  // Selected optional functions.
  void beginJob() override;
  void endJob() override;

private:
  CalibAnaParams fParams;

  void AddEventBranches(TTree * t);
  void ResetEventBranches();
  TTree * fEventTree;
  bool fIsRealData;
  int et_run, et_subrun, et_event;
  uint32_t et_evttime;
  int et_ntracks;
  int et_ncellhits;
  int et_npchits_tricell; 
  int et_npchits_traj;

  void AddTrackBranches(TTree * t);
  void ResetTrackBranches();
  TTree * fTrackTree;
  float tt_dcosz, tt_dcosy, tt_dcosx;
  float tt_totlength;
  float tt_totdisplace;
  int   tt_pdg;
  int   tt_ncellhits;
  int   tt_npchits_tricell; 
  int   tt_npchits_avg; 
  int   tt_npchits_z; 
  int   tt_npchits_traj;
  TVector3 tt_startpos;
  TVector3 tt_endpos;

  void AddHitBranches(TTree * t, std::string type);
  void ResetHitBranches();
  std::map< std::string, TTree * > fHitTree;
  bool  ht_istricell;
  bool  ht_istraj;
  bool  ht_isavg;
  bool  ht_isz;
  int   ht_diblock;
  int   ht_dcm;
  int   ht_apd;
  int   ht_pixel;
  int   ht_plane;  
  int   ht_cell;
  int   ht_view;
  float ht_w;
  float ht_readdist;
  float ht_xapprox, ht_yapprox, ht_zapprox;
  float ht_pe;
  float ht_pecorr;
  float ht_pecorrtomev;
  float ht_recoMeV;
  float ht_path;
  float ht_flightLen;
  float ht_cmFromEnd;
  float ht_tns;
  float ht_truepath;
  float ht_truepe;
  float ht_truew;
  float ht_trueMeV;
  float ht_poissonlambda;
  int   ht_brightbin;
  int   ht_tdc; // In CellHits,
  int   ht_adc; //  but not PCHits.
  void FillWithPCHit( caldp::PCHit const &pchit, 
                      bool is_trajectory, bool is_tricell,
                      bool is_avg, bool is_z,
                      float trk_len );
  void FillWithCellHit( art::Ptr< rb::CellHit > cellhit,
                      bool is_trajectory, bool is_tricell,
                      bool is_avg, bool is_z,
                      float trk_len );  // less info than PCHits

  // Use W from PCHit to get PECorr from calibrator
  float GetPECorrToMeV( caldp::PCHit const &pchit );
  float GetPECorr( caldp::PCHit const &pchit);

  void PrintHit ( caldp::PCHit const &hit );
  bool HitInList( caldp::PCHit const &hit,
                  std::vector<art::Ptr<caldp::PCHit> > list );
  bool HitInList( art::Ptr< rb::CellHit > hit,
                  std::vector<art::Ptr<caldp::PCHit> > list );
  bool MatchHits( caldp::PCHit const &h1,
                  caldp::PCHit const &h2 );
  bool MatchHits( art::Ptr< rb::CellHit > h1,
                  caldp::PCHit const &h2 );
  void ApproximateHitPos( caldp::PCHit const &hit,
                          float & x, float & y, float & z );
  void ApproximateHitPos( art::Ptr< rb::CellHit > cellhit,
                          float & x, float & y, float & z );
};


calib::CalibAna::CalibAna(const Parameters& params)
  :
  EDAnalyzer(params),
  fParams(params())
{}



// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void calib::CalibAna::analyze(art::Event const & e)
{
  fIsRealData = e.isRealData();
  ResetEventBranches();
  et_run = (int)(e.run());
  et_subrun = (int)(e.subRun());
  et_event = (int)(e.id().event());
  et_evttime = e.time().timeHigh();
  
  art::Handle<std::vector<rb::Track> > tracks;
  e.getByLabel(fParams.TrackLabel(), tracks);
  et_ntracks = tracks->size();


  // Annoyingly need all since (XY + Z + Avg) = all hits in track
  // FindTricellHits gets the main hits to be put in tree
  art::FindManyP<caldp::PCHit> FindTricellHits(tracks, e, 
                                               art::InputTag(fParams.PCHitsLabel(),
                                                             "pathQualXY")); 
  art::FindManyP<caldp::PCHit> FindTrajHits(tracks, e, 
                                            art::InputTag(fParams.PCHitsLabel(),
                                                          "pathQualTraj")); 
  art::FindManyP<caldp::PCHit> FindAvgHits( tracks, e, 
                                            art::InputTag(fParams.PCHitsLabel(),
                                                          "pathQualAvg")); 
  art::FindManyP<caldp::PCHit> FindZHits(   tracks, e, 
                                            art::InputTag(fParams.PCHitsLabel(),
                                                          "pathQualZ")); 
  if( !FindTricellHits.isValid() ){
    mf::LogInfo("CalibAna") << "No pathQualXY hits."; return; }
  if( !FindTrajHits.isValid() ){
    mf::LogInfo("CalibAna") << "No pathQualTraj hits."; return; }
  if( !FindAvgHits.isValid() ){
    mf::LogInfo("CalibAna") << "No pathQualAvg hits."; return; }
  if( !FindZHits.isValid() ){
    mf::LogInfo("CalibAna") << "No pathQualZ hits."; return; }
  

  for( unsigned int i=0; i<tracks->size(); ++i ){
    ResetTrackBranches();
    rb::Track trk = (*tracks)[i];
    tt_dcosx = trk.Dir().X()/trk.Dir().Mag();
    tt_dcosy = trk.Dir().Y()/trk.Dir().Mag();
    tt_dcosz = trk.Dir().Z()/trk.Dir().Mag();  

    tt_startpos = trk.Start();
    tt_endpos   = trk.Stop();
    tt_totdisplace = (tt_endpos-tt_startpos).Mag(); // displacement
    tt_totlength = trk.TotalLength();               // vs. distance

    // xy and z conditions are mutually exclusive,
    // avg hits are just the rest in the track
    //   so xy + z + avg is all the hits in the track
    std::vector<art::Ptr<caldp::PCHit> > trihits  = FindTricellHits.at(i);
    std::vector<art::Ptr<caldp::PCHit> > trajhits = FindTrajHits.at(i);
    std::vector<art::Ptr<caldp::PCHit> > avghits  = FindAvgHits.at(i);
    std::vector<art::Ptr<caldp::PCHit> > zhits    = FindZHits.at(i);
    //int npchits = trihits.size() + avghits.size() + zhits.size();
    tt_ncellhits = trk.NCell();
    //if(tt_ncellhits != npchits)
    //mf::LogInfo("CalibAna") << tt_ncellhits << " CellHits, but " << npchits << " PCHits";
    tt_npchits_tricell = trihits.size();
    tt_npchits_avg = avghits.size();
    tt_npchits_z = zhits.size();
    tt_npchits_traj = trajhits.size();
    //if(tt_ncellhits != tt_npchits_traj)
    //mf::LogInfo("CalibAna") << tt_ncellhits << " CellHits, but " << tt_npchits_traj << " Trajectory PCHits";    
    et_ncellhits += tt_ncellhits;
    et_npchits_tricell += tt_npchits_tricell;
    et_npchits_traj += tt_npchits_traj;

    /*
    std::cout << npchits << " PCHits (" 
              << tt_npchits_tricell << " Tricell, "
              << tt_npchits_avg << " Avg, "
              << zhits.size() << " Z), "
              << tt_npchits_traj << " Trajectory. "
              << trk.NCell() << " Cell Hits." << std::endl;
    */
 

    //daqchannelmap::DAQChannelMap* cmap = daqchannelmap::DAQChannelMap::getInstance(fDS->DetId());

    if(fHitTree.count("cellhit")==1){
      for(size_t c = 0; c < trk.NCell(); ++c){
        art::Ptr< rb::CellHit > cellhit = trk.Cell(c);

        // Check what hit qualities this cellhit has
        bool is_avg  = HitInList( cellhit, avghits  );
        bool is_z    = HitInList( cellhit, zhits    );
        bool is_traj = HitInList( cellhit, trajhits );
        bool is_tri  = HitInList( cellhit, trihits  );

        FillWithCellHit(cellhit, 
                        is_traj, is_tri,
                        is_avg, is_z,
                        tt_totdisplace); 
        fHitTree["cellhit"]->Fill();
      }
    }

 
    if(fHitTree.count("tricell")==1){
      for(size_t h = 0; h < trihits.size(); ++h){
        caldp::PCHit const & trihit = *trihits[h];

        // Check what other hit qualities this trajectory hit has
        bool is_avg  = HitInList( trihit, avghits  );
        bool is_z    = HitInList( trihit, zhits    );
        bool is_traj = HitInList( trihit, trajhits );
        
        FillWithPCHit( trihit,
                       is_traj, true,
                       is_avg, is_z,
                       tt_totdisplace); 
        fHitTree["tricell"]->Fill();
      }
    }

    if(fHitTree.count("trajectory")==1){
      size_t ntricell=0;
      for(size_t h = 0; h < trajhits.size(); ++h){
        caldp::PCHit const & trajhit = *trajhits[h];

        // Check what other hit qualities this trajectory hit has
        bool is_avg = HitInList( trajhit, avghits );
        bool is_z   = HitInList( trajhit, zhits   );
        bool is_tri = HitInList( trajhit, trihits );
        if( is_tri ) ntricell++; 

        FillWithPCHit( trajhit,
                       true, is_tri,
                       is_avg, is_z,
                       tt_totdisplace); 
        fHitTree["trajectory"]->Fill();
      }

      if(tt_dcosz > 0.95 && ntricell!=0 ){
        std::cout << "^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^\n\n" << std::endl;
        std::cout << "   Forward Track with " << ntricell << " tricell hits.\n\n" << std::endl;
        for(size_t h = 0; h < trajhits.size(); ++h){
          caldp::PCHit const & thit = *trajhits[h];
          if( HitInList( thit, trihits ) ) std::cout << "Tricell: " << std::endl;
          PrintHit(thit);
        }
      }

    }
     

    fTrackTree->Fill();
  }

  fEventTree->Fill();
}


// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void calib::CalibAna::PrintHit(caldp::PCHit const &hit){
  std::cout << std::right << std::setw(6) << std::setfill(' ')
            << "P:"  << hit.Plane()   << std::right << std::setw(6) << std::setfill(' ')
            << "C:"  << hit.Cell()    << std::right << std::setw(6) << std::setfill(' ')
            << "W:"  << hit.W()       << std::right << std::setw(6) << std::setfill(' ')    
            << "T:"  << hit.TNS()     << std::right << std::setw(6) << std::setfill(' ')
            << std::endl;
}



// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
bool calib::CalibAna::HitInList(caldp::PCHit const &hit,
                                std::vector<art::Ptr<caldp::PCHit> > list)
{
  for(size_t i = 0; i < list.size(); ++i)
    if( MatchHits( hit, *list[i] ) )
      return true;
  return false;
}
bool calib::CalibAna::HitInList(art::Ptr< rb::CellHit > cellhit, 
                                std::vector<art::Ptr<caldp::PCHit> > list)
{
  for(size_t i = 0; i < list.size(); ++i)
    if( MatchHits( cellhit, *list[i] ) )
      return true;
  return false;
}
bool calib::CalibAna::MatchHits(caldp::PCHit const &h1, 
                                caldp::PCHit const &h2)
{
  return (h1.Cell()  == h2.Cell()  &&
          h1.Plane() == h2.Plane() &&
          h1.TNS()   == h2.TNS());
}
bool calib::CalibAna::MatchHits(art::Ptr< rb::CellHit > h1, 
                                caldp::PCHit const &h2)
{
  return (h1->Cell()  == h2.Cell()  &&
          h1->Plane() == h2.Plane() &&
          h1->TNS()   == h2.TNS());
}


// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void calib::CalibAna::beginJob()
{
  art::ServiceHandle<art::TFileService> tfs;
  fEventTree = tfs->make<TTree>("event_tree","tree of variables per art event");
  AddEventBranches(fEventTree);

  fTrackTree = tfs->make<TTree>("track_tree","tree of variables per track");
  AddEventBranches(fTrackTree);
  AddTrackBranches(fTrackTree);

  std::vector< std::string > types = {"tricell","trajectory","cellhit"};
  for(auto& type: types){
    if(fParams.WhichHitTrees().find(type)==std::string::npos) continue;

    fHitTree[type] = tfs->make<TTree>((type+"_tree").c_str(),
                                      ("tree of variables per "+type+" hit").c_str());
    AddEventBranches(fHitTree[type]);
    AddTrackBranches(fHitTree[type]);
    AddHitBranches(fHitTree[type],type);
  }

}


// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void calib::CalibAna::FillWithPCHit(caldp::PCHit const &pc,
                                    bool is_trajectory, bool is_tricell,
                                    bool is_avg, bool is_z,
                                    float trk_length)
{
  ResetHitBranches();
  // straight from the hit
  ht_diblock  = pc.Diblock();
  ht_dcm      = pc.DCM();
  ht_apd      = pc.APD();
  ht_pixel    = pc.Pixel();
  ht_plane    = pc.Plane();
  ht_cell     = pc.Cell();
  ht_view     = pc.View();
  ht_w        = pc.W();
  ht_readdist = pc.ReadDist();
  ht_pe       = pc.PE();
  ht_path     = pc.Path();
  ht_tns      = pc.TNS();
  ht_truepath = pc.TruePath();
  ht_truepe   = pc.TruePE();
  ht_truew    = pc.TrueW();
  ht_trueMeV  = pc.TrueMeV();
  ht_flightLen = pc.FlightLen(); // comes from trk->DistanceFromStart( world coordinate );
  ht_cmFromEnd = trk_length - pc.FlightLen();
  ht_poissonlambda = pc.PoissonLambda();
  // use calibrator service
  ht_pecorr   = GetPECorr(pc);
  ht_pecorrtomev = GetPECorrToMeV(pc);
  ht_recoMeV  = ht_pecorr * ht_pecorrtomev;
  art::ServiceHandle<photrans::FiberBrightness> fb;
  ht_brightbin = fb->getBrightnessIndex(ht_plane, ht_cell);

  size_t ntypes = 0;
  if( is_tricell ) ntypes++;
  if( is_avg     ) ntypes++;
  if( is_z       ) ntypes++;
  if( ntypes > 1 ) 
    mf::LogError("CalibAna") 
      << "Tricell, Avg, and Z selections should be mutually exclusive, "
      << "but this hit is: " 
      << (is_tricell ? "Tricell, " : "")
      << (is_avg ? "Avg, " : "")
      << (is_z ? "Z. " : ".");
 
  ht_istraj    = is_trajectory;
  ht_istricell = is_tricell;
  ht_isavg     = is_avg;
  ht_isz       = is_z;

  ApproximateHitPos( pc, ht_xapprox, ht_yapprox, ht_zapprox );
}
// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void calib::CalibAna::FillWithCellHit(art::Ptr< rb::CellHit > ch,
                                      bool is_trajectory, bool is_tricell,
                                      bool is_avg, bool is_z,
                                      float trk_length)
{
  ResetHitBranches();
  // straight from the hit

  ht_plane    = ch->Plane();
  ht_cell     = ch->Cell();
  /*
  uint32_t lchan   = cmap->encodeLChan(fGeom->DetId(), chit->Plane(), chit->Cell());
  uint32_t dchan   = cmap->encodeDChan(lchan);
  uint32_t pixel   = cmap->getPixel(dchan);
  uint32_t apd     = cmap->getFEB(dchan);
  uint32_t dcm     = cmap->getDCM(dchan);
  uint32_t diblock = cmap->getDiBlock(dchan);
  */ 
  //ht_diblock  = ch->Diblock();
  //ht_dcm      = ch->DCM();
  //ht_apd      = ch->APD();
  //ht_pixel    = ch->Pixel();

  ht_view     = ch->View();
  //ht_w        = ch->W();
  ht_pe       = ch->PE();
  //ht_path     = ch->Path();
  ht_tns      = ch->TNS();
  //ht_truepath = ch->TruePath();
  //ht_truepe   = ch->TruePE();
  //ht_truew    = ch->TrueW();
  //ht_trueMeV  = ch->TrueMeV();
  //ht_x        = trk_length - ch->FlightLen();
  ht_tdc      = ch->TDC();
  ht_adc      = ch->ADC();
  // use calibrator service
  //ht_pecorr   = GetPECorr(ch); // need W from PCHits
  //ht_pecorrtomev = GetPECorrToMeV(ch);
  //ht_recoMeV  = ht_pecorr * ht_pecorrtomev;
  art::ServiceHandle<photrans::FiberBrightness> fb;
  ht_brightbin = fb->getBrightnessIndex(ht_plane, ht_cell);

  ht_istraj    = is_trajectory;
  ht_istricell = is_tricell;
  ht_isavg     = is_avg;
  ht_isz       = is_z;

  ApproximateHitPos( ch, ht_xapprox, ht_yapprox, ht_zapprox );
}


// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
float calib::CalibAna::GetPECorr(caldp::PCHit const &pchit){
  //Make mock cellhit
  rb::CellHit cellhit;
  cellhit.SetPlane(pchit.Plane());
  cellhit.SetCell(pchit.Cell());
  cellhit.SetView(pchit.View());
  cellhit.SetTNS(pchit.TNS(),pchit.GoodTime());
  if(!fIsRealData) cellhit.SetMC();
  cellhit.SetPE(pchit.PE());
  art::ServiceHandle<calib::Calibrator> calibrator;
  return calibrator->GetPECorr(cellhit, pchit.W());
}
float calib::CalibAna::GetPECorrToMeV(caldp::PCHit const &pchit){
  rb::CellHit cellhit;
  if(!fIsRealData) cellhit.SetMC();
  art::ServiceHandle<calib::Calibrator> calibrator;
  return 1000.*calibrator->GetPECorrToGeVScale(cellhit);
}


// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void calib::CalibAna::ApproximateHitPos( caldp::PCHit const &hit,
                                         float & x, float & y, float & z ){
  art::ServiceHandle<geo::Geometry> geo;
  const geo::PlaneGeo* planeg = geo->Plane(hit.Plane());
  const geo::CellGeo*  cellg  = planeg->Cell(hit.Cell());
  double cell_center[3];
  cellg->GetCenter(cell_center);

  z = cell_center[2];
  if( hit.View()==geo::kX ){
    x = cell_center[0];
    y = hit.W();
  }
  else if( hit.View()==geo::kY ){
    x = hit.W();
    y = cell_center[1];
  }
  else{
    mf::LogWarning("CalibAna") << "Bad hit.View()";
  }
}
void calib::CalibAna::ApproximateHitPos( art::Ptr< rb::CellHit > cellhit,
                                         float & x, float & y, float & z ){
  art::ServiceHandle<geo::Geometry> geo;
  const geo::CellGeo* cellg = geo->Plane(cellhit->Plane())->Cell(cellhit->Cell());
  double cell_center[3];
  cellg->GetCenter(cell_center);

  z = cell_center[2];
  if( cellhit->View()==geo::kX ){
    x = cell_center[0];
    y = -999.; // estimate W?
  }
  else if( cellhit->View()==geo::kY ){
    x = -999.; // estimate W?
    y = cell_center[1];
  }
  else{
    mf::LogWarning("CalibAna") << "Bad hit.View()";
  }
}


// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void calib::CalibAna::endJob()
{
  art::ServiceHandle<calib::Calibrator> calibrator;
  CalibratorParams cal_par = calibrator->GetCalibratorParams();

  std::string baserel = getenv("SRT_PUBLIC_CONTEXT");
  std::string testrel = getenv("SRT_PRIVATE_CONTEXT");
  std::string csvpath_atten  = calibrator->GetAttenCSVPath();
  std::string csvpath_abs    = calibrator->GetAbsConstsCSVPath();
  std::string csvpath_timing = calibrator->GetTimingConstsCSVPath();
  std::string csvpath_drift  = calibrator->GetDriftConstsCSVPath();
  std::string calibtag       = cal_par.Tag();
  bool        UseAbsEpochs    = cal_par.UseAbsEpochs();
  std::string AbsEpochTag     = cal_par.AbsEpochTag();
  std::string TimingTag       = cal_par.TimingTag();
  bool        UseTimingEpochs = cal_par.UseTimingEpochs();
  std::string TimingEpochTag  = cal_par.TimingEpochTag();
  std::string DriftCSV        = cal_par.DriftCSV();
  bool        UseDriftEpochs  = cal_par.UseDriftEpochs();
  std::string DriftEpochTag   = cal_par.DriftEpochTag();
  
  art::ServiceHandle<art::TFileService> tfs;  
  TTree* FileTree = tfs->make<TTree>("file_tree","tree of info per file");
  FileTree->Branch( "baserel",        &baserel        ); // $SRT_PUBLIC_CONTEXT 
  FileTree->Branch( "testrel",        &testrel        ); // $SRT_PRIVATE_CONTEXT
  FileTree->Branch( "csvpath_atten",  &csvpath_atten  ); // $CALIBCSVS_CSV_PATH or fcl provided
  FileTree->Branch( "csvpath_abs",    &csvpath_abs    ); 
  FileTree->Branch( "csvpath_timing", &csvpath_timing ); 
  FileTree->Branch( "csvpath_drift",  &csvpath_drift  ); 
  FileTree->Branch( "calibtag",       &calibtag       );
  
  // Other potentially important configurations
  //   TODO: Loop through parameter table and use key
  //         as branch name to do a full config dump.
  FileTree->Branch( "UseAbsEpochs", &UseAbsEpochs );
  FileTree->Branch( "AbsEpochTag", &AbsEpochTag );
  FileTree->Branch( "TimingTag", &TimingTag );
  FileTree->Branch( "UseTimingEpochs", &UseTimingEpochs );
  FileTree->Branch( "TimingEpochTag", &TimingEpochTag );
  FileTree->Branch( "DriftCSV", &DriftCSV );
  FileTree->Branch( "UseDriftEpochs", &UseDriftEpochs );
  FileTree->Branch( "DriftEpochTag", &DriftEpochTag );
  //FileTree->Branch("brightnessmodel",&brightnessmodel); // 


  FileTree->Fill();
}



// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void calib::CalibAna::AddEventBranches(TTree * t)
{
  t->Branch("run",&et_run);
  t->Branch("subrun",&et_subrun);
  t->Branch("event",&et_event);
  t->Branch("evt_ntracks",&et_ntracks);
  t->Branch("evt_ncellhits",&et_ncellhits);
  t->Branch("evt_nhits_tricell",&et_npchits_tricell);
  t->Branch("evt_nhits_traj",&et_npchits_traj);
  t->Branch("evt_time",&et_evttime);
}
void calib::CalibAna::ResetEventBranches()
{
  et_run=-5; et_subrun=-5; et_event=-5;
  et_ntracks=-5; et_ncellhits=-5; 
  et_npchits_tricell=-5; et_npchits_traj=5;
  et_evttime=-5;
}



// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void calib::CalibAna::AddTrackBranches(TTree * t)
{
  t->Branch("trk_nhits_tricell",&tt_npchits_tricell);
  t->Branch("trk_nhits_avg",&tt_npchits_avg);
  t->Branch("trk_nhits_z",&tt_npchits_z);
  t->Branch("trk_nhits_traj",&tt_npchits_traj);
  t->Branch("trk_ncellhits",&tt_ncellhits);
  //t->Branch("npts",&tt_npts);
  t->Branch("trk_dcosx",&tt_dcosx);
  t->Branch("trk_dcosy",&tt_dcosy);
  t->Branch("trk_dcosz",&tt_dcosz);
  t->Branch("trk_totlength",&tt_totlength);
  t->Branch("trk_totdisplace",&tt_totlength);
  t->Branch("trk_startpos",&tt_startpos);
  t->Branch("trk_endpos",&tt_endpos);
  t->Branch("trk_pdg",&tt_pdg);
}
void calib::CalibAna::ResetTrackBranches()
{
  tt_npchits_tricell=0;
  tt_npchits_avg=0;
  tt_npchits_z=0;
  tt_npchits_traj=0;
  tt_ncellhits=0;
  //tt_npts=0;
  tt_pdg=0;
  tt_dcosx=-5; tt_dcosy=-5; tt_dcosz=-5;
  tt_totlength=-5; tt_totdisplace=-5;
  tt_startpos.SetXYZ(0.,0.,0.);
  tt_endpos.SetXYZ(0.,0.,0.);
}


// ^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^~^
void calib::CalibAna::AddHitBranches(TTree * t, 
                                     std::string type)
{
  // Hit vars
  t->Branch("istricell",&ht_istricell);
  t->Branch("istraj",&ht_istraj);
  t->Branch("isavg",&ht_isavg);
  t->Branch("isz",&ht_isz);
  t->Branch("diblock",&ht_diblock);
  t->Branch("dcm",&ht_dcm);
  t->Branch("apd",&ht_apd);
  t->Branch("pixel",&ht_pixel);
  t->Branch("plane",&ht_plane);  
  t->Branch("cell",&ht_cell);
  t->Branch("w",&ht_w);
  t->Branch("readdist",&ht_readdist);
  t->Branch("xapprox",&ht_xapprox);
  t->Branch("yapprox",&ht_yapprox);
  t->Branch("zapprox",&ht_zapprox);
  t->Branch("pe",&ht_pe);
  t->Branch("pecorr",&ht_pecorr);
  t->Branch("pecorrtomev",&ht_pecorrtomev); 
  t->Branch("recoMeV",&ht_recoMeV);
  t->Branch("path",&ht_path);
  t->Branch("cmFromStart",&ht_flightLen); // cm from track start
  t->Branch("cmFromEnd",&ht_cmFromEnd); // cm from track endpt
  t->Branch("view",&ht_view);
  t->Branch("tns",&ht_tns);
  t->Branch("truepath",&ht_truepath);
  t->Branch("truepe",&ht_truepe);
  t->Branch("truew",&ht_truew);
  t->Branch("trueMeV",&ht_trueMeV);
  t->Branch("poissonLambda",&ht_poissonlambda);
  t->Branch("brightbin",&ht_brightbin);
  //if(type=="cellhit"){
    t->Branch("tdc",&ht_tdc); // time of last baseline sample
    t->Branch("adc",&ht_adc); // best estimate of hit charge
  //}
}
void calib::CalibAna::ResetHitBranches()
{
  ht_istricell=false; ht_isz=false;
  ht_istraj=false; ht_isavg=false;
  ht_diblock=-5; ht_dcm=-5; ht_apd=-5;
  ht_pixel=-5; ht_plane=-5; ht_cell=-5;
  ht_w=-999; ht_pe=-5; ht_pecorr=-5;
  ht_xapprox=-999; ht_yapprox=-999; ht_zapprox=-999;
  ht_pecorrtomev=5; ht_recoMeV=-5; ht_path=-5;
  ht_cmFromEnd=-5; ht_flightLen=-5; 
  ht_view=-5; ht_tns=-5;
  ht_truepath=-5; ht_truepe=-5; ht_truew=-5; 
  ht_trueMeV=-5; ht_brightbin=-5; 
  ht_poissonlambda=-5; ht_readdist=-5;
  ht_tdc=-5; ht_adc=-5;
}


DEFINE_ART_MODULE(calib::CalibAna)