Inheritance diagram for gar::rec::tpccathodestitch:

Public Member Functions
	tpccathodestitch (fhicl::ParameterSet const &p)

	tpccathodestitch (tpccathodestitch const &)=delete

	tpccathodestitch (tpccathodestitch &&)=delete

tpccathodestitch &	operator= (tpccathodestitch const &)=delete

tpccathodestitch &	operator= (tpccathodestitch &&)=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 Member Functions
bool	cathodematch (const gar::rec::Track &trackA, const gar::rec::Track &trackB, int &fwdflagA, int &fwdflagB, float &deltaX)

Private Attributes
std::string	fInputTrackLabel
	input tracks and associations with TPCClusters. Get vertices associated with these tracks from the assns More...

std::string	fInputVertexLabel
	input vertices and associations with tracks More...

std::string	fInputTPCClusterLabel
	input TPC Clusters More...

std::string	fInputHitLabel
	needed to make TPCCluster - hit associations More...

int	fPrintLevel
	debug printout: 0: none, 1: just track parameters and residuals, 2: all More...

float	fDistCut
	cut in distance between best-matching points More...

float	fCTCut
	cut on cosine of angle matching More...

float	fMaxDX
	cut on maximum translation of dX on one side. More...

float	fMinDX
	cut on minimum translation of dX on one side. 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 61 of file tpccathodestitch_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 100 of file tpccathodestitch_module.cc.

                                                                  : EDProducer{p}
     {
       fInputTrackLabel      = p.get<std::string>("InputTrackLabel","track");
       fInputVertexLabel     = p.get<std::string>("InputVertexLabel","vertex");
       fInputTPCClusterLabel = p.get<std::string>("InputTPCClusterLabel","tpccluster");
       fInputHitLabel        = p.get<std::string>("InputHitLabel","hit");
       fPrintLevel           = p.get<int>("PrintLevel",0);
       fDistCut              = p.get<float>("DistCut",3);
       fCTCut                = p.get<float>("CTCut",0.99);
       fMaxDX                = p.get<float>("MaxDX",50.0);
       fMinDX                = p.get<float>("MinDX",-50.0);
 
       art::InputTag inputTrackTag(fInputTrackLabel);
       consumes< std::vector<rec::Track> >(inputTrackTag);
       consumes< art::Assns<rec::TPCCluster, rec::Track> >(inputTrackTag);
       art::InputTag inputVertexTag(fInputVertexLabel);
       consumes< art::Assns<rec::Vertex, rec::Track> >(inputVertexTag);
       consumes<std::vector<rec::TrackIoniz>>(inputTrackTag);
       consumes<art::Assns<rec::TrackIoniz, rec::Track>>(inputTrackTag);
       art::InputTag inputHitTag(fInputHitLabel);
       consumes<art::Assns<rec::Hit, rec::TPCCluster>>(inputHitTag);
       art::InputTag inputTPCClusterTag(fInputTPCClusterLabel);
       consumes< std::vector<rec::TPCCluster>>(inputTPCClusterTag);
 
       // probably don't need the vector hits at this point if we have the TPCClusters
       //consumes< std::vector<rec::VecHit> >(patrecTag);
       //consumes< art::Assns<rec::VecHit, rec::Track> >(patrecTag);
 
       produces< std::vector<rec::Track> >();
       produces< art::Assns<rec::TPCCluster, rec::Track> >();
       produces<std::vector<rec::TrackIoniz>>();
       produces<art::Assns<rec::TrackIoniz, rec::Track>>();
       produces< std::vector<rec::Vertex> >();
       produces< art::Assns<rec::Track, rec::Vertex, rec::TrackEnd> >();
       produces< std::vector<gar::rec::TPCCluster> >();
       produces< art::Assns<gar::rec::Hit, gar::rec::TPCCluster> >();
     }

gar::rec::tpccathodestitch::tpccathodestitch ( tpccathodestitch const & )

delete

gar::rec::tpccathodestitch::tpccathodestitch ( tpccathodestitch && )

delete

Member Function Documentation

