Inheritance diagram for gar::rec::TPCECALAssociation:

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

	TPCECALAssociation (TPCECALAssociation const &)=delete

	TPCECALAssociation (TPCECALAssociation &&)=delete

TPCECALAssociation &	operator= (TPCECALAssociation const &)=delete

TPCECALAssociation &	operator= (TPCECALAssociation &&)=delete

void	beginJob () override

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	fTrackLabel
	label to find the reco tracks More...

std::string	fClusterLabel
	label to find the right reco caloclusters More...

std::string	fInstanceName

int	fVerbosity

float	fTrackEndXCut
	Extrapolate only track ends outside central drift. More...

float	fTrackEndRCut
	Extrapolate only track ends outside central region. More...

float	fPerpRCut
	Max dist cluster center to circle of track (z,y) only. More...

float	fBarrelXCut
	Max dist cluster center (drift time compensated) track in x. More...

float	fEndcapRphiCut
	Max dist cluster center (drift time compensated) track in r*phi. More...

float	fClusterDirNhitCut
	Do not cut on cluster direction unless you have this many hits or more. More...

float	fClusterDirCut
	Min cos angle between extrapolated track and cluster direction. More...

const detinfo::DetectorProperties *	fDetProp
	detector properties More...

const detinfo::DetectorClocks *	fClocks
	detector clock information More...

const geo::GeometryCore *	fGeo
	pointer to the geometry More...

float	maxXdisplacement

float	ItsInTulsa [3]

TH1F *	radClusTrack
	Cluster to track in transverse plane. More...

TH2F *	xClusTrack
	Cluster to track in barrel, x only, vs x of trackend. More...

TH1F *	phiClusTrack
	Cluster to track in endcap. More...

TH1F *	rPhiClusTrack

TH2F *	dotClustTrack
	cos(cluster axis to extrapolated track) More...

TH1F *	nEntryVsCut
	Number passing each cut. 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 39 of file TPCECALAssociation_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 88 of file TPCECALAssociation_module.cc.

                                                                           : EDProducer{p} {
 
             fTrackLabel        = p.get<std::string>("TrackLabel", "track");
             fClusterLabel      = p.get<std::string>("ClusterLabel","calocluster");
             fInstanceName      = p.get<std::string>("InstanceName","");
             fVerbosity         = p.get<int>("Verbosity", 0);
             // Needs to be computed in produce; so will be 0.0 to do that 
             // otherwise of course the fcl file value appears
             fTrackEndXCut      = p.get<float>("TrackEndXCut",     215.0);
             fTrackEndRCut      = p.get<float>("TrackEndRCut",     230.0);
             fPerpRCut          = p.get<float>("PerpRCut",          10.0);
             fBarrelXCut        = p.get<float>("BarrelXCut",        10.0);
             fEndcapRphiCut     = p.get<float>("EndXCut",           10.0);
             fClusterDirNhitCut = p.get<int>  ("ClusterDirNhitCut",    5);
             fClusterDirCut     = p.get<float>("ClusterDirCut",     0.40);
             fDetProp = gar::providerFrom<detinfo::DetectorPropertiesService>();
             fClocks  = gar::providerFrom<detinfo::DetectorClocksServiceGAr>();
             fGeo     = gar::providerFrom<geo::GeometryGAr>();
 
             produces< art::Assns<gar::rec::Cluster, gar::rec::Track, gar::rec::TrackEnd > >(fInstanceName);
         }

gar::rec::TPCECALAssociation::TPCECALAssociation ( TPCECALAssociation const & )

delete

gar::rec::TPCECALAssociation::TPCECALAssociation ( TPCECALAssociation && )

delete

Member Function Documentation

void gar::rec::TPCECALAssociation::beginJob ( )

overridevirtual

Reimplemented from art::EDProducer.

