BeamFlashTrackMatchTaggerAlg.cxx

Go to the documentation of this file.

#include "BeamFlashTrackMatchTaggerAlg.h"
#include "larcorealg/Geometry/OpDetGeo.h"
#include <limits>

cosmic::BeamFlashTrackMatchTaggerAlg::BeamFlashTrackMatchTaggerAlg(fhicl::ParameterSet const& p) 
  : COSMIC_TYPE_FLASHMATCH(anab::CosmicTagID_t::kFlash_BeamIncompatible),
    COSMIC_TYPE_OUTSIDEDRIFT(anab::CosmicTagID_t::kOutsideDrift_Partial),
    DEBUG_FLAG(p.get<bool>("RunDebugMode",false))
{
  this->reconfigure(p);
}

void cosmic::BeamFlashTrackMatchTaggerAlg::reconfigure(fhicl::ParameterSet const& p){
  fMinTrackLength = p.get<float>("MinTrackLength");
  fMinOpHitPE     = p.get<float>("MinOpHitPE",0.1);
  fMIPdQdx        = p.get<float>("MIPdQdx",2.1);
  fOpDetSaturation = p.get<float>("OpDetSaturation",200.);

  fSingleChannelCut           = p.get<float>("SingleChannelCut");
  fCumulativeChannelThreshold = p.get<float>("CumulativeChannelThreshold");
  fCumulativeChannelCut       = p.get<unsigned int>("CumulativeChannelCut");
  fIntegralCut                = p.get<float>("IntegralCut");
  
  fMakeOutsideDriftTags       = p.get<bool>("MakeOutsideDriftTags",false);
  fNormalizeHypothesisToFlash = p.get<bool>("NormalizeHypothesisToFlash");
}

void cosmic::BeamFlashTrackMatchTaggerAlg::SetHypothesisComparisonTree(TTree* tree, 
                                                                       TH1F* hist_flash, TH1F* hist_hyp){
  cTree = tree;

  cOpDetHist_flash = hist_flash;
  cOpDetHist_flash->SetNameTitle("opdet_hist_flash","Optical Detector Occupancy, Flash");

  cOpDetHist_hyp = hist_hyp;
  cOpDetHist_hyp->SetNameTitle("opdet_hist_hyp","Optical Detector Occupancy, Hypothesis");

  cTree->Branch("fcp",&cFlashComparison_p,cFlashComparison_p.leaf_structure.c_str());
  cTree->Branch("opdet_hyp",&cOpDetVector_hyp);
  cTree->Branch("opdet_flash",&cOpDetVector_flash);
  cTree->Branch("opdet_hist_flash","TH1F",cOpDetHist_flash);
  cTree->Branch("opdet_hist_hyp","TH1F",cOpDetHist_hyp);
}

