slid::ParticleIDAlg Class Reference

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-01-16/ShowerLID/ParticleIDAlg.h"

Public Member Functions
	ParticleIDAlg (fhicl::ParameterSet const &pset, novadaq::cnv::DetId detId, NuEEnergyAlg *NuEEnergyAlgIn, bool loadLibs=true)
	~ParticleIDAlg ()
void	reconfigure (const fhicl::ParameterSet &pset, novadaq::cnv::DetId detId)
	Reinitialize parameters read from FCL file. More...
double	DedxTransLL (const slid::DedxParticleType partHyp, const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
	Calculate transverse dE/dx log-likelihood for specific particle hypothesis. More...
double	DedxLongLL (const slid::DedxParticleType partHyp, const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
	Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis. More...
double	DedxInverseLongLL (const slid::DedxParticleType partHyp, const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
	Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis assuming that the particle direction if opposite to the reconstructed direction. More...
void	SetShower (const rb::Shower *vShower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
	Set rb::Prong to be analyzed. This must be set before any calculations can be done. More...
double	GetGapVertexToShowerStart (const rb::Shower *vShower, TVector3 evtVtx, const art::Event &evt)
	Calculate the gap between the vertex and the start of the shower. More...
double	Pi0Mass (const std::vector< const rb::Shower * > showerCollection, unsigned int iShower, const art::Event &evt)
	Calculate the mass of candidate and any other shower in the slice that is closest to the pi0 mass. More...
int	Pi0SecondPhoton (const std::vector< const rb::Shower * > showerCollection, unsigned int iShower, const art::Event &evt)
	Index in showerCollection of the shower for which the reconstructed invariant mass was closest to pi0 mass for the shower at index iShower. More...
double	Radius (const std::vector< const rb::Shower * > showerCollection, unsigned int iShower, const art::Event &evt)
	Calculate the shower radius. More...
double	PlaneLongDedx (unsigned int pIdx)
	return longitudinal dedx for a specified plane index More...
double	PlaneTransDedx (unsigned int tpIdx)
	return transverse dedx for specified transverse plane More...
double	CellPlaneDedx (int tpIdx, unsigned int pIdx)
	return transverse dedx for specified transverse plane More...
double	PlaneRadius (unsigned int pIdx)
	return shower radius for a plane index More...
TVector3	PlaneHitXYZ (unsigned int pIdx)
	Return point of intersection between shower core and plane. More...
double	PlaneHitCell (unsigned int pIdx)
	Return point of intersection cell between shower core and plane. More...
unsigned int	PlaneToGlobalPlane (unsigned int pIdx)
	Return the gloabel plane index for an index in the shower coordinate. More...
double	PlaneCentroid (unsigned int pIdx)
	Return the energy weighted mean shower positon for the plane in cm. More...
double	ShowerAsym () const
	Return the asymmetry of the shower. More...
TVector3	InertiaXYZ () const
	Return the inertia calculation. More...
double	PointDoca (TVector3 vtx, TVector3 start, TVector3 stop)
	Get distance of closest approach between shower and vertex. More...
int	GetMIPPlanes (const rb::Shower *shower, const std::vector< const rb::Shower * > showercol, const art::Event &evt)
	Number of MIP planes within the first 20 planes of the start of the shower. More...

Public Attributes
std::map< int, float >	fPartLongLL
	Map of the longitudinal ll by paricle type. More...
std::map< int, float >	fPartTransLL
	Map of the transverse ll by particle type. More...

Private Member Functions
double	CalcPlaneLongitudinalDedx (unsigned int plane)
	Calculate the longitudinal dE/dx in a given plane. Input is the actuall plane number in the detector. More...
double	CalcPlaneTransverseDedx (unsigned int iTransversePlane)
	Calculate the transverse dE/dx in a given tranverse slice. Input is the transverse slice index. More...
double	CalcCellPlaneTransverseDedx (int iTransversePlane, unsigned int iPlane)
	Calculate the transverse dE/dx in a given tranverse cell and plane. Input is the transverse slice index, including plus and minus to the shower core, and plane number. More...
std::map< int, std::map< int, int > >	CalcPlaneHits (const rb::Shower *vShower)
	Return a map off all hits in a given longitudinal plane indexed by the longitudinalplane for shower already read into object by SetShower. More...
std::map< int, double >	CalcPlaneCentroid (const rb::Shower *vShower)
	Return a map of mean cell posisiton for each plane for shower already read into object by SetShower. More...
double	CalcEnergyInLongitudinalPlane (int iPlane)
	Sum the energy in a longitudinal plane. More...
void	CalcDetectorXYRegionIndex ()
	Determine detector XY region. More...
double	CalcPlaneLongDedxProb (DedxParticleType particleType, int igev, unsigned int iPlane, unsigned int plane)
	Get longitudinal dE/dx probability for a given longitudinal plane, energy bin, and particle hypothesis. More...
double	CalcPlaneTransDedxProb (DedxParticleType particleType, int igev, unsigned int iTransversePlane)
	Get transverse dE/dx probability for a given transverse plane, energy bin, and particle hypothesis. More...
double	CalcInterPlaneDedxProb (DedxParticleType type, unsigned int iLongPlane, unsigned int plane)
	From Jianming's code, still working out changes. More...
double	CalcInterCellTransDedxProb (const DedxParticleType type, unsigned int iTransPlane)
	From Jianming's code, still working out changes. More...
bool	LoadDedxHistogramFiles (const fhicl::ParameterSet &pset, novadaq::cnv::DetId geom)
	Load all of the histograms from files into arrays of TH1D* objects. The details of handling the histograms are taken care of within object of the helper class slid::DedxDistribution. More...
double	CalcTrkHitPath (unsigned int plane)
	Calculate path length of shower through this plane. More...
void	CalculateEnergyBinIndex ()
	Calculate energy bin indices and energy values from the histograms corresponding energy of shower to be used to interpolate between bins in energy. More...
void	CalcAsymIneria ()
	Calculate the asymmetry and inertia of the shower. More...
TVector3	GetCellNodePos (unsigned int plane1, unsigned int cell1, unsigned int plane2, unsigned int cell2)
	Helper function for interia and asymmetry calculation. More...

Private Attributes
art::EventID	eventID
	data/MC flag More...
double	fClosestPi0Mass
	Closest mass to pi0 mass formed by this shower and any other shower in the slice. More...
std::vector< slid::DedxDistribution >	fDedxDistributions
	dE/dx distributions (histogram files and data) for each detector type and particle category More...
const rb::Shower *	fShower
	rb::Shower object for which particle ID is to be calculated. This is cached. More...
std::vector< const rb::Shower * >	fShowerCol
double	fShowerEnergyGev
	Shower energy. More...
art::ServiceHandle< geo::Geometry >	fGeom
	geometry More...
NuEEnergyAlg *	fNuEEnergyAlg
	Nue energy algorithms. More...
std::map< int, double >	fDedxLongByPlane
std::map< int, double >	fDedxTransByPlane
std::map< int, std::map< int, int > >	fPlaneHits
std::map< unsigned int, unsigned int >	fMapIdxToGlobalPlane
std::map< int, double >	fPlaneCentroid
std::map< int, double >	fTransEnergy
	Map of transverse energy for each transverse plane. More...
std::map< int, double >	fTransPathLength
	Map of path length for each transverse plane. More...
int	fPi0Index
int	iReg
	Detector sub-region bin. Calculated once and cached for a given shower/event combo. More...
TVector3	pos
	Shower start position. More...
int	igev0
	Shower adjacent energy bins. Calculated once and cached for a given shower/event combo. More...
int	igev1
double	e0
	Energy values of energy bins adjacent to shower. Calculated once and cached for a given shower/event combo. More...
double	e1
TVector3	fIxyz
	shower moment of ineria More...
double	fAsym
	shower asymmetry More...
double	fMipRangeHigh
	Upper bound of MIP Range. More...
double	fMipRangeLow
	Lower bound of MIP Range. More...
novadaq::cnv::DetId	fDetId
	Detector ID. More...

Detailed Description

Definition at line 50 of file ParticleIDAlg.h.

Constructor & Destructor Documentation

slid::ParticleIDAlg::ParticleIDAlg ( fhicl::ParameterSet const &  pset, novadaq::cnv::DetId  detId, NuEEnergyAlg *  NuEEnergyAlgIn, bool  loadLibs = true  )

explicit

Definition at line 30 of file ParticleIDAlg.cxx.

References LoadDedxHistogramFiles().

    : fClosestPi0Mass(-9999),
      fDedxDistributions(slid::kDEDXPARTICLETYPESIZE),
      fShower(0),
      fNuEEnergyAlg(NuEEnergyAlgIn),
      fAsym(1.0),
      fMipRangeHigh(0.005),
      fMipRangeLow(0.0005),
      fDetId(detId)
  {
    if (loadLibs){
      LoadDedxHistogramFiles(pset, detId);
    }
  }

slid::ParticleIDAlg::~ParticleIDAlg

(

)

Definition at line 49 of file ParticleIDAlg.cxx.

  {
  }

Member Function Documentation

void slid::ParticleIDAlg::CalcAsymIneria ( )

private

Calculate the asymmetry and inertia of the shower.

Definition at line 1037 of file ParticleIDAlg.cxx.

References slid::NuEEnergyAlg::CellEnergy(), rb::Prong::Dir(), e, eventID, stan::math::fabs(), fAsym, fIxyz, fMapIdxToGlobalPlane, fNuEEnergyAlg, fPlaneHits, fShower, fShowerCol, geom(), GetCellNodePos(), Mag(), file_size_ana::node, geo::GeometryBase::Plane(), PlaneHitXYZ(), cet::pow(), r(), std::sqrt(), rb::Prong::Start(), and geo::PlaneGeo::View().

Referenced by SetShower().

  {
    art::ServiceHandle<geo::Geometry> geom;
    double I11 = 0;
    double I12 = 0;
    double I21 = 0;
    double I22 = 0;
  
    // Using 3-D nodes to calculate shower balance
    //     TVector3 core(shower.Start()[0],shower.Start()[1],shower.Start()[2]); // original point of shower frame
    TVector3 iz = fShower->Dir().Unit(); // z in shower frame
    TVector3 ix0(1.0,0.0,0);
    TVector3 iy0(0.0,1.0,0);
    TVector3 ix = iy0.Cross(iz).Unit(); // x axis in shower frame
    TVector3 iy = iz.Cross(ix).Unit(); // y axis in shower frame
  
    std::map<unsigned int, unsigned int>::iterator itr2 = fMapIdxToGlobalPlane.begin();
    for(std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.begin(); itr != fMapIdxToGlobalPlane.end(); ++itr){
      itr2++;
      std::map<int, int> planecells1 = fPlaneHits[itr->second];
      if (itr2 == fMapIdxToGlobalPlane.end()) break;
      std::map<int, int> planecells2 = fPlaneHits[itr2->second];
      if(planecells1.size()==0)continue;
      unsigned int plane1 = itr->second;
      if(planecells2.size()==0)continue;
      unsigned int plane2 = itr2->second;
      geo::View_t view1 = geom->Plane(plane1)->View();
      geo::View_t view2 = geom->Plane(plane2)->View();
      if(view1==view2)continue;
      TVector3 core = PlaneHitXYZ(itr->first);
      TVector3 pcore = PlaneHitXYZ(itr2->first);
      std::vector<TVector3> planenodecol;
      std::vector<double> planenodeenergycol;
      planenodecol.clear();
      planenodeenergycol.clear();
      for (std::map<int, int>::iterator citr = planecells1.begin(); citr != planecells1.end(); ++citr) {
        unsigned int cell1 = citr->first;
        for(std::map<int, int>::iterator citr2 = planecells2.begin(); citr2 != planecells2.end(); ++citr2){
          unsigned int cell2 = citr2->first;
          TVector3 node = GetCellNodePos(plane1, cell1, plane2, cell2);// define 3-D node
          double nodez = (pcore - fShower->Start()).Mag();
          TVector3 nodetopcore = node - pcore;
          nodetopcore.SetZ(nodez);
          TVector3 nodetocore = node - core;
          nodetocore.SetZ(nodez);

          TVector3 xyz = node - core;//transform xy to shower frame 

          double p = fShower->Dir()[0];
          double q = fShower->Dir()[1];
          double r = fShower->Dir()[2];
          double x0 = xyz[0];
          double y0 = xyz[1];
          double z0 = xyz[2];
          TVector3 vt(q*z0-r*y0,r*x0-p*z0,p*y0-q*x0);

          double nodeenergy = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, citr->second, eventID) +
            fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, citr2->second, eventID);;

          I11 += nodeenergy*(nodetocore.Y()*nodetocore.Y());
          I12 += -nodeenergy*(nodetocore.X()*nodetocore.Y());
          I21 += -nodeenergy*(nodetocore.X()*nodetocore.Y());
          I22 += nodeenergy*(nodetocore.X()*nodetocore.X());
  
        }
      }
    }
    // calculate eiganvalues and diagonalize inertial tensor
    double Ixx = 0.5*(I11+I22+sqrt(pow(I11+I22,2)-4*(I11*I22-I12*I21)));
    double Iyy = 0.5*(I11+I22-sqrt(pow(I11+I22,2)-4*(I11*I22-I12*I21)));
    double asym = 1;
    if(fabs(Ixx+Iyy)>1e-10)asym=(Ixx-Iyy)/(Ixx+Iyy);

    TVector3 ixyz(Ixx,Iyy,0);
    fIxyz = ixyz;
    fAsym = asym;
  }

double slid::ParticleIDAlg::CalcCellPlaneTransverseDedx ( int  iTransversePlane, unsigned int  iPlane  )

private

Calculate the transverse dE/dx in a given tranverse cell and plane. Input is the transverse slice index, including plus and minus to the shower core, and plane number.

Definition at line 607 of file ParticleIDAlg.cxx.

References abs(), CalcTrkHitPath(), slid::NuEEnergyAlg::CellEnergy(), dir, rb::Prong::Dir(), eventID, stan::math::fabs(), fGeom, fMapIdxToGlobalPlane, fNuEEnergyAlg, fPlaneCentroid, fPlaneHits, fShower, fShowerCol, makeTrainCVSamples::int, std::isnan(), geo::PlaneGeo::Ncells(), path, geo::GeometryBase::Plane(), and PlaneHitCell().

Referenced by CellPlaneDedx().

  {
    double dedx = 0;
    double transcellE = 0; 
    //loop over hits in longitudinal planes to calculate transverse plane quantities
    unsigned int gplane = 9999;
    std::map<unsigned int, unsigned int>::iterator itr0 = fMapIdxToGlobalPlane.find(iPlane);
    if (itr0 != fMapIdxToGlobalPlane.end()){
//      std::cout<<iPlane<<" is global "<<"itr0->second "<<itr0->second<<std::endl;
      gplane = itr0->second;
    }else{
//      std::cout<<iPlane<<" "<<"out fMapIdxToGlobalPlane.end() total "<<fMapIdxToGlobalPlane.size()<<std::endl;
      return dedx;
    }
    int pcount=-1;
    for(std::map<int, std::map<int, int> >::iterator itr=fPlaneHits.begin(); itr!= fPlaneHits.end(); ++itr){
      pcount++;
      if((unsigned)itr->first!=gplane)continue;
      int cplane = itr->first;
      std::map<int, int> planecells = itr->second;
      if (planecells.size() == 0) break;
      if (fPlaneCentroid[itr->first] == -999.0) break;
//      double avgPos0 = fPlaneCentroid[itr->first];
      unsigned int planeidx = 0;
      if (fShower->Dir()[2] >= 0.0){ //forward going shower
        planeidx = (unsigned int)pcount;
      }else{
        planeidx = (unsigned int)(fPlaneHits.size()-pcount-1);
      }
      double avgPos = PlaneHitCell(planeidx); 
//      std::cout<<"avg Pos0, hit Pos "<<avgPos0<<" "<<avgPos<<std::endl;
      int meanCell = (int)avgPos;
      double wm = avgPos - (double)meanCell;//energy splitting for the centroid cell
      double wp = 1-wm; //energy splitting for the centroid cell
      if(meanCell<0)meanCell=0;
      if(meanCell>(int)fGeom->Plane(cplane)->Ncells()-1)meanCell=(int)fGeom->Plane(cplane)->Ncells()-1;
      double dir = 1.;
      if(fabs(fShower->Dir()[2])>0.000001)dir=fabs(fShower->Dir()[2]);
      int ic = abs((int)iTransversePlane) + 1;
      double dph = ((double)ic/dir)+wm;//upper limit for cells included in a transverse cell index
      double dpl = ((double)(ic-1)/dir)+wm;//lower limit for cells included in a transverse cell index
      //if the mean cell is inside the detector

      if(iTransversePlane>=0){
        if(meanCell+(int)dph<(int)fGeom->Plane(cplane)->Ncells()){
          double pfHigh = dph-(int)dph; //fraction of last cell on high end of transverse plane
          double pfLow = (int)dpl+1-dpl; //fraction of last cell on low end of transverse plane
          int cellHigh = meanCell + (int)dph;
          int cellLow  = meanCell + (int)dpl;
          std::map<int, int> cellsInTransPlane;
          cellsInTransPlane.insert(planecells.lower_bound(cellLow),
                                   planecells.upper_bound(cellHigh));
          for (std::map<int, int>::iterator itrC = cellsInTransPlane.begin(); itrC != cellsInTransPlane.end(); ++itrC){
            double cellE = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, itrC->second, eventID, false);
            if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
            if (itrC->first == cellLow)  transcellE += pfLow*cellE;
            else if (itrC->first == cellHigh) transcellE += pfHigh*cellE;
            else transcellE += cellE;
          }
        }
      } else {

        double dmh = ((double)ic/dir)+wp;//upper limit for cells included in a transeverse cell index
        double dml = ((double)(ic-1)/dir)+wp;//lower limit for cells included in a transeverse cell index
        if(meanCell-(int)dmh>=0){
          double pfHigh = dmh-(int)dmh; //fraction of last cell on high end of transverse plane
          double pfLow = (int)dml+1-dml; //fraction of last cell on low end of transverse plane
          int cellHigh = meanCell - (int)dmh;
          int cellLow  = meanCell - (int)dml;
          std::map<int, int> cellsInTransPlane;
          cellsInTransPlane.insert(planecells.lower_bound(cellHigh),
                                   planecells.upper_bound(cellLow));
          for (std::map<int, int>::iterator itrC = cellsInTransPlane.begin(); itrC != cellsInTransPlane.end(); ++itrC){
            double cellE = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, itrC->second, eventID, false);
            if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
            if (itrC->first == cellLow)  transcellE += pfLow*cellE;
            else if (itrC->first == cellHigh) transcellE += pfHigh*cellE;
            else transcellE += cellE;
          }
        }
      }
      break;
    }
    double path = CalcTrkHitPath(gplane); 
    if(path>0.00001)dedx = transcellE/path;
    return dedx;
  }