Definition at line 112 of file TPCECALAssociation_module.cc.

                                           {
 
             ItsInTulsa[0] = fGeo->TPCXCent();
             ItsInTulsa[1] = fGeo->TPCYCent();
             ItsInTulsa[2] = fGeo->TPCZCent();
 
             if (fVerbosity>0) {
                 art::ServiceHandle< art::TFileService > tfs;
                 radClusTrack     = tfs->make<TH1F>("radClusTrack",
                     "(z,y) distance from cluster to track", 100, 0.0,60.0);
                 xClusTrack       = tfs->make<TH2F>("xClusTrack",
                     "x of track minus x of cluster vs x position of track end in barrel",
                     100,-250.0,+250.0, 100,-200.0,+200.0);
                 phiClusTrack    = tfs->make<TH1F>("phiClusTrack",
                     "angular distance from cluster to track in endcap", 90,-M_PI,+M_PI);
                 rPhiClusTrack    = tfs->make<TH1F>("rPhiClusTrack",
                     "distance from cluster to track in endcap", 100,-100.0,+100.0);
                 dotClustTrack    = tfs->make<TH2F>("dotClustTrack",
                     "No hits in cluster vs dot product of cluster axis to extrapolated track direction",
                     100,-1.0,+1.0, 100,0.5,100.5);
                 nEntryVsCut      = tfs->make<TH1F>("nEntryVsCut",
                     "Number track end & cluster combos vs sequentially applied cut", 36, -0.5, 35.5);
             }
         }

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

delete

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

delete

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

overridevirtual

Implements art::EDProducer.

