MCUtil Namespace Reference

A collection of functions useful for MC processing. More...

Functions
int	GetMCTracks (const RBKTrack &trk, std::vector< int > &trkList)
	Fill a vector of MC track numbers which generated this track.

---

Detailed Description

A collection of functions useful for MC processing.

---

Function Documentation

int MCUtil::GetMCTracks	(	const RBKTrack &	trk,
		std::vector< int > &	trkList
	)

Fill a vector of MC track numbers which generated this track.

The first entry in the vector will be the most "significant" MCTrack. Use that for MCTrack-TPCRTrack association.

Returns:

The difference in the number of hits in the first track and the second most-significant track -- use it as a measure of pileup.

Definition at line 37 of file MCUtil.cxx.

References RBKTrack::NSpacePoints(), RBKTrack::NViewHits(), RBKTrack::SpacePoint(), tl, and RBKTrack::ViewHit().

Referenced by RICHReco::EvFill().

{
 // Loop over all hits in the track
  int nhit = trk.NSpacePoints();
  trkList.clear();

  map<int, int, TrkNumComp> trkHitCount;

  for (int i = 0; i < nhit; ++i) {
    const TPCRHit* hit =
      dynamic_cast<const TPCRHit*>(&trk.SpacePoint(i));
    if (!hit) continue;

    // Loop over all the digits that make up this hit
    int ndig = hit->Cluster()->NDigits();
    for (int j = 0; j < ndig; ++j) {
      const TPCDigit* dig = hit->Cluster()->Digit(j);
      if (!dig) continue;

      const MCCTrkList* tl = dig->TrkList();
      if (!tl) continue;

      // Loop over all the MCCHit's that make up this digit and
      // count the number of times each MC track happens
      for (int k = 0; k < tl->NTrack(); ++k) 
    ++(trkHitCount[tl->Track(k)]);
    }
  }

  // Now loop over wire clusters
  for (int i = 0; i < trk.NViewHits(); ++i) {
    const WireClust* hit = 
      dynamic_cast<const WireClust*>(&trk.ViewHit(i));
    if (!hit) continue;

    // And loop over wires that make up the cluster
    for (int j = 0; j < hit->NWire(); ++j) {
      const Wire* w = hit->GetWire(j);
      if (!w) continue;

      const MCCTrkList* tl = w->MCTracks();
      if (!tl) continue;

      for (int k = 0; k < tl->NTrack(); ++k) 
    ++(trkHitCount[tl->Track(k)]);
    }
  }
 
  if (trkHitCount.empty()) return 0;

  // Find the most populous track and the second most populous track
  int maxNmcTrk = 0;
  int maxIndx = 0;
  int maxNmcTrk2 = 0;
  int maxIndx2 = 0;
  map<int, int, TrkNumComp>::iterator itr(trkHitCount.begin());  
  map<int, int, TrkNumComp>::iterator itrEnd(trkHitCount.end());  
  for ( ; itr != itrEnd; ++itr) {
    int n = (*itr).second;
    if (n > maxNmcTrk) {
      maxNmcTrk2 = maxNmcTrk;
      maxIndx2 = maxIndx;

      maxNmcTrk = n;
      maxIndx = (*itr).first;
    }
    else if (n > maxNmcTrk2) {
      maxNmcTrk2 = n;
      maxIndx2 = (*itr).first;
    }
  }

  // Insert the "most visible" track first -- this is what most
  // everyone should use for TPCRTrack-MCTrack association
  trkList.push_back(maxIndx);
  if (trkHitCount.size() == 1) return 0;
  trkList.push_back(maxIndx2);

  itr = trkHitCount.begin();  
  for ( ; itr != itrEnd; ++itr) {
    int trkNum = (*itr).first;
    if (trkNum == maxIndx || trkNum == maxIndx2) continue;
    trkList.push_back(trkNum);
  }

  return maxNmcTrk - maxNmcTrk2;

}

---