#include <EgammaPCAHelper.h>

Public Types
typedef ROOT::Math::Transform3D::Point	Point

typedef ROOT::Math::Transform3D	Transform3D

Public Member Functions
const math::XYZVector &	axis () const

const math::XYZPoint &	barycenter () const

void	clear ()

float	clusterDepthCompatibility (const LongDeps &, float &measuredDepth, float &expectedDepth, float &expectedSigma)

void	computePCA (float radius, bool withHalo=true)

void	computeShowerWidth (float radius, bool withHalo=true)

	EGammaPCAHelper ()

const TVectorD &	eigenValues () const

LongDeps	energyPerLayer (float radius, bool withHalo=true)

const std::unordered_map < DetId, const HGCRecHit * > *	getHitMap ()

void	pcaInitialComputation ()

const TPrincipal &	pcaResult ()

void	printHits (float radius) const

void	setdEdXWeights (const std::vector< double > &dEdX)

void	setHitMap (const std::unordered_map< DetId, const HGCRecHit * > *hitMap)
	to set once per event More...

void	setRecHitTools (const hgcal::RecHitTools *recHitTools)

double	sigmaEE () const

double	sigmaPP () const

const TVectorD &	sigmas () const

double	sigmaUU () const

double	sigmaVV () const

void	storeRecHits (const reco::CaloCluster &theCluster)

void	storeRecHits (const reco::HGCalMultiCluster &cluster)

Private Member Functions
bool	checkIteration () const

float	findZFirstLayer (const LongDeps &) const

void	storeRecHits (const std::vector< std::pair< DetId, float >> &hf)

Private Attributes
math::XYZVector	axis_

math::XYZPoint	barycenter_

bool	debug_

std::vector< double >	dEdXWeights_

const std::unordered_map < DetId, const HGCRecHit * > *	hitMap_

std::vector< double >	invThicknessCorrection_

unsigned int	maxlayer_

std::unique_ptr< TPrincipal >	pca_

int	pcaIteration_

bool	recHitsStored_

const hgcal::RecHitTools *	recHitTools_

ShowerDepth	showerDepth_

double	sige_

double	sigp_

double	sigu_

double	sigv_

const reco::CaloCluster *	theCluster_

std::vector< Spot >	theSpots_

Transform3D	trans_

Detailed Description

Definition at line 33 of file EgammaPCAHelper.h.

Member Typedef Documentation

typedef ROOT::Math::Transform3D::Point hgcal::EGammaPCAHelper::Point

Definition at line 36 of file EgammaPCAHelper.h.

typedef ROOT::Math::Transform3D hgcal::EGammaPCAHelper::Transform3D

Definition at line 35 of file EgammaPCAHelper.h.

Constructor & Destructor Documentation

EGammaPCAHelper::EGammaPCAHelper ( )

Definition at line 17 of file EgammaPCAHelper.cc.

     :  // Thickness correction to dEdx weights
       // (100um, 200um, 300um silicon)
       // See RecoLocalCalo.HGCalRecProducers.HGCalRecHit_cfi
       invThicknessCorrection_({1. / 1.132, 1. / 1.092, 1. / 1.084}),

Member Function Documentation

const math::XYZVector& hgcal::EGammaPCAHelper::axis ( ) const

inline

Definition at line 57 of file EgammaPCAHelper.h.

References axis_.

Referenced by HGCalEgammaIDHelper::axis().

57 { return axis_; }

hgcal::EGammaPCAHelper::axis_

math::XYZVector axis_

Definition: EgammaPCAHelper.h:95

const math::XYZPoint& hgcal::EGammaPCAHelper::barycenter ( ) const

inline

Definition at line 56 of file EgammaPCAHelper.h.

References barycenter_.

Referenced by HGCalEgammaIDHelper::barycenter().

56 { return barycenter_; }

hgcal::EGammaPCAHelper::barycenter_

math::XYZPoint barycenter_

Definition: EgammaPCAHelper.h:94

bool EGammaPCAHelper::checkIteration ( ) const

private

Definition at line 215 of file EgammaPCAHelper.cc.

References gather_cfg::cout, debug_, and pcaIteration_.