Definition at line 139 of file TPCECALAssociation_module.cc.

                                                      {
         
             // Get tracks and clusters.  If either is missing, just skip this event
             // processing.  That's not an exception
 
             art::InputTag itag(fClusterLabel, fInstanceName);
             auto ClusterHandle = e.getHandle< std::vector<gar::rec::Cluster> >(itag);
             if (!ClusterHandle) return;
 
             auto TrackHandle = e.getHandle< std::vector<gar::rec::Track> >(fTrackLabel);
             if (!TrackHandle) return;
 
             // fClocks service provides this info on event-by-event basis not at beginJob.
             maxXdisplacement = fDetProp->DriftVelocity(fDetProp->Efield(),fDetProp->Temperature())
                               *fClocks->SpillLength();
 
             // Here are the associations
             std::unique_ptr<art::Assns<gar::rec::Cluster, gar::rec::Track, gar::rec::TrackEnd>>
                             ClusterTrackAssns(new art::Assns<gar::rec::Cluster,gar::rec::Track, gar::rec::TrackEnd>);
             auto const clusterPtrMaker = art::PtrMaker<rec::Cluster>(e, ClusterHandle.id());
             auto const   trackPtrMaker = art::PtrMaker<rec::Track>  (e, TrackHandle.id());
             // Here are some intermediate info needed because both ends of a track might
             // be associated to a cluster and that ain't rite!
             std::vector<art::Ptr<gar::rec::Track>> candTrack;
             std::vector<gar::rec::TrackEnd> candEnd;
             std::vector<float> candDist;
 
 
             for (size_t iCluster=0; iCluster<ClusterHandle->size(); ++iCluster) {
                 gar::rec::Cluster cluster = (*ClusterHandle)[iCluster];
                 TVector3 clusterCenter(cluster.Position());
                 bool inECALBarrel = fGeo->PointInECALBarrel(clusterCenter);
                 // fGeo uses one coordinate system, this code uses another.
                 clusterCenter -=ItsInTulsa;
                 float yClus = clusterCenter[1];
                 float zClus = clusterCenter[2];
                 float rClus = std::hypot(zClus,yClus);  
                 float xClus = clusterCenter[0];
 
                 for (size_t iTrack=0; iTrack<TrackHandle->size(); ++iTrack) {
                     gar::rec::Track track = (*TrackHandle)[iTrack];
 
                     // Which if any ends of the track are near the outer edges of the TPC?
                     bool outside[2];    outside[0] = outside[1] = false;
                     candTrack.clear();  candEnd.clear();   candDist.clear();
 
                     if ( abs(track.Vertex()[0] -ItsInTulsa[0]) > fTrackEndXCut ||
                         std::hypot(track.Vertex()[1] -ItsInTulsa[1],
                                    track.Vertex()[2] -ItsInTulsa[2]) > fTrackEndRCut ) {
                         outside[TrackEndBeg] = true;
                     }
                     if ( abs(track.End()[0] -ItsInTulsa[0])    > fTrackEndXCut ||
                         std::hypot(track.End()[1] -ItsInTulsa[1],
                                    track.End()[2] -ItsInTulsa[2])    > fTrackEndRCut ) {
                         outside[TrackEndEnd] = true;
                     }
 
                     // Plot and cut on the distance of cluster to helix in transverse plane.
                     for (gar::rec::TrackEnd iEnd = TrackEndBeg; iEnd >= TrackEndEnd; --iEnd) {
                         if (fVerbosity>0) nEntryVsCut->Fill(0);
 
                         if (outside[iEnd]) {
                             if (fVerbosity>0) nEntryVsCut->Fill(1);
 
                             float trackPar[5];    float trackEnd[3];
                             if (iEnd==TrackEndBeg) {
                                 for (int i=0; i<5; ++i) trackPar[i] = track.TrackParBeg()[i];
                                 for (int i=0; i<3; ++i) trackEnd[i] = track.Vertex()[i];
                             } else {
                                 for (int i=0; i<5; ++i) trackPar[i] = track.TrackParEnd()[i];
                                 for (int i=0; i<3; ++i) trackEnd[i] = track.End()[i];
                             }
                             // Translate to the center-of-MPD coordinate system
                             trackPar[0] -=ItsInTulsa[1];        trackPar[1] -=ItsInTulsa[2];
                             for (int i=0; i<3; ++i) trackEnd[i] -= ItsInTulsa[i];
                             float radius = 1.0/trackPar[2];
                             float zCent = trackPar[1] - radius*sin(trackPar[3]);
                             float yCent = trackPar[0] + radius*cos(trackPar[3]);
 
                             float distRadially = std::hypot(zClus-zCent,yClus-yCent) -abs(radius);
                             distRadially = abs(distRadially);
                             if (fVerbosity>0) radClusTrack->Fill( distRadially );
                             
                             if ( distRadially > fPerpRCut ) {
                                 // This track-cluster match fails
                                 continue;
                             }
                             if (fVerbosity>0) nEntryVsCut->Fill(2);
 
                             // Require plausible extrapolation in x as well
                             float retXYZ1[3];    float retXYZ2[3];
                             TVector3 trackXYZ;
                             if (inECALBarrel) {
                                 if (fVerbosity>0) nEntryVsCut->Fill(10);
 
                                 // Extrapolate track to that radius.  Using MPD center coords.
                                 int errcode = util::TrackPropagator::PropagateToCylinder(
                                     trackPar,trackEnd,rClus, 0.0, 0.0, retXYZ1,retXYZ2);
                                 if ( errcode!=0 ) {
                                     // This track-cluster match fails.  Error code 1, there is
                                     // no intersection at all, is possible
                                     continue;
                                 }
                                 if (fVerbosity>0) nEntryVsCut->Fill(11);
                                 float extrapXerr;
                                 float transDist1 = std::hypot(retXYZ1[2]-zClus,retXYZ1[1]-yClus);
                                 float transDist2 = std::hypot(retXYZ2[2]-zClus,retXYZ2[1]-yClus);
                                 bool looneyExtrap;
                                 if (transDist1<transDist2) {
                                     trackXYZ.SetXYZ(retXYZ1[0],retXYZ1[1],retXYZ1[2]);
                                 } else {
                                     trackXYZ.SetXYZ(retXYZ2[0],retXYZ2[1],retXYZ2[2]);
                                 }
                                 trackXYZ += TVector3(ItsInTulsa);
                                 looneyExtrap =  !fGeo->PointInECALBarrel(trackXYZ)
                                              && !fGeo->PointInECALEndcap(trackXYZ);
                                 if (looneyExtrap) continue;
                                 if (fVerbosity>0) nEntryVsCut->Fill(12);
                                 trackXYZ -= TVector3(ItsInTulsa);
 
                                 extrapXerr = trackXYZ.X() -xClus;
                                 if (fVerbosity>0) xClusTrack->Fill(trackEnd[0], extrapXerr);
                                 // extrapXerr is roughly maxXdisplacement/2 for trackend at
                                 // x < -25cm and -maxXdisplacement/2 for trackend at x > 25cm.
                                 float expected_mean = 0;
                                 if (trackEnd[0]<-25) expected_mean = +maxXdisplacement/2.0;
                                 if (trackEnd[0]>+25) expected_mean = -maxXdisplacement/2.0;
                                 float cutQuantity = extrapXerr -expected_mean;
                                 if ( abs(cutQuantity) > maxXdisplacement/2.0+fBarrelXCut ) {
                                     continue;
                                 }
                                 if (fVerbosity>0) nEntryVsCut->Fill(13);
 
                             } else {
                                 if (fVerbosity>0) nEntryVsCut->Fill(20);
 
                                 // In an endcap.  How many radians in a maxXdisplacement?
                                 float radiansInDrift = trackPar[2]*maxXdisplacement
                                                       / tan(trackPar[4]);
                                 if ( abs(radiansInDrift) >= 2.0*M_PI ) {
                                     // Drat!  No distinguishing power here.  So we
                                     // count it as associated
                                     goto angleCut;
                                 }
                                 if (fVerbosity>0) nEntryVsCut->Fill(21);
 
                                 int errcode = util::TrackPropagator::PropagateToX(
                                     trackPar,trackEnd, xClus, retXYZ1);
                                 if ( errcode!=0 ) {
                                     // This track-cluster match fails.
                                     continue;
                                 }
                                 if (fVerbosity>0) nEntryVsCut->Fill(22);
                                 trackXYZ.SetXYZ(retXYZ1[0],retXYZ1[1],retXYZ1[2]);
                                 trackXYZ += TVector3(ItsInTulsa);
                                 bool looneyExtrap =  !fGeo->PointInECALBarrel(trackXYZ)
                                                   && !fGeo->PointInECALEndcap(trackXYZ);
                                 if (looneyExtrap) continue;
                                 if (fVerbosity>0) nEntryVsCut->Fill(23);
                                 trackXYZ -= TVector3(ItsInTulsa);
 
                                 // Find how many radians the closest point on the propagated
                                 // track is from the cluster along the helix, projected onto
                                 // to the perpendicular plane.  (bound to -pi,+pi range)
                                 float angClus  = std::atan2(yClus-yCent,zClus-zCent);
                                 float angXtrap = std::atan2(retXYZ1[1] -yCent,retXYZ1[2] -zCent) -angClus;
                                 // angXtrap can indeed be outside of -PI to +PI
                                 if (angXtrap > +M_PI) angXtrap -= 2.0*M_PI;
                                 if (angXtrap < -M_PI) angXtrap += 2.0*M_PI;
                                 if (fVerbosity>0) phiClusTrack->Fill(angXtrap);
                                 float extrapRphiErr = abs(radius)*angXtrap;
                                 if (fVerbosity>0) rPhiClusTrack->Fill(extrapRphiErr);
                                 if (abs(extrapRphiErr) > fEndcapRphiCut) {
                                     continue;
                                 }
                                 if (fVerbosity>0) nEntryVsCut->Fill(24);
                             }
  
                             angleCut:
                             float trackDir[3];
                             float xEval = trackXYZ.x();  // Ya not the cluster X
                             int errcode = util::TrackPropagator::DirectionX(
                                 trackPar,trackEnd, xEval,trackDir);
                             if (fVerbosity>0) nEntryVsCut->Fill(30);
                             if (errcode!=0) {
                                 // Unlikely but possible condition
                                 continue;
                             }
                             if (fVerbosity>0) nEntryVsCut->Fill(31);
                             TVector3 clusterDir;
                             clusterDir.SetXYZ(cluster.EigenVectors()[0],
                                 cluster.EigenVectors()[1],cluster.EigenVectors()[2]);
                             float dotSee = clusterDir.Dot(trackDir);
                             int nCells = cluster.CalorimeterHits().size();
                             dotClustTrack->Fill(dotSee,nCells);
                             if (nCells>=fClusterDirNhitCut && dotSee<fClusterDirCut) continue;
                             if (fVerbosity>0) nEntryVsCut->Fill(32);
 
 
 
                              // Save the cluster-track association candidate
                             art::Ptr<gar::rec::Track> const trackPtr = trackPtrMaker(iTrack);
                             candTrack.push_back(trackPtr);
                             candEnd.push_back(iEnd);
                             float howFar = std::hypot(trackXYZ.x() -xClus, trackXYZ.y() -yClus, trackXYZ.z() -zClus);
                             candDist.push_back(howFar);
                         }
                     } // end loop over 2 ends of track
 
                     // Pick which, if either, end of the track to match
                     if (candEnd.size()==0) continue;
                     art::Ptr<gar::rec::Track> trackPtrBest = candTrack[0];
                     gar::rec::TrackEnd        candEndBest  = candEnd[0];
                     if (candEnd.size()==2) {
                         if (candDist[0]>candDist[1]) {
                             trackPtrBest = candTrack[1];
                             candEndBest  = candEnd[1];
                         }
                     }
                     art::Ptr<gar::rec::Cluster> const clusterPtr = clusterPtrMaker(iCluster);
                     ClusterTrackAssns->addSingle(clusterPtr,trackPtrBest,candEndBest);
                 }
             } // end loop over clusters
 
             e.put(std::move(ClusterTrackAssns), fInstanceName);
             return;
         }