Inheritance diagram for t0::MCTruthT0Matching:

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

	MCTruthT0Matching (MCTruthT0Matching const &)=delete

	MCTruthT0Matching (MCTruthT0Matching &&)=delete

MCTruthT0Matching &	operator= (MCTruthT0Matching const &)=delete

MCTruthT0Matching &	operator= (MCTruthT0Matching &&)=delete

void	produce (art::Event &e) override

void	beginJob () override

Private Attributes
art::InputTag	fTrackModuleLabel

art::InputTag	fShowerModuleLabel

art::InputTag	fPFParticleModuleLabel

bool	fMakeT0Assns

bool	fMakePFParticleAssns

art::InputTag	fHitModuleLabel

bool	fMakeHitAssns

bool	fOverrideRealData

TTree *	fTree

int	TrackID = 0

int	TrueTrackID = 0

int	TrueTriggerType = 0

double	TrueTrackT0 = 0

int	ShowerID = 0

int	ShowerMatchID = 0

int	ShowerTriggerType = 0

double	ShowerT0 = 0

Detailed Description

Definition at line 81 of file MCTruthT0Matching_module.cc.

Constructor & Destructor Documentation

t0::MCTruthT0Matching::MCTruthT0Matching ( fhicl::ParameterSet const & p )

explicit

Definition at line 125 of file MCTruthT0Matching_module.cc.

                                                                  : EDProducer{p}
 {
   // Call appropriate produces<>() functions here.
   fTrackModuleLabel = (p.get<art::InputTag>("TrackModuleLabel"));
   fShowerModuleLabel = (p.get<art::InputTag>("ShowerModuleLabel"));
   fPFParticleModuleLabel = (p.get<art::InputTag>("PFParticleModuleLabel", "pandoraNu"));
   fMakeT0Assns = (p.get<bool>("makeT0Assns", true));
   fMakePFParticleAssns = (p.get<bool>("makePFParticleAssns", false));
 
   fMakeHitAssns = p.get<bool>("makeHitAssns", true);
   if (fMakeHitAssns) fHitModuleLabel = p.get<art::InputTag>("HitModuleLabel");
   fOverrideRealData = p.get<bool>("OverrideRealData", false);
 
   if (
     fMakeT0Assns) { // T0 assns are deprecated - this allows one to use deprecated funcionality. Added 2017-08-15. Should not be kept around forever
     std::cout << "WARNING - You are using deprecated functionality\n";
     std::cout << "MCTruthT0Matching T0 assns will be removed soon\n";
     std::cout << "set your fcl parameter makeT0Assns to false and use MCParticle direct "
                  "associations instead"
               << std::endl;
     produces<std::vector<anab::T0>>();
     produces<art::Assns<recob::Track, anab::T0>>();
     produces<art::Assns<recob::Shower, anab::T0>>();
     if (fMakePFParticleAssns) {
       produces<art::Assns<recob::PFParticle, anab::T0>>();
     } // only do PFParticles if desired by the user
   }
 
   produces<art::Assns<recob::Track, simb::MCParticle, anab::BackTrackerMatchingData>>();
   produces<art::Assns<recob::Shower, simb::MCParticle, anab::BackTrackerMatchingData>>();
   if (fMakePFParticleAssns) { // only do PFParticles if desired by the user
     produces<art::Assns<recob::PFParticle, simb::MCParticle, anab::BackTrackerMatchingData>>();
   }
 
   if (fMakeHitAssns)
     produces<art::Assns<recob::Hit, simb::MCParticle, anab::BackTrackerHitMatchingData>>();
 }

t0::MCTruthT0Matching::MCTruthT0Matching ( MCTruthT0Matching const & )

delete

t0::MCTruthT0Matching::MCTruthT0Matching ( MCTruthT0Matching && )

delete

Member Function Documentation

void t0::MCTruthT0Matching::beginJob ( )

override

Definition at line 164 of file MCTruthT0Matching_module.cc.

 {
   // Implementation of optional member function here.
   art::ServiceHandle<art::TFileService const> tfs;
   fTree = tfs->make<TTree>("MCTruthT0Matching", "MCTruthT0");
   fTree->Branch("TrueTrackT0", &TrueTrackT0, "TrueTrackT0/D");
   fTree->Branch("TrueTrackID", &TrueTrackID, "TrueTrackID/D");
 }

MCTruthT0Matching& t0::MCTruthT0Matching::operator= ( MCTruthT0Matching const & )

