Inheritance diagram for sbnd::StoppingCosmicIdAna:

Classes
struct	Config

Public Types
using	Parameters = art::EDAnalyzer::Table< Config >

typedef art::Handle < std::vector < recob::PFParticle > >	PFParticleHandle

typedef std::map< size_t, art::Ptr< recob::PFParticle > >	PFParticleIdMap

Public Member Functions
	StoppingCosmicIdAna (Parameters const &config)

virtual void	beginJob () override

virtual void	analyze (const art::Event &event) override

virtual void	endJob () override

Private Member Functions
void	GetPFParticleIdMap (const PFParticleHandle &pfParticleHandle, PFParticleIdMap &pfParticleMap)

Private Attributes
art::InputTag	fSimModuleLabel
	name of detsim producer More...

art::InputTag	fTpcTrackModuleLabel
	name of TPC track producer More...

art::InputTag	fCaloModuleLabel
	name of calorimetry producer More...

art::InputTag	fPandoraLabel

bool	fVerbose
	print information about what's going on More...

StoppingParticleCosmicIdAlg	spTag

TPCGeoAlg	fTPCGeo

std::vector< std::string >	types {"NuMuTrack", "CrTrack", "DirtTrack", "NuTrack", "NuMuPfp", "CrPfp", "DirtPfp", "NuPfp"}

std::map< std::string, TH1D * >	hRatioTotal

std::map< std::string, TH1D * >	hRatioTag

std::map< std::string, TH1D * >	hStopLength

std::map< std::string, TH1D * >	hNoStopLength

std::map< std::string, TH1D * >	hStopChiSq

std::map< std::string, TH1D * >	hNoStopChiSq

Detailed Description

Definition at line 52 of file StoppingCosmicIdAna_module.cc.

Member Typedef Documentation

using sbnd::StoppingCosmicIdAna::Parameters = art::EDAnalyzer::Table<Config>

Definition at line 92 of file StoppingCosmicIdAna_module.cc.

typedef art::Handle< std::vector<recob::PFParticle> > sbnd::StoppingCosmicIdAna::PFParticleHandle

Definition at line 106 of file StoppingCosmicIdAna_module.cc.

typedef std::map< size_t, art::Ptr<recob::PFParticle> > sbnd::StoppingCosmicIdAna::PFParticleIdMap

Definition at line 107 of file StoppingCosmicIdAna_module.cc.

Constructor & Destructor Documentation

sbnd::StoppingCosmicIdAna::StoppingCosmicIdAna ( Parameters const & config )

explicit

Definition at line 137 of file StoppingCosmicIdAna_module.cc.

     : EDAnalyzer(config)
     , fSimModuleLabel       (config().SimModuleLabel())
     , fTpcTrackModuleLabel  (config().TpcTrackModuleLabel())
     , fCaloModuleLabel      (config().CaloModuleLabel())
     , fPandoraLabel        (config().PandoraLabel())
     , fVerbose              (config().Verbose())
     , spTag                 (config().SPTagAlg())
   {
 
 
   } // StoppingCosmicIdAna()

Member Function Documentation

void sbnd::StoppingCosmicIdAna::analyze ( const art::Event & event )

overridevirtual

Definition at line 170 of file StoppingCosmicIdAna_module.cc.

   {
     art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
 
     // Fetch basic event info
     if(fVerbose){
       std::cout<<"============================================"<<std::endl
                <<"Run = "<<event.run()<<", SubRun = "<<event.subRun()<<", Event = "<<event.id().event()<<std::endl
                <<"============================================"<<std::endl;
     }
 
     //----------------------------------------------------------------------------------------------------------
     //                                          GETTING PRODUCTS
     //----------------------------------------------------------------------------------------------------------
     
     // Retrieve all the truth info in the events
     auto particleHandle = event.getValidHandle<std::vector<simb::MCParticle>>(fSimModuleLabel);
 
     // Retrieve the TPC tracks
     auto tpcTrackHandle = event.getValidHandle<std::vector<recob::Track>>(fTpcTrackModuleLabel);
 
     // Get track to hit associations
     art::FindManyP<recob::Hit> findManyHits(tpcTrackHandle, event, fTpcTrackModuleLabel);
     art::FindManyP<anab::Calorimetry> findManyCalo(tpcTrackHandle, event, fCaloModuleLabel);
 
     // Get PFParticles from pandora
     PFParticleHandle pfParticleHandle;
     event.getByLabel(fPandoraLabel, pfParticleHandle);
     if( !pfParticleHandle.isValid() ){
       if(fVerbose) std::cout<<"Failed to find the PFParticles."<<std::endl;
       return;
     }
     PFParticleIdMap pfParticleMap;
     this->GetPFParticleIdMap(pfParticleHandle, pfParticleMap);
     // Get PFParticle to track associations
     art::FindManyP< recob::Track > pfPartToTrackAssoc(pfParticleHandle, event, fTpcTrackModuleLabel);
 
     //----------------------------------------------------------------------------------------------------------
     //                                          TRUTH MATCHING
     //----------------------------------------------------------------------------------------------------------
     
     std::map<int, simb::MCParticle> particles;
     std::vector<int> nuParticleIds;
     std::vector<int> lepParticleIds;
     std::vector<int> dirtParticleIds;
     std::vector<int> crParticleIds;
     // Loop over the true particles
     for (auto const& particle: (*particleHandle)){
       
       // Make map with ID
       int partID = particle.TrackId();
       particles[partID] = particle;
 
       // Get MCTruth
       art::Ptr<simb::MCTruth> truth = pi_serv->TrackIdToMCTruth_P(partID);
       int pdg = std::abs(particle.PdgCode());
 
       // If origin is a neutrino
       if(truth->Origin() == simb::kBeamNeutrino){
         geo::Point_t vtx;
         vtx.SetX(truth->GetNeutrino().Nu().Vx()); vtx.SetY(truth->GetNeutrino().Nu().Vy()); vtx.SetZ(truth->GetNeutrino().Nu().Vz());
         // If neutrino vertex is not inside the TPC then call it a dirt particle
         if(!fTPCGeo.InFiducial(vtx, 0, 0)){ 
           dirtParticleIds.push_back(partID);
         }
         // If it's a primary muon
         else if(pdg==13 && particle.Mother()==0){ 
           lepParticleIds.push_back(partID);
         }
         // Other nu particles
         else{
           nuParticleIds.push_back(partID);
         }
       }
 
       // If origin is a cosmic ray
       else if(truth->Origin() == simb::kCosmicRay){
         crParticleIds.push_back(partID);
       }
       
     }
 
     //----------------------------------------------------------------------------------------------------------
     //                                    STOPPING CHI2 ANALYSIS
     //----------------------------------------------------------------------------------------------------------
 
     auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
 
     for(auto const& tpcTrack : (*tpcTrackHandle)){
 
       // Match to the true particle
       std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
       int trueId = RecoUtils::TrueParticleIDFromTotalRecoHits(clockData, hits, false);
       std::string type = "none";
       if(std::find(lepParticleIds.begin(), lepParticleIds.end(), trueId) != lepParticleIds.end()) type = "NuMuTrack";
       if(std::find(nuParticleIds.begin(), nuParticleIds.end(), trueId) != nuParticleIds.end()) type = "NuTrack";
       if(std::find(crParticleIds.begin(), crParticleIds.end(), trueId) != crParticleIds.end()) type = "CrTrack";
       if(std::find(dirtParticleIds.begin(), dirtParticleIds.end(), trueId) != dirtParticleIds.end()) type = "DirtTrack";
       if(type == "none") continue;
 
       std::vector<art::Ptr<anab::Calorimetry>> calos = findManyCalo.at(tpcTrack.ID());
       if(calos.size()==0) continue;
 
       // Only focus on muon tracks or it'll be too hard
       if(std::abs(particles[trueId].PdgCode()) != 13) continue;
 
       geo::Point_t trueEnd {particles[trueId].EndX(), particles[trueId].EndY(), particles[trueId].EndZ()};
       bool stops = false;
 
       if(fTPCGeo.InFiducial(trueEnd, 0.)) stops = true;
 
       if(stops) hStopLength[type]->Fill(tpcTrack.Length());
       else hNoStopLength[type]->Fill(tpcTrack.Length());
 
       TVector3 trueEndVec (particles[trueId].EndX(), particles[trueId].EndY(), particles[trueId].EndZ());
       TVector3 start = tpcTrack.Vertex<TVector3>();
       TVector3 end = tpcTrack.End<TVector3>();
       geo::Point_t recoEnd = tpcTrack.End();
       if((start-trueEndVec).Mag() < (end-trueEndVec).Mag()) recoEnd = tpcTrack.Vertex();
 
       double chi2 = spTag.StoppingChiSq(recoEnd, calos);
       if(stops) hStopChiSq[type]->Fill(chi2);
       else hNoStopChiSq[type]->Fill(chi2);
 
       double startChi2 = spTag.StoppingChiSq(tpcTrack.Vertex(), calos);
       double endChi2 = spTag.StoppingChiSq(tpcTrack.End(), calos);
 
       int nbins = hRatioTotal.begin()->second->GetNbinsX();
       for(int i = 0; i < nbins; i++){
         double ratioCut = hRatioTotal.begin()->second->GetBinCenter(i);
 
         hRatioTotal[type]->Fill(ratioCut);
 
         // If closest hit is below limit and track matches any hits then fill efficiency
         if((startChi2 > ratioCut && !fTPCGeo.InFiducial(tpcTrack.End(), 5.) && tpcTrack.End().Y() > 0)
            || (endChi2 > ratioCut && !fTPCGeo.InFiducial(tpcTrack.Vertex(), 5.) && tpcTrack.Vertex().Y() > 0)){
           hRatioTag[type]->Fill(ratioCut);
         }
 
       }
 
     }
 
     //Loop over the pfparticle map
     for (PFParticleIdMap::const_iterator it = pfParticleMap.begin(); it != pfParticleMap.end(); ++it){
 
       const art::Ptr<recob::PFParticle> pParticle(it->second);
       // Only look for primary particles
       if (!pParticle->IsPrimary()) continue;
       // Check if this particle is identified as the neutrino
       const int pdg(pParticle->PdgCode());
       const bool isNeutrino(std::abs(pdg) == pandora::NU_E || std::abs(pdg) == pandora::NU_MU || std::abs(pdg) == pandora::NU_TAU);
       //Find neutrino pfparticle
       if(!isNeutrino) continue;
 
       std::string type = "none";
       std::vector<recob::Track> nuTracks;
 
       // Loop over daughters of pfparticle
       for (const size_t daughterId : pParticle->Daughters()){
 
         // Get tracks associated with daughter
         art::Ptr<recob::PFParticle> pDaughter = pfParticleMap.at(daughterId);
         const std::vector< art::Ptr<recob::Track> > associatedTracks(pfPartToTrackAssoc.at(pDaughter.key()));
         if(associatedTracks.size() != 1) continue;
 
         // Get the first associated track
         recob::Track tpcTrack = *associatedTracks.front();
         nuTracks.push_back(tpcTrack);
 
         // Truth match muon tracks and pfps
         std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
         int trueId = RecoUtils::TrueParticleIDFromTotalRecoHits(clockData, hits, false);
         if(std::find(lepParticleIds.begin(), lepParticleIds.end(), trueId) != lepParticleIds.end()){ 
           type = "NuMuPfp";
         }
         else if(std::find(nuParticleIds.begin(), nuParticleIds.end(), trueId) != nuParticleIds.end()){ 
           if(type != "NuMuPfp") type = "NuPfp";
         }
         else if(std::find(dirtParticleIds.begin(), dirtParticleIds.end(), trueId) != dirtParticleIds.end()){ 
           if(type != "NuMuPfp" && type != "NuPfp") type = "DirtPfp";
         }
         else if(std::find(crParticleIds.begin(), crParticleIds.end(), trueId) != crParticleIds.end()){
           if(type != "NuMuPfp" && type != "NuPfp" && type != "DirtPfp") type = "CrPfp";
         }
       }
 
       if(nuTracks.size() == 0) continue;
 
       std::sort(nuTracks.begin(), nuTracks.end(), [](auto& left, auto& right){
                 return left.Length() > right.Length();});
 
       recob::Track tpcTrack = nuTracks[0];
       std::vector<art::Ptr<recob::Hit>> hits = findManyHits.at(tpcTrack.ID());
       int trueId = RecoUtils::TrueParticleIDFromTotalRecoHits(clockData, hits, false);
 
       std::vector<art::Ptr<anab::Calorimetry>> calos = findManyCalo.at(tpcTrack.ID());
       if(calos.size()==0) continue;
 
       // Only focus on muon tracks or it'll be too hard
       if(std::abs(particles[trueId].PdgCode()) != 13) continue;
 
       geo::Point_t trueEnd {particles[trueId].EndX(), particles[trueId].EndY(), particles[trueId].EndZ()};
       bool stops = false;
 
       if(fTPCGeo.InFiducial(trueEnd, 0.)) stops = true;
 
       if(stops) hStopLength[type]->Fill(tpcTrack.Length());
       else hNoStopLength[type]->Fill(tpcTrack.Length());
 
       TVector3 trueEndVec (particles[trueId].EndX(), particles[trueId].EndY(), particles[trueId].EndZ());
       TVector3 start = tpcTrack.Vertex<TVector3>();
       TVector3 end = tpcTrack.End<TVector3>();
       geo::Point_t recoEnd = tpcTrack.End();
       if((start-trueEndVec).Mag() < (end-trueEndVec).Mag()) recoEnd = tpcTrack.Vertex();
 
       double chi2 = spTag.StoppingChiSq(recoEnd, calos);
       if(stops) hStopChiSq[type]->Fill(chi2);
       else hNoStopChiSq[type]->Fill(chi2);
 
       double startChi2 = spTag.StoppingChiSq(tpcTrack.Vertex(), calos);
       double endChi2 = spTag.StoppingChiSq(tpcTrack.End(), calos);
 
       int nbins = hRatioTotal.begin()->second->GetNbinsX();
       for(int i = 0; i < nbins; i++){
         double ratioCut = hRatioTotal.begin()->second->GetBinCenter(i);
 
         hRatioTotal[type]->Fill(ratioCut);
 
         // If closest hit is below limit and track matches any hits then fill efficiency
         if((startChi2 > ratioCut && !fTPCGeo.InFiducial(tpcTrack.End(), 5.) && tpcTrack.End().Y() > 0)
            || (endChi2 > ratioCut && !fTPCGeo.InFiducial(tpcTrack.Vertex(), 5.) && tpcTrack.Vertex().Y() > 0)){
           hRatioTag[type]->Fill(ratioCut);
         }
 
       }
 
     }
 
 
   } // StoppingCosmicIdAna::analyze()

void sbnd::StoppingCosmicIdAna::beginJob ( )

overridevirtual

Definition at line 151 of file StoppingCosmicIdAna_module.cc.

   {
     // Access tfileservice to handle creating and writing histograms
     art::ServiceHandle<art::TFileService> tfs;
     for(auto type : types){
       hRatioTotal[type] = tfs->make<TH1D>(Form("%sRatioTotal", type.c_str()), "", 20, 1, 3);
       hRatioTag[type] = tfs->make<TH1D>(Form("%sRatioTag", type.c_str()), "", 20, 1, 3);
       hStopLength[type] = tfs->make<TH1D>(Form("%sStopLength", type.c_str()), "", 50, 0, 500);
       hNoStopLength[type] = tfs->make<TH1D>(Form("%sNoStopLength", type.c_str()), "", 50, 0, 500);
       hStopChiSq[type] = tfs->make<TH1D>(Form("%sStopChiSq", type.c_str()), "", 50, 0, 5);
       hNoStopChiSq[type] = tfs->make<TH1D>(Form("%sNoStopChiSq", type.c_str()), "", 50, 0, 5);
     }
 
     // Initial output
     if(fVerbose) std::cout<<"----------------- Stopping Particle Cosmic ID Ana Module -------------------"<<std::endl;
 
   }// StoppingCosmicIdAna::beginJob()

void sbnd::StoppingCosmicIdAna::endJob ( )

overridevirtual

Definition at line 413 of file StoppingCosmicIdAna_module.cc.

                                   {
 
   } // StoppingCosmicIdAna::endJob()

void sbnd::StoppingCosmicIdAna::GetPFParticleIdMap	(	const PFParticleHandle &	pfParticleHandle,
		PFParticleIdMap &	pfParticleMap
	)

private

Definition at line 418 of file StoppingCosmicIdAna_module.cc.

                                                                                                                       {
       for (unsigned int i = 0; i < pfParticleHandle->size(); ++i){
           const art::Ptr<recob::PFParticle> pParticle(pfParticleHandle, i);
           if (!pfParticleMap.insert(PFParticleIdMap::value_type(pParticle->Self(), pParticle)).second){
               std::cout << "  Unable to get PFParticle ID map, the input PFParticle collection has repeat IDs!" <<"\n";
           }
       }
   }