gar::rec::TPCHitCluster Class Reference

Inheritance diagram for gar::rec::TPCHitCluster:

Public Member Functions
	TPCHitCluster (fhicl::ParameterSet const &p)
	TPCHitCluster (TPCHitCluster const &)=delete
	TPCHitCluster (TPCHitCluster &&)=delete
TPCHitCluster &	operator= (TPCHitCluster const &)=delete
TPCHitCluster &	operator= (TPCHitCluster &&)=delete
void	produce (art::Event &e) override
Public Member Functions inherited from art::EDProducer
	EDProducer (fhicl::ParameterSet const &pset)
template<typename Config >
	EDProducer (Table< Config > const &config)
std::string	workerType () const
Public Member Functions inherited from art::detail::Producer
virtual	~Producer () noexcept
	Producer (fhicl::ParameterSet const &)
	Producer (Producer const &)=delete
	Producer (Producer &&)=delete
Producer &	operator= (Producer const &)=delete
Producer &	operator= (Producer &&)=delete
void	doBeginJob (SharedResources const &resources)
void	doEndJob ()
void	doRespondToOpenInputFile (FileBlock const &fb)
void	doRespondToCloseInputFile (FileBlock const &fb)
void	doRespondToOpenOutputFiles (FileBlock const &fb)
void	doRespondToCloseOutputFiles (FileBlock const &fb)
bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)
bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
Public Member Functions inherited from art::Modifier
	~Modifier () noexcept
	Modifier ()
	Modifier (Modifier const &)=delete
	Modifier (Modifier &&)=delete
Modifier &	operator= (Modifier const &)=delete
Modifier &	operator= (Modifier &&)=delete
Public Member Functions inherited from art::ModuleBase
virtual	~ModuleBase () noexcept
	ModuleBase ()
ModuleDescription const &	moduleDescription () const
void	setModuleDescription (ModuleDescription const &)
std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const
void	sortConsumables (std::string const &current_process_name)
template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)
template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Private Attributes
std::string	fHitLabel
	label of module to get the hits from More...
int	fClusterHits
	hit clustering algorithm number More...
float	fHitClusterDx
	range in cm to look for hits to cluster in x More...
float	fHitClusterDyDz
	range in cm to look for hits to cluster in y and z More...
float	fHitClusterDr_I
	range in R to look for hits to cluster for IROC, in cm More...
float	fHitClusterDfi_I
	range in R*dPhi to look for hits to cluster for IROC, in cm More...
float	fHitClusterDrIO
	range in R to look for hits to cluster for IOROC, in cm More...
float	fHitClusterDfiIO
	range in R*dPhi to look for hits to cluster for IROC, in cm More...
float	fHitClusterDrOO
	range in R to look for hits to cluster for OOROC, in cm More...
float	fHitClusterDfiOO
	range in R*dPhi to look for hits to cluster for IROC, in cm More...
int	fPrintLevel
	debug printout: 0: none, 1: selected, 2: all More...

Additional Inherited Members
Public Types inherited from art::EDProducer
using	ModuleType = EDProducer
using	WorkerType = WorkerT< EDProducer >
Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using	Table = Modifier::Table< UserConfig, KeysToIgnore >
Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using	Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
Static Public Member Functions inherited from art::EDProducer
static void	commitEvent (EventPrincipal &ep, Event &e)
Protected Member Functions inherited from art::ModuleBase
ConsumesCollector &	consumesCollector ()
template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)
template<typename T , BranchType = InEvent>
void	consumesMany ()
template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)
template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)
template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Detailed Description

Definition at line 35 of file TPCHitCluster_module.cc.

Constructor & Destructor Documentation

gar::rec::TPCHitCluster::TPCHitCluster ( fhicl::ParameterSet const &  p )

explicit

Definition at line 73 of file TPCHitCluster_module.cc.

      : EDProducer{p}
    {
      fHitLabel         = p.get<std::string>("HitLabel","hit");
      fPrintLevel       = p.get<int>("PrintLevel",0);
      fClusterHits      = p.get<int>("ClusterHits",1);
      fHitClusterDx     = p.get<float>("HitClusterDx",   1.00);
      fHitClusterDyDz   = p.get<float>("HitClusterDyDz", 1.20);
      fHitClusterDr_I   = p.get<float>("HitClusterDr_I", 1.50);
      fHitClusterDfi_I  = p.get<float>("HitClusterDfi_I",1.06);
      fHitClusterDrIO   = p.get<float>("HitClusterDrIO", 2.00);
      fHitClusterDfiIO  = p.get<float>("HitClusterDfiIO",1.55);
      fHitClusterDrOO   = p.get<float>("HitClusterDrOO", 3.00);
      fHitClusterDfiOO  = p.get<float>("HitClusterDfiO)",1.73);

      art::InputTag hitTag(fHitLabel);
      consumes< std::vector<gar::rec::Hit> >(hitTag);
      produces< std::vector<gar::rec::TPCCluster> >();
      produces< art::Assns<gar::rec::Hit, gar::rec::TPCCluster> >();
    }

gar::rec::TPCHitCluster::TPCHitCluster ( TPCHitCluster const &  )

delete

gar::rec::TPCHitCluster::TPCHitCluster ( TPCHitCluster &&  )

delete

Member Function Documentation

TPCHitCluster& gar::rec::TPCHitCluster::operator= ( TPCHitCluster const &  )

delete

TPCHitCluster& gar::rec::TPCHitCluster::operator= ( TPCHitCluster &&  )

delete

void gar::rec::TPCHitCluster::produce ( art::Event &  e )

overridevirtual

Implements art::EDProducer.

Definition at line 99 of file TPCHitCluster_module.cc.

    {
      art::ServiceHandle<geo::GeometryGAr> euclid;
      float xtpccent = euclid->TPCXCent();
      float ytpccent = euclid->TPCYCent();
      float ztpccent = euclid->TPCZCent();

      // input: hits

      auto hitHandle = e.getValidHandle< std::vector<gar::rec::Hit> >(fHitLabel);
      auto const& hits = *hitHandle;

      // output: TPCCluster and associations with hits

      std::unique_ptr<std::vector<gar::rec::TPCCluster> > TPCClusterCol (new std::vector<gar::rec::TPCCluster> );
      std::unique_ptr<art::Assns<gar::rec::Hit,gar::rec::TPCCluster> > hitClusterAssns(new ::art::Assns<gar::rec::Hit,gar::rec::TPCCluster>);

      auto const TPCClusterPtrMaker = art::PtrMaker<gar::rec::TPCCluster>(e);
      auto const hitPtrMaker = art::PtrMaker<rec::Hit>(e, hitHandle.id());

      size_t nhits = hits.size();
      std::vector<float> hptmp(nhits);
      std::vector<int> hsi(nhits);
      std::vector<int> used(nhits,0);

      // sort hits by position in x

      for (size_t ihit=0;ihit<nhits;++ihit)
        {
          hptmp[ihit] = hits.at(ihit).Position()[0];
        }
      TMath::Sort( (int) nhits, hptmp.data(), hsi.data() );

      // loop through hits in x, clustering as we go, marking hits as used.

      for (size_t ihitx=0; ihitx< nhits; ++ihitx)
        {
          size_t ihit = hsi[ihitx];  // unwound index to hit array
          if (used[ihit]) continue;

          const float *xyz_fromhit = hits.at(ihit).Position();
              float xyz[3] = {xyz_fromhit[0] -xtpccent, xyz_fromhit[1] -ytpccent, xyz_fromhit[2] -ztpccent};

          // start a cluster with just this hit
          std::vector<size_t> hitsinclus;  // keep track of hit indices in this cluster so we can make associations
          hitsinclus.push_back(ihit);
          used[ihit] = 1;

          int side=0;                                   // detector side of the seed hit.
          auto hitchan = hits.at(ihit).Channel();
          float chanpos[3] = {0,0,0};
          euclid->ChannelToPosition(hitchan, chanpos);
          if (chanpos[0] > xtpccent)
            {
              side = 1;
            }
          else
            {
              side = -1;
            }

          // Clusters can be displaced in time.  We may want to save which end of the
          // chamber the data came from however and not cluster together charge
          // that comes from different sides.

          if (fPrintLevel > 1)
            {
              std::cout << "Started a new TPC cluster.  Pos= " << xyz[0] << " " << xyz[1] << " " << xyz[2] << std::endl;
              std::cout << "Side: " << side << std::endl;
            }

          // Determine cluster position in endcap for cuts to accumulate more hits
          float rHit = std::hypot(xyz[1],xyz[2]);
          bool In_CROC =                                                  rHit < euclid->GetIROCInnerRadius();
          bool In_IROC = euclid->GetIROCInnerRadius() < rHit               && rHit < euclid->GetIROCOuterRadius();
          bool InIOROC = euclid->GetOROCInnerRadius() < rHit               && rHit < euclid->GetOROCPadHeightChangeRadius();
          bool InOOROC = euclid->GetOROCPadHeightChangeRadius() < rHit && rHit < euclid->GetOROCOuterRadius();

          for (size_t ix = ihitx+1; ix<nhits; ++ix) // look for candidate hits to cluster in with this one
            {
              size_t ihc = hsi[ix];  // candidate hit to add to the cluster if it's in range
              if (used[ihc]) continue;

              int sidetest = 0;
              auto hitchantest = hits.at(ihc).Channel();
              float chanpostest[3] = {0,0,0};
              euclid->ChannelToPosition(hitchantest, chanpostest);
              if (chanpostest[0] > xtpccent)
                {
                  sidetest = 1;
                }
              else
                {
                  sidetest = -1;
                }
              if (sidetest != side) continue;  // don't cluster hits that were detected on opposite TPC endplates.

              const float *xyz2_fromhit = hits.at(ihc).Position();
                  float xyz2[3] = {xyz2_fromhit[0] -xtpccent, xyz2_fromhit[1] -ytpccent, xyz2_fromhit[2] -ztpccent};

              if (fPrintLevel > 1)
                {
                  std::cout << " Testing a hit: " << xyz2[0] << " " << xyz2[1] << " " << xyz2[2] << " " << sidetest << " "
                            << used[ihc] << " " << ihc << " " << ix << std::endl;
                }

              // Check if hit ix is close to ihit, first in drift direction, then
              // in transverse plane depending on which ROC has ihit
              if ( fabs(xyz2[0] -xyz[0]) > fHitClusterDx ) break;

              if (In_CROC) {
                if ( fabs(xyz2[1] -xyz[1]) > fHitClusterDyDz ) continue;
                if ( fabs(xyz2[2] -xyz[2]) > fHitClusterDyDz ) continue;
              }
              // Project out the radial and rotational components of xyz2 -xyz.
              float RhatY = xyz[1]/std::hypot(xyz[1],xyz[2]);
              float RhatZ = xyz[2]/std::hypot(xyz[1],xyz[2]);
              float clustDist = std::hypot(xyz2[1] -xyz[1],xyz2[2] -xyz[2]);
              float dR  = TMath::Abs((xyz2[1] -xyz[1])*RhatY +(xyz2[2] -xyz[2])*RhatZ);
              float Rdfi = sqrt(clustDist*clustDist -dR*dR);
              if (In_IROC) {
                if (dR   > fHitClusterDr_I)  continue;
                if (Rdfi > fHitClusterDfi_I) continue;
              }
              if (InIOROC) {
                if (dR   > fHitClusterDrIO)  continue;
                if (Rdfi > fHitClusterDfiIO) continue;
              }
              if (InOOROC) {
                if (dR   > fHitClusterDrOO)  continue;
                if (Rdfi > fHitClusterDfiOO) continue;
              }

              // add hit to cluster
              hitsinclus.push_back(ihc);
              used[ihc] = 1;
              if (fPrintLevel > 1)
                {
                  std::cout << "Added hit with pos: " << xyz2[0] << " " << xyz2[1] << " " << xyz2[2] << std::endl;
                }

            }

          // calculate cluster charge, centroid, min and max times, and covariance

          double cpos[3] = {0,0,0};  // charge-weighted cluster position
          double csig = 0;           // cluster signal total
          double cstime = 0;         // cluster start time
          double cetime = 0;         // cluster end time
          double crms = 0;           // cluster RMS in X
          double ctime = 0;          // cluster time centroid
          double cov[3][3] = { {0,0,0}, {0,0,0}, {0,0,0} };

          for (size_t ix = 0; ix < hitsinclus.size(); ++ix)
            {
              size_t ihx = hitsinclus.at(ix);
              double hsig = hits.at(ihx).Signal();
              csig += hsig;
              for (size_t idim=0; idim<3; ++idim)
                {
                  cpos[idim] += hsig * hits.at(ihx).Position()[idim];
                }
              ctime += hsig * hits.at(ihx).Time();
              if (ix == 0)
                {
                  cstime = hits.at(ihx).StartTime();
                  cetime = hits.at(ihx).EndTime();
                }
              else
                {
                  cstime = TMath::Min(cstime, (double) hits.at(ihx).StartTime());
                  cetime = TMath::Max(cetime, (double) hits.at(ihx).EndTime());
                }
            }
          if (csig != 0)
            {
              ctime /= csig;
              for (size_t idim=0; idim<3; ++idim)
                {
                  cpos[idim] /= csig;
                }

              // calculate covariance.  Do a numerical sum over hit widths assuming Gaussians.

              for (size_t ix = 0; ix < hitsinclus.size(); ++ix)
                {
                  size_t ihx = hitsinclus.at(ix);
                  double cfrac = hits.at(ihx).Signal()/csig;
                  const float *hpos = hits.at(ihx).Position();
                  for (size_t idim=0; idim<3; ++idim)
                    {
                      for (size_t jdim=0; jdim<3; ++jdim)
                        {
                          cov[idim][jdim] += cfrac*(hpos[idim]-cpos[idim])*(hpos[jdim]-cpos[jdim]);
                        }
                    }
                  cov[0][0] += cfrac*TMath::Sq(hits.at(ihx).RMS());
                }
            }
          crms = TMath::Sqrt(cov[0][0]);  // to mean what it meant before -- RMS along the drift direction

          float fcpos[3] = {0,0,0};
          for (int i=0;i<3;++i)
            {
              fcpos[i] = cpos[i];
            }
          float fccov[6] = {(float) cov[0][0], (float) cov[1][0], (float) cov[2][0],
                            (float) cov[1][1], (float) cov[1][2], (float) cov[2][2]};

          TPCClusterCol->emplace_back(csig,
                                      fcpos,
                                      cstime,
                                      cetime,
                                      ctime,
                                      crms,
                                      fccov);

          if (fPrintLevel > 0)
            {
              std::cout << "Made a TPC Cluster pos: " << fcpos[0] << " " << fcpos[1] << " " << fcpos[2] << 
                "  signal: " << csig << "  nHits: " << hitsinclus.size() << std::endl;
            }
          auto const TPCClusterPointer = TPCClusterPtrMaker(TPCClusterCol->size()-1);
          for (size_t i=0; i<hitsinclus.size(); ++i)
            {
              auto const HitPointer = hitPtrMaker(hitsinclus.at(i));
              hitClusterAssns->addSingle(HitPointer,TPCClusterPointer);
            }

        }

      e.put(std::move(TPCClusterCol));
      e.put(std::move(hitClusterAssns));

    }

Member Data Documentation

int gar::rec::TPCHitCluster::fClusterHits

private

hit clustering algorithm number

Definition at line 53 of file TPCHitCluster_module.cc.

float gar::rec::TPCHitCluster::fHitClusterDfi_I

private

range in R*dPhi to look for hits to cluster for IROC, in cm

Definition at line 58 of file TPCHitCluster_module.cc.

float gar::rec::TPCHitCluster::fHitClusterDfiIO

private

range in R*dPhi to look for hits to cluster for IROC, in cm

Definition at line 60 of file TPCHitCluster_module.cc.

float gar::rec::TPCHitCluster::fHitClusterDfiOO

private

range in R*dPhi to look for hits to cluster for IROC, in cm

Definition at line 62 of file TPCHitCluster_module.cc.

float gar::rec::TPCHitCluster::fHitClusterDr_I

private

range in R to look for hits to cluster for IROC, in cm

Definition at line 57 of file TPCHitCluster_module.cc.

float gar::rec::TPCHitCluster::fHitClusterDrIO

private

range in R to look for hits to cluster for IOROC, in cm

Definition at line 59 of file TPCHitCluster_module.cc.

float gar::rec::TPCHitCluster::fHitClusterDrOO

private

range in R to look for hits to cluster for OOROC, in cm

Definition at line 61 of file TPCHitCluster_module.cc.

float gar::rec::TPCHitCluster::fHitClusterDx

private

range in cm to look for hits to cluster in x

Definition at line 54 of file TPCHitCluster_module.cc.

float gar::rec::TPCHitCluster::fHitClusterDyDz

private

range in cm to look for hits to cluster in y and z

Definition at line 55 of file TPCHitCluster_module.cc.

std::string gar::rec::TPCHitCluster::fHitLabel

private

label of module to get the hits from

Definition at line 52 of file TPCHitCluster_module.cc.

int gar::rec::TPCHitCluster::fPrintLevel

private

debug printout: 0: none, 1: selected, 2: all

Definition at line 64 of file TPCHitCluster_module.cc.

---

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

• garsoft/Reco/TPCHitCluster_module.cc