delete

MCTruthT0Matching& t0::MCTruthT0Matching::operator= ( MCTruthT0Matching && )

delete

void t0::MCTruthT0Matching::produce ( art::Event & e )

override

Definition at line 174 of file MCTruthT0Matching_module.cc.

 {
   if (evt.isRealData() && !fOverrideRealData) return;
 
   // Access art services...
   art::ServiceHandle<geo::Geometry const> geom;
   art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
   art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
   auto const clockData =
     art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
 
   //TrackList handle
   art::Handle<std::vector<recob::Track>> trackListHandle;
   std::vector<art::Ptr<recob::Track>> tracklist;
   if (evt.getByLabel(fTrackModuleLabel, trackListHandle))
     art::fill_ptr_vector(tracklist, trackListHandle);
 
   //ShowerList handle
   art::Handle<std::vector<recob::Shower>> showerListHandle;
   std::vector<art::Ptr<recob::Shower>> showerlist;
   if (evt.getByLabel(fShowerModuleLabel, showerListHandle))
     art::fill_ptr_vector(showerlist, showerListHandle);
 
   //PFParticleList handle
   art::Handle<std::vector<recob::PFParticle>> pfparticleListHandle;
   std::vector<art::Ptr<recob::PFParticle>> pfparticlelist;
   if (evt.getByLabel(fPFParticleModuleLabel, pfparticleListHandle))
     art::fill_ptr_vector(pfparticlelist, pfparticleListHandle);
 
   auto mcpartHandle = evt.getValidHandle<std::vector<simb::MCParticle>>("largeant");
   //  simb::MCParticle const *firstParticle = &mcpartHandle->front();
   // Create anab::T0 objects and make association with recob::Track, recob::Shower, and recob::PFParticle objects
   std::unique_ptr<std::vector<anab::T0>> T0col(new std::vector<anab::T0>);
   //  if (fMakeT0Assns){ // T0 assns are deprecated - this allows one to use deprecated funcionality. Added 2017-08-15. Should not be kept around forever
   std::unique_ptr<art::Assns<recob::Track, anab::T0>> Trackassn(
     new art::Assns<recob::Track, anab::T0>);
   std::unique_ptr<art::Assns<recob::Shower, anab::T0>> Showerassn(
     new art::Assns<recob::Shower, anab::T0>);
   std::unique_ptr<art::Assns<recob::PFParticle, anab::T0>> PFParticleassn(
     new art::Assns<recob::PFParticle, anab::T0>);
 
   // Create associations directly between MCParticle and either recob::Track, recob::Shower, or recob::PFParticle
   // These associations contain metadata summarising the quality of the match
   std::unique_ptr<art::Assns<recob::Track, simb::MCParticle, anab::BackTrackerMatchingData>>
     MCPartTrackassn(new art::Assns<recob::Track, simb::MCParticle, anab::BackTrackerMatchingData>);
   std::unique_ptr<art::Assns<recob::Shower, simb::MCParticle, anab::BackTrackerMatchingData>>
     MCPartShowerassn(
       new art::Assns<recob::Shower, simb::MCParticle, anab::BackTrackerMatchingData>);
   std::unique_ptr<art::Assns<recob::PFParticle, simb::MCParticle, anab::BackTrackerMatchingData>>
     MCPartPFParticleassn(
       new art::Assns<recob::PFParticle, simb::MCParticle, anab::BackTrackerMatchingData>);
   // Association block for the hits<-->MCParticles
   std::unique_ptr<art::Assns<recob::Hit, simb::MCParticle, anab::BackTrackerHitMatchingData>>
     MCPartHitassn(new art::Assns<recob::Hit, simb::MCParticle, anab::BackTrackerHitMatchingData>);
 
   double maxe = -1;
   double tote = 0;
   // int    trkid = -1;
   int maxtrkid = -1;
   double maxn = -1;
   double totn = 0;
   int maxntrkid = -1;
   anab::BackTrackerHitMatchingData bthmd;
 
   std::unordered_map<int, int> trkid_lookup; //indexed by geant4trkid, delivers MC particle location
 
   //if we want to make per-hit assns
   if (fMakeHitAssns) {
 
     art::Handle<std::vector<recob::Hit>> hitListHandle;
     evt.getByLabel(fHitModuleLabel, hitListHandle);
 
     if (hitListHandle.isValid()) {
 
       auto const& hitList(*hitListHandle);
       auto const& mcpartList(*mcpartHandle);
       for (size_t i_h = 0; i_h < hitList.size(); ++i_h) {
         art::Ptr<recob::Hit> hitPtr(hitListHandle, i_h);
         auto trkide_list = bt_serv->HitToTrackIDEs(clockData, hitPtr);
         struct TrackIDEinfo {
           float E;
           float NumElectrons;
         };
         std::map<int, TrackIDEinfo> trkide_collector;
         maxe = -1;
         tote = 0;
         maxtrkid = -1;
         maxn = -1;
         totn = 0;
         maxntrkid = -1;
         for (auto const& t : trkide_list) {
           trkide_collector[t.trackID].E += t.energy;
           tote += t.energy;
           if (trkide_collector[t.trackID].E > maxe) {
             maxe = trkide_collector[t.trackID].E;
             maxtrkid = t.trackID;
           }
           trkide_collector[t.trackID].NumElectrons += t.numElectrons;
           totn += t.numElectrons;
           if (trkide_collector[t.trackID].NumElectrons > maxn) {
             maxn = trkide_collector[t.trackID].NumElectrons;
             maxntrkid = t.trackID;
           }
 
           //if not found, find mc particle...
           if (trkid_lookup.find(t.trackID) == trkid_lookup.end()) {
             size_t i_p = 0;
             while (i_p < mcpartList.size()) {
               if (mcpartList[i_p].TrackId() == t.trackID) {
                 trkid_lookup[t.trackID] = (int)i_p;
                 break;
               }
               ++i_p;
             }
             if (i_p == mcpartList.size()) trkid_lookup[t.trackID] = -1;
           }
 
         } //end loop on TrackIDs
 
         //now find the mcparticle and loop back through ...
         for (auto const& t : trkide_collector) {
           int mcpart_i = trkid_lookup[t.first];
           if (mcpart_i == -1) continue; //no mcparticle here
           art::Ptr<simb::MCParticle> mcpartPtr(mcpartHandle, mcpart_i);
           bthmd.ideFraction = t.second.E / tote;
           bthmd.isMaxIDE = (t.first == maxtrkid);
           bthmd.ideNFraction = t.second.NumElectrons / totn;
           bthmd.isMaxIDEN = (t.first == maxntrkid);
           MCPartHitassn->addSingle(hitPtr, mcpartPtr, bthmd);
         }
 
       } //end loop on hits
 
     } //end if handle valid
 
   } //end if make hit/mcpart assns
 
   if (trackListHandle.isValid()) {
     //Access tracks and hits
     art::FindManyP<recob::Hit> fmtht(trackListHandle, evt, fTrackModuleLabel);
 
     size_t NTracks = tracklist.size();
 
     // Now to access MCTruth for each track...
     for (size_t iTrk = 0; iTrk < NTracks; ++iTrk) {
       TrueTrackT0 = 0;
       TrackID = 0;
       TrueTrackID = 0;
       anab::BackTrackerMatchingData btdata;
       std::vector<art::Ptr<recob::Hit>> allHits = fmtht.at(iTrk);
 
       std::map<int, double> trkide;
       for (size_t h = 0; h < allHits.size(); ++h) {
         art::Ptr<recob::Hit> hit = allHits[h];
         std::vector<sim::IDE> ides;
         std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
 
         for (size_t e = 0; e < TrackIDs.size(); ++e) {
           trkide[TrackIDs[e].trackID] += TrackIDs[e].energy;
         }
       }
       // Work out which IDE despoited the most charge in the hit if there was more than one.
       maxe = -1;
       tote = 0;
       for (std::map<int, double>::iterator ii = trkide.begin(); ii != trkide.end(); ++ii) {
         tote += ii->second;
         if ((ii->second) > maxe) {
           maxe = ii->second;
           TrackID = ii->first;
         }
       }
       btdata.cleanliness = maxe / tote;
 
       // Now have trackID, so get PdG code and T0 etc.
       const simb::MCParticle* tmpParticle = pi_serv->TrackIdToParticle_P(TrackID);
       if (!tmpParticle)
         continue; // Retain this check that the BackTracker can find the right particle
       // Now, loop through the MCParticle's myself to find the correct match
       int mcpart_i(-1);
       for (auto const particle : *mcpartHandle) {
         mcpart_i++;
         if (TrackID == particle.TrackId()) { break; }
       }
       const simb::MCParticle particle = mcpartHandle.product()->at(mcpart_i);
       TrueTrackT0 = particle.T();
       TrueTrackID = particle.TrackId();
       TrueTriggerType = 2; // Using MCTruth as trigger, so tigger type is 2.
 
       T0col->push_back(anab::T0(TrueTrackT0, TrueTriggerType, TrueTrackID, (*T0col).size()));
       //auto diff = particle - firstParticle;
       auto diff = mcpart_i; // check I have a sensible value for this counter
       if (diff >= (int)mcpartHandle->size()) {
         std::cout << "Error, the backtracker is doing weird things to your pointers!" << std::endl;
         throw std::exception();
       }
 
       art::Ptr<simb::MCParticle> mcpartPtr(mcpartHandle, mcpart_i);
       MCPartTrackassn->addSingle(tracklist[iTrk], mcpartPtr, btdata);
       if (fMakeT0Assns) { util::CreateAssn(*this, evt, *T0col, tracklist[iTrk], *Trackassn); }
       fTree->Fill();
 
     } // Loop over tracks
   }
 
   if (showerListHandle.isValid()) {
     art::FindManyP<recob::Hit> fmsht(showerListHandle, evt, fShowerModuleLabel);
     // Now Loop over showers....
     size_t NShowers = showerlist.size();
     for (size_t Shower = 0; Shower < NShowers; ++Shower) {
       ShowerMatchID = 0;
       ShowerID = 0;
       ShowerT0 = 0;
       std::vector<art::Ptr<recob::Hit>> allHits = fmsht.at(Shower);
       anab::BackTrackerMatchingData btdata;
 
       std::map<int, double> showeride;
       for (size_t h = 0; h < allHits.size(); ++h) {
         art::Ptr<recob::Hit> hit = allHits[h];
         std::vector<sim::IDE> ides;
         std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
 
         for (size_t e = 0; e < TrackIDs.size(); ++e) {
           showeride[TrackIDs[e].trackID] += TrackIDs[e].energy;
         }
       }
       // Work out which IDE despoited the most charge in the hit if there was more than one.
       maxe = -1;
       tote = 0;
       for (std::map<int, double>::iterator ii = showeride.begin(); ii != showeride.end(); ++ii) {
         tote += ii->second;
         if ((ii->second) > maxe) {
           maxe = ii->second;
           ShowerID = ii->first;
         }
       }
       btdata.cleanliness = maxe / tote;
 
       // Now have MCParticle trackID corresponding to shower, so get PdG code and T0 etc.
       const simb::MCParticle* tmpParticle = pi_serv->TrackIdToParticle_P(ShowerID);
       if (!tmpParticle)
         continue; // Retain this check that the BackTracker can find the right particle
       // Now, loop through the MCParticle's myself to find the correct match
       int mcpart_i(-1);
       for (auto const particle : *mcpartHandle) {
         mcpart_i++;
         if (ShowerID == particle.TrackId()) { break; }
       }
       const simb::MCParticle particle = mcpartHandle.product()->at(mcpart_i);
       ShowerT0 = particle.T();
       ShowerID = particle.TrackId();
       ShowerTriggerType = 2; // Using MCTruth as trigger, so tigger type is 2.
       T0col->push_back(anab::T0(ShowerT0, ShowerTriggerType, ShowerID, (*T0col).size()));
 
       auto diff = mcpart_i; // check I have a sensible value for this counter
       if (diff >= (int)mcpartHandle->size()) {
         std::cout << "Error, the backtracker is doing weird things to your pointers!" << std::endl;
         throw std::exception();
       }
       art::Ptr<simb::MCParticle> mcpartPtr(mcpartHandle, mcpart_i);
       if (fMakeT0Assns) { util::CreateAssn(*this, evt, *T0col, showerlist[Shower], *Showerassn); }
       MCPartShowerassn->addSingle(showerlist[Shower], mcpartPtr, btdata);
 
     } // Loop over showers
   }
 
   if (pfparticleListHandle.isValid()) {
     //Access pfparticles and hits
     art::FindManyP<recob::Cluster> fmcp(pfparticleListHandle, evt, fPFParticleModuleLabel);
 
     size_t NPfparticles = pfparticlelist.size();
 
     // Now to access MCTruth for each pfparticle...
     for (size_t iPfp = 0; iPfp < NPfparticles; ++iPfp) {
       TrueTrackT0 = 0;
       TrackID = 0;
       TrueTrackID = 0;
       anab::BackTrackerMatchingData btdata;
 
       std::vector<art::Ptr<recob::Hit>> allHits;
       //Get all hits through associated clusters
       std::vector<art::Ptr<recob::Cluster>> allClusters = fmcp.at(iPfp);
       art::FindManyP<recob::Hit> fmhcl(allClusters, evt, fPFParticleModuleLabel);
       for (size_t iclu = 0; iclu < allClusters.size(); ++iclu) {
         std::vector<art::Ptr<recob::Hit>> hits = fmhcl.at(iclu);
         allHits.insert(allHits.end(), hits.begin(), hits.end());
       }
 
       std::map<int, double> trkide;
       for (size_t h = 0; h < allHits.size(); ++h) {
         art::Ptr<recob::Hit> hit = allHits[h];
         std::vector<sim::IDE> ides;
         std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
 
         for (size_t e = 0; e < TrackIDs.size(); ++e) {
           trkide[TrackIDs[e].trackID] += TrackIDs[e].energy;
         }
       }
       // Work out which IDE despoited the most charge in the hit if there was more than one.
       double maxe = -1;
       double tote = 0;
       for (std::map<int, double>::iterator ii = trkide.begin(); ii != trkide.end(); ++ii) {
         tote += ii->second;
         if ((ii->second) > maxe) {
           maxe = ii->second;
           TrackID = ii->first;
         }
       }
       btdata.cleanliness = maxe / tote;
 
       // Now have trackID, so get PdG code and T0 etc.
       const simb::MCParticle* tmpParticle = pi_serv->TrackIdToParticle_P(TrackID);
       if (!tmpParticle)
         continue; // Retain this check that the BackTracker can find the right particle
       // Now, loop through the MCParticle's myself to find the correct match
       int mcpart_i(-1);
       for (auto const particle : *mcpartHandle) {
         mcpart_i++;
         if (TrackID == particle.TrackId()) { break; }
       }
       const simb::MCParticle particle = mcpartHandle.product()->at(mcpart_i);
       TrueTrackT0 = particle.T();
       TrueTrackID = particle.TrackId();
       TrueTriggerType = 2; // Using MCTruth as trigger, so tigger type is 2.
       //std::cout << "Got particle, PDG = " << particle.PdgCode() << std::endl;
 
       //std::cout << "Filling T0col with " << TrueTrackT0 << " " << TrueTriggerType << " " << TrueTrackID << " " << (*T0col).size() << std::endl;
 
       T0col->push_back(anab::T0(TrueTrackT0, TrueTriggerType, TrueTrackID, (*T0col).size()));
       //auto diff = particle - firstParticle;
       auto diff = mcpart_i; // check I have a sensible value for this counter
       if (diff >= (int)mcpartHandle->size()) {
         std::cout << "Error, the backtracker is doing weird things to your pointers!" << std::endl;
         throw std::exception();
       }
       //      art::Ptr<simb::MCParticle> mcpartPtr(mcpartHandle, particle - firstParticle);
       art::Ptr<simb::MCParticle> mcpartPtr(mcpartHandle, mcpart_i);
       //std::cout << "made MCParticle Ptr" << std::endl;
       //      const std::vertex< std::pair<double, std::string> > cleanlinessCompleteness;
       //      const std::pair<double, double> tmp(0.5, 0.5);
       //      MCPartassn->addSingle(tracklist[iPfp], mcpartPtr, tmp);
       if (fMakePFParticleAssns) {
         if (fMakeT0Assns) {
           util::CreateAssn(*this, evt, *T0col, pfparticlelist[iPfp], *PFParticleassn);
         }
         MCPartPFParticleassn->addSingle(pfparticlelist[iPfp], mcpartPtr, btdata);
       }
       fTree->Fill();
     } // Loop over tracks
   }
 
   if (fMakeT0Assns) {
     evt.put(std::move(T0col));
     evt.put(std::move(Trackassn));
     evt.put(std::move(Showerassn));
     if (fMakePFParticleAssns) { evt.put(std::move(PFParticleassn)); }
   }
   evt.put(std::move(MCPartTrackassn));
   evt.put(std::move(MCPartShowerassn));
   if (fMakePFParticleAssns) { evt.put(std::move(MCPartPFParticleassn)); }
   if (fMakeHitAssns) evt.put(std::move(MCPartHitassn));
 } // Produce