void slid::ParticleIDAlg::CalcDetectorXYRegionIndex ( )

private

Determine detector XY region.

Definition at line 416 of file ParticleIDAlg.cxx.

References fDetId, fShower, iReg, cmf::kFARDET, pos, rb::Prong::Start(), and rb::Shower::Stop().

Referenced by SetShower().

  {

    TVector3 pos((fShower->Start()[0] + fShower->Stop()[0]) / 2.0,
                 (fShower->Start()[1] + fShower->Stop()[1]) / 2.0,
                 (fShower->Start()[2] + fShower->Stop()[2]) / 2.0);
    iReg = 0;
    //only switch quadrants if this is the Far Detector, for ND whole detector is ~quadrant 0
    // (near the readout in both views)
    if (fDetId == novadaq::cnv::DetId::kFARDET){
      if (pos[0] >= 0 && pos[1] >= 0)iReg = 0;
      if (pos[0] >= 0 && pos[1] < 0)iReg = 1;
      if (pos[0] < 0 && pos[1] >= 0)iReg = 2;
      if (pos[0] <= 0 && pos[1] <= 0)iReg = 3;
    }
    return;
  }

double slid::ParticleIDAlg::CalcEnergyInLongitudinalPlane ( int  iPlane )

private

Sum the energy in a longitudinal plane.

Definition at line 492 of file ParticleIDAlg.cxx.

References slid::NuEEnergyAlg::CellEnergy(), slid::NuEEnergyAlg::ECorrMC(), eventID, fNuEEnergyAlg, fPlaneHits, fShower, fShowerCol, and std::isnan().

Referenced by CalcPlaneLongitudinalDedx().

                                                                {
    double energySum = 0;
    std::map<int, int> planecells = fPlaneHits[iPlane];
    for (std::map<int, int>::iterator itr=planecells.begin(); itr!= planecells.end(); ++itr) {
      double cellE = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, itr->second, eventID);
      if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
      energySum += cellE;
    }
    double ecor = fNuEEnergyAlg->ECorrMC(energySum);
    return energySum/ecor;
  }

double slid::ParticleIDAlg::CalcInterCellTransDedxProb ( const DedxParticleType  type, unsigned int  iTransPlane  )

private

From Jianming's code, still working out changes.

Definition at line 756 of file ParticleIDAlg.cxx.

References CalcPlaneTransDedxProb(), e0, e1, igev0, and igev1.

Referenced by DedxTransLL().

  {
    double probPXPY = 0;
    double probe0 = CalcPlaneTransDedxProb(type, igev0, iTransPlane);
    double probe1 = CalcPlaneTransDedxProb(type, igev1, iTransPlane);
    if (igev0 != igev1)probPXPY = (e1 * probe0 + e0 * probe1) / (e0 + e1);
    if (igev0 == igev1)probPXPY = probe0;
    return probPXPY;
  }

double slid::ParticleIDAlg::CalcInterPlaneDedxProb ( DedxParticleType  type, unsigned int  iLongPlane, unsigned int  plane  )

private

From Jianming's code, still working out changes.

Definition at line 729 of file ParticleIDAlg.cxx.

References CalcPlaneLongDedxProb(), std::cos(), rb::Prong::Dir(), e0, e1, stan::math::fabs(), fShower, igev0, igev1, and makeTrainCVSamples::int.

Referenced by DedxLongLL().

  {

    
    double probPXPY = 0;
    double cos = fShower->Dir()[2];
    unsigned int plane0 = int( int(iLongPlane) / fabs(cos));
    unsigned int plane1 = int( int(iLongPlane) / fabs(cos)) + 1;
    double r0 = double(iLongPlane / fabs(cos)) - double(plane0);
    double r1 = 1 - r0;
    double probe0r0 = 1;
    double probe0r1 = 1;
    double probe1r0 = 1;
    double probe1r1 = 1;
    probe0r0 = CalcPlaneLongDedxProb(type, igev0, plane0, plane);
    probe0r1 = CalcPlaneLongDedxProb(type, igev0, plane1, plane);
    probe1r0 = CalcPlaneLongDedxProb(type, igev1, plane0, plane);
    probe1r1 = CalcPlaneLongDedxProb(type, igev1, plane1, plane);

    double probe0 = (r1 * probe0r0 + r0 * probe0r1);
    double probe1 = (r1 * probe1r0 + r0 * probe1r1);
    if (igev0 != igev1)probPXPY = (e1 * probe0 + e0 * probe1) / (e0 + e1);
    if (igev0 == igev1)probPXPY = probe0;
    return probPXPY;
  }

std::map< int, double > slid::ParticleIDAlg::CalcPlaneCentroid ( const rb::Shower *  vShower )

private

Return a map of mean cell posisiton for each plane for shower already read into object by SetShower.

Definition at line 278 of file ParticleIDAlg.cxx.

References slid::NuEEnergyAlg::CellEnergy(), cellPos(), eventID, fNuEEnergyAlg, fPlaneHits, fShowerCol, and std::isnan().

Referenced by SetShower().

  {

    std::map< int, double > mapT;

    //loop plane by plane to calculate weighted cell center and then transverse cells
    for(std::map<int, std::map<int, int> >::iterator itr=fPlaneHits.begin(); itr!=fPlaneHits.end(); ++itr){
      std::map<int, int> planecells = itr->second;
      //loop over hits in plane to determine energy weighted center
      double sumPos = 0.0;
      double planeGev = 0.0;
      for(std::map<int, int>::iterator itrC = planecells.begin(); itrC!=planecells.end(); ++itrC){
        double cellPos = itrC->first + 0.5;
        double cellE   = fNuEEnergyAlg->CellEnergy(vShower, fShowerCol, itrC->second, eventID );
        if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
        sumPos += cellPos*cellE;
        planeGev += cellE;
      }
      double avgPos = -999.0; //non physical cell position as the default
      if (planeGev > 0.00001) avgPos = sumPos/planeGev;
      mapT[itr->first] = avgPos;
    }

    return mapT;
  }

std::map< int, std::map< int, int > > slid::ParticleIDAlg::CalcPlaneHits ( const rb::Shower *  vShower )

private

Return a map off all hits in a given longitudinal plane indexed by the longitudinalplane for shower already read into object by SetShower.

Definition at line 241 of file ParticleIDAlg.cxx.

References rb::CellHit::Cell(), rb::Cluster::Cell(), rb::Prong::Dir(), fMapIdxToGlobalPlane, MECModelEnuComparisons::i, rb::Cluster::MaxPlane(), rb::Cluster::MinPlane(), rb::Cluster::NCell(), ncells, and rb::CellHit::Plane().

Referenced by SetShower().

  {
    std::map< int, std::map< int, int> > mapP;

    std::map<int, int> dummymap;

    int startPlane = vShower->MinPlane();
    int stopPlane  = vShower->MaxPlane();
    for (int i=startPlane; i<=stopPlane; ++i){
      mapP[i] = dummymap;
    }
    
    int ncells = vShower->NCell();
    for (int i = 0; i < ncells; ++i) {

      mapP[ vShower->Cell(i)->Plane() ][vShower->Cell(i)->Cell()] = i;
    }// end loop over cells

    //now map the plane index with the global plane number
    unsigned int planecount = 0;
    if (vShower->Dir()[2] >= 0.0){ //forward going shower
      for (std::map<int, std::map<int, int> >::iterator itr = mapP.begin(); itr != mapP.end(); ++itr){
        fMapIdxToGlobalPlane[planecount] = itr->first;
        planecount++;
      }
    }
    else{ //backward going shower
      for (std::map<int, std::map<int, int> >::reverse_iterator itr = mapP.rbegin(); itr != mapP.rend(); ++itr){
        fMapIdxToGlobalPlane[planecount] = itr->first;
        planecount++;
      }
    }

    return mapP;
  }// end of CalcPlaneHits

double slid::ParticleIDAlg::CalcPlaneLongDedxProb ( DedxParticleType  particleType, int  igev, unsigned int  iPlane, unsigned int  plane  )

private

Get longitudinal dE/dx probability for a given longitudinal plane, energy bin, and particle hypothesis.

Definition at line 696 of file ParticleIDAlg.cxx.

References CalcPlaneLongitudinalDedx(), fDedxDistributions, GetEntries(), HandyFuncs::Interpolate(), iReg, slid::kNumEnergyBin, slid::kNumLongitudinalPlane, and slid::kNumXYRegion.

Referenced by CalcInterPlaneDedxProb().

  {
    double prob = 0;

    if (iReg < 0 || iReg > kNumXYRegion - 1) return 1;
    if (iGev >= kNumEnergyBin || iGev < 0) return 1;
    if (iPlane >= kNumLongitudinalPlane) return 1;
    double dedx = CalcPlaneLongitudinalDedx(plane);
    double normHist = double(fDedxDistributions[particleType].GetLongitudinalPlaneHitDedxHist(iReg, iGev, iPlane)->GetEntries());
    double dedxHistValue = double(fDedxDistributions[particleType].GetLongitudinalPlaneHitDedxHist(iReg, iGev, iPlane)->Interpolate(dedx));
    double nHistBins = fDedxDistributions[particleType].GetLongitudinalPlaneHitDedxHist(iReg, iGev, iPlane)->GetNbinsX();
    if (dedxHistValue < 0.0001) dedxHistValue = 0.0001;
    if (normHist >= 1) prob = dedxHistValue * nHistBins / normHist;
    return prob;
  }

double slid::ParticleIDAlg::CalcPlaneLongitudinalDedx ( unsigned int  plane )

private

Calculate the longitudinal dE/dx in a given plane. Input is the actuall plane number in the detector.

Definition at line 506 of file ParticleIDAlg.cxx.

References CalcEnergyInLongitudinalPlane(), CalcTrkHitPath(), and fDedxLongByPlane.

Referenced by CalcPlaneLongDedxProb(), GetMIPPlanes(), and PlaneLongDedx().

  {
    double dedx = 0.0;
    //
    //  Lazy caching-- only calculate if fDedxLongByPlane[iPlane] not yet set.
    //
    std::map< int, double>::iterator planeIter = fDedxLongByPlane.find(iPlane);
    //std::cout<<"Asked to calculate long dedx for plane "<<iPlane<<std::endl;
    if (planeIter == fDedxLongByPlane.end()) {
      double energyInPlane = CalcEnergyInLongitudinalPlane(iPlane);
      double pathLength = CalcTrkHitPath(iPlane);

      if (pathLength > 0.0) dedx = energyInPlane / pathLength;
      
      fDedxLongByPlane.insert(planeIter, std::pair< int, double>(iPlane, dedx));
      //std::cout<<"plane idx: "<<iPlane<<" path: "<<pathLength<<" ep: "<<energyInPlane<<" dedx: "<<dedx<<std::endl;
    }
    else dedx = planeIter->second;
    return dedx;
  }

double slid::ParticleIDAlg::CalcPlaneTransDedxProb ( DedxParticleType  particleType, int  igev, unsigned int  iTransversePlane  )

private

Get transverse dE/dx probability for a given transverse plane, energy bin, and particle hypothesis.

Definition at line 713 of file ParticleIDAlg.cxx.

References CalcPlaneTransverseDedx(), fDedxDistributions, GetEntries(), HandyFuncs::Interpolate(), iReg, slid::kNumEnergyBin, slid::kNumTransversePlane, and slid::kNumXYRegion.

Referenced by CalcInterCellTransDedxProb().

  {
    double prob = 0;
    if (iReg < 0 || iReg > kNumXYRegion - 1)return 1;
    if (iGev >= kNumEnergyBin || iGev < 0)return 1;
    if (iTransversePlane >= kNumTransversePlane)return 1;
    double dedx = CalcPlaneTransverseDedx(iTransversePlane);
    double normHist = double(fDedxDistributions[particleType].GetTransversePlaneCellDedxHist(iReg, iGev, iTransversePlane)->GetEntries());
    double dedxHistValue = double(fDedxDistributions[particleType].GetTransversePlaneCellDedxHist(iReg, iGev, iTransversePlane)->Interpolate(dedx));
    double nHistBins = fDedxDistributions[particleType].GetTransversePlaneCellDedxHist(iReg, iGev, iTransversePlane)->GetNbinsX();
    if (dedxHistValue < 0.0001)dedxHistValue = 0.0001;
    if (normHist >= 1) prob = dedxHistValue * nHistBins / normHist;
    return prob;
  }

double slid::ParticleIDAlg::CalcPlaneTransverseDedx ( unsigned int  iTransversePlane )

private

Calculate the transverse dE/dx in a given tranverse slice. Input is the transverse slice index.

Definition at line 530 of file ParticleIDAlg.cxx.

References CalcTrkHitPath(), slid::NuEEnergyAlg::CellEnergy(), dir, rb::Prong::Dir(), eventID, stan::math::fabs(), fDedxTransByPlane, fGeom, fNuEEnergyAlg, fPlaneCentroid, fPlaneHits, fShower, fShowerCol, fTransEnergy, fTransPathLength, makeTrainCVSamples::int, std::isnan(), geo::PlaneGeo::Ncells(), and geo::GeometryBase::Plane().

Referenced by CalcPlaneTransDedxProb(), DedxTransLL(), and PlaneTransDedx().

  {
    double dedx = -1.0;
    //
    //  Lazy caching-- only calculate if fDedxTransByPlane[iPlane] not yet set.
    //
    std::map< int, double>::iterator transPlaneIter = fDedxTransByPlane.find(iTransversePlane);
    if (transPlaneIter == fDedxTransByPlane.end()) {
      //loop over hits in longitudinal planes to calculate transverse plane quantities
      for(std::map<int, std::map<int, int> >::iterator itr=fPlaneHits.begin(); itr!= fPlaneHits.end(); ++itr){
        int cplane = itr->first;
        std::map<int, int> planecells = itr->second;
        if (planecells.size() == 0) continue;
        if (fPlaneCentroid[itr->first] == -999.0) continue;
        double avgPos = fPlaneCentroid[itr->first];
        int meanCell = (int)avgPos;
        double wm = avgPos - (double)meanCell;//energy splitting for the centroid cell
        double wp = 1-wm; //energy splitting for the centroid cell
        if(meanCell<0)meanCell=0;
        if(meanCell>(int)fGeom->Plane(cplane)->Ncells()-1)meanCell=(int)fGeom->Plane(cplane)->Ncells()-1;
        double dir = 1.;
        if(fabs(fShower->Dir()[2])>0.000001)dir=fabs(fShower->Dir()[2]);
        int ic = (int)iTransversePlane + 1;
        double dph = ((double)ic/dir)+wm;//upper limit for cells included in a transverse cell index
        double dpl = ((double)(ic-1)/dir)+wm;//lower limit for cells included in a transverse cell index
        //if the mean cell is inside the detector
        if(meanCell+(int)dph<(int)fGeom->Plane(cplane)->Ncells()){
          double pfHigh = dph-(int)dph; //fraction of last cell on high end of transverse plane
          double pfLow = (int)dpl+1-dpl; //fraction of last cell on low end of transverse plane
          int cellHigh = meanCell + (int)dph;
          int cellLow  = meanCell + (int)dpl;
          std::map<int, int> cellsInTransPlane;
          cellsInTransPlane.insert(planecells.lower_bound(cellLow),
                                   planecells.upper_bound(cellHigh));
          for (std::map<int, int>::iterator itrC = cellsInTransPlane.begin(); itrC != cellsInTransPlane.end(); ++itrC){
            double cellE = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, itrC->second, eventID);
            if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
            if (itrC->first == cellLow)  fTransEnergy[iTransversePlane] += pfLow*cellE;
            else if (itrC->first == cellHigh) fTransEnergy[iTransversePlane] += pfHigh*cellE;
            else fTransEnergy[iTransversePlane] += cellE;
          }

          fTransPathLength[iTransversePlane] += CalcTrkHitPath(cplane);
        }
        double dmh = ((double)ic/dir)+wp;//upper limit for cells included in a transeverse cell index
        double dml = ((double)(ic-1)/dir)+wp;//lower limit for cells included in a transeverse cell index

        if(meanCell-(int)dmh>=0){
          double pfHigh = dmh-(int)dmh; //fraction of last cell on high end of transverse plane
          double pfLow = (int)dml+1-dml; //fraction of last cell on low end of transverse plane
          int cellHigh = meanCell - (int)dmh;
          int cellLow  = meanCell - (int)dml;
          std::map<int, int> cellsInTransPlane;
          cellsInTransPlane.insert(planecells.lower_bound(cellHigh),
                                   planecells.upper_bound(cellLow));
          for (std::map<int, int>::iterator itrC = cellsInTransPlane.begin(); itrC != cellsInTransPlane.end(); ++itrC){
            double cellE = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, itrC->second, eventID);
            if ((cellE == -5) || std::isnan(cellE)) cellE = 0.0;
            if (itrC->first == cellLow)  fTransEnergy[iTransversePlane] += pfLow*cellE;
            else if (itrC->first == cellHigh) fTransEnergy[iTransversePlane] += pfHigh*cellE;
            else fTransEnergy[iTransversePlane] += cellE;
          }
          fTransPathLength[iTransversePlane] += CalcTrkHitPath(cplane);
        }
      }
      fTransPathLength[iTransversePlane] /= 2.0;
      if(fTransPathLength[iTransversePlane]>0.00001)dedx = fTransEnergy[iTransversePlane]/fTransPathLength[iTransversePlane];
      fDedxTransByPlane[iTransversePlane] = dedx;
    }
    else dedx = transPlaneIter->second;


    return dedx;
  }

double slid::ParticleIDAlg::CalcTrkHitPath ( unsigned int  plane )

private

Calculate path length of shower through this plane.

Definition at line 473 of file ParticleIDAlg.cxx.

References geo::PlaneGeo::Cell(), rb::Prong::Dir(), e, stan::math::fabs(), fGeom, fPlaneHits, fShower, geo::CellGeo::HalfD(), Mag(), geo::GeometryBase::Plane(), util::pythag(), rb::Prong::Start(), and rb::Shower::Stop().

Referenced by CalcCellPlaneTransverseDedx(), CalcPlaneLongitudinalDedx(), and CalcPlaneTransverseDedx().

  {

    TVector3 trkDirection = fShower->Dir();
    TVector3 trkstart = fShower->Start();
    TVector3 trkstop = fShower->Stop();
    double trklength = fabs((trkstart - trkstop).Mag());
    double cellDepth = 2 * fGeom->Plane(plane)->Cell(0)->HalfD();

    // if track is only one cell in length or short, return the calculated full track length
    if (fabs(trkDirection.z()) < 1e-6 || fPlaneHits.size() == 1) {
      return trklength;
    }
    // Calculate path length through cell for track with given direction
    double trkhitpath = cellDepth * (util::pythag(trkDirection.x(),trkDirection.y(),trkDirection.z()) / fabs(trkDirection.z()));
    return trkhitpath;
  }

void slid::ParticleIDAlg::CalculateEnergyBinIndex ( )

private

Calculate energy bin indices and energy values from the histograms corresponding energy of shower to be used to interpolate between bins in energy.

Definition at line 440 of file ParticleIDAlg.cxx.

References e0, e1, fShowerEnergyGev, igev0, igev1, slid::kEnergyBins, and slid::kNumEnergyBin.

Referenced by SetShower().

  {
    //        int& igev0, int& igev1,
    //            double& e0, double& e1

    double energyBinBounds[slid::kNumEnergyBin];
    for (int ig = 0; ig < slid::kNumEnergyBin; ig++) {
      energyBinBounds[ig] = (slid::kEnergyBins[ig] + kEnergyBins[ig + 1]) / 2.;
    }

    igev0 = 0;
    igev1 = 0;
    if (fShowerEnergyGev > energyBinBounds[slid::kNumEnergyBin - 1]) {
      igev0 = slid::kNumEnergyBin - 1;
      igev1 = slid::kNumEnergyBin - 1;
    } else if (fShowerEnergyGev < energyBinBounds[0]) {
      igev0 = 0;
      igev1 = 0;
    } else {
      for (int ig = 0; ig < slid::kNumEnergyBin - 1; ig++) {
        if (fShowerEnergyGev > energyBinBounds[ig] && fShowerEnergyGev < energyBinBounds[ig + 1]) {
          igev0 = ig;
          igev1 = ig + 1;
        }
      }
    }
    if (igev0 > slid::kNumEnergyBin - 1) igev0 = slid::kNumEnergyBin - 1;
    if (igev1 > slid::kNumEnergyBin - 1) igev1 = slid::kNumEnergyBin - 1;
    e0 = fShowerEnergyGev - energyBinBounds[igev0];
    e1 = energyBinBounds[igev1] - fShowerEnergyGev;
  }

double slid::ParticleIDAlg::CellPlaneDedx	(	int	tpIdx,
		unsigned int	pIdx
	)

return transverse dedx for specified transverse plane

Definition at line 891 of file ParticleIDAlg.cxx.

References CalcCellPlaneTransverseDedx(), om::cout, allTimeWatchdog::endl, and slid::kNumTransversePlane.

Referenced by slid::LIDTraining::analyze(), and slid::LIDBuilder::produce().

  {
    double dedx = -1.0;
    if (tpIdx < slid::kNumTransversePlane) dedx = CalcCellPlaneTransverseDedx(tpIdx, pIdx);
    else std::cout<<"Requested transverse cell plane index exceeds valid range, maximum plane is: "<<slid::kNumTransversePlane<<", returning dedx=-1"<<std::endl;

    return dedx;
  }

double slid::ParticleIDAlg::DedxInverseLongLL	(	const slid::DedxParticleType	partHyp,
		const rb::Shower *	vShower,
		const std::vector< const rb::Shower * >	showercol,
		const art::Event &	evt
	)

Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis assuming that the particle direction if opposite to the reconstructed direction.

Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis.

Definition at line 150 of file ParticleIDAlg.cxx.

References DedxLongLL(), rb::Prong::Dir(), rb::Prong::SetDir(), rb::Prong::SetStart(), and rb::Shower::Stop().

Referenced by slid::LIDBuilder::produce().

  {
    rb::Shower newShower = *vShower;
    newShower.SetStart( vShower->Stop());
    newShower.SetDir( -(vShower->Dir()) );
    return DedxLongLL(partHyp, &(newShower), showercol, evt);
  }

double slid::ParticleIDAlg::DedxLongLL	(	const slid::DedxParticleType	partHyp,
		const rb::Shower *	vShower,
		const std::vector< const rb::Shower * >	showercol,
		const art::Event &	evt
	)

Calculate longitudinal dE/dx log-likelihood for specific particle hypothesis.

Definition at line 101 of file ParticleIDAlg.cxx.

References CalcInterPlaneDedxProb(), rb::Prong::Dir(), fPartLongLL, fPlaneHits, test_ParserArtEvents::log, and SetShower().

Referenced by DedxInverseLongLL(), slid::LIDBuilder::produce(), and slid::SPIDBuilder::produce().

  {
    // Lazy setting of shower
    SetShower(vShower, showercol, evt);

    std::map< int, float>::iterator llIter = fPartLongLL.find(partHyp);
    if (llIter == fPartLongLL.end()) {
      double sumLogLikelihood = 0;
      unsigned int totplane = fPlaneHits.size();
      //unsigned int nonzeroplane = 0;
      unsigned int planecount = 0;

      if (vShower->Dir()[2] >= 0.0){ //forward going shower
        for (std::map<int, std::map<int, int> >::iterator itr = fPlaneHits.begin(); itr != fPlaneHits.end(); ++itr){
          double prob = 0;
          prob = CalcInterPlaneDedxProb(partHyp, planecount, itr->first);
          if (prob < 0.0001)prob = 0.0001;
          double lgprob = log(prob);
          //nonzeroplane++;
          sumLogLikelihood += lgprob;
          planecount++;
        }
      }
      else{ //backward going shower
        for (std::map<int, std::map<int, int> >::reverse_iterator itr = fPlaneHits.rbegin(); itr != fPlaneHits.rend(); ++itr){
          double prob = 0;
          prob = CalcInterPlaneDedxProb(partHyp, planecount, itr->first);
        
          if (prob < 0.0001)prob = 0.0001;
          double lgprob = log(prob);
          //nonzeroplane++;
          sumLogLikelihood += lgprob;
          planecount++;
        }
      }
      if (totplane > 0) sumLogLikelihood = sumLogLikelihood / (double) totplane;
      if (sumLogLikelihood > 999) sumLogLikelihood = 999;
      if (sumLogLikelihood<-999) sumLogLikelihood = -999;
      fPartLongLL[partHyp] = sumLogLikelihood;
      return sumLogLikelihood;
    }
    return llIter->second;
  }

double slid::ParticleIDAlg::DedxTransLL	(	const slid::DedxParticleType	partHyp,
		const rb::Shower *	vShower,
		const std::vector< const rb::Shower * >	showercol,
		const art::Event &	evt
	)

Calculate transverse dE/dx log-likelihood for specific particle hypothesis.

Definition at line 68 of file ParticleIDAlg.cxx.

References CalcInterCellTransDedxProb(), CalcPlaneTransverseDedx(), fPartTransLL, slid::kNumTransversePlane, test_ParserArtEvents::log, and SetShower().

Referenced by slid::LIDBuilder::produce(), and slid::SPIDBuilder::produce().

  {
    // Lazy setting of shower
    SetShower(vShower, showercol, evt);
    std::map< int, float>::iterator llIter = fPartTransLL.find(partHyp);
    if (llIter == fPartTransLL.end()) {
      double sumLogProb = 0;
      int transPlane = 0;
      for (int iTransPlane = 0; iTransPlane < kNumTransversePlane; iTransPlane++) {
        double dedx = CalcPlaneTransverseDedx(iTransPlane);
        if (dedx < -0.5) continue;
        double prob = 0;
        prob = CalcInterCellTransDedxProb(partHyp, iTransPlane);
        if (prob < 0.0001) prob = 0.0001;
        double logProb = log(prob);
        sumLogProb += logProb;
        transPlane++;
      }
      if (transPlane >0 ) sumLogProb = sumLogProb / (double)transPlane;
      if (sumLogProb > 999)sumLogProb = 999;
      if (sumLogProb<-999)sumLogProb = -999;

      fPartTransLL[partHyp] = sumLogProb;
      return sumLogProb;
    }
    return llIter->second;
  }

TVector3 slid::ParticleIDAlg::GetCellNodePos ( unsigned int  plane1, unsigned int  cell1, unsigned int  plane2, unsigned int  cell2  )

private

Helper function for interia and asymmetry calculation.

Definition at line 1116 of file ParticleIDAlg.cxx.

References geo::PlaneGeo::Cell(), geom(), geo::CellGeo::GetCenter(), geo::kX, geo::kY, geo::GeometryBase::Plane(), geo::PlaneGeo::View(), submit_syst::x, submit_syst::y, and test::z.

Referenced by CalcAsymIneria().

  {
    double cellXYZ1[3]={0,0,0};
    double cellXYZ2[3]={0,0,0};
    art::ServiceHandle<geo::Geometry> geom;
    geom->Plane(plane1)->Cell(cell1)->GetCenter(cellXYZ1);
    geom->Plane(plane2)->Cell(cell2)->GetCenter(cellXYZ2);
    geo::View_t view1 = geom->Plane(plane1)->View();
    geo::View_t view2 = geom->Plane(plane2)->View();
    TVector3 nullpos(-9999,-9999,-9999);
    if(view1==view2)return nullpos;
    double x=0,y=0,z=0;
    if(view1 == geo::kX && view2 == geo::kY){
      x = cellXYZ1[0];
      y = cellXYZ2[1];
    }else if(view1 == geo::kY && view2 == geo::kX){
      x = cellXYZ2[0];
      y = cellXYZ1[1];
    }
    z = (cellXYZ1[2]+cellXYZ2[2])/2.;
    TVector3 nodepos(x,y,z);
    return nodepos;
  }

double slid::ParticleIDAlg::GetGapVertexToShowerStart	(	const rb::Shower *	vShower,
		TVector3	evtVtx,
		const art::Event &	evt
	)

Calculate the gap between the vertex and the start of the shower.

Definition at line 767 of file ParticleIDAlg.cxx.

References Mag(), PointDoca(), std::sqrt(), rb::Prong::Start(), and rb::Shower::Stop().

Referenced by slid::LIDBuilder::produce().

  {

    //initialize event vertex to start of track and first plane in track
    //TVector3 evtVtx(vShower->Start()[0], vShower->Start()[1], vShower->Start()[2]);
    double vtxDoca = PointDoca(evtVtx, vShower->Start(), vShower->Stop());
    double vtxToStart = (evtVtx - vShower->Start()).Mag();
    double vtxIntToStart = sqrt(vtxToStart * vtxToStart - vtxDoca * vtxDoca);
    //(Jianming's comments) record the gap from vertex to the start of the track
    //if the vertex is inside a track record distance of closest approach
    //TODO: For fuzzyK there can be an internal gap, where the start of the track is misleading
    //Ruth has looked into this, will try to incorporate internal gaps into this algorithm
    double shGap = vtxIntToStart;
    if (vtxIntToStart >= 0 && vtxIntToStart < vtxDoca)shGap = vtxDoca;
    if (vtxIntToStart < 0)shGap = vtxDoca;
    return shGap;
  }

int slid::ParticleIDAlg::GetMIPPlanes	(	const rb::Shower *	shower,
		const std::vector< const rb::Shower * >	showercol,
		const art::Event &	evt
	)

Number of MIP planes within the first 20 planes of the start of the shower.

Definition at line 1141 of file ParticleIDAlg.cxx.

References CalcPlaneLongitudinalDedx(), confusionMatrixTree::count, fMipRangeHigh, fMipRangeLow, fPlaneHits, and SetShower().

Referenced by slid::LIDBuilder::produce().

  {
    if(shower == NULL){
      return 0;
    }
    SetShower( shower, showercol, evt);
    
    int nMipPlanes = 0;
    int nPlanes = 20;
    if( fPlaneHits.size() < 20 ) 
      nPlanes = fPlaneHits.size();
    if(nPlanes == 0)
      return 0;
   
    int count = 0;

    for( auto const & iPlaneHit : fPlaneHits ){
      if(count > nPlanes)
        break;
      ++count;

      double dedx = CalcPlaneLongitudinalDedx( iPlaneHit.first );

      if( dedx > fMipRangeHigh || dedx < fMipRangeLow )
        continue;

      nMipPlanes++;
    }
    return nMipPlanes;
  }// end of GetMIPPlanes

TVector3 slid::ParticleIDAlg::InertiaXYZ ( ) const

inline

Return the inertia calculation.

Definition at line 134 of file ParticleIDAlg.h.

References febshutoff_auto::start.

{return fIxyz;}

bool slid::ParticleIDAlg::LoadDedxHistogramFiles ( const fhicl::ParameterSet &  pset, novadaq::cnv::DetId  geom  )

private

Load all of the histograms from files into arrays of TH1D* objects. The details of handling the histograms are taken care of within object of the helper class slid::DedxDistribution.

FD dE/dx input histogram files

NDOS dE/dx input histograms

ND dE/dx input histograms

Definition at line 306 of file ParticleIDAlg.cxx.

References om::cout, allTimeWatchdog::endl, util::EnvExpansion(), fDedxDistributions, fhicl::ParameterSet::get(), slid::kELECTRON, slid::kELECTRONCCCOH, slid::kELECTRONCCDIS, slid::kELECTRONCCQE, slid::kELECTRONCCRES, cmf::kFARDET, slid::kMUON, novadaq::cnv::kNDOS, cmf::kNEARDET, slid::kNEUTRON, slid::kPHOTON, slid::kPI0, slid::kPION, slid::kPROTON, and string.