bool gar::rec::tpccathodestitch::cathodematch	(	const gar::rec::Track &	trackA,
		const gar::rec::Track &	trackB,
		int &	fwdflagA,
		int &	fwdflagB,
		float &	deltaX
	)

private

Definition at line 642 of file tpccathodestitch_module.cc.

     {
       afw = 0;
       bfw = 0;
       bool matchable = false;
 
       // may not need this yet -- but magnetic field will be a function of position someday
       // auto const *magFieldService = gar::providerFrom<mag::MagneticFieldService>();
       // G4ThreeVector zerovec(0,0,0);
       // G4ThreeVector magfield = magFieldService->FieldAtPoint(zerovec);
 
       art::ServiceHandle<geo::GeometryGAr> geo;
       double xtpccent = geo->TPCXCent();
       TVector3 xhat(1,0,0);
       TVector3 xcentv = xhat*xtpccent;
 
       std::vector<TVector3> apt;  // endpoint
       std::vector<TVector3> bpt;
       std::vector<TVector3> adir; // direction
       std::vector<TVector3> bdir;
 
       apt.emplace_back(atrack.Vertex());
       apt.emplace_back(atrack.End());
       adir.emplace_back(atrack.VtxDir());
       adir.emplace_back(atrack.EndDir());
       bpt.emplace_back(btrack.Vertex());
       bpt.emplace_back(btrack.End());
       bdir.emplace_back(btrack.VtxDir());
       bdir.emplace_back(btrack.EndDir());
 
       // center the postions on the center of the detector
 
       for (size_t i=0; i<apt.size(); ++i) apt.at(i) -= xcentv;
       for (size_t i=0; i<bpt.size(); ++i) bpt.at(i) -= xcentv;
 
       std::vector<int> fwflag = {2, 1};
 
       for (size_t ia=0; ia<apt.size(); ++ia)
         {
           for (size_t ib=0; ib<bpt.size(); ++ib)
             {
               // directions should point along track pieces to stitch, so they should be back to back
               TVector3 v=adir.at(ia);
               TVector3 u=bdir.at(ib);
 
               if (v.Dot(u) > -fCTCut) continue;
 
               // simple comparisons of vertex positions.  Ignores the fact that we may be missing
               // part of one of the tracks near the endpoint due to the gaps or patrec failure
               //if ( ((apt.at(ia)-bpt.at(ib)).Cross(xhat)).Mag() > dYZCut) continue;
               //if (TMath::Abs(apt.at(ia).X() + bpt.at(ib).X()) > fdXCut) continue;
 
               // solved for X shift called tau
 
               float tdenom = 1.0 - TMath::Sq(xhat.Dot(v));
               if (tdenom == 0) continue;
               TVector3 d = apt.at(ia) - bpt.at(ib);
               float tau = (0.5/tdenom) * d.Dot( (v.Dot(xhat))*v - xhat );
 
               if (tau > fMaxDX) continue;
               if (tau < fMinDX) continue;
 
               // solve for displacement along adir.
 
               float gamma = -v.Dot(2*tau*xhat + d);
               TVector3 chivec = d + gamma*v + 2.0*tau*xhat;
               if (chivec.Mag() > fDistCut) continue;
 
               // we have a match
               afw = fwflag.at(ia);
               bfw = fwflag.at(1-ib);
               deltaX = tau;
               matchable = true;
 
               break;   // just find the first match -- possible to match the same tracks with different endpoints?
             }
           if (matchable) break;
         }
       return matchable;
     }

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

delete

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

delete

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

overridevirtual

Implements art::EDProducer.

