LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
LArG4ParticleFilter_module.cc
Go to the documentation of this file.
1 //Author: Dom Brailsford
2 //An ART filter that selects events which contain 'interesting' simb::MCParticles. Interesting is user-defined.
3 //The filter parameters are FHICL configurable but can select based on, momentum, particle flavour, TPC trajectory length and whether the particle stops/starts in the TPC.
4 //The user can ask for events which contain multiple types of interesting particles, each with their own requirements.
5 
6 //ART
12 
13 //LArSoft
17 
18 // ROOT
19 #include "TLorentzVector.h"
20 #include "TVector3.h"
21 
22 namespace filt {
23 
25  public:
26  explicit LArG4ParticleFilter(fhicl::ParameterSet const& pset);
27  virtual bool filter(art::Event& e);
28 
29  private:
31  const unsigned int index) const;
32  bool PDGCheck(const art::Ptr<simb::MCParticle>& particle, const unsigned int index) const;
33  bool IsPrimaryCheck(const art::Ptr<simb::MCParticle>& particle, const unsigned int index) const;
34  bool MinMomentumCheck(const art::Ptr<simb::MCParticle>& particle,
35  const unsigned int index) const;
36  bool MaxMomentumCheck(const art::Ptr<simb::MCParticle>& particle,
37  const unsigned int index) const;
38  bool StartInTPCCheck(const art::Ptr<simb::MCParticle>& particle,
39  const unsigned int index) const;
40  bool StopInTPCCheck(const art::Ptr<simb::MCParticle>& particle, const unsigned int index) const;
42  const unsigned int index) const;
43 
44  double CalculateLength(const std::vector<TVector3>& position_segment) const;
45 
46  void VerifyDataMembers() const;
47  void VerifyElementRequest(const unsigned int index) const;
48 
50 
51  const std::string fLArG4ModuleLabel;
52  const std::vector<int> fInterestingPDGs;
53  const std::vector<int> fIsPrimary;
54  const std::vector<double> fParticleMinMomentum;
55  const std::vector<double> fParticleMaxMomentum;
56  const std::vector<int> fStartInTPC;
57  const std::vector<int> fStopInTPC;
58  const std::vector<double> fParticleMinTPCLength;
60 
61  std::vector<bool> fFoundInterestingParticles;
62  };
63 
64  LArG4ParticleFilter::LArG4ParticleFilter::LArG4ParticleFilter(fhicl::ParameterSet const& pset)
65  : EDFilter{pset}
66  , fLArG4ModuleLabel(pset.get<std::string>("LArG4ModuleLabel"))
67  , fInterestingPDGs(pset.get<std::vector<int>>("InterestingPDGs"))
68  , fIsPrimary(pset.get<std::vector<int>>("IsPrimary"))
69  , fParticleMinMomentum(pset.get<std::vector<double>>("ParticleMinMomentum"))
70  , fParticleMaxMomentum(pset.get<std::vector<double>>("ParticleMaxMomentum"))
71  , fStartInTPC(pset.get<std::vector<int>>("StartInTPC"))
72  , fStopInTPC(pset.get<std::vector<int>>("StopInTPC"))
73  , fParticleMinTPCLength(pset.get<std::vector<double>>("ParticleMinTPCLength"))
74  , fRequireAllInterestingParticles(pset.get<bool>("RequireAllInterestingParticles"))
76  {
78  }
79 
81  {
82  //Reset the found vector
83  for (unsigned int i = 0; i < fFoundInterestingParticles.size(); i++) {
84  fFoundInterestingParticles[i] = false;
85  }
86 
87  //Get the vector of particles
89  e.getByLabel(fLArG4ModuleLabel, particles);
90  //Loop through the particles
91  for (unsigned int part_i = 0; part_i < particles->size(); part_i++) {
92  //Get the particle
93  const art::Ptr<simb::MCParticle> particle(particles, part_i);
94  //Loop over the list of particles we want and compare it with the particle we are looking it
95  for (unsigned int interest_i = 0; interest_i < fInterestingPDGs.size(); interest_i++) {
96  if (IsInterestingParticle(particle, interest_i)) {
98  return true; //If we only require at least one of the particles be found then we have already done our job
99  fFoundInterestingParticles[interest_i] = true;
100  bool foundThemAll = true;
101  for (unsigned int found_i = 0; found_i < fFoundInterestingParticles.size(); found_i++) {
102  if (fFoundInterestingParticles[found_i] == false) {
103  foundThemAll = false;
104  break;
105  }
106  }
107  if (foundThemAll) { return true; }
108  }
109  }
110  }
111  //Assume that the event is not worth saving
112  return false;
113  }
114 
116  const unsigned int index) const
117  {
118  VerifyElementRequest(index);
119  //Run the checks
120  if (!PDGCheck(particle, index)) return false;
121  if (!IsPrimaryCheck(particle, index)) return false;
122  if (!MinMomentumCheck(particle, index)) return false;
123  if (!MaxMomentumCheck(particle, index)) return false;
124  if (!StartInTPCCheck(particle, index)) return false;
125  if (!StopInTPCCheck(particle, index)) return false;
126  if (!TPCTrajLengthCheck(particle, index)) return false;
127 
128  return true;
129  }
130 
132  const unsigned int index) const
133  {
134  int pdg = fInterestingPDGs[index];
135  if (pdg == 0) return true; //User does not care what the PDG is
136  if (particle->PdgCode() != pdg) return false;
137  return true;
138  }
139 
141  const unsigned int index) const
142  {
143  const int isPrimaryCondition(fIsPrimary[index]);
144  if (isPrimaryCondition == -1) return true;
145  bool particlePrimaryStatus(particle->Mother() > 0 ? false : true);
146  if (particlePrimaryStatus != isPrimaryCondition) return false;
147  return true;
148  }
149 
151  const unsigned int index) const
152  {
153  if (fParticleMinMomentum[index] > 0 &&
154  particle->Momentum(0).Vect().Mag() < fParticleMinMomentum[index])
155  return false;
156  return true;
157  }
158 
160  const unsigned int index) const
161  {
162  if (fParticleMaxMomentum[index] > 0 &&
163  particle->Momentum(0).Vect().Mag() > fParticleMaxMomentum[index])
164  return false;
165  return true;
166  }
167 
169  const unsigned int index) const
170  {
171  //Firstly check if we even care if the particle starts in the TPC or not
172  int demand = fStartInTPC[index];
173  if (demand == 0)
174  return true; //We don't care if the particle starts in the TPC or not so pass the check
175  //Get TPC corresponding to starting position of particle
176  geo::TPCID tpcid = fGeom->FindTPCAtPosition(geo::vect::toPoint(particle->Position(0).Vect()));
177  //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
178  if (tpcid.isValid) {
179  //The particle DID start in a TPC. Now, did we WANT this to happen
180  return demand == 1;
181  }
182  else {
183  //The particle did NOT start in a TPC. Did we WANT this to happen?
184  return demand == 2;
185  }
186 
187  //Assume true by default
188  return true;
189  }
190 
192  const unsigned int index) const
193  {
194  //Firstly check if we even care if the particle stops in the TPC or not
195  int demand = fStopInTPC[index];
196  if (demand == 0)
197  return true; //We don't care if the particle stops in the TPC or not so pass the check
198  //Get final position of particle
199  auto const final_position = particle->Position(particle->NumberTrajectoryPoints() - 1).Vect();
200 
201  geo::TPCID tpcid = fGeom->FindTPCAtPosition(geo::vect::toPoint(final_position));
202  //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
203  if (tpcid.isValid) {
204  //The particle DID stop in a TPC. Now, did we WANT this to happen
205  return demand == 1;
206  }
207  else {
208  //The particle did NOT stop in a TPC. Did we WANT this to happen?
209  return demand == 2;
210  }
211 
212  //Assume true by default
213  return true;
214  }
215 
217  const unsigned int index) const
218  {
219  double min_traj_length = fParticleMinTPCLength[index];
220 
221  //Firstly, if we don't care about the TPC trajectory length then pass the check
222  if (min_traj_length < 0) return true;
223 
224  //Now the hard bit
225  //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
226  //It is easiest to store the positions in a vector of vectors
227  bool OK = false;
228  std::vector<std::vector<TVector3>> position_segments;
229  //We are also going to need an empty vector to store in the above vector
230  std::vector<TVector3> position_segment;
231  //Loop through the trajectory points
232  for (unsigned int i = 0; i < particle->NumberTrajectoryPoints(); i++) {
233  //Extract the current position of the particle
234  geo::TPCID curr_tpcid =
235  fGeom->FindTPCAtPosition(geo::vect::toPoint(particle->Position(i).Vect()));
236  //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
237  //If it is currently in the TPC
238  if (curr_tpcid.isValid) position_segment.push_back(particle->Position(i).Vect());
239  //It has just exited the TPC
240  else if (position_segment.size() > 0) {
241  //Store the segment
242  position_segments.push_back(position_segment);
243  //Now reset the segment
244  position_segment.clear();
245  }
246  //There is nothing to do because the particle has remained outside of the TPC
247  }
248  //We need to check once more if the position_segment vector has been filled
249  if (position_segment.size() > 0) {
250  position_segments.push_back(position_segment);
251  position_segment.clear();
252  }
253  //Now lets check the length of each segment
254  //Firstly, if we didn't store a segment then the particle fails the check
255  if (position_segments.size() == 0) return false;
256  //Now loop through the segments and check if they are above threshold
257  for (unsigned int i = 0; i < position_segments.size(); i++) {
258  double segment_length = CalculateLength(position_segments[i]);
259  if (segment_length > min_traj_length) {
260  //We found a track segment in the TPC which passes the length threshold so don't flag as bad
261  OK = true;
262  break;
263  }
264  }
265  if (!OK) return false;
266 
267  return true;
268  }
269 
270  double LArG4ParticleFilter::CalculateLength(const std::vector<TVector3>& position_segment) const
271  {
272  double length = 0;
273  //Check how many points we have in the segment. If it is one or less, there is nothing to calculate so return 0
274  if (position_segment.size() <= 1) return length;
275 
276  //Now we need to compare every adjacent pair of positions to work out the length of this segment
277  for (unsigned int i = 1; i < position_segment.size(); i++) {
278  length += (position_segment[i] - position_segment[i - 1]).Mag();
279  }
280 
281  return length;
282  }
283 
285  {
286  if (fInterestingPDGs.size() != fIsPrimary.size())
288  << "Config error: InterestingPDGs and IsPrimary FCL vectors are different sizes\n";
289  if (fInterestingPDGs.size() != fParticleMinMomentum.size())
291  << "Config error: InterestingPDGs and ParticleMinMomentum FCL vectors are different "
292  "sizes\n";
293  if (fInterestingPDGs.size() != fParticleMaxMomentum.size())
295  << "Config error: InterestingPDGs and ParticleMaxMomentum FCL vectors are different "
296  "sizes\n";
297  if (fInterestingPDGs.size() != fStartInTPC.size())
299  << "Config error: InterestingPDGs and StartInTPC FCL vectors are different sizes\n";
300  if (fInterestingPDGs.size() != fStopInTPC.size())
302  << "Config error: InterestingPDGs and StopInTPC FCL vectors are different sizes\n";
303  if (fInterestingPDGs.size() != fParticleMinTPCLength.size())
305  << "Config error: InterestingPDGs and ParticleMinTPCLength FCL vectors are different "
306  "sizes\n";
307 
308  for (unsigned int iIsPrimary = 0; iIsPrimary < fIsPrimary.size(); iIsPrimary++)
309  if (!(-1 == fIsPrimary[iIsPrimary] || 0 == fIsPrimary[iIsPrimary] ||
310  1 == fIsPrimary[iIsPrimary]))
312  << "Config error: Element " << iIsPrimary << " of the IsPrimary vector equals "
313  << fIsPrimary[iIsPrimary]
314  << ". Valid options are -1 (not set), 0 (not primary) or 1 (is primary)" << std::endl;
315 
316  for (unsigned int iStartInTPC = 0; iStartInTPC < fStartInTPC.size(); iStartInTPC++)
317  if (!(0 == fStartInTPC[iStartInTPC] || 1 == fStartInTPC[iStartInTPC] ||
318  2 == fStartInTPC[iStartInTPC]))
320  << "Config error: Element " << iStartInTPC << " of the StartInTPC vector equals "
321  << fStartInTPC[iStartInTPC]
322  << ". Valid options are 0 (not set), 1 (start in TPC) or 2 (do not start in TPC)"
323  << std::endl;
324 
325  for (unsigned int iStopInTPC = 0; iStopInTPC < fStopInTPC.size(); iStopInTPC++)
326  if (!(0 == fStopInTPC[iStopInTPC] || 1 == fStopInTPC[iStopInTPC] ||
327  2 == fStopInTPC[iStopInTPC]))
329  << "Config error: Element " << iStopInTPC << " of the StopInTPC vector equals "
330  << fStopInTPC[iStopInTPC]
331  << ". Valid options are 0 (not set), 1 (stop in TPC) or 2 (do not stop in TPC)"
332  << std::endl;
333 
334  return;
335  }
336 
337  void LArG4ParticleFilter::VerifyElementRequest(const unsigned int index) const
338  {
339  if (index > fInterestingPDGs.size())
341  << "Bounds error: Requested element " << index
342  << " from InterestingPDGs vector which is size " << fInterestingPDGs.size() << "\n";
343  if (index > fIsPrimary.size())
345  << "Bounds error: Requested element " << index << " from IsPrimary vector which is size "
346  << fIsPrimary.size() << "\n";
347  if (index > fParticleMinMomentum.size())
349  << "Bounds error: Requested element " << index
350  << " from ParticleMinMomentum vector which is size " << fParticleMinMomentum.size() << "\n";
351  if (index > fParticleMaxMomentum.size())
353  << "Bounds error: Requested element " << index
354  << " from ParticleMaxMomentum vector which is size " << fParticleMaxMomentum.size() << "\n";
355  if (index > fStartInTPC.size())
357  << "Bounds error: Requested element " << index << " from StartInTPC vector which is size "
358  << fStartInTPC.size() << "\n";
359  if (index > fStopInTPC.size())
361  << "Bounds error: Requested element " << index << " from StopInTPC vector which is size "
362  << fStopInTPC.size() << "\n";
363  if (index > fParticleMinTPCLength.size())
365  << "Bounds error: Requested element " << index
366  << " from ParticleMinTPCLength vector which is size " << fParticleMinTPCLength.size()
367  << "\n";
368  return;
369  }
371 }
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:219
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:220
int PdgCode() const
Definition: MCParticle.h:213
std::vector< bool > fFoundInterestingParticles
bool IsInterestingParticle(const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const
int Mother() const
Definition: MCParticle.h:214
art::ServiceHandle< geo::Geometry const > fGeom
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:210
double CalculateLength(const std::vector< TVector3 > &position_segment) const
::geo::Point_t toPoint(Point const &p)
Convert the specified point into a geo::Point_t.
const std::vector< int > fInterestingPDGs
bool IsPrimaryCheck(const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const
const std::vector< double > fParticleMaxMomentum
bool MaxMomentumCheck(const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const
Particle class.
void VerifyElementRequest(const unsigned int index) const
bool StartInTPCCheck(const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const
TPCID FindTPCAtPosition(Point_t const &point) const
Returns the ID of the TPC at specified location.
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:65
bool MinMomentumCheck(const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const
bool TPCTrajLengthCheck(const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
const std::vector< int > fStopInTPC
Definition of data types for geometry description.
cet::coded_exception< errors::ErrorCodes, ExceptionDetail::translate > Exception
Definition: Exception.h:66
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
const std::vector< double > fParticleMinMomentum
const std::vector< int > fIsPrimary
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:221
const std::vector< double > fParticleMinTPCLength
EDFilter(fhicl::ParameterSet const &pset)
Definition: EDFilter.cc:6
virtual bool filter(art::Event &e)
const std::vector< int > fStartInTPC
Float_t e
Definition: plot.C:35
Definition: fwd.h:26
art framework interface to geometry description
bool PDGCheck(const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const
bool StopInTPCCheck(const art::Ptr< simb::MCParticle > &particle, const unsigned int index) const
LArG4ParticleFilter(fhicl::ParameterSet const &pset)