Referenced by ParticleIDAlg().

  {

    std::string libPath = util::EnvExpansion(pset.get< std::string >("LibraryPath"));
    std::string thistname = pset.get< std::string >("HistFDDedxElectron");
    std::cout << "Loading LID library file from " << libPath  << std::endl;
    if (detectorID == novadaq::cnv::DetId::kFARDET) {

      /// FD dE/dx input histogram files
      std::cout<<"---- dE/dx histogram file for electron\n";
      fDedxDistributions[slid::kELECTRON].Initialize(pset, "HistFDDedxElectron");
      std::cout<<"---- dE/dx histogram file for photon\n";
      fDedxDistributions[slid::kPHOTON].Initialize(pset, "HistFDDedxPhoton");
      std::cout<<"---- dE/dx histogram file for muon\n";
      fDedxDistributions[slid::kMUON].Initialize(pset, "HistFDDedxMuon");
      std::cout<<"---- dE/dx histogram file for pi0\n";
      fDedxDistributions[slid::kPI0].Initialize(pset, "HistFDDedxPi0");
      std::cout<<"---- dE/dx histogram file for proton\n";
      fDedxDistributions[slid::kPROTON].Initialize(pset, "HistFDDedxProton");
      std::cout<<"---- dE/dx histogram file for neutron\n";
      fDedxDistributions[slid::kNEUTRON].Initialize(pset, "HistFDDedxNeutron");
      std::cout<<"---- dE/dx histogram file for charged pion\n";
      fDedxDistributions[slid::kPION].Initialize(pset, "HistFDDedxPion");
      std::cout<<"---- dE/dx histogram file for nue CCQE\n";
      fDedxDistributions[slid::kELECTRONCCQE].Initialize(pset, "HistFDDedxElectronCCQE");
      std::cout<<"---- dE/dx histogram file for nue CCRes\n";
      fDedxDistributions[slid::kELECTRONCCRES].Initialize(pset, "HistFDDedxElectronCCRES");
      std::cout<<"---- dE/dx histogram file for nue CC DIS\n";
      fDedxDistributions[slid::kELECTRONCCDIS].Initialize(pset, "HistFDDedxElectronCCDIS");
      std::cout<<"---- dE/dx histogram file for nue Coherent\n";
      fDedxDistributions[slid::kELECTRONCCCOH].Initialize(pset, "HistFDDedxElectronCCCOH");

    } else if (detectorID == novadaq::cnv::DetId::kNDOS) {

      /// NDOS dE/dx input histograms
      
      std::cout<<"---- dE/dx histogram file for electron\n";
      fDedxDistributions[slid::kELECTRON].Initialize(pset, "HistNDOSDedxElectron");
      std::cout<<"---- dE/dx histogram file for photon\n";
      fDedxDistributions[slid::kPHOTON].Initialize(pset, "HistNDOSDedxPhoton");
      std::cout<<"---- dE/dx histogram file for muon\n";
      fDedxDistributions[slid::kMUON].Initialize(pset, "HistNDOSDedxMuon");
      std::cout<<"---- dE/dx histogram file for pi0\n";
      fDedxDistributions[slid::kPI0].Initialize(pset, "HistNDOSDedxPi0");
      std::cout<<"---- dE/dx histogram file for proton\n";
      fDedxDistributions[slid::kPROTON].Initialize(pset, "HistNDOSDedxProton");
      std::cout<<"---- dE/dx histogram file for neutron\n";
      fDedxDistributions[slid::kNEUTRON].Initialize(pset, "HistNDOSDedxNeutron");
      std::cout<<"---- dE/dx histogram file for charged pion\n";
      fDedxDistributions[slid::kPION].Initialize(pset, "HistNDOSDedxPion");
      std::cout<<"---- dE/dx histogram file for nue CCQE\n";
      fDedxDistributions[slid::kELECTRONCCQE].Initialize(pset, "HistNDOSDedxElectronCCQE");
      std::cout<<"---- dE/dx histogram file for nue CCRes\n";
      fDedxDistributions[slid::kELECTRONCCRES].Initialize(pset, "HistNDOSDedxElectronCCRES");
      std::cout<<"---- dE/dx histogram file for nue CC DIS\n";
      fDedxDistributions[slid::kELECTRONCCDIS].Initialize(pset, "HistNDOSDedxElectronCCDIS");
      std::cout<<"---- dE/dx histogram file for nue Coherent\n";
      fDedxDistributions[slid::kELECTRONCCCOH].Initialize(pset, "HistNDOSDedxElectronCCCOH");

    } else if (detectorID == novadaq::cnv::DetId::kNEARDET) {

      /// ND dE/dx input histograms
      
      std::cout<<"---- dE/dx histogram file for electron\n";
      fDedxDistributions[slid::kELECTRON].Initialize(pset, "HistNDDedxElectron");
      std::cout<<"---- dE/dx histogram file for photon\n";
      fDedxDistributions[slid::kPHOTON].Initialize(pset, "HistNDDedxPhoton");
      std::cout<<"---- dE/dx histogram file for muon\n";
      fDedxDistributions[slid::kMUON].Initialize(pset, "HistNDDedxMuon");
      std::cout<<"---- dE/dx histogram file for pi0\n";
      fDedxDistributions[slid::kPI0].Initialize(pset, "HistNDDedxPi0");
      std::cout<<"---- dE/dx histogram file for proton\n";
      fDedxDistributions[slid::kPROTON].Initialize(pset, "HistNDDedxProton");
      std::cout<<"---- dE/dx histogram file for neutron\n";
      fDedxDistributions[slid::kNEUTRON].Initialize(pset, "HistNDDedxNeutron");
      std::cout<<"---- dE/dx histogram file for charged pion\n";
      fDedxDistributions[slid::kPION].Initialize(pset, "HistNDDedxPion");
      std::cout<<"---- dE/dx histogram file for nue CCQE\n";
      fDedxDistributions[slid::kELECTRONCCQE].Initialize(pset, "HistNDDedxElectronCCQE");
      std::cout<<"---- dE/dx histogram file for nue CC Res\n";
      fDedxDistributions[slid::kELECTRONCCRES].Initialize(pset, "HistNDDedxElectronCCRES");
      std::cout<<"---- dE/dx histogram file for nue CC DIS\n";
      fDedxDistributions[slid::kELECTRONCCDIS].Initialize(pset, "HistNDDedxElectronCCDIS");
      std::cout<<"---- dE/dx histogram file for nue Coherent\n";
      fDedxDistributions[slid::kELECTRONCCCOH].Initialize(pset, "HistNDDedxElectronCCCOH");
    } else {
      throw cet::exception("DetectorTypeNotDefined")
        << "Detector type " << detectorID
        << " not recognized.  See LIDBuilder_module.cc for list of defined algorithms.  "
        << "Only novadaq::cnv::DetId::kNEARDET, kFARDET, kNDOS defined.";

    }
    return true;
  }

double slid::ParticleIDAlg::Pi0Mass	(	const std::vector< const rb::Shower * >	showerCollection,
		unsigned int	iShower,
		const art::Event &	evt
	)

Calculate the mass of candidate and any other shower in the slice that is closest to the pi0 mass.

Definition at line 799 of file ParticleIDAlg.cxx.

References e, evt, stan::math::fabs(), fClosestPi0Mass, fNuEEnergyAlg, fPi0Index, std::isnan(), calib::j, SetShower(), and slid::NuEEnergyAlg::ShowerEnergy().

Referenced by Pi0SecondPhoton(), and slid::LIDBuilder::produce().

  {
    // Load shower if this is a new one.
    fClosestPi0Mass = -999.;
    SetShower(showerCollection.at(iShower), showerCollection, evt);
    
    double thisShowerERaw = fNuEEnergyAlg->ShowerEnergy(showerCollection[iShower], showerCollection, evt);
    if ((thisShowerERaw < 1e-10) || std::isnan(thisShowerERaw)) thisShowerERaw = 0.0;
    TLorentzVector thisShower4Vec(showerCollection.at(iShower)->Dir()[0] * thisShowerERaw,
                                  showerCollection.at(iShower)->Dir()[1] * thisShowerERaw,
                                  showerCollection.at(iShower)->Dir()[2] * thisShowerERaw, thisShowerERaw);
    
    for (unsigned int j = 0; j < showerCollection.size(); j++) { // pi0 mass calculation
      if (j == iShower)continue;
      // NOT NECESSARY?                if (showerCollection.at(j).ISlice() != showerCollection.at(iShower).ISlice())continue;
      double showereraw = fNuEEnergyAlg->ShowerEnergy(showerCollection[j], showerCollection, evt);
      if ((showereraw < 1e-10) || std::isnan(showereraw)) showereraw = 0.0;
      TLorentzVector otherShower4Vec(showerCollection.at(j)->Dir()[0] * showereraw,
                                     showerCollection.at(j)->Dir()[1] * showereraw,
                                     showerCollection.at(j)->Dir()[2] * showereraw, showereraw);
      TLorentzVector p2g = thisShower4Vec + otherShower4Vec;
      double mgg = p2g.M();
      if (fabs(fClosestPi0Mass - 0.15) > fabs(mgg - 0.15)) {
        fClosestPi0Mass = mgg;
        fPi0Index = j;
      }
    }
    if (fClosestPi0Mass < 0)fClosestPi0Mass = 0.;
    return fClosestPi0Mass;
  }

int slid::ParticleIDAlg::Pi0SecondPhoton	(	const std::vector< const rb::Shower * >	showerCollection,
		unsigned int	iShower,
		const art::Event &	evt
	)

Index in showerCollection of the shower for which the reconstructed invariant mass was closest to pi0 mass for the shower at index iShower.

Definition at line 835 of file ParticleIDAlg.cxx.

References fPi0Index, fShower, and Pi0Mass().

Referenced by slid::LIDBuilder::produce().

  {
    if( fShower == NULL)
      Pi0Mass( showerCollection, iShower, evt);
    
    return fPi0Index;
          
  }// end of Pi0SecondPhoton

double slid::ParticleIDAlg::PlaneCentroid

(

unsigned int

pIdx

)

Return the energy weighted mean shower positon for the plane in cm.

Definition at line 1024 of file ParticleIDAlg.cxx.

References fMapIdxToGlobalPlane, fPlaneCentroid, and it.

  {
    double centroid = 9999.0;
    std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
    if (itr != fMapIdxToGlobalPlane.end()){
      std::map<int, double>::iterator it = fPlaneCentroid.find(itr->second);
      if (it != fPlaneCentroid.end()) centroid = it->second;
    }

    return centroid;
  }

double slid::ParticleIDAlg::PlaneHitCell

(

unsigned int

pIdx

)

Return point of intersection cell between shower core and plane.

Definition at line 968 of file ParticleIDAlg.cxx.

References geo::PlaneGeo::Cell(), om::cout, rb::Prong::Dir(), allTimeWatchdog::endl, stan::math::fabs(), fMapIdxToGlobalPlane, fShower, geom(), geo::CellGeo::GetCenter(), geo::CellGeo::HalfW(), geo::kX, geo::kY, geo::PlaneGeo::Ncells(), geo::GeometryBase::Plane(), NDAPDHVSetting::plane, rb::Prong::Start(), rb::Shower::Stop(), POTSpillRate::view, and geo::PlaneGeo::View().

Referenced by CalcCellPlaneTransverseDedx().

  {
    double trkcellId = 9999.0;
    std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
    if (itr != fMapIdxToGlobalPlane.end()){
      unsigned int plane = itr->second;
      Double_t cellXYZ[3]={0,0,0};
      Double_t cellXYZ1[3];
      Double_t cellXYZmax[3];
      art::ServiceHandle<geo::Geometry> geom;
      geom->Plane(plane)->Cell(0)->GetCenter(cellXYZ);
      geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZ1);
      geom->Plane(plane)->Cell(geom->Plane(plane)->Ncells()-1)->GetCenter(cellXYZmax);
      double cellW = 2*geom->Plane(plane)->Cell(0)->HalfW();
      double cellwx = (cellXYZ1[0] - cellXYZ[0] + cellW)/geom->Plane(plane)->Ncells();
      double cellwy = (cellXYZ1[1] - cellXYZ[1] + cellW)/geom->Plane(plane)->Ncells();
      TVector3 trkdir=fShower->Dir();
      TVector3 trkstart=fShower->Start();
      TVector3 trkstop=fShower->Stop();
      double trkplanex=99999,trkplaney=99999;//, trkplanez=0;
  //    double trkplanex0=99999,trkplaney0=99999;
  //    double trkplanex1=99999,trkplaney1=99999;
      if(fabs(trkdir.z())>0.0001){
        trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
        trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
  //      trkplanex1 = cellXYZmax[0];
  //      trkplaney1 = cellXYZmax[1];
      }else{
        trkplanex = (trkstart.x()+ trkstop.x())/2.0;
        trkplaney = (trkstart.y()+ trkstop.y())/2.0;
      }
  //    trkplanez = cellXYZ[2];
      geo::View_t view = geom->Plane(plane)->View();
  /*
      if(view == geo::kX){
       std::cout<<"x view plane"<<std::endl;
      std::cout<<"x1-x0,cellW "<<" "<<cellXYZ1[0] - cellXYZ[0]<<" "<<cellW<<std::endl;
      }else if(view == geo::kY){
       std::cout<<"y view plane"<<std::endl;
      std::cout<<"y1-y0,cellW "<<" "<<cellXYZ1[1] - cellXYZ[1]<<" "<<cellW<<std::endl;
      }
  */
  //    int trkcellId0 = -9999;
  //    int trkcellId1 = -9999;
  //    int trkcellIdmax = geom->Plane(plane)->Ncells()-1;
      if(view == geo::kX)trkcellId = (trkplanex+cellwx/2-cellXYZ[0])/cellwx;
      if(view == geo::kY)trkcellId = (trkplaney+cellwy/2-cellXYZ[1])/cellwy;
  //    if(view == geo::kX)trkcellId0 = int((trkplanex0+cellwx/2-cellXYZ[0])/cellwx);
  //    if(view == geo::kY)trkcellId0 = int((trkplaney0+cellwy/2-cellXYZ[1])/cellwy);
  //    if(view == geo::kX)trkcellId1 = int((trkplanex1+cellwx/2-cellXYZ[0])/cellwx);
  //    if(view == geo::kY)trkcellId1 = int((trkplaney1+cellwy/2-cellXYZ[1])/cellwy);
    }else{std::cout<<pIdx<<" end !!! "<<std::endl;}
    return trkcellId;
  }