Referenced by clusterDepthCompatibility(), energyPerLayer(), sigmaEE(), sigmaPP(), sigmaUU(), and sigmaVV().

                                            {
   if (pcaIteration_ == 0) {
     if (debug_)
       std::cout << " The PCA has not been run yet " << std::endl;
     return false;
   } else if (pcaIteration_ == 1) {
     if (debug_)
       std::cout << " The PCA has been run only once - careful " << std::endl;
     return false;
   } else if (pcaIteration_ == -1) {
     if (debug_)
       std::cout << " Not enough layers to perform PCA " << std::endl;
     return false;
   }
   return true;
 }

void EGammaPCAHelper::clear ( void )

Definition at line 232 of file EgammaPCAHelper.cc.

References pcaIteration_, sige_, sigp_, sigu_, sigv_, and theSpots_.

Referenced by HGCalEgammaIDHelper::computeHGCAL().

                             {
   theSpots_.clear();
   pcaIteration_ = 0;
   sigu_ = 0.;
   sigv_ = 0.;
   sigp_ = 0.;
   sige_ = 0.;
 }

float EGammaPCAHelper::clusterDepthCompatibility	(	const LongDeps &	ld,
		float &	measuredDepth,
		float &	expectedDepth,
		float &	expectedSigma
	)

Definition at line 307 of file EgammaPCAHelper.cc.

References funct::abs(), axis_, barycenter_, checkIteration(), DeadROC_duringRun::dir, hgcal::LongDeps::energyEE(), findZFirstLayer(), hgcal::ShowerDepth::getClusterDepthCompatibility(), and showerDepth_.

Referenced by HGCalEgammaIDHelper::clusterDepthCompatibility().

                                                                        {
   expectedDepth = -999.;
   expectedSigma = -999.;
   measuredDepth = -999.;
   if (!checkIteration())
     return -999.;
 
   float z = findZFirstLayer(ld);
   math::XYZVector dir = axis_.unit();
   measuredDepth = std::abs((z - std::abs(barycenter_.z())) / dir.z());
   return showerDepth_.getClusterDepthCompatibility(measuredDepth, ld.energyEE(), expectedDepth, expectedSigma);
 }

void EGammaPCAHelper::computePCA	(	float	radius,
		bool	withHalo = `true`
	)

Definition at line 106 of file EgammaPCAHelper.cc.

References axis_, barycenter_, gather_cfg::cout, debug_, mps_fire::i, hgcal::RecHitTools::lastLayerEE(), LayerTriplets::layers(), compareTotals::means, pca_, pcaIteration_, CosmicsPD_Skims::radius, recHitTools_, theSpots_, and trans_.

Referenced by HGCalEgammaIDHelper::computeHGCAL(), and pcaInitialComputation().

                                                             {
   // very important - to reset
   pca_ = std::make_unique<TPrincipal>(3, "D");
   bool initialCalculation = radius < 0;
   if (debug_)
     std::cout << " Initial calculation " << initialCalculation << std::endl;
   if (initialCalculation && withHalo) {
     edm::LogWarning("EGammaPCAHelper") << "Warning - in the first iteration, the halo hits are excluded " << std::endl;
     withHalo = false;
   }
 
   float radius2 = radius * radius;
   if (!initialCalculation) {
     math::XYZVector mainAxis(axis_);
     mainAxis.unit();
     math::XYZVector phiAxis(barycenter_.x(), barycenter_.y(), 0);
     math::XYZVector udir(mainAxis.Cross(phiAxis));
     udir = udir.unit();
     trans_ = Transform3D(Point(barycenter_),
                          Point(barycenter_ + axis_),
                          Point(barycenter_ + udir),
                          Point(0, 0, 0),
                          Point(0., 0., 1.),
                          Point(1., 0., 0.));
   }
 
   std::set<int> layers;
   for (const auto& spot : theSpots_) {
     if (spot.layer() > recHitTools_->lastLayerEE())
       continue;
     if (!withHalo && (!spot.isCore()))
       continue;
     if (initialCalculation) {
       // initial calculation, take only core hits
       if (!spot.isCore())
         continue;
       layers.insert(spot.layer());
       for (int i = 0; i < spot.multiplicity(); ++i)
         pca_->AddRow(spot.row());
     } else {
       // use a cylinder, include all hits
       math::XYZPoint local = trans_(Point(spot.row()[0], spot.row()[1], spot.row()[2]));
       if (local.Perp2() > radius2)
         continue;
       layers.insert(spot.layer());
       for (int i = 0; i < spot.multiplicity(); ++i)
         pca_->AddRow(spot.row());
     }
   }
   if (debug_)
     std::cout << " Nlayers " << layers.size() << std::endl;
   if (layers.size() < 3) {
     pcaIteration_ = -1;
     return;
   }
   pca_->MakePrincipals();
   ++pcaIteration_;
   const TVectorD& means = *(pca_->GetMeanValues());
   const TMatrixD& eigens = *(pca_->GetEigenVectors());
 
   barycenter_ = math::XYZPoint(means[0], means[1], means[2]);
   axis_ = math::XYZVector(eigens(0, 0), eigens(1, 0), eigens(2, 0));
   if (axis_.z() * barycenter_.z() < 0.0) {
     axis_ = -1. * axis_;
   }
 }

void EGammaPCAHelper::computeShowerWidth	(	float	radius,
		bool	withHalo = `true`
	)

Definition at line 173 of file EgammaPCAHelper.cc.

References srCondWrite_cfg::deltaPhi, reco::CaloCluster::eta(), reco::CaloCluster::phi(), CosmicsPD_Skims::radius, sige_, sigp_, sigu_, sigv_, mathSSE::sqrt(), theCluster_, theSpots_, and trans_.

Referenced by HGCalEgammaIDHelper::computeHGCAL().

                                                                     {
   sigu_ = 0.;
   sigv_ = 0.;
   sigp_ = 0.;
   sige_ = 0.;
   double cyl_ene = 0.;
 
   float radius2 = radius * radius;
   for (const auto& spot : theSpots_) {
     Point globalPoint(spot.row()[0], spot.row()[1], spot.row()[2]);
     math::XYZPoint local = trans_(globalPoint);
     if (local.Perp2() > radius2)
       continue;
 
     // Select halo hits or not
     if (withHalo && spot.fraction() < 0)
       continue;
     if (!withHalo && !(spot.isCore()))
       continue;
 
     sige_ += (globalPoint.eta() - theCluster_->eta()) * (globalPoint.eta() - theCluster_->eta()) * spot.energy();
     sigp_ += deltaPhi(globalPoint.phi(), theCluster_->phi()) * deltaPhi(globalPoint.phi(), theCluster_->phi()) *
              spot.energy();
 
     sigu_ += local.x() * local.x() * spot.energy();
     sigv_ += local.y() * local.y() * spot.energy();
     cyl_ene += spot.energy();
   }
 
   if (cyl_ene > 0.) {
     const double inv_cyl_ene = 1. / cyl_ene;
     sigu_ = sigu_ * inv_cyl_ene;
     sigv_ = sigv_ * inv_cyl_ene;
     sigp_ = sigp_ * inv_cyl_ene;
     sige_ = sige_ * inv_cyl_ene;
   }
   sigu_ = std::sqrt(sigu_);
   sigv_ = std::sqrt(sigv_);
   sigp_ = std::sqrt(sigp_);
   sige_ = std::sqrt(sige_);
 }

const TVectorD& hgcal::EGammaPCAHelper::eigenValues ( ) const

inline

Definition at line 66 of file EgammaPCAHelper.h.

References pca_.

Referenced by HGCalEgammaIDHelper::eigenValues().

66 { return *pca_->GetEigenValues(); }

hgcal::EGammaPCAHelper::pca_

std::unique_ptr< TPrincipal > pca_

Definition: EgammaPCAHelper.h:101

LongDeps EGammaPCAHelper::energyPerLayer	(	float	radius,
		bool	withHalo = `true`
	)

Definition at line 241 of file EgammaPCAHelper.cc.

References axis_, barycenter_, checkIteration(), debug_, validate-o2o-wbm::f, DetId::HGCalEE, DetId::HGCalHSc, DetId::HGCalHSi, HGCEE, HGCHEB, HGCHEF, LayerTriplets::layers(), maxlayer_, or, CosmicsPD_Skims::radius, theSpots_, and trans_.

Referenced by HGCalEgammaIDHelper::energyPerLayer().

                                                                     {
   if (debug_)
     checkIteration();
   std::set<int> layers;
   float radius2 = radius * radius;
   std::vector<float> energyPerLayer(maxlayer_ + 1, 0.f);
   math::XYZVector mainAxis(axis_);
   mainAxis.unit();
   math::XYZVector phiAxis(barycenter_.x(), barycenter_.y(), 0);
   math::XYZVector udir(mainAxis.Cross(phiAxis));
   udir = udir.unit();
   trans_ = Transform3D(Point(barycenter_),
                        Point(barycenter_ + axis_),
                        Point(barycenter_ + udir),
                        Point(0, 0, 0),
                        Point(0., 0., 1.),
                        Point(1., 0., 0.));
   float energyEE = 0.;
   float energyFH = 0.;
   float energyBH = 0.;
 
   for (const auto& spot : theSpots_) {
     if (!withHalo && !spot.isCore())
       continue;
     math::XYZPoint local = trans_(Point(spot.row()[0], spot.row()[1], spot.row()[2]));
     if (local.Perp2() > radius2)
       continue;
     energyPerLayer[spot.layer()] += spot.energy();
     layers.insert(spot.layer());
     if (spot.detId().det() == DetId::HGCalEE or spot.subdet() == HGCEE) {
       energyEE += spot.energy();
     } else if (spot.detId().det() == DetId::HGCalHSi or spot.subdet() == HGCHEF) {
       energyFH += spot.energy();
     } else if (spot.detId().det() == DetId::HGCalHSc or spot.subdet() == HGCHEB) {
       energyBH += spot.energy();
     }
   }
   return LongDeps(radius, energyPerLayer, energyEE, energyFH, energyBH, layers);
 }

float EGammaPCAHelper::findZFirstLayer ( const LongDeps & ld ) const

private

Definition at line 295 of file EgammaPCAHelper.cc.

References hgcal::LongDeps::energyPerLayer(), hgcal::RecHitTools::getPositionLayer(), maxlayer_, recHitTools_, and PV3DBase< T, PVType, FrameType >::z().

Referenced by clusterDepthCompatibility().

                                                                {
   unsigned int firstLayer = 0;
   for (unsigned il = 1; il <= maxlayer_; ++il) {
     if (ld.energyPerLayer()[il] > 0.) {
       firstLayer = il;
       break;
     }
   }
   // Make dummy DetId to get abs(z) for layer
   return recHitTools_->getPositionLayer(firstLayer).z();
 }

const std::unordered_map<DetId, const HGCRecHit *>* hgcal::EGammaPCAHelper::getHitMap ( )

inline

Definition at line 47 of file EgammaPCAHelper.h.

References hitMap_.

47 { return hitMap_; }

hgcal::EGammaPCAHelper::hitMap_

const std::unordered_map< DetId, const HGCRecHit * > * hitMap_

Definition: EgammaPCAHelper.h:88

void hgcal::EGammaPCAHelper::pcaInitialComputation ( )

inline

Definition at line 53 of file EgammaPCAHelper.h.

References computePCA().

Referenced by HGCalEgammaIDHelper::computeHGCAL().

53 { computePCA(-1., false); }

hgcal::EGammaPCAHelper::computePCA

void computePCA(float radius, bool withHalo=true)

Definition: EgammaPCAHelper.cc:106

const TPrincipal& hgcal::EGammaPCAHelper::pcaResult ( )

void EGammaPCAHelper::printHits ( float radius ) const

Definition at line 281 of file EgammaPCAHelper.cc.

References gather_cfg::cout, mps_fire::i, CosmicsPD_Skims::radius, hgcal::Spot::row(), mathSSE::sqrt(), theSpots_, and trans_.

Referenced by ntuplePrintersDiff.TrackingParticlePrinter::__call__(), ntuplePrintersDiff.SeedPrinter::diff(), ntuplePrintersDiff.TrackPrinter::diff(), ntuplePrintersDiff.TrackingParticlePrinter::diff(), HGCalEgammaIDHelper::printHits(), ntuplePrintersDiff.SeedPrinter::printSeed(), and ntuplePrintersDiff.TrackPrinter::printTrack().

                                                   {
   unsigned nSpots = theSpots_.size();
   float radius2 = radius * radius;
   for (unsigned i = 0; i < nSpots; ++i) {
     Spot spot(theSpots_[i]);
     math::XYZPoint local = trans_(Point(spot.row()[0], spot.row()[1], spot.row()[2]));
     if (local.Perp2() < radius2) {
       std::cout << i << "  " << theSpots_[i].detId().rawId() << " " << theSpots_[i].layer() << " "
                 << theSpots_[i].energy() << " " << theSpots_[i].isCore();
       std::cout << " " << std::sqrt(local.Perp2()) << std::endl;
     }
   }
 }

void hgcal::EGammaPCAHelper::setdEdXWeights ( const std::vector< double > & dEdX )

inline

Definition at line 51 of file EgammaPCAHelper.h.

References HGCalRecHit_cfi::dEdX, and dEdXWeights_.

Referenced by HGCalEgammaIDHelper::HGCalEgammaIDHelper().

51 { dEdXWeights_ = dEdX; }

hgcal::EGammaPCAHelper::dEdXWeights_

std::vector< double > dEdXWeights_

Definition: EgammaPCAHelper.h:84

HGCalRecHit_cfi.dEdX

tuple dEdX

Definition: HGCalRecHit_cfi.py:62

void EGammaPCAHelper::setHitMap ( const std::unordered_map< DetId, const HGCRecHit * > * hitMap )

to set once per event

Definition at line 27 of file EgammaPCAHelper.cc.

References hitMap_, and pcaIteration_.

Referenced by HGCalEgammaIDHelper::eventInit().

                                                                                        {
   hitMap_ = hitMap;
   pcaIteration_ = 0;
 }

void EGammaPCAHelper::setRecHitTools ( const hgcal::RecHitTools * recHitTools )

Definition at line 32 of file EgammaPCAHelper.cc.

References hgcal::RecHitTools::lastLayerBH(), maxlayer_, and recHitTools_.

Referenced by HGCalEgammaIDHelper::eventInit().

                                                                         {
   recHitTools_ = recHitTools;
   maxlayer_ = recHitTools_->lastLayerBH();
 }

double hgcal::EGammaPCAHelper::sigmaEE ( ) const

inline

Definition at line 63 of file EgammaPCAHelper.h.

References checkIteration(), and sige_.

Referenced by HGCalEgammaIDHelper::sigmaEE().

63 { return checkIteration() ? sige_ : -1.; }

hgcal::EGammaPCAHelper::checkIteration

bool checkIteration() const

Definition: EgammaPCAHelper.cc:215

hgcal::EGammaPCAHelper::sige_

double sige_

Definition: EgammaPCAHelper.h:98

double hgcal::EGammaPCAHelper::sigmaPP ( ) const

inline

Definition at line 64 of file EgammaPCAHelper.h.

References checkIteration(), and sigp_.

Referenced by HGCalEgammaIDHelper::sigmaPP().

64 { return checkIteration() ? sigp_ : -1.; }

hgcal::EGammaPCAHelper::sigp_

double sigp_

Definition: EgammaPCAHelper.h:98

hgcal::EGammaPCAHelper::checkIteration

bool checkIteration() const

Definition: EgammaPCAHelper.cc:215

const TVectorD& hgcal::EGammaPCAHelper::sigmas ( ) const

inline

Definition at line 67 of file EgammaPCAHelper.h.

References pca_.

Referenced by HGCalEgammaIDHelper::sigmas().

67 { return *pca_->GetSigmas(); }

hgcal::EGammaPCAHelper::pca_

std::unique_ptr< TPrincipal > pca_

Definition: EgammaPCAHelper.h:101

double hgcal::EGammaPCAHelper::sigmaUU ( ) const

inline

Definition at line 61 of file EgammaPCAHelper.h.

References checkIteration(), and sigu_.

Referenced by HGCalEgammaIDHelper::sigmaUU().

61 { return checkIteration() ? sigu_ : -1.; }

hgcal::EGammaPCAHelper::sigu_

double sigu_

Definition: EgammaPCAHelper.h:98

hgcal::EGammaPCAHelper::checkIteration

bool checkIteration() const

Definition: EgammaPCAHelper.cc:215

double hgcal::EGammaPCAHelper::sigmaVV ( ) const

inline

Definition at line 62 of file EgammaPCAHelper.h.

References checkIteration(), and sigv_.

Referenced by HGCalEgammaIDHelper::sigmaVV().

62 { return checkIteration() ? sigv_ : -1.; }

hgcal::EGammaPCAHelper::sigv_

double sigv_

Definition: EgammaPCAHelper.h:98

hgcal::EGammaPCAHelper::checkIteration

bool checkIteration() const

Definition: EgammaPCAHelper.cc:215

void EGammaPCAHelper::storeRecHits ( const reco::CaloCluster & theCluster )

Definition at line 47 of file EgammaPCAHelper.cc.

References reco::CaloCluster::hitsAndFractions(), and theCluster_.

Referenced by HGCalEgammaIDHelper::computeHGCAL(), and storeRecHits().

                                                                  {
   theCluster_ = &cluster;
   storeRecHits(cluster.hitsAndFractions());
 }

void EGammaPCAHelper::storeRecHits ( const reco::HGCalMultiCluster & cluster )

Definition at line 37 of file EgammaPCAHelper.cc.

References reco::HGCalMultiCluster::begin(), reco::HGCalMultiCluster::end(), mps_fire::result, storeRecHits(), and theCluster_.

                                                                        {
   theCluster_ = &cluster;
   std::vector<std::pair<DetId, float>> result;
   for (reco::HGCalMultiCluster::component_iterator it = cluster.begin(); it != cluster.end(); it++) {
     const std::vector<std::pair<DetId, float>>& hf = (*it)->hitsAndFractions();
     result.insert(result.end(), hf.begin(), hf.end());
   }
   storeRecHits(result);
 }

void EGammaPCAHelper::storeRecHits ( const std::vector< std::pair< DetId, float >> & hf )

private

Definition at line 52 of file EgammaPCAHelper.cc.

References gather_cfg::cout, debug_, dEdXWeights_, DetId::det(), first, HLT_FULL_cff::fraction, hgcal::RecHitTools::getLayerWithOffset(), hgcal::RecHitTools::getPosition(), hgcal::RecHitTools::getSiThickIndex(), hitMap_, invThicknessCorrection_, dqmiolumiharvest::j, phase1PixelTopology::layer, pcaIteration_, DetId::rawId(), recHitTools_, sige_, sigp_, sigu_, sigv_, theSpots_, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

                                                                              {
   std::vector<double> pcavars;
   pcavars.resize(3, 0.);
   theSpots_.clear();
   pcaIteration_ = 0;
 
   sigu_ = 0.;
   sigv_ = 0.;
   sigp_ = 0.;
   sige_ = 0.;
 
   unsigned hfsize = hf.size();
   if (debug_)
     std::cout << "The seed cluster constains " << hfsize << " hits " << std::endl;
 
   if (hfsize == 0)
     return;
 
   for (unsigned int j = 0; j < hfsize; j++) {
     unsigned int layer = recHitTools_->getLayerWithOffset(hf[j].first);
 
     const DetId rh_detid = hf[j].first;
     std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap_->find(rh_detid);
     if (itcheck == hitMap_->end()) {
       edm::LogWarning("EgammaPCAHelper") << " Big problem, unable to find a hit " << rh_detid.rawId() << " "
                                          << rh_detid.det() << " " << HGCalDetId(rh_detid) << std::endl;
       continue;
     }
     if (debug_) {
       std::cout << "DetId " << rh_detid.rawId() << " " << layer << " " << itcheck->second->energy() << std::endl;
       std::cout << " Hit " << itcheck->second << " " << itcheck->second->energy() << std::endl;
     }
     float fraction = hf[j].second;
 
     int thickIndex = recHitTools_->getSiThickIndex(rh_detid);
     double mip = dEdXWeights_[layer] * 0.001;  // convert in GeV
     if (thickIndex > -1 and thickIndex < 3)
       mip *= invThicknessCorrection_[thickIndex];
 
     pcavars[0] = recHitTools_->getPosition(rh_detid).x();
     pcavars[1] = recHitTools_->getPosition(rh_detid).y();
     pcavars[2] = recHitTools_->getPosition(rh_detid).z();
     if (pcavars[2] == 0.)
       edm::LogWarning("EgammaPCAHelper") << " Problem, hit with z =0 ";
     else {
       Spot mySpot(rh_detid, itcheck->second->energy(), pcavars, layer, fraction, mip);
       theSpots_.push_back(mySpot);
     }
   }
   if (debug_) {
     std::cout << " Stored " << theSpots_.size() << " hits " << std::endl;
   }
 }