Inheritance diagram for filt::LArG4ParticleFilter:

Public Member Functions
	LArG4ParticleFilter (fhicl::ParameterSet const &pset)

virtual bool	filter (art::Event &e)

Private Member Functions
bool	IsInterestingParticle (const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const

bool	PDGCheck (const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const

bool	IsPrimaryCheck (const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const

bool	MinMomentumCheck (const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const

bool	MaxMomentumCheck (const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const

bool	StartInTPCCheck (const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const

bool	StopInTPCCheck (const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const

bool	TPCTrajLengthCheck (const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const

double	CalculateLength (const std::vector< TVector3 > &position_segment) const

void	VerifyDataMembers () const

void	VerifyElementRequest (const unsigned int index) const

Private Attributes
art::ServiceHandle < geo::Geometry const >	fGeom

const std::string	fLArG4ModuleLabel

const std::vector< int >	fInterestingPDGs

const std::vector< int >	fIsPrimary

const std::vector< double >	fParticleMinMomentum

const std::vector< double >	fParticleMaxMomentum

const std::vector< int >	fStartInTPC

const std::vector< int >	fStopInTPC

const std::vector< double >	fParticleMinTPCLength

const bool	fRequireAllInterestingParticles

std::vector< bool >	fFoundInterestingParticles

Detailed Description

Definition at line 23 of file LArG4ParticleFilter_module.cc.

Constructor & Destructor Documentation

filt::LArG4ParticleFilter::LArG4ParticleFilter ( fhicl::ParameterSet const & pset )

explicit

Member Function Documentation

double filt::LArG4ParticleFilter::CalculateLength ( const std::vector< TVector3 > & position_segment ) const

private

Definition at line 273 of file LArG4ParticleFilter_module.cc.

   {
     double length = 0;
     //Check how many points we have in the segment.  If it is one or less, there is nothing to calculate so return 0
     if (position_segment.size() <= 1) return length;
 
     //Now we need to compare every adjacent pair of positions to work out the length of this segment
     for (unsigned int i = 1; i < position_segment.size(); i++){
       length += (position_segment[i] - position_segment[i-1]).Mag();
     }
 
     return length;
   }

bool filt::LArG4ParticleFilter::filter ( art::Event & e )

virtual

Definition at line 76 of file LArG4ParticleFilter_module.cc.

   {
     //Reset the found vector
     for (unsigned int i = 0; i < fFoundInterestingParticles.size(); i++){
       fFoundInterestingParticles[i] = false;
     }
 
     //Get the vector of particles
     art::Handle<std::vector<simb::MCParticle> > particles;
     e.getByLabel(fLArG4ModuleLabel, particles);
     //Loop through the particles
     for (unsigned int part_i = 0; part_i < particles->size(); part_i++){
       //Get the particle
       const art::Ptr<simb::MCParticle> particle(particles,part_i);
       //Loop over the list of particles we want and compare it with the particle we are looking it
       for (unsigned int interest_i = 0; interest_i < fInterestingPDGs.size(); interest_i++){
         if (IsInterestingParticle(particle,interest_i)) {
           if (!fRequireAllInterestingParticles) return true; //If we only require at least one of the particles be found then we have already done our job
           fFoundInterestingParticles[interest_i] = true;
           bool foundThemAll = true;
           for (unsigned int found_i = 0; found_i < fFoundInterestingParticles.size(); found_i++){
             if (fFoundInterestingParticles[found_i] == false){
               foundThemAll = false;
               break;
             }
           }
           if (foundThemAll){
             return true;
           }
         }
       }
     }
     //Assume that the event is not worth saving
     return false;
   }

bool filt::LArG4ParticleFilter::IsInterestingParticle	(	const art::Ptr< simb::MCParticle > &	particle,
		const unsigned int	index
	)		const

private

Definition at line 112 of file LArG4ParticleFilter_module.cc.

   {
       VerifyElementRequest(index);
     //Run the checks
     if (!PDGCheck(particle,index)) return false;
     if (!IsPrimaryCheck(particle,index)) return false;
     if (!MinMomentumCheck(particle,index)) return false;
     if (!MaxMomentumCheck(particle,index)) return false;
     if (!StartInTPCCheck(particle,index)) return false;
     if (!StopInTPCCheck(particle,index)) return false;
     if (!TPCTrajLengthCheck(particle,index)) return false;
 
     return true;
   }

bool filt::LArG4ParticleFilter::IsPrimaryCheck	(	const art::Ptr< simb::MCParticle > &	particle,
		const unsigned int	index
	)		const

private

Definition at line 135 of file LArG4ParticleFilter_module.cc.

   {
       const int isPrimaryCondition(fIsPrimary[index]);
       if (isPrimaryCondition==-1) return true;
       bool particlePrimaryStatus(particle->Mother() > 0 ? false: true);
       if (particlePrimaryStatus != isPrimaryCondition) return false;
       return true;
   }

bool filt::LArG4ParticleFilter::MaxMomentumCheck	(	const art::Ptr< simb::MCParticle > &	particle,
		const unsigned int	index
	)		const

private

Definition at line 150 of file LArG4ParticleFilter_module.cc.

   {
     if (fParticleMaxMomentum[index] > 0 && particle->Momentum(0).Vect().Mag() > fParticleMaxMomentum[index]) return false;
     return true;
   }

bool filt::LArG4ParticleFilter::MinMomentumCheck	(	const art::Ptr< simb::MCParticle > &	particle,
		const unsigned int	index
	)		const

private

Definition at line 144 of file LArG4ParticleFilter_module.cc.

   {
     if (fParticleMinMomentum[index] > 0 && particle->Momentum(0).Vect().Mag() < fParticleMinMomentum[index]) return false;
     return true;
   }

bool filt::LArG4ParticleFilter::PDGCheck	(	const art::Ptr< simb::MCParticle > &	particle,
		const unsigned int	index
	)		const

private

Definition at line 127 of file LArG4ParticleFilter_module.cc.

   {
     int pdg = fInterestingPDGs[index];
     if (pdg == 0) return true; //User does not care what the PDG is
     if (particle->PdgCode() != pdg) return false;
     return true;
   }

bool filt::LArG4ParticleFilter::StartInTPCCheck	(	const art::Ptr< simb::MCParticle > &	particle,
		const unsigned int	index
	)		const

private

Definition at line 156 of file LArG4ParticleFilter_module.cc.

   {
     //Firstly check if we even care if the particle starts in the TPC or not
     int demand = fStartInTPC[index];
     if (demand == 0) return true; //We don't care if the particle starts in the TPC or not so pass the check
     //Get starting position of particle
     TLorentzVector starting_position_4vect = particle->Position(0);
     double starting_position[3];
     starting_position[0] = starting_position_4vect.X();
     starting_position[1] = starting_position_4vect.Y();
     starting_position[2] = starting_position_4vect.Z();
 
     geo::TPCID tpcid = fGeom->FindTPCAtPosition(starting_position);
     bool validtpc = tpcid.isValid;
     //Now we need to compare if we have a TPC that we started in with whether we wanted to start in a TPC at all
     if (validtpc){
       //The particle DID start in a TPC.  Now, did we WANT this to happen
       if (demand == 1) return true; //We DID want this to happen
       else return false;
     }
     else{
       //The particle did NOT start in a TPC.  Did we WANT this to happen?
       if (demand == 2) return true; //We DID want this to happen
       else return false;
     }
 
     //Assume true by default
     return true;
   }

bool filt::LArG4ParticleFilter::StopInTPCCheck	(	const art::Ptr< simb::MCParticle > &	particle,
		const unsigned int	index
	)		const

private

Definition at line 187 of file LArG4ParticleFilter_module.cc.

   {
     //Firstly check if we even care if the particle stops in the TPC or not
     int demand = fStopInTPC[index];
     if (demand == 0) return true; //We don't care if the particle stops in the TPC or not so pass the check
     //Get final position of particle
     TLorentzVector final_position_4vect = particle->Position(particle->NumberTrajectoryPoints()-1);
     double final_position[3];
     final_position[0] = final_position_4vect.X();
     final_position[1] = final_position_4vect.Y();
     final_position[2] = final_position_4vect.Z();
 
     geo::TPCID tpcid = fGeom->FindTPCAtPosition(final_position);
     bool validtpc = tpcid.isValid;
     //Now we need to compare if we have a TPC that we stopped in with whether we wanted to stop in a TPC at all
     if (validtpc){
       //The particle DID stop in a TPC.  Now, did we WANT this to happen
       if (demand == 1) return true; //We DID want this to happen
       else return false;
     }
     else{
       //The particle did NOT stop in a TPC.  Did we WANT this to happen?
       if (demand == 2) return true; //We DID want this to happen
       else return false;
     }
 
     //Assume true by default
     return true;
   }

bool filt::LArG4ParticleFilter::TPCTrajLengthCheck	(	const art::Ptr< simb::MCParticle > &	particle,
		const unsigned int	index
	)		const

private

Definition at line 217 of file LArG4ParticleFilter_module.cc.

   {
     double min_traj_length = fParticleMinTPCLength[index];
 
     //Firstly, if we don't care about the TPC trajectory length then pass the check
     if (min_traj_length < 0) return true;
 
     //Now the hard bit
     //To do this, we need to collect the sequential particle positions which are contained in the TPC into segments.  The reason for doing this is that the particle may enter a TPC, leave it and enter it again and if this isn't taken into account, the length might be grossly overestimated
     //It is easiest to store the positions in a vector of vectors
     bool OK = false;
     std::vector< std::vector<TVector3> > position_segments;
     //We are also going to need an empty vector to store in the above vector
     std::vector<TVector3> position_segment;
     //Loop through the trajectory points
     for (unsigned int i = 0; i < particle->NumberTrajectoryPoints(); i++){
       //Extract the current position of the particle
       double curr_pos[3];
       curr_pos[0] = particle->Position(i).X();
       curr_pos[1] = particle->Position(i).Y();
       curr_pos[2] = particle->Position(i).Z();
       geo::TPCID curr_tpcid = fGeom->FindTPCAtPosition(curr_pos);
       //There are a couple of things to check here.  If the particle is currently in the TPC, then we need to store that particular position.  If it is NOT in the TPC, then its either exited the TPC or has not yet entered.  If it has just exited, then the position_segment should have some positions stored in it, it which case we now need to store this segment.  If it has not yet entered the TPC, then we don't need to do anything
       //If it is currently in the TPC
       if (curr_tpcid.isValid) position_segment.push_back(particle->Position(i).Vect());
       //It has just exited the TPC
       else if (position_segment.size() > 0){
         //Store the segment
         position_segments.push_back(position_segment);
         //Now reset the segment
         position_segment.clear();
       }
       //There is nothing to do because the particle has remained outside of the TPC
     }
     //We need to check once more if the position_segment vector has been filled
     if (position_segment.size() > 0){
       position_segments.push_back(position_segment);
       position_segment.clear();
     }
     //Now lets check the length of each segment
     //Firstly, if we didn't store a segment then the particle fails the check
     if (position_segments.size() == 0) return false;
     //Now loop through the segments and check if they are above threshold
     for (unsigned int i = 0; i < position_segments.size(); i++){
       double segment_length = CalculateLength(position_segments[i]);
       if (segment_length > min_traj_length){
         //We found a track segment in the TPC which passes the length threshold so don't flag as bad
         OK = true;
         break;
       }
     }
     if (!OK) return false;
 
     return true;
   }

void filt::LArG4ParticleFilter::VerifyDataMembers ( ) const

private

Definition at line 287 of file LArG4ParticleFilter_module.cc.

   {
       if (fInterestingPDGs.size() != fIsPrimary.size())
            throw art::Exception(art::errors::Configuration) << "Config error: InterestingPDGs and IsPrimary FCL vectors are different sizes\n";
       if (fInterestingPDGs.size() != fParticleMinMomentum.size())
            throw art::Exception(art::errors::Configuration) << "Config error: InterestingPDGs and ParticleMinMomentum FCL vectors are different sizes\n";
       if (fInterestingPDGs.size() != fParticleMaxMomentum.size())
            throw art::Exception(art::errors::Configuration) << "Config error: InterestingPDGs and ParticleMaxMomentum FCL vectors are different sizes\n";
       if (fInterestingPDGs.size() != fStartInTPC.size())
            throw art::Exception(art::errors::Configuration) << "Config error: InterestingPDGs and StartInTPC FCL vectors are different sizes\n";
       if (fInterestingPDGs.size() != fStopInTPC.size())
            throw art::Exception(art::errors::Configuration) << "Config error: InterestingPDGs and StopInTPC FCL vectors are different sizes\n";
       if (fInterestingPDGs.size() != fParticleMinTPCLength.size())
            throw art::Exception(art::errors::Configuration) << "Config error: InterestingPDGs and ParticleMinTPCLength FCL vectors are different sizes\n";
 
       for (unsigned int iIsPrimary = 0; iIsPrimary < fIsPrimary.size(); iIsPrimary++)
           if (!(-1 == fIsPrimary[iIsPrimary] || 0 == fIsPrimary[iIsPrimary] || 1 == fIsPrimary[iIsPrimary]))
               throw art::Exception(art::errors::Configuration) << "Config error: Element " << iIsPrimary << " of the IsPrimary vector equals "
               << fIsPrimary[iIsPrimary] << ".  Valid options are -1 (not set), 0 (not primary) or 1 (is primary)" << std::endl;
 
       for (unsigned int iStartInTPC = 0; iStartInTPC < fStartInTPC.size(); iStartInTPC++)
           if (!(0 == fStartInTPC[iStartInTPC] || 1 == fStartInTPC[iStartInTPC] || 2 == fStartInTPC[iStartInTPC]))
               throw art::Exception(art::errors::Configuration) << "Config error: Element " << iStartInTPC << " of the StartInTPC vector equals "
               << fStartInTPC[iStartInTPC] << ".  Valid options are 0 (not set), 1 (start in TPC) or 2 (do not start in TPC)" << std::endl;
 
       for (unsigned int iStopInTPC = 0; iStopInTPC < fStopInTPC.size(); iStopInTPC++)
           if (!(0 == fStopInTPC[iStopInTPC] || 1 == fStopInTPC[iStopInTPC] || 2 == fStopInTPC[iStopInTPC]))
               throw art::Exception(art::errors::Configuration) << "Config error: Element " << iStopInTPC << " of the StopInTPC vector equals "
               << fStopInTPC[iStopInTPC] << ".  Valid options are 0 (not set), 1 (stop in TPC) or 2 (do not stop in TPC)" << std::endl;
 
       return;
   }

void filt::LArG4ParticleFilter::VerifyElementRequest ( const unsigned int index ) const

private

Definition at line 320 of file LArG4ParticleFilter_module.cc.

   {
       if (index > fInterestingPDGs.size())
               throw art::Exception(art::errors::InvalidNumber) << "Bounds error:  Requested element " << index <<" from InterestingPDGs vector which is size "
               << fInterestingPDGs.size() << "\n";
       if (index > fIsPrimary.size())
               throw art::Exception(art::errors::InvalidNumber) << "Bounds error:  Requested element " << index <<" from IsPrimary vector which is size "
               << fIsPrimary.size() << "\n";
       if (index > fParticleMinMomentum.size())
               throw art::Exception(art::errors::InvalidNumber) << "Bounds error:  Requested element " << index <<" from ParticleMinMomentum vector which is size "
               << fParticleMinMomentum.size() << "\n";
       if (index > fParticleMaxMomentum.size())
               throw art::Exception(art::errors::InvalidNumber) << "Bounds error:  Requested element " << index <<" from ParticleMaxMomentum vector which is size "
               << fParticleMaxMomentum.size() << "\n";
       if (index > fStartInTPC.size())
               throw art::Exception(art::errors::InvalidNumber) << "Bounds error:  Requested element " << index <<" from StartInTPC vector which is size "
               << fStartInTPC.size() << "\n";
       if (index > fStopInTPC.size())
               throw art::Exception(art::errors::InvalidNumber) << "Bounds error:  Requested element " << index <<" from StopInTPC vector which is size "
               << fStopInTPC.size() << "\n";
       if (index > fParticleMinTPCLength.size())
               throw art::Exception(art::errors::InvalidNumber) << "Bounds error:  Requested element " << index <<" from ParticleMinTPCLength vector which is size "
               << fParticleMinTPCLength.size() << "\n";
       return;
   }