TVector3 slid::ParticleIDAlg::PlaneHitXYZ

(

unsigned int

pIdx

)

Return point of intersection between shower core and plane.

Definition at line 938 of file ParticleIDAlg.cxx.

References geo::PlaneGeo::Cell(), rb::Prong::Dir(), stan::math::fabs(), fMapIdxToGlobalPlane, fShower, geom(), geo::CellGeo::GetCenter(), geo::GeometryBase::Plane(), rb::Prong::Start(), and rb::Shower::Stop().

Referenced by slid::LIDTraining::analyze(), and CalcAsymIneria().

  {
    TVector3 trkplanepos;
    std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
    if (itr != fMapIdxToGlobalPlane.end()){
      Double_t cellXYZ[3]={0,0,0};
      art::ServiceHandle<geo::Geometry> geom;
      geom->Plane(itr->second)->Cell(0)->GetCenter(cellXYZ);
      TVector3 trkdir=fShower->Dir();
      TVector3 trkstart=fShower->Start();
      TVector3 trkstop=fShower->Stop();
      double trkplanex=99999,trkplaney=99999, trkplanez=0;
      if(fabs(trkdir.z())>0.0001){
        trkplanex = trkstart.x()+(cellXYZ[2]-trkstart.z())*(trkdir.x()/trkdir.z());
        trkplaney = trkstart.y()+(cellXYZ[2]-trkstart.z())*(trkdir.y()/trkdir.z());
      }
      else{
        trkplanex = (trkstart.x()+trkstop.x())/2.0; 
        trkplaney = (trkstart.y()+trkstop.y())/2.0;
      }
      trkplanez = cellXYZ[2];
      trkplanepos.SetX(trkplanex);
      trkplanepos.SetY(trkplaney);
      trkplanepos.SetZ(trkplanez);
    }

    return trkplanepos;
  }

double slid::ParticleIDAlg::PlaneLongDedx

(

unsigned int

pIdx

)

return longitudinal dedx for a specified plane index

Definition at line 870 of file ParticleIDAlg.cxx.

References CalcPlaneLongitudinalDedx(), om::cout, allTimeWatchdog::endl, and fMapIdxToGlobalPlane.

Referenced by slid::LIDTraining::analyze(), slid::LIDBuilder::produce(), and slid::SPIDBuilder::produce().

  {
    std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
    double dedx = -1.0;
    if (itr != fMapIdxToGlobalPlane.end()) dedx = CalcPlaneLongitudinalDedx(itr->second);
    else  std::cout<<"Requested longitudinal plane index not found, returning dedx=-1"<<std::endl;

    return dedx;
  }

double slid::ParticleIDAlg::PlaneRadius

(

unsigned int

pIdx

)

return shower radius for a plane index

Definition at line 901 of file ParticleIDAlg.cxx.

References rb::CellHit::Cell(), rb::Cluster::Cell(), slid::NuEEnergyAlg::CellEnergy(), geo::ClosestApproach(), om::cout, dir, rb::Prong::Dir(), dist, allTimeWatchdog::endl, eventID, fGeom, fMapIdxToGlobalPlane, fNuEEnergyAlg, fPlaneHits, fShower, fShowerCol, it, geo::kX, rb::CellHit::Plane(), geo::GeometryBase::Plane(), Munits::rad, std::sqrt(), febshutoff_auto::start, rb::Prong::Start(), and rb::CellHit::View().

Referenced by slid::LIDTraining::analyze().

  {
    std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
    double rad = -1.0;
    double gev = 0;
    double dist, di;
    if (itr != fMapIdxToGlobalPlane.end()){
      std::map<int, int> planecells = fPlaneHits[itr->second];
      rad = 0.0;
      for (std::map<int, int>::iterator it=planecells.begin(); it!= planecells.end(); ++it) {
        art::Ptr<rb::CellHit>  cHit = fShower->Cell(it->second);
        double cellEnergy = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, it->second, eventID);
        TVector3 start  = fShower->Start();
        TVector3 dir    = fShower->Dir();
        TVector3 xyz;
        fGeom->Plane(cHit->Plane())->Cell(cHit->Cell())->GetCenter(xyz);
        const TVector3 dxyz = (cHit->View() == geo::kX) ? TVector3(0, 1, 0): TVector3(1, 0, 0);
        dist = sqrt(geo::ClosestApproach(start, start+dir, xyz, xyz+dxyz, &di));
        rad += dist*cellEnergy;
        gev += cellEnergy;
      }
    }
    else  std::cout<<"Requested longitudinal plane index not found, returning dedx=-1"<<std::endl;

    if(gev>0.001)rad = rad/gev;
    return rad;
  }

unsigned int slid::ParticleIDAlg::PlaneToGlobalPlane

(

unsigned int

pIdx

)

Return the gloabel plane index for an index in the shower coordinate.

Definition at line 930 of file ParticleIDAlg.cxx.

References fMapIdxToGlobalPlane.

Referenced by slid::LIDTraining::analyze().

  {
    std::map<unsigned int, unsigned int>::iterator itr = fMapIdxToGlobalPlane.find(pIdx);
    if(itr != fMapIdxToGlobalPlane.end()) return itr->second;
    return 99999;
  }

double slid::ParticleIDAlg::PlaneTransDedx

(

unsigned int

tpIdx

)

return transverse dedx for specified transverse plane

Definition at line 881 of file ParticleIDAlg.cxx.

References CalcPlaneTransverseDedx(), om::cout, allTimeWatchdog::endl, and slid::kNumTransversePlane.

Referenced by slid::LIDTraining::analyze().

  {
    double dedx = -1.0;
    if (tpIdx < slid::kNumTransversePlane) dedx = CalcPlaneTransverseDedx(tpIdx);
    else std::cout<<"Requested transverse plane index exceeds valid range, maximum plane is: "<<slid::kNumTransversePlane<<", returning dedx=-1"<<std::endl;

    return dedx;
  }

double slid::ParticleIDAlg::PointDoca	(	TVector3	vtx,
		TVector3	start,
		TVector3	stop
	)

Get distance of closest approach between shower and vertex.

Definition at line 786 of file ParticleIDAlg.cxx.

References d, Mag(), and febshutoff_auto::start.

Referenced by GetGapVertexToShowerStart(), slid::LIDBuilder::produce(), slid::SPIDBuilder::produce(), and cosrej::MakeNueCosRej::produce().

  {
    double d = 9999;
    double denorm = (start - stop).Mag();
    if (denorm < 0.00001) {
      d = (vtx - start).Mag();
    } else {
      d = ((vtx - start).Cross(vtx - stop)).Mag() / denorm;
    }
    return d;
  }

double slid::ParticleIDAlg::Radius	(	const std::vector< const rb::Shower * >	showerCollection,
		unsigned int	iShower,
		const art::Event &	evt
	)

Calculate the shower radius.

Definition at line 846 of file ParticleIDAlg.cxx.

References rb::CellHit::Cell(), rb::Cluster::Cell(), slid::NuEEnergyAlg::CellEnergy(), geo::ClosestApproach(), dir, rb::Prong::Dir(), dist, eventID, evt, fGeom, fNuEEnergyAlg, fShower, fShowerCol, geo::kX, rb::Cluster::NCell(), rb::CellHit::Plane(), geo::GeometryBase::Plane(), Munits::rad, SetShower(), std::sqrt(), febshutoff_auto::start, rb::Prong::Start(), and rb::CellHit::View().

Referenced by slid::LIDBuilder::produce().

  {
    SetShower(showerCollection.at(iShower), showerCollection, evt);
    double rad = 0;
    double gev = 0;
    double dist, di;
    for (unsigned int nh=0;nh<fShower->NCell();nh++) {
      art::Ptr<rb::CellHit>  cHit = fShower->Cell(nh);
      double cellEnergy = fNuEEnergyAlg->CellEnergy(fShower, fShowerCol, nh, eventID);;
      TVector3 start  = fShower->Start();
      TVector3 dir    = fShower->Dir();
      TVector3 xyz;
      fGeom->Plane(cHit->Plane())->Cell(cHit->Cell())->GetCenter(xyz);
      const TVector3 dxyz = (cHit->View() == geo::kX) ? TVector3(0, 1, 0): TVector3(1, 0, 0);
      dist = sqrt(geo::ClosestApproach(start, start+dir, xyz, xyz+dxyz, &di));
      rad += dist*cellEnergy;
      gev += cellEnergy;
    }
    if(gev>0.001)rad = rad/gev;
    return rad;
  }

void slid::ParticleIDAlg::reconfigure	(	const fhicl::ParameterSet &	pset,
		novadaq::cnv::DetId	detId
	)

Reinitialize parameters read from FCL file.

Definition at line 54 of file ParticleIDAlg.cxx.

  {
  }

void slid::ParticleIDAlg::SetShower	(	const rb::Shower *	vShower,
		const std::vector< const rb::Shower * >	showercol,
		const art::Event &	evt
	)

Set rb::Prong to be analyzed. This must be set before any calculations can be done.

Set rb::Shower object to be analyzed. Calculate its total energy, which will be needed for calculations.

Parameters

vShower
evt

Definition at line 169 of file ParticleIDAlg.cxx.

References CalcAsymIneria(), CalcDetectorXYRegionIndex(), CalcPlaneCentroid(), CalcPlaneHits(), CalculateEnergyBinIndex(), e, gov::fnal::cd::rms::provider::equal(), eventID, fAsym, fDedxLongByPlane, fDedxTransByPlane, fIxyz, fMapIdxToGlobalPlane, fNuEEnergyAlg, fPartLongLL, fPartTransLL, fPi0Index, fPlaneCentroid, fPlaneHits, fShower, fShowerCol, fShowerEnergyGev, fTransEnergy, fTransPathLength, art::Event::id(), rb::Cluster::ID(), std::isnan(), and slid::NuEEnergyAlg::ShowerEnergy().

