Calorimetry.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
//
// \brief Definition of Calorimetry analysis object
//
// \author brebel@fnal.gov, tjyang@fnal.gov
////////////////////////////////////////////////////////////////////////

#include "lardataobj/AnalysisBase/Calorimetry.h"
#include "cetlib_except/exception.h"

#include <iostream>
#include <iomanip>

namespace anab{

  //----------------------------------------------------------------------
  Calorimetry::Calorimetry()
    : fKineticEnergy(0.)
    , fRange(0.)
  {
    fdEdx.clear();
    fdQdx.clear();
    fResidualRange.clear();
    fDeadWireResR.clear();
    fTrkPitch.clear();
    fXYZ.clear();
    fTpIndices.clear();
  }


 //----------------------------------------------------------------------
  Calorimetry::Calorimetry(float KineticEnergy,
                           std::vector<float> const& dEdx,
                           std::vector<float> const& dQdx,
                           std::vector<float> const& resRange,
                           std::vector<float> const& deadwire,
                           float Range,
                           float TrkPitch,
                           geo::PlaneID planeID)
  {

    fKineticEnergy = KineticEnergy;
    fRange = Range;
    for(size_t i=0; i!=dQdx.size(); ++i){
      fTrkPitch.push_back(TrkPitch);
      fXYZ.push_back({-999.,-999.,-999.});
    }
    if(dEdx.size() != resRange.size())
      throw cet::exception("anab::Calorimetry") << "dE/dx and residual range vectors "
                                                << "have different sizes, this is a problem.\n";
    fdEdx.resize(dEdx.size());
    fdQdx.resize(dQdx.size());
    fResidualRange.resize(resRange.size());
    for(size_t i = 0; i < dEdx.size(); ++i){
      fdEdx[i]          = dEdx[i];
      fdQdx[i]          = dQdx[i];
      fResidualRange[i] = resRange[i];
    }

    fDeadWireResR.resize(deadwire.size());
    for(size_t i = 0; i<deadwire.size(); ++i){
      fDeadWireResR[i] = deadwire[i];
    }

    fPlaneID = planeID;
  }


  //----------------------------------------------------------------------
  Calorimetry::Calorimetry(float KineticEnergy,
                           std::vector<float> const& dEdx,
                           std::vector<float> const& dQdx,
                           std::vector<float> const& resRange,
                           std::vector<float> const& deadwire,
                           float Range,
                           std::vector<float> const& TrkPitch,
                           geo::PlaneID planeID)
  {

    fPlaneID = planeID;
    fKineticEnergy = KineticEnergy;
    fRange = Range;
    fTrkPitch = TrkPitch;
    if(dEdx.size() != resRange.size())
      throw cet::exception("anab::Calorimetry") << "dE/dx and residual range vectors "
                                                << "have different sizes, this is a problem.\n";
    for(size_t i=0; i!=dQdx.size(); ++i){
      fXYZ.push_back({-999.,-999.,-999.});
    }
    fdEdx.resize(dEdx.size());
    fdQdx.resize(dQdx.size());
    fResidualRange.resize(resRange.size());
    for(size_t i = 0; i < dEdx.size(); ++i){
      fdEdx[i]          = dEdx[i];
      fdQdx[i]          = dQdx[i];
      fResidualRange[i] = resRange[i];
    }

    fDeadWireResR.resize(deadwire.size());
    for(size_t i = 0; i<deadwire.size(); ++i){
      fDeadWireResR[i] = deadwire[i];
    }

  }

  //----------------------------------------------------------------------
  Calorimetry::Calorimetry(float KineticEnergy,
                           std::vector<float> const& dEdx,
                           std::vector<float> const& dQdx,
                           std::vector<float> const& resRange,
                           std::vector<float> const& deadwire,
                           float Range,
                           std::vector<float> const& TrkPitch,
                           std::vector<anab::Point_t> const& XYZ,
                           geo::PlaneID planeID)
    : Calorimetry(KineticEnergy, dEdx, dQdx, resRange, deadwire, Range, TrkPitch, XYZ, std::vector<size_t>(),planeID) { }
  //----------------------------------------------------------------------
  Calorimetry::Calorimetry(float KineticEnergy,
                           std::vector<float> const& dEdx,
                           std::vector<float> const& dQdx,
                           std::vector<float> const& resRange,
                           std::vector<float> const& deadwire,
                           float Range,
                           std::vector<float> const& TrkPitch,
                           std::vector<anab::Point_t> const& XYZ,
                           std::vector<size_t> const& TpIndices,
                           geo::PlaneID planeID)
  {

    if( dEdx.size() != resRange.size() ||
        dEdx.size() != dQdx.size()     ||
        dEdx.size() != TrkPitch.size() ||
        dEdx.size() != XYZ.size()      ||
        (TpIndices.size()>0 && dEdx.size() != TpIndices.size()) )
      throw cet::exception("anab::Calorimetry") << "Input vectors "
                                                << "have different sizes, this is a problem.\n";
    fPlaneID = planeID;
    fKineticEnergy = KineticEnergy;
    fRange = Range;
    fTrkPitch = TrkPitch;
    fdEdx = dEdx;
    fdQdx = dQdx;
    fResidualRange = resRange;
    fXYZ = XYZ;
    fTpIndices = TpIndices;
    fDeadWireResR = deadwire;
  }

  //----------------------------------------------------------------------
  // ostream operator.
  //
  std::ostream& operator<< (std::ostream & o, Calorimetry const& a)
  {
    o << "Kinetic Energy: " << a.fKineticEnergy
      << "\n Range: "         << a.fRange << std::endl;

    for(size_t n = 0; n < a.fdEdx.size(); ++n)
      o << "dE/dx="           << a.fdEdx[n]
        << " Residual range=" << a.fResidualRange[n]
        << " dQ/dx=" << a.fdQdx[n]
        << " (x,y,z)=(" << a.fXYZ[n].X() << "," << a.fXYZ[n].Y() << "," << a.fXYZ[n].Z() << ")"
        << " pitch=" << a.fTrkPitch[n]
        << " planeID=(" << a.fPlaneID.Cryostat << "," << a.fPlaneID.TPC << "," << a.fPlaneID.Plane << ")"
        << std::endl;

    return o;
  }

}