Definition at line 140 of file tpccathodestitch_module.cc.

     {
 
       // some useful structs
 
       typedef struct locVtx
       {
         float fPosition[3];
         float fCovMat[9];
         double fTime;
       } locVtx_t;
 
       // for keeping a list of track endpoints associated with vertices.  Keep dx in here so we can average
       // it at a later step
 
       typedef struct trkiend
       {
         size_t trkindex;
         rec::TrackEnd trkend;
         float dx;
       } trkiend_t;
 
       art::ServiceHandle<geo::GeometryGAr> geo;
       float xtpccent = geo->TPCXCent();
 
       // get the distance corresponding to one ADC tick so we can report the time in ticks
 
       auto detProp = gar::providerFrom<detinfo::DetectorPropertiesService>();
       double DriftVelocity = detProp->DriftVelocity(detProp->Efield(),detProp->Temperature());       // in cm per microsecond
       //auto clockService = gar::providerFrom<detinfo::DetectorClocksServiceGAr>();
       //double distonetick = DriftVelocity * (clockService->TPCTick2Time(1) - clockService->TPCTick2Time(0)) ;
 
       // output collections
 
       std::unique_ptr< std::vector<rec::Track> >                            trkCol(new std::vector<rec::Track>);
       std::unique_ptr<std::vector<gar::rec::TPCCluster> >                   TPCClusterCol (new std::vector<gar::rec::TPCCluster> );
       std::unique_ptr< art::Assns<rec::TPCCluster,rec::Track> >             TPCClusterTrkAssns(new art::Assns<rec::TPCCluster,rec::Track>);
       std::unique_ptr< std::vector<rec::TrackIoniz> >                       ionCol(new std::vector<rec::TrackIoniz>);
       std::unique_ptr< art::Assns<rec::TrackIoniz,rec::Track> >             ionTrkAssns(new art::Assns<rec::TrackIoniz,rec::Track>);
       std::unique_ptr< std::vector<rec::Vertex> >                           vtxCol(new std::vector<rec::Vertex>);
       std::unique_ptr< art::Assns<rec::Track, rec::Vertex, rec::TrackEnd> > trkVtxAssns(new art::Assns<rec::Track, rec::Vertex, rec::TrackEnd>);
       std::unique_ptr< art::Assns<rec::Hit,rec::TPCCluster> >               HitTPCClusterAssns(new art::Assns<rec::Hit,rec::TPCCluster>);
 
       // inputs
 
       auto inputTrackHandle = e.getValidHandle< std::vector<rec::Track> >(fInputTrackLabel);
       auto const& inputTracks = *inputTrackHandle;
       auto inputTPCClusterHandle = e.getValidHandle< std::vector<rec::TPCCluster> >(fInputTPCClusterLabel);
       auto const& inputTPCClusters = *inputTPCClusterHandle;
 
       const art::FindManyP<rec::TPCCluster>              TPCClustersFromInputTracks(inputTrackHandle,e,fInputTrackLabel);
       const art::FindManyP<rec::TrackIoniz>              TrackIonizFromInputTracks(inputTrackHandle,e,fInputTrackLabel);
       const art::FindManyP<rec::Vertex, rec::TrackEnd>   VerticesFromInputTracks(inputTrackHandle,e,fInputVertexLabel);
       const art::FindManyP<rec::Hit>                     HitsFromInputTPCClusters(inputTPCClusterHandle,e,fInputTPCClusterLabel);
 
       // we get the input ionization data products from the associations
       //auto inputIonHandle = e.getValidHandle< std::vector<rec::TrackIoniz> >(fInputTrackLabel);
       //auto const& inputIoniz = *inputIonHandle;
 
       // for making output associations
 
       auto const trackPtrMaker =       art::PtrMaker<rec::Track>(e);
       auto const ionizPtrMaker =       art::PtrMaker<rec::TrackIoniz>(e);
       auto const TPCClusterPtrMaker =  art::PtrMaker<rec::TPCCluster>(e);
       auto const vtxPtrMaker   =       art::PtrMaker<rec::Vertex>(e);
 
 
       // and a map of TPC Cluster ID numbers to locations in the hit-cluster assocation FindManyP
 
       std::map<rec::IDNumber,size_t> hcim;
       for (size_t iclus = 0; iclus<inputTPCClusters.size(); ++iclus)
         {
           rec::IDNumber clusid = inputTPCClusters.at(iclus).getIDNumber();
           hcim[clusid] = iclus;
         }
 
       // make a list of which side each track is on by checking the majority of hits
       // this is used so we do not stitch tracks on the same side of the cathode
 
       std::vector<int> trackside;
       float chanpos[3] = {0,0,0};
       for (size_t itrack = 0; itrack < inputTracks.size(); ++itrack)
         {
           size_t nplus = 0;
           size_t nminus = 0;
           for (size_t iTPCCluster=0; iTPCCluster<TPCClustersFromInputTracks.at(itrack).size(); ++iTPCCluster)
             {
               size_t icindex = hcim[TPCClustersFromInputTracks.at(itrack).at(iTPCCluster)->getIDNumber()];
               for (size_t ihit=0; ihit < HitsFromInputTPCClusters.at(icindex).size(); ++ihit)
                 {
                   auto hitchan = HitsFromInputTPCClusters.at(icindex).at(ihit)->Channel();
                   geo->ChannelToPosition(hitchan, chanpos);
                   if (chanpos[0] > xtpccent)
                     {
                       nplus++;
                     }
                   else
                     {
                       nminus++;
                     }
                 }
             }
           if (nplus > nminus)
             {
               trackside.push_back(1);
             }
           else
             {
               trackside.push_back(-1);
             }
         }
 
       // vector of merged flags contains a list of which pairs of tracks are merged
       // A more general track stitcher may want to stitch kinked tracks together and thus be able to stitch
       // more than a pair of tracks, but this module so far just stitches across the cathode.
 
       size_t nti = inputTracks.size();
       std::vector<int> mergedflag(nti,-1); // save indexes here; negative means unmerged
       std::vector<int> fwdflag(nti,0);     // 1 for forwards, 2 for backwards
       std::vector<float> dx(nti,0);       // shift in X needed
 
       for (size_t itrack = 0; itrack < inputTracks.size(); ++itrack)
         {
           if (mergedflag.at(itrack) >= 0) continue;  // this track is already part of an earlier merge
           for (size_t jtrack = itrack+1; jtrack < inputTracks.size(); ++jtrack)
             {
               if (mergedflag.at(jtrack)>=0) continue;  // already merged, don't merge another one
               if (trackside.at(itrack) == trackside.at(jtrack)) continue;  // don't merge tracks on the same side
               if (cathodematch(inputTracks.at(itrack),inputTracks.at(jtrack),fwdflag.at(itrack),fwdflag.at(jtrack),dx.at(itrack)))
                 {
                   if (fPrintLevel >0)
                     {
                       std::cout << "found a merge.  dx= " <<
                         dx.at(itrack) <<
                         " flip flags: " << fwdflag.at(itrack) << " " << fwdflag.at(jtrack) <<  std::endl;
                     }
                   mergedflag.at(itrack) = jtrack;
                   mergedflag.at(jtrack) = itrack;
                   dx.at(jtrack) = -dx.at(itrack);
                   break;
                 }
             }
         }
 
       // we have lists of tracks to merge and how much to shift them by, and which ones to flip.  Look for associated vertices
       // and follow the chain of associated tracks, proposing new shifts to them.
 
       // make a local copy of vertex information so we can change the X locations
 
       std::map<rec::IDNumber,locVtx_t> vtxmap;  // index of vertices by Leo's IDNumber
 
       // a map of vertex to track associations, the inverse of the track to vertex assocation list
       std::map<rec::IDNumber,std::vector<trkiend_t>>  vttrackmap;
 
       // keep track of which vertices we have already shifted
       std::map<rec::IDNumber,bool>  vtxshifted;
 
       // fill out the list of vertices.  A track may belong to more than one vertex (either end, and a track end
       // itself may belong to more than one vertex.  fill in the dx's as needed.
 
       for (size_t itrack = 0; itrack < inputTracks.size(); ++itrack)
         {
           trkiend_t tmpti;
           tmpti.trkindex = itrack;
           tmpti.dx = dx.at(itrack);
 
           for (size_t ivtx = 0; ivtx < VerticesFromInputTracks.at(itrack).size(); ++ivtx)
             {
               auto vtxp = VerticesFromInputTracks.at(itrack).at(ivtx);
               rec::IDNumber vtxid = vtxp->getIDNumber();
               auto vtxiter = vtxmap.find(vtxid);
               if (vtxiter == vtxmap.end())
                 {
                   locVtx_t tmpvtxdat;
                   tmpvtxdat.fTime = vtxp->Time();
                   for (int i=0; i<3; ++i)
                     {
                       tmpvtxdat.fPosition[i] = vtxp->Position()[i];
                     }
                   for (int i=0; i<9; ++i)
                     {
                       tmpvtxdat.fCovMat[i] = vtxp->CovMat()[i];
                     }
                   vtxmap[vtxid] = tmpvtxdat;
                   vtxshifted[vtxid] = false;
 
                   std::vector<trkiend_t> trkindexvector;
                   tmpti.trkend = *(VerticesFromInputTracks.data(itrack).at(ivtx));
                   trkindexvector.push_back(tmpti);
                   vttrackmap[vtxid] = trkindexvector;
                 }
               else
                 {
                   tmpti.trkend = *(VerticesFromInputTracks.data(itrack).at(ivtx));
                   vttrackmap[vtxid].push_back(tmpti);
                 }
             }
         }
 
       // Find new vertex shift positions by averaging the shifts of associated tracks that are on the
       // cathode-crosser list.  To do -- carry uncertainties around and include in the average.
       // This while loop searches for chains of tracks and vertices to shift.
 
       bool moretodo = true;
       while (moretodo)
         {
           moretodo = false;
           for (const auto &vtxi : vttrackmap)
             {
               if (vtxshifted[vtxi.first]) continue;  // skip the ones we've already shifted
 
               float avgdx = 0;
               size_t numdx = 0;
               for (size_t i=0; i< vtxi.second.size(); ++i)
                 {
                   size_t itrack = vtxi.second.at(i).trkindex;
                   if (mergedflag.at(itrack) >= 0 || dx.at(itrack) != 0)
                     {
                       avgdx += dx.at(itrack);
                       ++numdx;
                     }
                 }
               if (numdx > 0)
                 {
                   avgdx /= (float) numdx;
                   vtxshifted[vtxi.first] = true;  // we are shifting a vertex.
                   moretodo = true;        // we'll shift its tracks and thus need to go around looking for more vertices
                 }
               vtxmap[vtxi.first].fPosition[0] += avgdx;
 
               double ts = vtxmap[vtxi.first].fTime;
               double deltat = avgdx/DriftVelocity;
               ts += deltat;
               vtxmap[vtxi.first].fTime = ts;
 
               // shift the as-yet-unshifted tracks (or rather mark them for shifting)
               for (size_t i=0; i< vtxi.second.size(); ++i)
                 {
                   size_t itrack = vtxi.second.at(i).trkindex;
                   if (mergedflag.at(itrack) < 0 && dx.at(itrack) == 0)
                     {
                       dx.at(itrack) = avgdx;
                     }
                 }
             }
         }
 
       // put the new vertices into the output collection, and remember where we put them so we can
       // use it in associations.
 
       std::map<rec::IDNumber,size_t> vtxoutmap;
       for (const auto &vtxr : vtxmap)
         {
           vtxCol->emplace_back(vtxr.second.fPosition,vtxr.second.fCovMat,vtxr.second.fTime);
           vtxoutmap[vtxr.first] = vtxCol->size() - 1;
         }
 
 
       // fill the output track collection, merging, shifting, and flipping as needed.
       // fill in the associations, updating the
       // trkend flags if the tracks have been flipped.
 
       for (size_t itrack = 0; itrack < inputTracks.size(); ++itrack)
         {
           if (mergedflag.at(itrack) < 0)   // not merged, but make a new shifted track if need be
             {
               TrackPar tpi(inputTracks.at(itrack));
               double ts = tpi.getTime();
               double deltat = dx.at(itrack)/DriftVelocity;
               ts += deltat;
               TrackPar tpm(tpi.getLengthForwards(),
                            tpi.getLengthBackwards(),
                            tpi.getNTPCClusters(),
                            tpi.getXBeg() + dx.at(itrack),
                            tpi.getTrackParametersBegin(),
                            tpi.getCovMatBeg(),
                            tpi.getChisqForwards(),
                            tpi.getXEnd() + dx.at(itrack),
                            tpi.getTrackParametersEnd(),
                            tpi.getCovMatEnd(),
                            tpi.getChisqBackwards(),
                            ts);
               trkCol->push_back(tpm.CreateTrack());
 
               //Make copies of the associated TPCClusters and shift them as need be.
 
               auto const trackpointer = trackPtrMaker(trkCol->size()-1);
               float cpos[3] = {0,0,0};
               for (size_t iTPCCluster=0; iTPCCluster<TPCClustersFromInputTracks.at(itrack).size(); ++iTPCCluster)
                 {
                   auto const &tc = TPCClustersFromInputTracks.at(itrack).at(iTPCCluster);
                   for (size_t i=0; i<3; ++i)
                     {
                       cpos[i] = tc->Position()[i];
                     }
                   cpos[0] += dx.at(itrack);
                   TPCClusterCol->emplace_back(tc->Signal(),
                                               cpos,
                                               tc->StartTime(),
                                               tc->EndTime(),
                                               tc->Time(),
                                               tc->RMS(),
                                               tc->CovMatPacked());
                   auto const TPCClusterPointer = TPCClusterPtrMaker(TPCClusterCol->size()-1);
                   TPCClusterTrkAssns->addSingle(TPCClusterPointer,trackpointer);
                   rec::IDNumber oldclusid = tc->getIDNumber();
                   size_t icindex = hcim[oldclusid];
                   for (size_t ihit=0; ihit < HitsFromInputTPCClusters.at(icindex).size(); ++ihit)
                     {
                       HitTPCClusterAssns->addSingle(HitsFromInputTPCClusters.at(icindex).at(ihit),TPCClusterPointer);
                     }
                 }
 
               // copy trkioniz and associations to the output collections.
               // this loop may be unnecessary as there is only one TrackIoniz per Track, but just in case
 
               for (size_t iTrkIoniz=0; iTrkIoniz<TrackIonizFromInputTracks.at(itrack).size(); ++iTrkIoniz)
                 {
                   ionCol->push_back(*TrackIonizFromInputTracks.at(itrack).at(iTrkIoniz));
                   auto const ionizpointer = ionizPtrMaker(ionCol->size()-1);
                   ionTrkAssns->addSingle(ionizpointer, trackpointer);
                 }
 
               // copy over associations with vertices to the new vertex collection.  Retain
               // track-end identification as these tracks are not flipped.
 
               for (size_t ivtx = 0; ivtx < VerticesFromInputTracks.at(itrack).size(); ++ivtx)
                 {
                    auto vtxp = VerticesFromInputTracks.at(itrack).at(ivtx);
                    auto trkend = *(VerticesFromInputTracks.data(itrack).at(ivtx));
                    rec::IDNumber vtxid = vtxp->getIDNumber();
                    size_t ivout = vtxoutmap[vtxid];
                    auto const vtxptr = vtxPtrMaker(ivout);
                    trkVtxAssns->addSingle(trackpointer,vtxptr,trkend);
                 }
             }
           else  // merge them -- adjust track parameters and put in the time
             {
               // assume directions are set so that itrack is the "beginning" and jtrack is the "end"
               if (mergedflag.at(itrack) < 0)
                 {
                   MF_LOG_WARNING("tpccathodestitch: inconsistent merge flags ") << itrack << " " << mergedflag.at(itrack);
                   continue;
                 }
               size_t jtrack = mergedflag.at(itrack);
               if (jtrack < itrack) continue;   // don't merge the same tracks twice
 
               // flip the tracks around according to the fwdflag returned by the cathode match method
               // fwdflag == 1 means don't flip
 
               TrackPar tpi(inputTracks.at(itrack), fwdflag.at(itrack) != 1);
               TrackPar tpj(inputTracks.at(jtrack), fwdflag.at(jtrack) != 1);
 
               tpi.setXBeg( tpi.getXBeg() + dx.at(itrack) );
               tpi.setXEnd( tpi.getXEnd() + dx.at(itrack) );
               tpj.setXBeg( tpj.getXBeg() + dx.at(jtrack) );
               tpj.setXEnd( tpj.getXEnd() + dx.at(jtrack) );
 
               // some checking due to the unsigned nature of the timestmap.  Assume in ticks.
               double ts = tpi.getTime();
               int deltat = dx.at(itrack)/DriftVelocity;
               if ( (int) ts + deltat >= 0)
                 {
                   ts += deltat;
                 }
 
               TrackPar tpm(tpi.getLengthForwards() + tpj.getLengthForwards(),
                            tpi.getLengthBackwards() + tpj.getLengthBackwards(),
                            tpi.getNTPCClusters() + tpj.getNTPCClusters(),
                            tpi.getXBeg(),
                            tpi.getTrackParametersBegin(),
                            tpi.getCovMatBeg(),
                            tpi.getChisqForwards() + tpj.getChisqForwards(),
                            tpj.getXEnd(),
                            tpj.getTrackParametersEnd(),
                            tpj.getCovMatEnd(),
                            tpi.getChisqBackwards() + tpj.getChisqBackwards(),
                            ts);
               trkCol->push_back(tpm.CreateTrack());
               auto const trackpointer = trackPtrMaker(trkCol->size()-1);
 
               // make new TPCClusters, shifted with the track dx.  To think about: the times.  They are arrival times of hit charge at the moment
 
               float cpos[3] = {0,0,0};
               for (size_t iTPCCluster=0; iTPCCluster<TPCClustersFromInputTracks.at(itrack).size(); ++iTPCCluster)
                 {
                   auto const &tc = TPCClustersFromInputTracks.at(itrack).at(iTPCCluster);
                   for (size_t i=0; i<3; ++i)
                     {
                       cpos[i] = tc->Position()[i];
                     }
                   cpos[0] += dx.at(itrack);
                   TPCClusterCol->emplace_back(tc->Signal(),
                                               cpos,
                                               tc->StartTime(),
                                               tc->EndTime(),
                                               tc->Time(),
                                               tc->RMS(),
                                               tc->CovMatPacked());
                   auto const TPCClusterPointer = TPCClusterPtrMaker(TPCClusterCol->size()-1);
                   TPCClusterTrkAssns->addSingle(TPCClusterPointer,trackpointer);
                   rec::IDNumber oldclusid = tc->getIDNumber();
                   size_t icindex = hcim[oldclusid];
                   for (size_t ihit=0; ihit < HitsFromInputTPCClusters.at(icindex).size(); ++ihit)
                     {
                       HitTPCClusterAssns->addSingle(HitsFromInputTPCClusters.at(icindex).at(ihit),TPCClusterPointer);
                     }
                 }
               for (size_t iTPCCluster=0; iTPCCluster<TPCClustersFromInputTracks.at(jtrack).size(); ++iTPCCluster)
                 {
                   auto const &tc = TPCClustersFromInputTracks.at(jtrack).at(iTPCCluster);
                   for (size_t i=0; i<3; ++i)
                     {
                       cpos[i] = tc->Position()[i];
                     }
                   cpos[0] += dx.at(jtrack);
                   TPCClusterCol->emplace_back(tc->Signal(),
                                               cpos,
                                               tc->StartTime(),
                                               tc->EndTime(),
                                               tc->Time(),
                                               tc->RMS(),
                                               tc->CovMatPacked());
                   auto const TPCClusterPointer = TPCClusterPtrMaker(TPCClusterCol->size()-1);
                   TPCClusterTrkAssns->addSingle(TPCClusterPointer,trackpointer);
                   rec::IDNumber oldclusid = tc->getIDNumber();
                   size_t icindex = hcim[oldclusid];
                   for (size_t ihit=0; ihit < HitsFromInputTPCClusters.at(icindex).size(); ++ihit)
                     {
                       HitTPCClusterAssns->addSingle(HitsFromInputTPCClusters.at(icindex).at(ihit),TPCClusterPointer);
                     }
                 }
 
               // merge the track ionization data -- assume just one TrackIoniz per track, with forward and backward vectors
               // reverse them if need be.
 
               const auto itif = TrackIonizFromInputTracks.at(itrack).at(0)->getFWD_dSigdXs();
               const auto itib = TrackIonizFromInputTracks.at(itrack).at(0)->getBAK_dSigdXs();
               const auto jtif = TrackIonizFromInputTracks.at(jtrack).at(0)->getFWD_dSigdXs();
               const auto jtib = TrackIonizFromInputTracks.at(jtrack).at(0)->getBAK_dSigdXs();
 
               auto mtif = itif;
               auto mtib = itib;
               if (fwdflag.at(itrack) != 1)      //  flip itrack's ionization vector if need be
                 {
                   mtif = itib;
                   mtib = itif;
                 }
               auto jmtif = jtif;
               auto jmtib = jtib;
               if (fwdflag.at(jtrack) != 1)
                 {
                   jmtif = jtib;
                   jmtib = jtif;
                 }
 
               for (size_t i=0; i<jmtif.size(); ++i)
                 {
                   mtif.push_back(jmtif.at(i));
                 }
               for (size_t i=0; i<jmtib.size(); ++i)
                 {
                   mtib.push_back(jmtib.at(i));
                 }
               TrackIoniz tim;
               tim.setData(mtif,mtib);
               ionCol->push_back(tim);
               auto const ionizpointer = ionizPtrMaker(ionCol->size()-1);
               ionTrkAssns->addSingle(ionizpointer, trackpointer);
 
               for (size_t ivtx = 0; ivtx < VerticesFromInputTracks.at(itrack).size(); ++ivtx)
                 {
                    auto vtxp = VerticesFromInputTracks.at(itrack).at(ivtx);
                    auto trkend = *(VerticesFromInputTracks.data(itrack).at(ivtx));
                    if (fwdflag.at(itrack) != 1)
                      {
                        trkend = 1 - trkend;
                      }
                    rec::IDNumber vtxid = vtxp->getIDNumber();
                    size_t ivout = vtxoutmap[vtxid];
                    auto const vtxptr = vtxPtrMaker(ivout);
                    trkVtxAssns->addSingle(trackpointer,vtxptr,trkend);
                 }
 
             }
         }
 
       e.put(std::move(trkCol));
       e.put(std::move(TPCClusterCol));
       e.put(std::move(TPCClusterTrkAssns));
       e.put(std::move(ionCol));
       e.put(std::move(ionTrkAssns));
       e.put(std::move(vtxCol));
       e.put(std::move(trkVtxAssns));
       e.put(std::move(HitTPCClusterAssns));
     }