Referenced by slid::LIDTraining::analyze(), DedxLongLL(), DedxTransLL(), GetMIPPlanes(), Pi0Mass(), slid::LIDBuilder::produce(), slid::SPIDBuilder::produce(), and Radius().

  {
    //
    //  Only recalculate shower quantities if a different event or shower is presented or if 
    //  fShower is not already set
    //

    if ( fShower  != NULL && 
         evt.id() == eventID &&
         fShowerCol.size() == showercol.size() &&
         std::equal( fShowerCol.begin(), fShowerCol.end(), showercol.begin()) &&
         fShower->ID() == vShower->ID()) 
      return;

    //  Save to compare next one
    fShower    = vShower;
    fShowerCol = showercol; 
    eventID    = evt.id();

    //clear map of plane index to global plane
    fMapIdxToGlobalPlane.clear();

    //  Clear cache of longitudinal dedx values by longitudinal and tranverse plane
    fDedxLongByPlane.clear();
    fDedxTransByPlane.clear();

    fPartLongLL.clear();
    fPartTransLL.clear();

    fTransEnergy.clear();
    fTransPathLength.clear();
    fPi0Index = -5;

    fIxyz.SetXYZ(0.0,0.0,0.0);
    fAsym = 1.0;

    //  Calculate and save the total energy of the shower in a member field
    fShowerEnergyGev = fNuEEnergyAlg->ShowerEnergy(fShower, showercol, evt);
    if ((fShowerEnergyGev<1e-10)|| std::isnan(fShowerEnergyGev)) fShowerEnergyGev = 0.0;


    // Fill a  map where the first element is
    // plane number and the second is vector of global
    // cell indices in a shower of cells that belong to
    // that plane.
    fPlaneHits = CalcPlaneHits( vShower);

    // Fill a map where the first index is the
    // number of planes away from the track core, and the
    // second is a vector of global cell indices in a
    // shower that are that many planes from the track
    // core. It collapses the transverse profile of the
    // shower to a single plane.
    fPlaneCentroid = CalcPlaneCentroid( vShower );

    //
    // Determine detector XY region index
    //
    CalcDetectorXYRegionIndex();
    //
    //  Determine energy bins
    //
    CalculateEnergyBinIndex();

    //calculate shower asymmetry and ineria
    CalcAsymIneria();

  }

double slid::ParticleIDAlg::ShowerAsym ( ) const

inline

Return the asymmetry of the shower.

Definition at line 131 of file ParticleIDAlg.h.

{return fAsym;}

Member Data Documentation

double slid::ParticleIDAlg::e0

private

Energy values of energy bins adjacent to shower. Calculated once and cached for a given shower/event combo.

Definition at line 280 of file ParticleIDAlg.h.

Referenced by CalcInterCellTransDedxProb(), CalcInterPlaneDedxProb(), and CalculateEnergyBinIndex().

double slid::ParticleIDAlg::e1

private

Definition at line 281 of file ParticleIDAlg.h.

Referenced by CalcInterCellTransDedxProb(), CalcInterPlaneDedxProb(), and CalculateEnergyBinIndex().

art::EventID slid::ParticleIDAlg::eventID

private

data/MC flag

event ID

Definition at line 210 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), CalcCellPlaneTransverseDedx(), CalcEnergyInLongitudinalPlane(), CalcPlaneCentroid(), CalcPlaneTransverseDedx(), PlaneRadius(), Radius(), and SetShower().

double slid::ParticleIDAlg::fAsym

private

shower asymmetry

Definition at line 286 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), and SetShower().

double slid::ParticleIDAlg::fClosestPi0Mass

private

Closest mass to pi0 mass formed by this shower and any other shower in the slice.

Definition at line 213 of file ParticleIDAlg.h.

Referenced by Pi0Mass().

std::vector<slid::DedxDistribution> slid::ParticleIDAlg::fDedxDistributions

private

dE/dx distributions (histogram files and data) for each detector type and particle category

Definition at line 216 of file ParticleIDAlg.h.

Referenced by CalcPlaneLongDedxProb(), CalcPlaneTransDedxProb(), and LoadDedxHistogramFiles().

std::map< int, double> slid::ParticleIDAlg::fDedxLongByPlane

private

Map of longitudinal dedx by plane. The first element is plane number and the second is the dedx for that plane. This is used for lazy caching of the dedx plane values in order to save recalculating for the same shower.

Definition at line 236 of file ParticleIDAlg.h.

Referenced by CalcPlaneLongitudinalDedx(), and SetShower().

std::map< int, double> slid::ParticleIDAlg::fDedxTransByPlane

private

Map of transverse dedx by plane. The first element is plane number and the second is the dedx for that plane. This is used for lazy caching of the dedx plane values in order to save recalculating for the same shower.

Definition at line 242 of file ParticleIDAlg.h.

Referenced by CalcPlaneTransverseDedx(), and SetShower().

novadaq::cnv::DetId slid::ParticleIDAlg::fDetId

private

Detector ID.

Definition at line 295 of file ParticleIDAlg.h.

Referenced by CalcDetectorXYRegionIndex().

art::ServiceHandle<geo::Geometry> slid::ParticleIDAlg::fGeom

private

geometry

Definition at line 226 of file ParticleIDAlg.h.

Referenced by CalcCellPlaneTransverseDedx(), CalcPlaneTransverseDedx(), CalcTrkHitPath(), PlaneRadius(), and Radius().

TVector3 slid::ParticleIDAlg::fIxyz

private

shower moment of ineria

Definition at line 284 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), and SetShower().

std::map<unsigned int, unsigned int> slid::ParticleIDAlg::fMapIdxToGlobalPlane

private

Map where the key is the index of the plane in the shower and the second is the global plane This map saves looping in some of the utility functions

Definition at line 253 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), CalcCellPlaneTransverseDedx(), CalcPlaneHits(), PlaneCentroid(), PlaneHitCell(), PlaneHitXYZ(), PlaneLongDedx(), PlaneRadius(), PlaneToGlobalPlane(), and SetShower().

double slid::ParticleIDAlg::fMipRangeHigh

private

Upper bound of MIP Range.

Definition at line 289 of file ParticleIDAlg.h.

Referenced by GetMIPPlanes().

double slid::ParticleIDAlg::fMipRangeLow

private

Lower bound of MIP Range.

Definition at line 292 of file ParticleIDAlg.h.

Referenced by GetMIPPlanes().

NuEEnergyAlg* slid::ParticleIDAlg::fNuEEnergyAlg

private

Nue energy algorithms.

Definition at line 229 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), CalcCellPlaneTransverseDedx(), CalcEnergyInLongitudinalPlane(), CalcPlaneCentroid(), CalcPlaneTransverseDedx(), Pi0Mass(), PlaneRadius(), Radius(), and SetShower().

std::map<int, float> slid::ParticleIDAlg::fPartLongLL

Map of the longitudinal ll by paricle type.

Definition at line 145 of file ParticleIDAlg.h.

Referenced by DedxLongLL(), slid::LIDBuilder::produce(), slid::SPIDBuilder::produce(), and SetShower().

std::map<int, float> slid::ParticleIDAlg::fPartTransLL

Map of the transverse ll by particle type.

Definition at line 148 of file ParticleIDAlg.h.

Referenced by DedxTransLL(), slid::LIDBuilder::produce(), slid::SPIDBuilder::produce(), and SetShower().

int slid::ParticleIDAlg::fPi0Index

private

Index of the other photon shower used for pi0 mass calculation for a given shower

Definition at line 267 of file ParticleIDAlg.h.

Referenced by Pi0Mass(), Pi0SecondPhoton(), and SetShower().

std::map< int, double > slid::ParticleIDAlg::fPlaneCentroid

private

Map where the first index is plane number and the second is the energy weighted center position for that plane

Definition at line 257 of file ParticleIDAlg.h.

Referenced by CalcCellPlaneTransverseDedx(), CalcPlaneTransverseDedx(), PlaneCentroid(), and SetShower().

std::map< int, std::map< int, int > > slid::ParticleIDAlg::fPlaneHits

private

Map where the first element is plane number and the second is a map with the key being the cell number in the plane and the second value is the global cell index in a shower

Definition at line 249 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), CalcCellPlaneTransverseDedx(), CalcEnergyInLongitudinalPlane(), CalcPlaneCentroid(), CalcPlaneTransverseDedx(), CalcTrkHitPath(), DedxLongLL(), GetMIPPlanes(), PlaneRadius(), and SetShower().

const rb::Shower* slid::ParticleIDAlg::fShower

private

rb::Shower object for which particle ID is to be calculated. This is cached.

Definition at line 219 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), CalcCellPlaneTransverseDedx(), CalcDetectorXYRegionIndex(), CalcEnergyInLongitudinalPlane(), CalcInterPlaneDedxProb(), CalcPlaneTransverseDedx(), CalcTrkHitPath(), Pi0SecondPhoton(), PlaneHitCell(), PlaneHitXYZ(), PlaneRadius(), Radius(), and SetShower().

std::vector< const rb::Shower* > slid::ParticleIDAlg::fShowerCol

private

Definition at line 220 of file ParticleIDAlg.h.

Referenced by CalcAsymIneria(), CalcCellPlaneTransverseDedx(), CalcEnergyInLongitudinalPlane(), CalcPlaneCentroid(), CalcPlaneTransverseDedx(), PlaneRadius(), Radius(), and SetShower().

double slid::ParticleIDAlg::fShowerEnergyGev

private

Shower energy.

Definition at line 223 of file ParticleIDAlg.h.

Referenced by CalculateEnergyBinIndex(), and SetShower().

std::map<int, double > slid::ParticleIDAlg::fTransEnergy

private

Map of transverse energy for each transverse plane.

Definition at line 260 of file ParticleIDAlg.h.

Referenced by CalcPlaneTransverseDedx(), and SetShower().

std::map<int, double > slid::ParticleIDAlg::fTransPathLength

private

Map of path length for each transverse plane.

Definition at line 263 of file ParticleIDAlg.h.

Referenced by CalcPlaneTransverseDedx(), and SetShower().

int slid::ParticleIDAlg::igev0

private

Shower adjacent energy bins. Calculated once and cached for a given shower/event combo.

Definition at line 276 of file ParticleIDAlg.h.

Referenced by CalcInterCellTransDedxProb(), CalcInterPlaneDedxProb(), and CalculateEnergyBinIndex().

int slid::ParticleIDAlg::igev1

private

Definition at line 277 of file ParticleIDAlg.h.

Referenced by CalcInterCellTransDedxProb(), CalcInterPlaneDedxProb(), and CalculateEnergyBinIndex().

int slid::ParticleIDAlg::iReg

private

Detector sub-region bin. Calculated once and cached for a given shower/event combo.

Definition at line 270 of file ParticleIDAlg.h.

Referenced by CalcDetectorXYRegionIndex(), CalcPlaneLongDedxProb(), and CalcPlaneTransDedxProb().

TVector3 slid::ParticleIDAlg::pos

private

Shower start position.

Definition at line 273 of file ParticleIDAlg.h.

Referenced by CalcDetectorXYRegionIndex().

---

The documentation for this class was generated from the following files:

• /cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-01-16/ShowerLID/ParticleIDAlg.h
• /cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N21-01-16/ShowerLID/ParticleIDAlg.cxx