void cosmic::BeamFlashTrackMatchTaggerAlg::RunCompatibilityCheck(std::vector<recob::OpFlash> const& flashVector,
                                                                 std::vector<recob::Track> const& trackVector,
                                                                 std::vector<anab::CosmicTag>& cosmicTagVector,
                                                                 std::vector<size_t>& assnTrackTagVector,
                                                                 Providers_t providers,
                                                                 phot::PhotonVisibilityService const& pvs,
                                                                 opdet::OpDigiProperties const& opdigip){

  auto const& geom = *(providers.get<geo::GeometryCore>());
  
  std::vector< const recob::OpFlash* > flashesOnBeamTime;
  for(auto const& flash : flashVector){
    if(!flash.OnBeamTime()) continue;
    flashesOnBeamTime.push_back(&flash);
  }

  //make sure this association vector is initialized properly
  assnTrackTagVector.resize(trackVector.size(),std::numeric_limits<size_t>::max());
  cosmicTagVector.reserve(trackVector.size());

  for(size_t track_i=0; track_i<trackVector.size(); track_i++){

    recob::Track const& track(trackVector[track_i]);

    if(track.Length() < fMinTrackLength) continue;

    //get the begin and end points of this track
    TVector3 const& pt_begin = track.LocationAtPoint(0);
    TVector3 const& pt_end = track.LocationAtPoint(track.NumberTrajectoryPoints()-1);
    std::vector<float> xyz_begin = { (float)pt_begin.x(), (float)pt_begin.y(), (float)pt_begin.z()};
    std::vector<float> xyz_end = {(float)pt_end.x(), (float)pt_end.y(), (float)pt_end.z()};

    //check if this track is outside the drift window, and if it is continue
    if(!InDriftWindow(pt_begin.x(),pt_end.x(),geom)) {
      if(fMakeOutsideDriftTags){
        cosmicTagVector.emplace_back(xyz_begin,xyz_end,1.,COSMIC_TYPE_OUTSIDEDRIFT);
        assnTrackTagVector[track_i] = cosmicTagVector.size()-1;
      }
      continue;
    }

    //get light hypothesis for track
    std::vector<float> lightHypothesis = GetMIPHypotheses(track,providers,pvs,opdigip);

    //check compatibility with beam flash
    bool compatible=false;
    for(const recob::OpFlash* flashPointer : flashesOnBeamTime){
      CompatibilityResultType result = CheckCompatibility(lightHypothesis,flashPointer, geom);
      if(result==CompatibilityResultType::kCompatible) compatible=true;
      if(DEBUG_FLAG){
        PrintTrackProperties(track);
        PrintFlashProperties(*flashPointer);
        PrintHypothesisFlashComparison(lightHypothesis,flashPointer,geom,result);
      }
    }

    //make tag
    float cosmicScore=1.;
    if(compatible) cosmicScore=0.;
    cosmicTagVector.emplace_back(xyz_begin,xyz_end,cosmicScore,COSMIC_TYPE_FLASHMATCH);
    assnTrackTagVector[track_i]=cosmicTagVector.size()-1;
  }

}