Member Data Documentation

float gar::rec::tpccathodestitch::fCTCut

private

cut on cosine of angle matching

Definition at line 86 of file tpccathodestitch_module.cc.

float gar::rec::tpccathodestitch::fDistCut

private

cut in distance between best-matching points

Definition at line 85 of file tpccathodestitch_module.cc.

std::string gar::rec::tpccathodestitch::fInputHitLabel

private

needed to make TPCCluster - hit associations

Definition at line 83 of file tpccathodestitch_module.cc.

std::string gar::rec::tpccathodestitch::fInputTPCClusterLabel

private

input TPC Clusters

Definition at line 82 of file tpccathodestitch_module.cc.

std::string gar::rec::tpccathodestitch::fInputTrackLabel

private

input tracks and associations with TPCClusters. Get vertices associated with these tracks from the assns

Definition at line 80 of file tpccathodestitch_module.cc.

std::string gar::rec::tpccathodestitch::fInputVertexLabel

private

input vertices and associations with tracks

Definition at line 81 of file tpccathodestitch_module.cc.

float gar::rec::tpccathodestitch::fMaxDX

private

cut on maximum translation of dX on one side.

Definition at line 87 of file tpccathodestitch_module.cc.

float gar::rec::tpccathodestitch::fMinDX

private

cut on minimum translation of dX on one side.

Definition at line 88 of file tpccathodestitch_module.cc.

int gar::rec::tpccathodestitch::fPrintLevel

private

debug printout: 0: none, 1: just track parameters and residuals, 2: all

Definition at line 84 of file tpccathodestitch_module.cc.

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

garsoft/Reco/tpccathodestitch_module.cc

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation