FLSHit.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   FLSHit class
/// \author  brebel@fnal.gov
/// \date
////////////////////////////////////////////////////////////////////////
#include "Simulation/FLSHit.h"
#include <cmath>

namespace sim{
  
  //---------------------------------------------------------------
  FLSHit::FLSHit():
    fId     (0)
    ,fPlaneId(0)
    ,fCellId (0)
    ,fPDG    (0)
    ,fTrackId(0)
  { }
  //__________________________________________________________
  FLSHit::~FLSHit() { }
  
  //__________________________________________________________
  bool FLSHit::operator==(const FLSHit& rhs) const
  {
    if (fId == rhs.fId && GetEntryT() == rhs.GetEntryT()) return true;
    return false;
  }
  
  //__________________________________________________________
  bool FLSHit::operator<(const FLSHit& rhs) const
  {
    if (fId<rhs.fId) return true;
    if (fId == rhs.fId)
    if (GetEntryT() < rhs.GetEntryT()) return true;
    return false;
  }

  //__________________________________________________________
  void FLSHit::Clear() {

    // Clear almost everything. We don't want to clear fId, though

    fPlaneId = 0;
    fCellId  = 0;
    fTrackId = 0;
    fPDG     = 0;

    fX.clear();
    fY.clear();
    fZ.clear();
    fT.clear();
    fEnergy.clear();
    fVectorEdep.clear();
    fVectorEdepBirks.clear();
    fVectorNCerenkov.clear();
  }

//----------------------------------------------------------------------------
// Get path length of a particular step
float FLSHit::GetPathLength(const int step) const{

    const double x = fX[step+1] - fX[step];
    const double y = fY[step+1] - fY[step];
    const double z = fZ[step+1] - fZ[step];

    return sqrt(x*x + y*y + z*z);
}
//----------------------------------------------------------------------------
// Get path length of all steps in FLSHit
float FLSHit::GetTotalPathLength() const{

    const int nsteps = GetNSteps();

    double path_length = 0.0;
    for(int i=0; i<nsteps; ++i){
        path_length += GetPathLength(i);
    }

    return path_length;
}

//----------------------------------------------------------------------------
// Get energy loss of a particular step
float FLSHit:: GetEnergyLoss(const int step) const{
    return fEnergy[step] - fEnergy[step+1];
}

//----------------------------------------------------------------------------
// Get total energy loss in the entire FLSHit
float FLSHit::GetTotalEnergyLoss() const{
    return GetEntryEnergy() - GetExitEnergy();
}

//----------------------------------------------------------------------------
float FLSHit::GetXYZAverage(const std::vector<float>& xyz, const int step) const{
    return (xyz[step] + xyz[step+1]) * 0.5;
}

//----------------------------------------------------------------------------
float FLSHit::GetXYZAverage(const std::vector<float>& xyz) const{
    const int nsteps = GetNSteps();

    double xyz_av = 0.0;
    double weight = 0.0;
    
    for(int i=0; i<nsteps; ++i){
        const double path_length = GetPathLength(i);

        weight += path_length;
        xyz_av += path_length * GetXYZAverage(xyz, i);
    }

    return xyz_av / weight;
}

  //----------------------------------------------------------------------------
  // Performs various checks on the size of the data and its content
  // Has an option to print everything out
  bool FLSHit::IsFLSHitReasonable(const bool    print_errors
                                 ,std::ostream& os) const
  {
      bool is_everything_ok = true;

      /// All sizes should be at least 2 (Entry+Exit)
      if(fX.size() < 2 )               {is_everything_ok = false; if(print_errors){os<<"Vector fX is empty\n";}}
      if(fY.size() < 2 )               {is_everything_ok = false; if(print_errors){os<<"Vector fY is empty\n";}}
      if(fZ.size() < 2 )               {is_everything_ok = false; if(print_errors){os<<"Vector fZ is empty\n";}}
      if(fT.size() < 2 )               {is_everything_ok = false; if(print_errors){os<<"Vector fT is empty\n";}}
      if(fEnergy.size() < 2 )          {is_everything_ok = false; if(print_errors){os<<"Vector fEnergy is empty\n";}}
      if(fVectorEdep.size() < 2 )      {is_everything_ok = false; if(print_errors){os<<"Vector fVectorEdep is empty\n";}}
      if(fVectorEdepBirks.size() < 2 ) {is_everything_ok = false; if(print_errors){os<<"Vector fVectorEdepBirks is empty\n";}}

      /// Compare the sizes. Should be identical
      if(fX.size() != fY.size() ||
         fX.size() != fZ.size() ||
         fX.size() != fT.size() ||
         fX.size() != fEnergy.size()
         ){ 
          is_everything_ok = false;
          if(print_errors){
              os<<"Vectors fX, fY, fZ, fT, fEnergy should have the same size, but they have sizes";
              os<<" "<<fX.size();
              os<<" "<<fY.size();
              os<<" "<<fZ.size();
              os<<" "<<fT.size();
              os<<" "<<fEnergy.size();
              os<<" respectively.\n";
         }// end of printing
      }// end of checking

      /// Compare the sizes. Should be identical
      if(fVectorEdep.size() != fVectorEdepBirks.size()){
          is_everything_ok = false;
          if(print_errors){
              os<<"Vectors fVectorEdep and fVectorEdepBirks, but they have sizes";
              os<<" "<<fVectorEdep.size();
              os<<" "<<fVectorEdepBirks.size();
              os<<" respectively.\n";
          }// end of printing
      }// end of checking
      
      /// Compare the sizes. fX should have one more entry than fVectorEdep
      if(fX.size() != fVectorEdep.size() + 1){
          is_everything_ok = false;
          if(print_errors){
              os<<"Vector fX should have one more entry than fVectorEdep, but they have sizes";
              os<<" "<<fX.size();
              os<<" "<<fVectorEdep.size();
              os<<" respectively.\n";
          }// end of printing
      }// end of checking

      /// Time should increase monotonically
      if(!DoesIncreaseMonotonically(fT)){
          is_everything_ok = false;
          if(print_errors){
              os<<"Vector fT did not increase monotonically\n";
          }// end of printing
      }// end of checking

      /// Energy should decrease monotonically
      if(!DoesDecreaseMonotonically(fEnergy)){
          is_everything_ok = false;
          if(print_errors){
              os<<"Vector fEnergy did not decrease monotonically\n";
          }// end of printing
      }// end of checking

      /// ID of the FLSHit should be non-zero
      if(fId == 0){
          is_everything_ok = false;
          if(print_errors){
              os<<"fID is zero\n";
          }// end of printing
      }// end of checking fID

      return is_everything_ok;
  }// end of IsFLSHitReasonable

  //----------------------------------------------------------------------
  // ostream operator.  
  //
  std::ostream& operator<< (std::ostream & o, const FLSHit & a)
  {
    o << std::setiosflags(std::ios::fixed) << std::setprecision(2);
    o << " CellUniqueId "       << std::setw(10)  << std::right << a.GetCellUniqueId()
      << " Entry (x,y,z): ("    << a.GetEntryX()  << "," << a.GetEntryY() << "," << a.GetEntryZ() << ")"
      << " Entry E: "           << std::setw(5)   << std::right << a.GetEntryEnergy()
      << " Exit (x,y,z): ("     << a.GetExitX()   << "," << a.GetExitY() << "," << a.GetExitZ() << ")"
      << " Exit E: "            << std::setw(5)   << std::right << a.GetExitEnergy()
      << " Entry Time: "        << std::setw(5)   << std::right << a.GetEntryT()
      << " Exit Time: "         << std::setw(5)   << std::right << a.GetExitT()
      << " E deposited: "       << std::setw(5)   << std::right << a.GetEdep()
      << " E Birks Weight: "    << std::setw(5)   << std::right << a.GetEdepBirks()
      << " N Cerenkov Photons " << std::setw(5)   << std::right << a.GetNCerenkov()
      << " Plane,Cell: "        << a.GetPlaneID() << "," << a.GetCellID()
      << " Track ID: "          << std::setw(5)   << std::right << a.GetTrackID();
    
    return o;
  }

}// end  namespace sim

////////////////////////////////////////////////////////////////////////