//this compares the hypothesis to the flash itself.
void cosmic::BeamFlashTrackMatchTaggerAlg::RunHypothesisComparison(unsigned int const run,
                                                                   unsigned int const event,
                                                                   std::vector<recob::OpFlash> const& flashVector,
                                                                   std::vector<recob::Track> const& trackVector,
                                                                   Providers_t providers,
                                                                   phot::PhotonVisibilityService const& pvs,
                                                                   opdet::OpDigiProperties const& opdigip){

  auto const& geom = *(providers.get<geo::GeometryCore>());
  
  cFlashComparison_p.run = run;
  cFlashComparison_p.event = event;

  std::vector< std::pair<unsigned int, const recob::OpFlash*> > flashesOnBeamTime;
  for(unsigned int i=0; i<flashVector.size(); i++){
    recob::OpFlash const& flash = flashVector[i];
    if(!flash.OnBeamTime()) continue;
    flashesOnBeamTime.push_back(std::make_pair(i,&flash));
  }
  
  for(size_t track_i=0; track_i<trackVector.size(); track_i++){

    recob::Track const& track(trackVector[track_i]);

    if(track.Length() < fMinTrackLength) continue;

    //get the begin and end points of this track
    TVector3 const& pt_begin = track.LocationAtPoint(0);
    TVector3 const& pt_end = track.LocationAtPoint(track.NumberTrajectoryPoints()-1);
    std::vector<float> xyz_begin = { (float)pt_begin.x(), (float)pt_begin.y(), (float)pt_begin.z()};
    std::vector<float> xyz_end = {(float)pt_end.x(), (float)pt_end.y(), (float)pt_end.z()};

    //check if this track is outside the drift window, and if it is continue
    if(!InDriftWindow(pt_begin.x(),pt_end.x(),geom)) continue; 

    cFlashComparison_p.trk_startx = pt_begin.x();
    cFlashComparison_p.trk_starty = pt_begin.y();
    cFlashComparison_p.trk_startz = pt_begin.z();
    cFlashComparison_p.trk_endx = pt_end.x();
    cFlashComparison_p.trk_endy = pt_end.y();
    cFlashComparison_p.trk_endz = pt_end.z();

    //get light hypothesis for track
    cOpDetVector_hyp = GetMIPHypotheses(track,providers,pvs,opdigip);

    cFlashComparison_p.hyp_index = track_i;
    FillFlashProperties(cOpDetVector_hyp,
                        cFlashComparison_p.hyp_totalPE,
                        cFlashComparison_p.hyp_y,cFlashComparison_p.hyp_sigmay,
                        cFlashComparison_p.hyp_z,cFlashComparison_p.hyp_sigmaz,
                        geom);
    
    for(auto flash : flashesOnBeamTime){
      cOpDetVector_flash = std::vector<float>(geom.NOpDets(),0);  
      cFlashComparison_p.flash_nOpDet = 0;
      for(size_t c=0; c<=geom.MaxOpChannel(); c++){
        if ( geom.IsValidOpChannel( c ) ) {
          unsigned int OpDet = geom.OpDetFromOpChannel(c);
          cOpDetVector_flash[OpDet] += flash.second->PE(c);
        }
      }
      for(size_t o=0; o<cOpDetVector_flash.size(); o++)
        if(cOpDetVector_flash[o] < fMinOpHitPE) cFlashComparison_p.flash_nOpDet++;

      cFlashComparison_p.flash_index = flash.first;
      cFlashComparison_p.flash_totalPE = flash.second->TotalPE();
      cFlashComparison_p.flash_y = flash.second->YCenter();
      cFlashComparison_p.flash_sigmay = flash.second->YWidth();
      cFlashComparison_p.flash_z = flash.second->ZCenter();
      cFlashComparison_p.flash_sigmaz = flash.second->ZWidth();

      cFlashComparison_p.chi2 = CalculateChi2(cOpDetVector_flash,cOpDetVector_hyp);

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


}

//this compares the hypothesis to the flash itself.
void cosmic::BeamFlashTrackMatchTaggerAlg::RunHypothesisComparison(unsigned int const run,
                                                                   unsigned int const event,
                                                                   std::vector<recob::OpFlash> const& flashVector,
                                                                   std::vector<simb::MCParticle> const& mcParticleVector,
                                                                   Providers_t providers,
                                                                   phot::PhotonVisibilityService const& pvs,
                                                                   opdet::OpDigiProperties const& opdigip){

  auto const& geom = *(providers.get<geo::GeometryCore>());
  
  cFlashComparison_p.run = run;
  cFlashComparison_p.event = event;

  std::vector< std::pair<unsigned int, const recob::OpFlash*> > flashesOnBeamTime;
  for(unsigned int i=0; i<flashVector.size(); i++){
    recob::OpFlash const& flash = flashVector[i];
    if(!flash.OnBeamTime()) continue;
    flashesOnBeamTime.push_back(std::make_pair(i,&flash));
  }
  
  for(size_t particle_i=0; particle_i<mcParticleVector.size(); particle_i++){

    simb::MCParticle const& particle(mcParticleVector[particle_i]);
    if(particle.Process().compare("primary")!=0) continue;
    if(particle.Trajectory().TotalLength() < fMinTrackLength) continue;

    //get the begin and end points of this track
    size_t start_i=0, end_i=particle.NumberTrajectoryPoints()-1;
    bool prev_inside=false;
    for(size_t pt_i=0; pt_i < particle.NumberTrajectoryPoints(); pt_i++){
      bool inside = InDetector(particle.Position(pt_i).Vect(),geom);
      if(inside && !prev_inside) start_i = pt_i;
      if(!inside && prev_inside) { end_i = pt_i-1; break; }
      prev_inside = inside;
    }
    TVector3 const& pt_begin = particle.Position(start_i).Vect();
    TVector3 const& pt_end = particle.Position(end_i).Vect();
    std::vector<float> xyz_begin = { (float)pt_begin.x(), (float)pt_begin.y(), (float)pt_begin.z()};
    std::vector<float> xyz_end = {(float)pt_end.x(), (float)pt_end.y(), (float)pt_end.z()};

    //check if this track is outside the drift window, and if it is continue
    if(!InDriftWindow(pt_begin.x(),pt_end.x(),geom)) continue; 

    cFlashComparison_p.trk_startx = pt_begin.x();
    cFlashComparison_p.trk_starty = pt_begin.y();
    cFlashComparison_p.trk_startz = pt_begin.z();
    cFlashComparison_p.trk_endx = pt_end.x();
    cFlashComparison_p.trk_endy = pt_end.y();
    cFlashComparison_p.trk_endz = pt_end.z();

    //get light hypothesis for track
    cOpDetVector_hyp = GetMIPHypotheses(particle,start_i,end_i,providers,pvs,opdigip);

    cFlashComparison_p.hyp_index = particle_i;
    FillFlashProperties(cOpDetVector_hyp,
                        cFlashComparison_p.hyp_totalPE,
                        cFlashComparison_p.hyp_y,cFlashComparison_p.hyp_sigmay,
                        cFlashComparison_p.hyp_z,cFlashComparison_p.hyp_sigmaz,
                        geom);
    
    for(auto flash : flashesOnBeamTime){
      cOpDetVector_flash = std::vector<float>(geom.NOpDets(),0);  
      cFlashComparison_p.flash_nOpDet = 0;
      //for(size_t c=0; c<geom.NOpChannels(); c++){
      for(size_t c=0; c<=geom.MaxOpChannel(); c++){
        if ( geom.IsValidOpChannel(c) ) {
          unsigned int OpDet = geom.OpDetFromOpChannel(c);
          cOpDetVector_flash[OpDet] += flash.second->PE(c);
        }
      }
      for(size_t o=0; o<cOpDetVector_flash.size(); o++)
        if(cOpDetVector_flash[o] < fMinOpHitPE) cFlashComparison_p.flash_nOpDet++;
      cFlashComparison_p.flash_index = flash.first;
      cFlashComparison_p.flash_totalPE = flash.second->TotalPE();
      cFlashComparison_p.flash_y = flash.second->YCenter();
      cFlashComparison_p.flash_sigmay = flash.second->YWidth();
      cFlashComparison_p.flash_z = flash.second->ZCenter();
      cFlashComparison_p.flash_sigmaz = flash.second->ZWidth();

      cFlashComparison_p.chi2 = CalculateChi2(cOpDetVector_flash,cOpDetVector_hyp);

      //cOpDetHist_flash->SetBins(cOpDetVector_flash.size(),0,cOpDetVector_flash.size());
      for(size_t i=0; i<cOpDetVector_flash.size(); i++)
        cOpDetHist_flash->SetBinContent(i+1,cOpDetVector_flash[i]);

      //cOpDetHist_hyp->SetBins(cOpDetVector_hyp.size(),0,cOpDetVector_hyp.size());
      for(size_t i=0; i<cOpDetVector_hyp.size(); i++)
        cOpDetHist_hyp->SetBinContent(i+1,cOpDetVector_hyp[i]);

      for(size_t i=0; i<cOpDetVector_flash.size(); i++){
        std::cout << "Flash/Hyp " << i << " : " 
                  << cOpDetHist_flash->GetBinContent(i+1) << " "
                  << cOpDetHist_hyp->GetBinContent(i+1) << std::endl;
      }

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


}

void cosmic::BeamFlashTrackMatchTaggerAlg::FillFlashProperties(std::vector<float> const& opdetVector,
                                                               float& sum,
                                                               float& y, float& sigmay,
                                                               float& z, float& sigmaz,
                                                               geo::GeometryCore const& geom){
  y=0; sigmay=0; z=0; sigmaz=0; sum=0;
  double xyz[3];
  for(unsigned int opdet=0; opdet<opdetVector.size(); opdet++){
    sum+=opdetVector[opdet];
    geom.Cryostat(0).OpDet(opdet).GetCenter(xyz);
    y += opdetVector[opdet]*xyz[1];
    z += opdetVector[opdet]*xyz[2];
  }

  y /= sum; z /= sum;

  for(unsigned int opdet=0; opdet<opdetVector.size(); opdet++){
    geom.Cryostat(0).OpDet(opdet).GetCenter(xyz);
    sigmay += (opdetVector[opdet]*xyz[1]-y)*(opdetVector[opdet]*xyz[1]-y);
    sigmaz += (opdetVector[opdet]*xyz[2]-y)*(opdetVector[opdet]*xyz[2]-y);
  }

  sigmay = std::sqrt(sigmay)/sum; 
  sigmaz = std::sqrt(sigmaz)/sum; 

}

bool cosmic::BeamFlashTrackMatchTaggerAlg::InDetector(TVector3 const& pt, geo::GeometryCore const& geom){
  if(pt.x() < 0 || pt.x() > 2*geom.DetHalfWidth()) return false;
  if(std::abs(pt.y()) > geom.DetHalfHeight()) return false;
  if(pt.z() < 0 || pt.z() > geom.DetLength()) return false;
  return true;
}

bool cosmic::BeamFlashTrackMatchTaggerAlg::InDriftWindow(double start_x, double end_x, geo::GeometryCore const& geom){
  if(start_x < 0. || end_x < 0.) return false;
  if(start_x > 2*geom.DetHalfWidth() || end_x > 2*geom.DetHalfWidth()) return false;
  return true;
}

void cosmic::BeamFlashTrackMatchTaggerAlg::AddLightFromSegment(TVector3 const& pt1,
                                                               TVector3 const& pt2,
                                                               std::vector<float> & lightHypothesis,
                                                               float & totalHypothesisPE,
                                                               geo::GeometryCore const& geom,
                                                               phot::PhotonVisibilityService const& pvs,
                                                               float const& PromptMIPScintYield,
                                                               float XOffset){

  double xyz_segment[3];
  xyz_segment[0] = 0.5*(pt2.x()+pt1.x()) + XOffset;
  xyz_segment[1] = 0.5*(pt2.y()+pt1.y());
  xyz_segment[2] = 0.5*(pt2.z()+pt1.z());
    
  //get the visibility vector
  const float* PointVisibility = pvs.GetAllVisibilities(xyz_segment);
  
  //check pointer, as it may be null if given a y/z outside some range
  if(!PointVisibility) return;
  
  //get the amount of light
  float LightAmount = PromptMIPScintYield*(pt2-pt1).Mag();
  
  for(size_t opdet_i=0; opdet_i<pvs.NOpChannels(); opdet_i++){
    lightHypothesis[opdet_i] += PointVisibility[opdet_i]*LightAmount;
    totalHypothesisPE += PointVisibility[opdet_i]*LightAmount;
   
    //apply saturation limit
    if(lightHypothesis[opdet_i]>fOpDetSaturation){
      totalHypothesisPE -= (lightHypothesis[opdet_i]-fOpDetSaturation);
      lightHypothesis[opdet_i] = fOpDetSaturation;
    }
  }
  
}//end AddLightFromSegment

void cosmic::BeamFlashTrackMatchTaggerAlg::NormalizeLightHypothesis(std::vector<float> & lightHypothesis,
                                                                    float const& totalHypothesisPE,
                                                                    geo::GeometryCore const& geom){
  for(size_t opdet_i=0; opdet_i<geom.NOpDets(); opdet_i++)
    lightHypothesis[opdet_i] /= totalHypothesisPE;
}
                                                                    

// Get a hypothesis for the light collected for a track
std::vector<float> cosmic::BeamFlashTrackMatchTaggerAlg::GetMIPHypotheses(recob::Track const& track, 
                                                                          Providers_t providers,
                                                                          phot::PhotonVisibilityService const& pvs,
                                                                          opdet::OpDigiProperties const& opdigip,
                                                                          float XOffset)
{
  auto const& geom = *(providers.get<geo::GeometryCore>());
  auto const& larp = *(providers.get<detinfo::LArProperties>());
  std::vector<float> lightHypothesis(geom.NOpDets(),0);  
  float totalHypothesisPE=0;
  const float PromptMIPScintYield = larp.ScintYield()*larp.ScintYieldRatio()*opdigip.QE()*fMIPdQdx;

  //get QE from ubChannelConfig, which gives per tube, so goes in AddLightFromSegment
  //VisibleEnergySeparation(step);
  
  for(size_t pt=1; pt<track.NumberTrajectoryPoints(); pt++)    
    AddLightFromSegment(track.LocationAtPoint(pt-1),track.LocationAtPoint(pt),
                        lightHypothesis,totalHypothesisPE,
                        geom,pvs,PromptMIPScintYield,
                        XOffset);

  if(fNormalizeHypothesisToFlash && totalHypothesisPE > std::numeric_limits<float>::epsilon())
    NormalizeLightHypothesis(lightHypothesis,totalHypothesisPE,geom);

  return lightHypothesis;

}//end GetMIPHypotheses


// Get a hypothesis for the light collected for a particle trajectory
std::vector<float> cosmic::BeamFlashTrackMatchTaggerAlg::GetMIPHypotheses(simb::MCParticle const& particle, 
                                                                          size_t start_i, size_t end_i,
                                                                          Providers_t providers,
                                                                          phot::PhotonVisibilityService const& pvs,
                                                                          opdet::OpDigiProperties const& opdigip,
                                                                          float XOffset)
{
  auto const& geom = *(providers.get<geo::GeometryCore>());
  auto const& larp = *(providers.get<detinfo::LArProperties>());
  std::vector<float> lightHypothesis(geom.NOpDets(),0);  
  float totalHypothesisPE=0;
  const float PromptMIPScintYield = larp.ScintYield()*larp.ScintYieldRatio()*opdigip.QE()*fMIPdQdx;

  for(size_t pt=start_i+1; pt<=end_i; pt++)
    AddLightFromSegment(particle.Position(pt-1).Vect(),particle.Position(pt).Vect(),
                        lightHypothesis,totalHypothesisPE,
                        geom,pvs,PromptMIPScintYield,
                        XOffset);

  if(fNormalizeHypothesisToFlash && totalHypothesisPE > std::numeric_limits<float>::epsilon())
    NormalizeLightHypothesis(lightHypothesis,totalHypothesisPE,geom);

  return lightHypothesis;

}//end GetMIPHypotheses


//---------------------------------------
//  Check whether a hypothesis can be accomodated in a flash
//   Flashes fail if 1 bin is far in excess of the observed signal 
//   or if the whole flash intensity is much too large for the hypothesis.
//  MIP dEdx is assumed for now.  Accounting for real dQdx will 
//   improve performance of this algorithm.
//---------------------------------------
cosmic::BeamFlashTrackMatchTaggerAlg::CompatibilityResultType 
cosmic::BeamFlashTrackMatchTaggerAlg::CheckCompatibility(std::vector<float> const& lightHypothesis, 
                                                         const recob::OpFlash* flashPointer,
                                                         geo::GeometryCore const& geom)
{
  float hypothesis_integral=0;
  float flash_integral=0;
  unsigned int cumulativeChannels=0;

  std::vector<double> PEbyOpDet(geom.NOpDets(),0);
  //for (unsigned int c = 0; c < geom.NOpChannels(); c++){
  for (unsigned int c = 0; c <= geom.MaxOpChannel(); c++){
    if ( geom.IsValidOpChannel(c) ) {
      unsigned int o = geom.OpDetFromOpChannel(c);
      PEbyOpDet[o] += flashPointer->PE(c);
    }
  }
  
  float hypothesis_scale=1.;
  if(fNormalizeHypothesisToFlash) hypothesis_scale = flashPointer->TotalPE();

  for(size_t pmt_i=0; pmt_i<lightHypothesis.size(); pmt_i++){

    flash_integral += PEbyOpDet[pmt_i];

    if(lightHypothesis[pmt_i] < std::numeric_limits<float>::epsilon() ) continue;
    hypothesis_integral += lightHypothesis[pmt_i]*hypothesis_scale;

    if(PEbyOpDet[pmt_i] < fMinOpHitPE) continue;

    float diff_scaled = (lightHypothesis[pmt_i]*hypothesis_scale - PEbyOpDet[pmt_i])/std::sqrt(lightHypothesis[pmt_i]*hypothesis_scale);

    if( diff_scaled > fSingleChannelCut ) return CompatibilityResultType::kSingleChannelCut;

    if( diff_scaled > fCumulativeChannelThreshold ) cumulativeChannels++;
    if(cumulativeChannels >= fCumulativeChannelCut) return CompatibilityResultType::kCumulativeChannelCut;

  }

  if( (hypothesis_integral - flash_integral)/std::sqrt(hypothesis_integral) 
      > fIntegralCut) return CompatibilityResultType::kIntegralCut;

  return CompatibilityResultType::kCompatible;
}


float cosmic::BeamFlashTrackMatchTaggerAlg::CalculateChi2(std::vector<float> const& light_flash,
                                                          std::vector<float> const& light_track){

  float chi2=0;
  for(size_t pmt_i=0; pmt_i<light_flash.size(); pmt_i++){

    if(light_flash[pmt_i] < fMinOpHitPE) continue;

    float err2 = 1;
    if(light_track[pmt_i] > 1) err2 = light_track[pmt_i];

    chi2 += (light_flash[pmt_i]-light_track[pmt_i])*(light_flash[pmt_i]-light_track[pmt_i]) / err2;
  }

  return chi2;
}
                                                           

void cosmic::BeamFlashTrackMatchTaggerAlg::PrintTrackProperties(recob::Track const& track, std::ostream* output)
{
  *output << "----------------------------------------------" << std::endl;
  *output << "Track properties: ";
  *output << "\n\tLength=" << track.Length();

  TVector3 const& pt_begin = track.LocationAtPoint(0);
  *output << "\n\tBegin Location (x,y,z)=(" << pt_begin.x() << "," << pt_begin.y() << "," << pt_begin.z() << ")";

  TVector3 const& pt_end = track.LocationAtPoint(track.NumberTrajectoryPoints()-1);
  *output << "\n\tEnd Location (x,y,z)=(" << pt_end.x() << "," << pt_end.y() << "," << pt_end.z() << ")";

  *output << "\n\tTrajectoryPoints=" << track.NumberTrajectoryPoints();
  *output << std::endl;
  *output << "----------------------------------------------" << std::endl;
}

void cosmic::BeamFlashTrackMatchTaggerAlg::PrintFlashProperties(recob::OpFlash const& flash, std::ostream* output)
{
  *output << "----------------------------------------------" << std::endl;
  *output << "Flash properties: ";

  *output << "\n\tTime=" << flash.Time();
  *output << "\n\tOnBeamTime=" << flash.OnBeamTime();
  *output << "\n\ty position (center,width)=(" << flash.YCenter() << "," << flash.YWidth() << ")";
  *output << "\n\tz position (center,width)=(" << flash.ZCenter() << "," << flash.ZWidth() << ")";
  *output << "\n\tTotal PE=" << flash.TotalPE();

  *output << std::endl;
  *output << "----------------------------------------------" << std::endl;

}

void cosmic::BeamFlashTrackMatchTaggerAlg::PrintHypothesisFlashComparison(std::vector<float> const& lightHypothesis,
                                                                          const recob::OpFlash* flashPointer,
                                                                          geo::GeometryCore const& geom,
                                                                          CompatibilityResultType result,
                                                                          std::ostream* output)
{

  *output << "----------------------------------------------" << std::endl;
  *output << "Hypothesis-flash comparison: ";

  float hypothesis_integral=0;
  float flash_integral=0;

  float hypothesis_scale=1.;
  if(fNormalizeHypothesisToFlash) hypothesis_scale = flashPointer->TotalPE();

  
  std::vector<double> PEbyOpDet(geom.NOpDets(),0);
  //for (unsigned int c = 0; c < geom.NOpChannels(); c++){
  for ( unsigned int c = 0; c <= geom.MaxOpChannel(); c++){
    if ( geom.IsValidOpChannel(c) ) {
      unsigned int o = geom.OpDetFromOpChannel(c);
      PEbyOpDet[o] += flashPointer->PE(c);
    }
  }

  for(size_t pmt_i=0; pmt_i<lightHypothesis.size(); pmt_i++){

    flash_integral += PEbyOpDet[pmt_i];

    *output << "\n\t pmt_i=" << pmt_i << ", (hypothesis,flash)=(" 
           << lightHypothesis[pmt_i]*hypothesis_scale << "," << PEbyOpDet[pmt_i] << ")";

    if(lightHypothesis[pmt_i] < std::numeric_limits<float>::epsilon() ) continue;

    *output << "  difference=" 
            << (lightHypothesis[pmt_i]*hypothesis_scale - PEbyOpDet[pmt_i])/std::sqrt(lightHypothesis[pmt_i]*hypothesis_scale);

    hypothesis_integral += lightHypothesis[pmt_i]*hypothesis_scale;
  }

  *output << "\n\t TOTAL (hypothesis,flash)=(" 
         << hypothesis_integral << "," << flash_integral << ")"
         << "  difference=" << (hypothesis_integral - flash_integral)/std::sqrt(hypothesis_integral);

  *output << std::endl;
  *output << "End result=" << result << std::endl;
  *output << "----------------------------------------------" << std::endl;

}