Inheritance diagram for gar::MatchingPerformance:

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

	MatchingPerformance (MatchingPerformance const &)=delete

	MatchingPerformance (MatchingPerformance &&)=delete

MatchingPerformance &	operator= (MatchingPerformance const &)=delete

MatchingPerformance &	operator= (MatchingPerformance &&)=delete

virtual void	beginJob () override

void	analyze (art::Event const &e) override

Public Member Functions inherited from art::EDAnalyzer
	EDAnalyzer (fhicl::ParameterSet const &pset)

template<typename Config >
	EDAnalyzer (Table< Config > const &config)

std::string	workerType () const

Public Member Functions inherited from art::detail::Analyzer
virtual	~Analyzer () noexcept

	Analyzer (fhicl::ParameterSet const &pset)

template<typename Config >
	Analyzer (Table< Config > const &config)

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::Observer
	~Observer () noexcept

	Observer (Observer const &)=delete

	Observer (Observer &&)=delete

Observer &	operator= (Observer const &)=delete

Observer &	operator= (Observer &&)=delete

void	registerProducts (ProductDescriptions &, ModuleDescription const &)

void	fillDescriptions (ModuleDescription const &)

fhicl::ParameterSetID	selectorConfig () const

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
void	ClearVectors ()

bool	FillVectors (art::Event const &e)

Private Attributes
float	ItsInTulsa [3]

std::string	fGeantLabel

std::string	fTrackLabel

std::string	fClusterLabel

std::string	fECALAssnLabel

std::string	fInstanceLabelECAL
	Instance name for ECAL. More...

TTree *	fTree

const gar::geo::GeometryCore *	fGeo

cheat::BackTrackerCore *	fBack

TDatabasePDG *	pdgInstance

const detinfo::DetectorProperties *	fDetProp

const detinfo::DetectorClocks *	fClocks

float	fMaxXdisplacement

int	fVerbosity

int	fClusterDirNhitCut
	Do not use cluster direction unless you have this many hits or more. More...

float	fMinPvalCut
	Do not analyse track fits with this pVal or less. More...

TH1F *	chargeFrac
	Fraction total ionization from theMCPart in any given track. More...

float	fMinIonFracCut
	Do not analyse track fits with low ionization fraction. More...

Int_t	fRun

Int_t	fSubRun

Int_t	fEvent

std::vector< ULong64_t >	fTrackIDNumber

std::vector< Float_t >	fTrackX

std::vector< Float_t >	fTrackY

std::vector< Float_t >	fTrackZ

std::vector< Float_t >	fTrackPX

std::vector< Float_t >	fTrackPY

std::vector< Float_t >	fTrackPZ

std::vector< Float_t >	fTrackPmag

std::vector< Int_t >	fTrackQ

std::vector< Float_t >	fTrackLen

std::vector< Float_t >	fTrackChi2

std::vector< Int_t >	fNTPCClustersOnTrack

std::vector< Float_t >	fTrackPval

std::vector< Double_t >	fTrackTime

std::vector< Int_t >	fMCPDG

std::vector< Int_t >	fMCPDGMother

std::vector< Float_t >	fMCPStartX

std::vector< Float_t >	fMCPStartY

std::vector< Float_t >	fMCPStartZ

std::vector< Float_t >	fMCPStartPX

std::vector< Float_t >	fMCPStartPY

std::vector< Float_t >	fMCPStartPZ

std::vector< Float_t >	fMCPStartE

std::vector< std::string >	fMCPProc

std::vector< std::string >	fMCPEndProc

std::vector< Float_t >	fMCPTime

std::vector< Float_t >	fEassocNoInTruth

std::vector< Int_t >	fNassocNoInTruth

std::vector< Float_t >	fEunassocInTruth

std::vector< Int_t >	fNunassocInTruth

std::vector< Float_t >	fEisassocInTruth

std::vector< Int_t >	fNisassocInTruth

std::vector< Float_t >	fRadClusTrackTru

std::vector< Float_t >	fRadClusTrackFal

std::vector< Float_t >	fXClusTrackOver25Tru

std::vector< Float_t >	fXClusTrackOver25Fal

std::vector< Float_t >	fXClusTrackUnder25Tru

std::vector< Float_t >	fXClusTrackUnder25Fal

std::vector< Float_t >	fRPhiClusTrackTru

std::vector< Float_t >	fRPhiClusTrackFal

std::vector< Float_t >	fDotClustTrackTru

std::vector< Float_t >	fDotClustTrackFal

std::vector< Float_t >	fDistClustTrackTru

std::vector< Float_t >	fDistClustTrackFal

Additional Inherited Members
Public Types inherited from art::EDAnalyzer
using	WorkerType = WorkerT< EDAnalyzer >

using	ModuleType = EDAnalyzer

Protected Member Functions inherited from art::Observer
std::string const &	processName () const

bool	wantAllEvents () const noexcept

bool	wantEvent (ScheduleID id, Event const &e) const

Handle< TriggerResults >	getTriggerResults (Event const &e) const

	Observer (fhicl::ParameterSet const &config)

	Observer (std::vector< std::string > const &select_paths, std::vector< std::string > const &reject_paths, fhicl::ParameterSet const &config)

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 68 of file MatchingPerformance_module.cc.

Constructor & Destructor Documentation

gar::MatchingPerformance::MatchingPerformance ( fhicl::ParameterSet const & p )

explicit

Definition at line 195 of file MatchingPerformance_module.cc.

                                                                        : EDAnalyzer(p) {
 
   fGeantLabel        = p.get<std::string>("GEANTLabel",       "geant");
   fTrackLabel        = p.get<std::string>("TrackLabel",       "track");
   fClusterLabel      = p.get<std::string>("ClusterLabel",     "calocluster");
   fECALAssnLabel     = p.get<std::string>("ECALAssnLabel",    "trkecalassn");
   fInstanceLabelECAL = p.get<std::string>("InstanceLabelCalo","ECAL");
   fVerbosity         = p.get<int>        ("Verbosity",         0);
   fClusterDirNhitCut = p.get<int>        ("ClusterDirNhitCut", 5);
   fMinPvalCut        = p.get<float>      ("MinPvalCut",        0.0);
   fMinIonFracCut     = p.get<float>      ("MinIonFracCut",     0.50);
 
   pdgInstance = TDatabasePDG::Instance();
 
   consumesMany<std::vector<simb::MCTruth> >();
   consumes<std::vector<simb::MCParticle> >(fGeantLabel);
   consumes<std::vector<rec::Track> >(fTrackLabel);
   consumes<std::vector<rec::Cluster> >(fClusterLabel);
   consumes<art::Assns<rec::Cluster, rec::Track>>(fECALAssnLabel);
 
   return;
 } // end constructor

gar::MatchingPerformance::MatchingPerformance ( MatchingPerformance const & )

delete

gar::MatchingPerformance::MatchingPerformance ( MatchingPerformance && )

delete

Member Function Documentation

void gar::MatchingPerformance::analyze ( art::Event const & e )

overridevirtual

Implements art::EDAnalyzer.

Definition at line 302 of file MatchingPerformance_module.cc.

                                                            {
 
   // Need this here because the clock service is invoked at preProcessEvent.
   fDetProp = gar::providerFrom<detinfo::DetectorPropertiesService>();
   fClocks  = gar::providerFrom<detinfo::DetectorClocksServiceGAr>();
   fMaxXdisplacement =
     fDetProp->DriftVelocity(fDetProp->Efield(),fDetProp->Temperature())
     *fClocks->SpillLength();
 
   // Is the backtracker good to go?  Need to fill fBack on per-event basis
   // not in beginJob()
   cheat::BackTrackerCore const* const_bt = gar::providerFrom<cheat::BackTracker>();
   fBack = const_cast<cheat::BackTrackerCore*>(const_bt);
   if ( !fBack->HasTracks() || !fBack->HasClusters() ) {
     throw cet::exception("MatchingPerformance") << " Not enough backtracker info"
                                                 << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
   }
 
   ClearVectors();
   if ( FillVectors(event) ) {
     fTree->Fill();
   }
 
   return;
 }

void gar::MatchingPerformance::beginJob ( )

overridevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 225 of file MatchingPerformance_module.cc.

                                       {
 
   fGeo = gar::providerFrom<geo::GeometryGAr>();
   ItsInTulsa[0] = fGeo->TPCXCent();        // 1 S Boston Ave
   ItsInTulsa[1] = fGeo->TPCYCent();
   ItsInTulsa[2] = fGeo->TPCZCent();
 
 
 
   art::ServiceHandle<art::TFileService> tfs;
   fTree = tfs->make<TTree>("GArPerfTree","GArPerfTree");
 
   if (fVerbosity>0) {
     chargeFrac = tfs->make<TH1F>("chargeFrac", "Fraction of MCPart's ionization in track",
                                  101,0.00,1.01);
   }
 
 
 
   fTree->Branch("Run",                    &fRun,              "Run/I");
   fTree->Branch("SubRun",                 &fSubRun,           "SubRun/I");
   fTree->Branch("Event",                  &fEvent,            "Event/I");
 
   fTree->Branch("TrackIDNumber",          &fTrackIDNumber);
   fTree->Branch("TrackX",                 &fTrackX);
   fTree->Branch("TrackY",                 &fTrackY);
   fTree->Branch("TrackZ",                 &fTrackZ);
   fTree->Branch("TrackPX",                &fTrackPX);
   fTree->Branch("TrackPY",                &fTrackPY);
   fTree->Branch("TrackPZ",                &fTrackPZ);
   fTree->Branch("TrackPmag",              &fTrackPmag);
   fTree->Branch("TrackQ",                 &fTrackQ);
   fTree->Branch("TrackLen",               &fTrackLen);
   fTree->Branch("TrackChi2",              &fTrackChi2);
   fTree->Branch("NTPCClustersOnTrack",    &fNTPCClustersOnTrack);
   fTree->Branch("TrackPval",              &fTrackPval);
   fTree->Branch("TrackTime",              &fTrackTime);
 
   fTree->Branch("PDG",                    &fMCPDG);
   fTree->Branch("PDGMother",              &fMCPDGMother);
   fTree->Branch("MCPStartX",              &fMCPStartX);
   fTree->Branch("MCPStartY",              &fMCPStartY);
   fTree->Branch("MCPStartZ",              &fMCPStartZ);
   fTree->Branch("MCPStartPX",             &fMCPStartPX);
   fTree->Branch("MCPStartPY",             &fMCPStartPY);
   fTree->Branch("MCPStartPZ",             &fMCPStartPZ);
   fTree->Branch("MCPTime",                &fMCPTime);
 
   fTree->Branch("EassocNoInTruth",        &fEassocNoInTruth);
   fTree->Branch("NassocNoInTruth",        &fNassocNoInTruth);
   fTree->Branch("EunassocInTruth",        &fEunassocInTruth);
   fTree->Branch("NunassocInTruth",        &fNunassocInTruth);
   fTree->Branch("EisassocInTruth",        &fEisassocInTruth);
   fTree->Branch("NisassocInTruth",        &fNisassocInTruth);
 
   fTree->Branch("RadClusTrackTru",        &fRadClusTrackTru);
   fTree->Branch("RadClusTrackFal",        &fRadClusTrackFal);
   fTree->Branch("XClusTrackOver25Tru",    &fXClusTrackOver25Tru);
   fTree->Branch("XClusTrackOver25Fal",    &fXClusTrackOver25Fal);
   fTree->Branch("XClusTrackUnder25Tru",   &fXClusTrackUnder25Tru);
   fTree->Branch("XClusTrackUnder25Fal",   &fXClusTrackUnder25Fal);
   fTree->Branch("RPhiClusTrackTru",       &fRPhiClusTrackTru);
   fTree->Branch("RPhiClusTrackFal",       &fRPhiClusTrackFal);
   fTree->Branch("DotClustTrackTru",       &fDotClustTrackTru);
   fTree->Branch("DotClustTrackFal",       &fDotClustTrackFal);
   fTree->Branch("DistClustTrackTru",      &fDistClustTrackTru);
   fTree->Branch("DistClustTrackFal",      &fDistClustTrackFal);
   return;
 }  // End of :MatchingPerformance::beginJob

void gar::MatchingPerformance::ClearVectors ( )

private

Definition at line 335 of file MatchingPerformance_module.cc.

                                           {
   // clear out all our vectors
 
   fTrackIDNumber.clear();
   fTrackX.clear();
   fTrackY.clear();
   fTrackZ.clear();
   fTrackPX.clear();
   fTrackPY.clear();
   fTrackPZ.clear();
   fTrackPmag.clear();
   fTrackQ.clear();
   fTrackLen.clear();
   fTrackChi2.clear();
   fNTPCClustersOnTrack.clear();
   fTrackPval.clear();
   fTrackTime.clear();
 
   fMCPDG.clear();
   fMCPDGMother.clear();
   fMCPStartX.clear();
   fMCPStartY.clear();
   fMCPStartZ.clear();
   fMCPStartPX.clear();
   fMCPStartPY.clear();
   fMCPStartPZ.clear();
   fMCPTime.clear();
 
   fEassocNoInTruth.clear();
   fNassocNoInTruth.clear();
   fEunassocInTruth.clear();
   fNunassocInTruth.clear();
   fEisassocInTruth.clear();
   fNisassocInTruth.clear();
 
   fRadClusTrackTru.clear();
   fRadClusTrackFal.clear();
   fXClusTrackOver25Tru.clear();
   fXClusTrackOver25Fal.clear();
   fXClusTrackUnder25Tru.clear();
   fXClusTrackUnder25Fal.clear();
   fRPhiClusTrackTru.clear();
   fRPhiClusTrackFal.clear();
   fDotClustTrackTru.clear();
   fDotClustTrackFal.clear();
   fDistClustTrackTru.clear();
   fDistClustTrackFal.clear();
 
   return;
 } // end :MatchingPerformance::ClearVectors

bool gar::MatchingPerformance::FillVectors ( art::Event const & e )

private

Want to be able to know what you've matched back to at MCTruth

Definition at line 393 of file MatchingPerformance_module.cc.

                                                               {
 
   // =============  Get art handles ==========================================
   // Get handles for Tracks, clusters, associations between them
   auto TrackHandle = event.getHandle< std::vector<rec::Track> >(fTrackLabel);
   if (!TrackHandle) {
     throw cet::exception("MatchingPerformance") << " No rec::Track"
                                                 << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
   }
 
   art::InputTag ecalclustertag(fClusterLabel, fInstanceLabelECAL);
   auto ClusterHandle = event.getHandle< std::vector<rec::Cluster> >(ecalclustertag);
   if (!ClusterHandle) {
     throw cet::exception("MatchingPerformance") << " No rec::Cluster branch."
                                                 << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
   }
 
   art::FindMany<rec::Cluster, rec::TrackEnd>* findManyTrackEndCAL = NULL;
   art::InputTag ecalassntag(fECALAssnLabel, fInstanceLabelECAL);
   findManyTrackEndCAL = new art::FindMany<rec::Cluster, rec::TrackEnd>
     (TrackHandle,event,ecalassntag);
   if ( !findManyTrackEndCAL->isValid() ) {
     throw cet::exception("MatchingPerformance") << " Bad TrackEnd-ECAL Assn."
                                                 << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
   }
  
   // Get handles for MCinfo, also good for MCPTrajectory.  Want to get all
   // MCTruths, regardless of generator label.
 
   auto mctruthHandles = event.getMany<std::vector<simb::MCTruth> >();
 
   if (mctruthHandles.size()!=1) {
     throw cet::exception("MatchingPerformance") << " Need just 1 simb::MCTruth"
                                                 << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
   }
 
   auto MCPHandle = event.getHandle< std::vector<simb::MCParticle> >(fGeantLabel);
   if (!MCPHandle) {
     throw cet::exception("MatchingPerformance") << " No simb::MCParticle"
                                                 << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
   }
 
 
 
   // Need something to look at!
   if ( TrackHandle->size()==0 || ClusterHandle->size()==0 ) return false;
 
   // Make the usual map for the MC info
   typedef int TrkId;
   std::unordered_map<TrkId, Int_t> TrackIdToIndex;
   Int_t index = 0;
   for ( auto const& mcp : (*MCPHandle) ) {
     int TrackId = mcp.TrackId();
     TrackIdToIndex[TrackId] = index++;
   }
 
   // Will need a vector of art::Ptrs to all the ECAL clusters for the
   // backtracker calls.
   std::vector<art::Ptr<rec::Cluster>> ClusterGrabber;
   art::PtrMaker<rec::Cluster> makeClusterPtr(event,ClusterHandle.id());
   for (size_t iCluster=0; iCluster<ClusterHandle->size(); ++iCluster) {
     art::Ptr<rec::Cluster> aPtr = makeClusterPtr(iCluster);
     ClusterGrabber.push_back(aPtr);
   }
   // Will need a vector<art::Ptr<rec::Track>> of tracks also
   std::vector<art::Ptr<rec::Track>> TrackGrabber;
   art::PtrMaker<rec::Track> makeTrackPtr(event,TrackHandle.id());
   for (size_t iTrack=0; iTrack<TrackHandle->size(); ++iTrack) {
     art::Ptr<rec::Track> aPtr = makeTrackPtr(iTrack);
     TrackGrabber.push_back(aPtr);
   }
 
 
 
 
 
   // =============  Pull art's handles =======================================
   fRun    = event.run();
   fSubRun = event.subRun();
   fEvent  = event.id().event();
 
 
 
   // Require tracks match to right kind of MC first
   struct niceNice {rec::Track recoTrk; simb::MCParticle simiTrk;};
   std::vector<niceNice> niceTracks;
   for ( auto const& track : (*TrackHandle) ) {
 
     /// Want to be able to know what you've matched back to at MCTruth
     std::vector<std::pair<simb::MCParticle*,float>> MCsfromTrack;
     gar::rec::Track* nonconst_track = const_cast<gar::rec::Track*>(&track);
     MCsfromTrack = fBack->TrackToMCParticles(nonconst_track);
     int nMCsfromTrack = MCsfromTrack.size();
     if (nMCsfromTrack!=1) continue;
 
     simb::MCParticle theMCPart  = *(MCsfromTrack[0].first);
     TParticlePDG* particle = pdgInstance->GetParticle(theMCPart.PdgCode());
     if (particle==NULL         ) continue;   // What causes this err? If anything.
     if (particle->Charge()==0.0) continue;
     if (particle->Stable()==0  ) continue;
 
     // Does theMCParticle go out to the ECAL?
     int nTraj = theMCPart.NumberTrajectoryPoints();
     TVector3 doink;            bool whackThatECAL = false;
     for (int iTraj=0; iTraj<nTraj; ++iTraj) {
       doink.SetXYZ(theMCPart.Vx(iTraj),theMCPart.Vy(iTraj),theMCPart.Vz(iTraj));
       if ( fGeo->PointInECALBarrel(doink) || fGeo->PointInECALEndcap(doink) ) {
         whackThatECAL = true;
         break;
       }
     }
     if (!whackThatECAL) continue;
 
     // This track is nice
     niceNice tmp = {track, theMCPart};
     niceTracks.push_back(tmp);
   }
 
 
 
   // Anti-stubbiness cut.
   std::vector<niceNice>::iterator iNiceTrk=niceTracks.begin();
   for (; iNiceTrk!=niceTracks.end(); ++iNiceTrk) {
     simb::MCParticle theMCPart = iNiceTrk->simiTrk;
     std::vector<art::Ptr<rec::Track>> tracksFromSameMCP =
       fBack->MCParticleToTracks(&theMCPart,TrackGrabber);
 
     // Get ionization fraction for every track with contribution from this MCP
     std::vector<float> ionzFromSameMCP;    ionzFromSameMCP.resize(tracksFromSameMCP.size());
     for (size_t iTrkFromSame=0; iTrkFromSame<tracksFromSameMCP.size(); ++iTrkFromSame) {
       rec::Track tmpTrk = *(tracksFromSameMCP[iTrkFromSame]);
       std::vector<art::Ptr<rec::Hit>> hitList = fBack->TrackToHits(&tmpTrk);
       float thisTracksI = 0;
       for (auto& ptrHit : hitList) thisTracksI += ptrHit->Signal();
       ionzFromSameMCP[iTrkFromSame] = thisTracksI;
     }
     float ionzTotSameMCP = 0;
     for (size_t iIonSame=0; iIonSame<ionzFromSameMCP.size(); ++iIonSame) 
       ionzTotSameMCP += ionzFromSameMCP[iIonSame];
     for (size_t iIonSame=0; iIonSame<ionzFromSameMCP.size(); ++iIonSame) {
       ionzFromSameMCP[iIonSame] /= ionzTotSameMCP;
       if (fVerbosity>0) chargeFrac->Fill(ionzFromSameMCP[iIonSame]);
     }
 
     // Cut on ion fraction in niceTracks
     for (size_t iTrkFromSame=0; iTrkFromSame<tracksFromSameMCP.size(); ++iTrkFromSame) {
       if ( *(tracksFromSameMCP[iTrkFromSame])==iNiceTrk->recoTrk
            && ionzFromSameMCP[iTrkFromSame] <fMinIonFracCut) {
         niceTracks.erase(iNiceTrk);
         --iNiceTrk;
         break;
       }
     }
   }
 
 
   iNiceTrk = niceTracks.begin();
   std::vector<niceNice>::iterator jNiceTrk = iNiceTrk +1;
   for (; iNiceTrk!=niceTracks.end(); ++iNiceTrk) {
     for (; jNiceTrk!=niceTracks.end(); ++jNiceTrk) {
       if ( iNiceTrk->simiTrk.TrackId() == jNiceTrk->simiTrk.TrackId() ){
         niceTracks.erase(jNiceTrk);
         --jNiceTrk;
       }
     }
   }
 
 
 
   for (size_t iNiceTrk_local=0; iNiceTrk_local<niceTracks.size(); ++iNiceTrk_local) {
     rec::Track       track      = niceTracks[iNiceTrk_local].recoTrk;
     simb::MCParticle theMCPart  = niceTracks[iNiceTrk_local].simiTrk;
 
 
 
     // Examine matched clusters on downstream ends of the track unless 
     // MC vertex in gas - don't test the matching quality then.
     // Note: rec::TrackEndBeg is at the front of the deque here!
     std::deque<rec::TrackEnd> endList = {rec::TrackEndBeg,rec::TrackEndEnd};
 
     TVector3 positionMCP = theMCPart.Position(0).Vect();
     if (fGeo->PointInGArTPC(positionMCP)) {
       // Pick the end where the particle ended.
       // Don't include the drift distance in matching.
       float distStart = std::hypot(track.Vertex()[1] -positionMCP[1],
                                    track.Vertex()[2] -positionMCP[2]);
       float distEnd   = std::hypot(track.End()[1]    -positionMCP[1],
                                    track.End()[2]    -positionMCP[2]);
       if ( distStart < distEnd ) {
         endList.pop_front();
       } else {
         endList.pop_back();
       }
     }
 
     // Having selected now the end we wil extrapolate,
     std::deque<rec::TrackEnd>::iterator iEnd = endList.begin();
     // Record info about the track at extrapolated end - need trackPar, 
     // TrackEnd later and they must be in MPD coordinates.
     fTrackIDNumber.push_back(track.getIDNumber());
     Float_t saveChi2;
     float trackPar[5];      float trackEnd[3];
 
     if ( (*iEnd)==rec::TrackEndBeg ) {
       for (size_t i=0; i<5; ++i) trackPar[i] = track.TrackParBeg()[i];
       for (size_t i=0; i<3; ++i) trackEnd[i] = track.Vertex()[i];
       fTrackX.push_back   ( track.Vertex()[0] -ItsInTulsa[0] );
       fTrackY.push_back   ( track.Vertex()[1] -ItsInTulsa[1] );
       fTrackZ.push_back   ( track.Vertex()[2] -ItsInTulsa[2] );
       fTrackPX.push_back  (-track.Momentum_beg()*track.VtxDir()[0] );
       fTrackPY.push_back  (-track.Momentum_beg()*track.VtxDir()[1] );
       fTrackPZ.push_back  (-track.Momentum_beg()*track.VtxDir()[2] );
       fTrackPmag.push_back( track.Momentum_beg() );
       fTrackQ.push_back   ( track.ChargeEnd() );
       fTrackLen.push_back ( track.LengthBackward() );
       saveChi2 = track.ChisqBackward();
       fTrackChi2.push_back( saveChi2 );
     } else {
       for (size_t i=0; i<5; ++i) trackPar[i] = track.TrackParEnd()[i];
       for (size_t i=0; i<3; ++i) trackEnd[i] = track.End()[i];
       fTrackX.push_back   ( track.End()[0] -ItsInTulsa[0] );
       fTrackY.push_back   ( track.End()[1] -ItsInTulsa[1] );
       fTrackZ.push_back   ( track.End()[2] -ItsInTulsa[2] );
       fTrackPX.push_back  (-track.Momentum_end()*track.EndDir()[0] );
       fTrackPY.push_back  (-track.Momentum_end()*track.EndDir()[1] );
       fTrackPZ.push_back  (-track.Momentum_end()*track.EndDir()[2] );
       fTrackPmag.push_back( track.Momentum_end() );
       fTrackQ.push_back   ( track.ChargeBeg() );
       fTrackLen.push_back ( track.LengthForward() );        
       saveChi2 = track.ChisqForward();
       fTrackChi2.push_back( saveChi2);
 
     }
 
     // Apply the pVal cut
     Int_t nHits = track.NHits();
     fNTPCClustersOnTrack.push_back(nHits);
     Float_t pVal = ROOT::Math::chisquared_cdf_c(saveChi2,nHits-5);
     if (pVal <= fMinPvalCut) return false;
     fTrackPval.push_back(pVal);
     fTrackTime.push_back(track.Time());
 
     trackPar[0] -=ItsInTulsa[1];        trackPar[1] -=ItsInTulsa[2];
     for (int i=0; i<3; ++i) trackEnd[i] -= ItsInTulsa[i];
 
     // Record some info about the matching MC
     fMCPDG.push_back(theMCPart.PdgCode());
     TrkId momTrkId = theMCPart.Mother();
     int momPDG = 0;
     if (momTrkId>0) {
       int momIndex = TrackIdToIndex[momTrkId];
       momPDG   = (*MCPHandle).at(momIndex).PdgCode();
     }
     fMCPDGMother.push_back(momPDG);
 
     fMCPStartX.push_back(positionMCP.X());
     fMCPStartY.push_back(positionMCP.Y());
     fMCPStartZ.push_back(positionMCP.Z());
     const TLorentzVector& momentumMCP = theMCPart.Momentum(0);
     fMCPStartPX.push_back(momentumMCP.Px());
     fMCPStartPY.push_back(momentumMCP.Py());
     fMCPStartPZ.push_back(momentumMCP.Pz());
     fMCPTime.push_back(theMCPart.T());    
 
 
 
     // Now, look for the ECAL-track matching info
     std::vector<rec::Cluster> clustersOnMCP;
     for (rec::Cluster cluster : *ClusterHandle) {
       if (fBack->MCParticleCreatedCluster(&theMCPart,&cluster)) {
         clustersOnMCP.push_back(cluster);
       }
       if (fBack->ClusterCreatedMCParticle(&theMCPart,&cluster)) {
         clustersOnMCP.push_back(cluster);
       }
     }
 
     // Because of how IsForebearOf works, the ionization in the cluster
     // which is from theMCPart makes theMCPart both a parent and descendent
     // of the MCParticles in the cluster.  So we must remove duplicate 
     // clustersOnMCP, which means sorting them, which means we need a 
     // lambda function, etc. 
     std::sort(clustersOnMCP.begin(),clustersOnMCP.end(),
               [](rec::Cluster a,rec::Cluster b){return a.getIDNumber() > b.getIDNumber();});
     std::vector<rec::Cluster>::iterator itr =
       std::unique(clustersOnMCP.begin(),clustersOnMCP.end());
     clustersOnMCP.resize(std::distance(clustersOnMCP.begin(),itr));
 
 
 
 
     float eAssocNoInTruth = 0.0;            int nAssocNoInTruth = 0;
     float eUnassocInTruth = 0.0;            int nUnassocInTruth = 0;
     float eIsassocInTruth = 0.0;            int nIsassocInTruth = 0;
     size_t   nCALedClusts =   0;
     std::vector<rec::Cluster  const*> clustersFromAssn;
     std::vector<rec::TrackEnd const*> trackEndsFromAssn;
 
     // The art::FindMany object findManyTrackEndCAL is indexed over 
     // all the tracks in *TrackHandle (which is bigger than niceTracks).
     size_t iTrack=0;
     for (; iTrack<TrackHandle->size(); ++iTrack) {
       if ( (*TrackHandle)[iTrack] == track ) break;
     }
 
     // nCALedClusts.size() == trackEndsFromAssn.size() and is the same
     // number regardless of std::deque<rec::TrackEnd>::iterator iEnd but
     // the code reads better here
     findManyTrackEndCAL->get(iTrack, clustersFromAssn,trackEndsFromAssn);
     nCALedClusts = clustersFromAssn.size();
 
     for ( auto const& cluster : (*ClusterHandle) ) {
 
       // recompute the matching variables, similar to (if not exactly
       // the same as) in Reco/TPCECALAssociation_module.cc  They will be:
       // inECALBarrel, distRadially, Over25, cutQuantity, dotSee
       float radius = 1.0/trackPar[2];
       float zCent = trackPar[1] - radius*sin(trackPar[3]);
       float yCent = trackPar[0] + radius*cos(trackPar[3]);
       float xClus = cluster.Position()[0] -ItsInTulsa[0];
       float yClus = cluster.Position()[1] -ItsInTulsa[1];
       float zClus = cluster.Position()[2] -ItsInTulsa[2];
       float rClus = std::hypot(zClus,yClus);
       bool  inECALBarrel = fGeo->PointInECALBarrel(cluster.Position());
 
       float distRadially = std::hypot(zClus-zCent,yClus-yCent) -abs(radius);
       distRadially = abs(distRadially);
 
       float retXYZ1[3];    float retXYZ2[3];    TVector3 trackXYZ;
       float cutQuantity = -1.0;      bool over25 = false;
 
       if (inECALBarrel) {
         int errcode = util::TrackPropagator::PropagateToCylinder(
                                                                  trackPar,trackEnd,rClus, 0.0, 0.0, retXYZ1,retXYZ2);
         if ( errcode==0 ) {
           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);
           trackXYZ -= TVector3(ItsInTulsa);
           if (!looneyExtrap) {
             extrapXerr = trackXYZ.X() -xClus;
             float expected_mean = 0;
             if (trackEnd[0]<-25) {
               expected_mean = +fMaxXdisplacement/2.0;
               over25 = true;
             }
             if (trackEnd[0]>+25) {
               expected_mean = -fMaxXdisplacement/2.0;
               over25 = true;
             }
             cutQuantity = abs(extrapXerr -expected_mean);
           }
         }
       } else {
         // In an endcap.  How many radians in a fMaxXdisplacement?
         float radiansInDrift = trackPar[2]*fMaxXdisplacement
           / tan(trackPar[4]);
         if ( abs(radiansInDrift) >= 2.0*M_PI ) goto angleCut;
         int errcode = util::TrackPropagator::PropagateToX(
                                                           trackPar,trackEnd, xClus, retXYZ1);
         if ( errcode==0 ) {
           trackXYZ.SetXYZ(retXYZ1[0],retXYZ1[1],retXYZ1[2]);
           trackXYZ += TVector3(ItsInTulsa);
           bool looneyExtrap =  !fGeo->PointInECALBarrel(trackXYZ)
             && !fGeo->PointInECALEndcap(trackXYZ);
           trackXYZ -= TVector3(ItsInTulsa);
           if (!looneyExtrap) {
                             
             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;
             cutQuantity = abs(radius*angXtrap);
           }
         }
       }
 
     angleCut:
       float trackDir[3];
       float dotSee = -2.0;
       int nCells = cluster.CalorimeterHits().size();
       if (nCells >= fClusterDirNhitCut) {
         float xEval = trackXYZ.x();  // Ya not the cluster X
         int errcode = util::TrackPropagator::DirectionX(
                                                         trackPar,trackEnd, xEval,trackDir);
         if (errcode==0) {
           TVector3 clusterDir;
           clusterDir.SetXYZ(cluster.EigenVectors()[0],
                             cluster.EigenVectors()[1],cluster.EigenVectors()[2]);
           dotSee = clusterDir.Dot(trackDir);
         }
       }
 
       // Look at distance from extrapolated track to cluster.
       float dist3dXtrap = std::hypot(
                                      trackXYZ.x()-xClus, trackXYZ.y()-yClus, trackXYZ.z()-zClus);
 
 
 
       bool clusterOnTrack = false;
       for (size_t iCALedClust=0; iCALedClust<nCALedClusts; ++iCALedClust) {
         if (*trackEndsFromAssn[iCALedClust] != (*iEnd) ) continue;
         if (*clustersFromAssn[iCALedClust] == cluster) {
           clusterOnTrack = true;
           break;
         }
       }
 
       bool clusterOnMCP = false;
       for (size_t iClusOnMCP = 0; iClusOnMCP<clustersOnMCP.size(); ++iClusOnMCP) {
         if ( clustersOnMCP[iClusOnMCP] == cluster ) {
           clusterOnMCP = true;
           break;
         }
       }
 
 
 
  
       if (!clusterOnMCP ) {
         fRadClusTrackFal.push_back(distRadially);
         if (inECALBarrel) {
           if (over25) {
             fXClusTrackOver25Fal.push_back(cutQuantity);
           } else {
             fXClusTrackUnder25Fal.push_back(cutQuantity);
           }
         } else {
           fRPhiClusTrackFal.push_back(cutQuantity);
         }
         fDotClustTrackFal.push_back(dotSee);
         fDistClustTrackFal.push_back(dist3dXtrap);
 
         if ( clusterOnTrack ) {
           eAssocNoInTruth += cluster.Energy();
           nAssocNoInTruth += cluster.CalorimeterHits().size();
         }
       } else {
         fRadClusTrackTru.push_back(distRadially);
         if (inECALBarrel) {
           if (over25) {
             fXClusTrackOver25Tru.push_back(cutQuantity);
           } else {
             fXClusTrackUnder25Tru.push_back(cutQuantity);
           }
         } else {
           fRPhiClusTrackTru.push_back(cutQuantity);
         }
         fDotClustTrackTru.push_back(dotSee);
         fDistClustTrackTru.push_back(dist3dXtrap);
 
         if (!clusterOnTrack && clusterOnMCP ) {
           eUnassocInTruth += cluster.Energy();
           nUnassocInTruth += cluster.CalorimeterHits().size();
         } else {
           eIsassocInTruth += cluster.Energy();
           nIsassocInTruth += cluster.CalorimeterHits().size();
         }
       }
     }
 
     fEassocNoInTruth.push_back(eAssocNoInTruth);
     fNassocNoInTruth.push_back(nAssocNoInTruth);
     fEunassocInTruth.push_back(eUnassocInTruth);
     fNunassocInTruth.push_back(nUnassocInTruth);
     fEisassocInTruth.push_back(eIsassocInTruth);
     fNisassocInTruth.push_back(nIsassocInTruth);
 
   } // end loop over TrackHandle
 
   return true;
 } // end MatchingPerformance::FillVectors