All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Member Functions | Private Attributes | List of all members
filt::LArG4CRTFilter Class Reference
Inheritance diagram for filt::LArG4CRTFilter:

Public Member Functions

 LArG4CRTFilter (fhicl::ParameterSet const &pset)
 
virtual bool filter (art::Event &e) override
 
void reconfigure (fhicl::ParameterSet const &pset)
 
virtual void beginJob () override
 

Private Member Functions

bool IsInterestingParticle (const art::Ptr< simb::MCParticle > particle)
 
void LoadCRTAuxDetIDs ()
 
bool UsesCRTAuxDets (const art::Ptr< simb::MCParticle > particle, const std::vector< unsigned int > &crt_auxdet_vector)
 
bool EntersTPC (const art::Ptr< simb::MCParticle > particle)
 
std::pair< double, double > XLimitsTPC (const art::Ptr< simb::MCParticle > particle)
 

Private Attributes

std::vector< unsigned int > fTopHighCRTAuxDetIDs
 
std::vector< unsigned int > fTopLowCRTAuxDetIDs
 
std::vector< unsigned int > fBottomCRTAuxDetIDs
 
std::vector< unsigned int > fFrontCRTAuxDetIDs
 
std::vector< unsigned int > fBackCRTAuxDetIDs
 
std::vector< unsigned int > fLeftCRTAuxDetIDs
 
std::vector< unsigned int > fRightCRTAuxDetIDs
 
bool fUseTopHighCRTs
 
bool fUseTopLowCRTs
 
bool fUseBottomCRTs
 
bool fUseFrontCRTs
 
bool fUseBackCRTs
 
bool fUseLeftCRTs
 
bool fUseRightCRTs
 
std::vector< int > fPDGs
 
std::vector< double > fMinMomentums
 
std::vector< double > fMaxMomentums
 
std::string fLArG4ModuleName
 
bool fUseReadoutWindow
 
bool fUseTightReadoutWindow
 
bool fUseTPC
 
geo::GeometryCore const * fGeometryService
 
double readoutWindow
 
double driftTime
 

Detailed Description

Definition at line 21 of file LArG4CRTFilter_module.cc.

Constructor & Destructor Documentation

filt::LArG4CRTFilter::LArG4CRTFilter ( fhicl::ParameterSet const &  pset)
explicit

Definition at line 65 of file LArG4CRTFilter_module.cc.

66  : EDFilter(pset)
67  {
68  this->reconfigure(pset);
69 
70  fGeometryService = lar::providerFrom<geo::Geometry>();
71  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
72  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob(clockData);
73  double rw_ticks = detProp.ReadOutWindowSize();
74  double inv_samp_freq = clockData.TPCClock().Frequency();
75  readoutWindow = rw_ticks/inv_samp_freq;
76  // the function TPCTick2Time does not do what you think it does!
77  // It converts ticks to absolute time where 0 is the trigger time and not the start of the readout
78  // window, so if there is a front porch, this value is larger than the actually readout window
79  // readoutWindow = clockData.TPCTick2Time(static_cast<double>(detProp.ReadOutWindowSize())); // [us]
80  driftTime = (2.*fGeometryService->DetHalfWidth())/detProp.DriftVelocity(); // [us]
81  std::cout << "readout window is " << readoutWindow << " us and drift time is " << driftTime << " us "<< std::endl;
82  }
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
void reconfigure(fhicl::ParameterSet const &pset)
geo::GeometryCore const * fGeometryService
BEGIN_PROLOG could also be cout
auto const detProp

Member Function Documentation

void filt::LArG4CRTFilter::beginJob ( )
overridevirtual

Definition at line 173 of file LArG4CRTFilter_module.cc.

173  {
175  }
bool filt::LArG4CRTFilter::EntersTPC ( const art::Ptr< simb::MCParticle >  particle)
private

Definition at line 304 of file LArG4CRTFilter_module.cc.

304  {
305  bool enters = false;
306  double xmin = -2.0 * fGeometryService->DetHalfWidth();
307  double xmax = 2.0 * fGeometryService->DetHalfWidth();
308  double ymin = -fGeometryService->DetHalfHeight();
309  double ymax = fGeometryService->DetHalfHeight();
310  double zmin = 0.;
311  double zmax = fGeometryService->DetLength();
312 
313  int nTrajPoints = particle->NumberTrajectoryPoints();
314  for (int traj_i = 0; traj_i < nTrajPoints; traj_i++){
315  TVector3 trajPoint(particle->Vx(traj_i), particle->Vy(traj_i), particle->Vz(traj_i));
316  // Check if point is within reconstructable volume
317  if (trajPoint[0] >= xmin && trajPoint[0] <= xmax && trajPoint[1] >= ymin && trajPoint[1] <= ymax && trajPoint[2] >= zmin && trajPoint[2] <= zmax){
318  enters = true;
319  }
320  }
321  return enters;
322  }
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
process_name pandoraGausCryo1 vertexChargeCryo1 vertexStubCryo1 xmin
geo::Length_t DetHalfHeight(geo::TPCID const &tpcid) const
Returns the half height of the active volume of the specified TPC.
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
geo::GeometryCore const * fGeometryService
bool filt::LArG4CRTFilter::filter ( art::Event &  e)
overridevirtual

Definition at line 103 of file LArG4CRTFilter_module.cc.

103  {
104  art::Handle<std::vector<simb::MCParticle> > particles;
105  e.getByLabel(fLArG4ModuleName,particles);
106 
107  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(e);
108 
109  for (unsigned int part_i = 0; part_i < particles->size(); part_i++){
110  const art::Ptr<simb::MCParticle> particle(particles,part_i);
111  if (!IsInterestingParticle(particle)) continue;
112  //std::cout<<"particlePDG: " << particle->PdgCode() << std::endl;
113  //particle->Position().Vect().Print();
114  //particle->Momentum().Vect().Print();
115  double time = particle->T() * 1e-3; //[us]
116  if (fUseReadoutWindow){
117  if (time <= -driftTime || time >= readoutWindow) continue;
118  // Get the minimum and maximum |x| position in the TPC
119  std::pair<double, double> xLimits = XLimitsTPC(particle);
120  // Calculate the expected time of arrival of those points
121  double minTime = time + (2.0 * fGeometryService->DetHalfWidth() - xLimits.second)/detProp.DriftVelocity();
122  double maxTime = time + (2.0 * fGeometryService->DetHalfWidth() - xLimits.first)/detProp.DriftVelocity();
123  // If both times are below or above the readout window time then skip
124  if((minTime < 0 && maxTime < 0) || (minTime > readoutWindow && maxTime > readoutWindow)) continue;
125  }
127  if (time<0 || time>(readoutWindow-driftTime)) continue;
128  }
129  if (fUseTPC && !EntersTPC(particle)) continue;
130  if (fUseTopHighCRTs){
131  bool OK = UsesCRTAuxDets(particle,fTopHighCRTAuxDetIDs);
132  if (!OK) continue;
133  //std::cout<<"TopHighCRTs: " << OK << std::endl;
134  }
135  if (fUseTopLowCRTs){
136  bool OK = UsesCRTAuxDets(particle,fTopLowCRTAuxDetIDs);
137  if (!OK) continue;
138  //std::cout<<"TopLowCRTs: " << OK << std::endl;
139  }
140  if (fUseBottomCRTs){
141  bool OK = UsesCRTAuxDets(particle,fBottomCRTAuxDetIDs);
142  if (!OK) continue;
143  //std::cout<<"BottomCRTs: " << OK << std::endl;
144  }
145  if (fUseFrontCRTs){
146  bool OK = UsesCRTAuxDets(particle,fFrontCRTAuxDetIDs);
147  if (!OK) continue;
148  //std::cout<<"FrontCRTs: " << OK << std::endl;
149  }
150  if (fUseBackCRTs){
151  bool OK = UsesCRTAuxDets(particle,fBackCRTAuxDetIDs);
152  if (!OK) continue;
153  //std::cout<<"BackCRTs: " << OK << std::endl;
154  }
155  if (fUseLeftCRTs){
156  bool OK = UsesCRTAuxDets(particle,fLeftCRTAuxDetIDs);
157  if (!OK) continue;
158  //std::cout<<"LeftCRTs: " << OK << std::endl;
159  }
160  if (fUseRightCRTs){
161  bool OK = UsesCRTAuxDets(particle,fRightCRTAuxDetIDs);
162  if (!OK) continue;
163  //std::cout<<"RightCRTs: " << OK << std::endl;
164  }
165  //std::cout<<"Particle hit all CRTs"<<std::endl;
166  return true;
167  }
168 
169  return false;
170  }
std::vector< unsigned int > fTopLowCRTAuxDetIDs
bool IsInterestingParticle(const art::Ptr< simb::MCParticle > particle)
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
std::vector< unsigned int > fRightCRTAuxDetIDs
std::vector< unsigned int > fLeftCRTAuxDetIDs
std::vector< unsigned int > fTopHighCRTAuxDetIDs
std::vector< unsigned int > fFrontCRTAuxDetIDs
bool EntersTPC(const art::Ptr< simb::MCParticle > particle)
bool UsesCRTAuxDets(const art::Ptr< simb::MCParticle > particle, const std::vector< unsigned int > &crt_auxdet_vector)
do i e
std::pair< double, double > XLimitsTPC(const art::Ptr< simb::MCParticle > particle)
geo::GeometryCore const * fGeometryService
auto const detProp
std::vector< unsigned int > fBottomCRTAuxDetIDs
std::vector< unsigned int > fBackCRTAuxDetIDs
bool filt::LArG4CRTFilter::IsInterestingParticle ( const art::Ptr< simb::MCParticle >  particle)
private

Definition at line 178 of file LArG4CRTFilter_module.cc.

178  {
179  double mom = particle->Momentum().Vect().Mag();
180  int pdg = particle->PdgCode();
181 
182  for (unsigned int particle_i = 0; particle_i < fPDGs.size(); particle_i++){
183  bool matched = true;
184  //Check minimum momentum
185  if (fMinMomentums[particle_i] > 0 && mom < fMinMomentums[particle_i]) matched = false;
186  //Check maximum momentum
187  if (fMaxMomentums[particle_i] > 0 && mom > fMaxMomentums[particle_i]) matched = false;
188  //Check PDG
189  if (fPDGs[particle_i]!=0 && pdg!=fPDGs[particle_i]) matched = false;
190  if(matched) return true;
191  }
192  return false;
193  }
var pdg
Definition: selectors.fcl:14
std::vector< double > fMinMomentums
std::vector< double > fMaxMomentums
std::vector< int > fPDGs
void filt::LArG4CRTFilter::LoadCRTAuxDetIDs ( )
private

Definition at line 196 of file LArG4CRTFilter_module.cc.

196  {
197  art::ServiceHandle<geo::Geometry> geom;
198 
199  for (unsigned int auxdet_i = 0; auxdet_i < geom->NAuxDets(); auxdet_i++){
200  geo::AuxDetGeo const& crt = geom->AuxDet(auxdet_i);
201  const TGeoVolume* volModule = crt.TotalVolume();
202  std::set<std::string> volNames = { volModule->GetName() };
203  std::vector<std::vector<TGeoNode const*> > paths = geom->FindAllVolumePaths(volNames);
204 
205  std::string path = "";
206  for (size_t inode=0; inode<paths.at(0).size(); inode++) {
207  path += paths.at(0).at(inode)->GetName();
208  if (inode < paths.at(0).size() - 1) {
209  path += "/";
210  }
211  }
212 
213  TGeoManager* manager = geom->ROOTGeoManager();
214  manager->cd(path.c_str());
215 
216  manager->GetCurrentNode();
217  manager->GetMother(1);
218  TGeoNode* nodeTagger = manager->GetMother(2);
219 
220  std::string taggerName = nodeTagger->GetName();
221 
222  //Now a bunch of if statements so that we can fill our CRT AuxDet ID vectors
223  //BLEH
224  if (taggerName.find("TopHigh")!=std::string::npos){
225  fTopHighCRTAuxDetIDs.push_back(auxdet_i);
226  }
227  else if (taggerName.find("TopLow")!=std::string::npos){
228  fTopLowCRTAuxDetIDs.push_back(auxdet_i);
229  }
230  else if (taggerName.find("Bot")!=std::string::npos){
231  fBottomCRTAuxDetIDs.push_back(auxdet_i);
232  }
233  else if (taggerName.find("South")!=std::string::npos){
234  fFrontCRTAuxDetIDs.push_back(auxdet_i);
235  }
236  else if (taggerName.find("North")!=std::string::npos){
237  fBackCRTAuxDetIDs.push_back(auxdet_i);
238  }
239  else if (taggerName.find("West")!=std::string::npos){
240  fLeftCRTAuxDetIDs.push_back(auxdet_i);
241  }
242  else if (taggerName.find("East")!=std::string::npos){
243  fRightCRTAuxDetIDs.push_back(auxdet_i);
244  }
245  else {
246  std::cout << "Tagger with name: " << taggerName
247  << " does not fit the logic. This should not happen!!!!!!" << std::endl;
248  }
249  }
250 
251  std::cout<< "No. top high CRT AuxDets found: " << fTopHighCRTAuxDetIDs.size() << std::endl;
252  std::cout<< "No. top low CRT AuxDets found: " << fTopLowCRTAuxDetIDs.size() << std::endl;
253  std::cout<< "No. bottom CRT AuxDets found: " << fBottomCRTAuxDetIDs.size() << std::endl;
254  std::cout<< "No. front CRT AuxDets found: " << fFrontCRTAuxDetIDs.size() << std::endl;
255  std::cout<< "No. back CRT AuxDets found: " << fBackCRTAuxDetIDs.size() << std::endl;
256  std::cout<< "No. left CRT AuxDets found: " << fLeftCRTAuxDetIDs.size() << std::endl;
257  std::cout<< "No. right CRT AuxDets found: " << fRightCRTAuxDetIDs.size() << std::endl;
258  return;
259  }
std::vector< unsigned int > fTopLowCRTAuxDetIDs
std::vector< unsigned int > fRightCRTAuxDetIDs
std::vector< unsigned int > fLeftCRTAuxDetIDs
std::vector< unsigned int > fTopHighCRTAuxDetIDs
std::vector< unsigned int > fFrontCRTAuxDetIDs
BEGIN_PROLOG triggeremu_data_config_icarus settings PMTADCthresholds sequence::icarus_stage0_multiTPC_TPC physics sequence::icarus_stage0_EastHits_TPC physics sequence::icarus_stage0_WestHits_TPC physics producers purityana0 caloskimCalorimetryCryoE physics caloskimCalorimetryCryoW physics path
const TGeoVolume * TotalVolume() const
Definition: AuxDetGeo.h:106
process_name crt
BEGIN_PROLOG could also be cout
std::vector< unsigned int > fBottomCRTAuxDetIDs
std::vector< unsigned int > fBackCRTAuxDetIDs
void filt::LArG4CRTFilter::reconfigure ( fhicl::ParameterSet const &  pset)

Definition at line 85 of file LArG4CRTFilter_module.cc.

85  {
86  fUseTopHighCRTs = pset.get<bool>("UseTopHighCRTs");
87  fUseTopLowCRTs = pset.get<bool>("UseTopLowCRTs");
88  fUseBottomCRTs = pset.get<bool>("UseBottomCRTs");
89  fUseFrontCRTs = pset.get<bool>("UseFrontCRTs");
90  fUseBackCRTs = pset.get<bool>("UseBackCRTs");
91  fUseLeftCRTs = pset.get<bool>("UseLeftCRTs");
92  fUseRightCRTs = pset.get<bool>("UseRightCRTs");
93  fPDGs = pset.get<std::vector<int> >("PDGs");
94  fMinMomentums = pset.get<std::vector<double> >("MinMomentums");
95  fMaxMomentums = pset.get<std::vector<double> >("MaxMomentums");
96  fLArG4ModuleName = pset.get<std::string>("LArG4ModuleName");
97  fUseReadoutWindow = pset.get<bool>("UseReadoutWindow");
98  fUseTightReadoutWindow = pset.get<bool>("UseTightReadoutWindow");
99  fUseTPC = pset.get<bool>("UseTPC");
100  }
std::vector< double > fMinMomentums
std::vector< double > fMaxMomentums
std::vector< int > fPDGs
bool filt::LArG4CRTFilter::UsesCRTAuxDets ( const art::Ptr< simb::MCParticle >  particle,
const std::vector< unsigned int > &  crt_auxdet_vector 
)
private

Definition at line 262 of file LArG4CRTFilter_module.cc.

262  {
263  //Loop over the aux dets, extract each one and then perform the test
264  art::ServiceHandle<geo::Geometry> geom;
265  //art:: ServiceHandle <geo:: AuxDetGeometry > adGeoService;
266  //const geo:: AuxDetGeometry* adG = &(* adGeoService); // !!
267  //const geo:: AuxDetGeometryCore* adGeoCore = adG ->GetProviderPtr ();
268 
269  //Loop over the trajectory points made by this particle
270  for (unsigned int pt_i = 0; pt_i < particle->NumberTrajectoryPoints(); pt_i++){
271  //Get the position of the particle
272  TLorentzVector position_lvector = particle->Position(pt_i);
273  //We are going to use a function in the geometry service to see if there is a CRT at this particular position. Said function requires an array, lets make one
274  double position[3];
275  position[0] = position_lvector.X();
276  position[1] = position_lvector.Y();
277  position[2] = position_lvector.Z();
278  //The find the auxdet function throws a wobbler (an exception) if it can't find an auxdet. Wrap what we want to do in a try catch statement pair
279  try{
280  unsigned int crt_id = geom->FindAuxDetAtPosition(position);
281  //size_t adID , svID;
282  //adGeoCore ->PositionToAuxDetChannel(position , adID , svID);
283  //So we got an ID. Lets compare it to all of the CRT IDs we are interested in
284  for (unsigned int crt_i = 0; crt_i < crt_auxdet_vector.size(); crt_i++){
285  if (crt_id == crt_auxdet_vector[crt_i]){
286  //if (adID == crt_auxdet_vector[crt_i]){
287  /*
288  std::cout<<"Hit a CRT we wanted: DUMPING" <<std::endl;
289  for (unsigned int i = 0; i < particle->NumberTrajectoryPoints(); i++){
290  particle->Position(i).Print();
291  }
292  */
293  //We found a CRT that we are interested in at this poisition!!!
294  //We can leave now
295  return true;
296  }
297  }
298  }
299  catch(...){}; //no CRT here, move along
300  }
301  return false;
302  }
std::pair< double, double > filt::LArG4CRTFilter::XLimitsTPC ( const art::Ptr< simb::MCParticle >  particle)
private

Definition at line 324 of file LArG4CRTFilter_module.cc.

324  {
325  double xmin = -2.0 * fGeometryService->DetHalfWidth();
326  double xmax = 2.0 * fGeometryService->DetHalfWidth();
327  double ymin = -fGeometryService->DetHalfHeight();
328  double ymax = fGeometryService->DetHalfHeight();
329  double zmin = 0.;
330  double zmax = fGeometryService->DetLength();
331 
332  double minimum = 99999;
333  double maximum = -99999;
334 
335  int nTrajPoints = particle->NumberTrajectoryPoints();
336  for (int traj_i = 0; traj_i < nTrajPoints; traj_i++){
337  TVector3 trajPoint(particle->Vx(traj_i), particle->Vy(traj_i), particle->Vz(traj_i));
338  // Check if point is within reconstructable volume
339  if (trajPoint[0] >= xmin && trajPoint[0] <= xmax && trajPoint[1] >= ymin && trajPoint[1] <= ymax && trajPoint[2] >= zmin && trajPoint[2] <= zmax){
340  if(std::abs(trajPoint[0]) < minimum) minimum = std::abs(trajPoint[0]);
341  if(std::abs(trajPoint[0]) > maximum) maximum = std::abs(trajPoint[0]);
342  }
343  }
344  return std::make_pair(minimum, maximum);
345  }
geo::Length_t DetHalfWidth(geo::TPCID const &tpcid) const
Returns the half width of the active volume of the specified TPC.
process_name pandoraGausCryo1 vertexChargeCryo1 vertexStubCryo1 xmin
T abs(T value)
geo::Length_t DetHalfHeight(geo::TPCID const &tpcid) const
Returns the half height of the active volume of the specified TPC.
geo::Length_t DetLength(geo::TPCID const &tpcid) const
Returns the length of the active volume of the specified TPC.
geo::GeometryCore const * fGeometryService

Member Data Documentation

double filt::LArG4CRTFilter::driftTime
private

Definition at line 55 of file LArG4CRTFilter_module.cc.

std::vector<unsigned int> filt::LArG4CRTFilter::fBackCRTAuxDetIDs
private

Definition at line 34 of file LArG4CRTFilter_module.cc.

std::vector<unsigned int> filt::LArG4CRTFilter::fBottomCRTAuxDetIDs
private

Definition at line 32 of file LArG4CRTFilter_module.cc.

std::vector<unsigned int> filt::LArG4CRTFilter::fFrontCRTAuxDetIDs
private

Definition at line 33 of file LArG4CRTFilter_module.cc.

geo::GeometryCore const* filt::LArG4CRTFilter::fGeometryService
private

Definition at line 53 of file LArG4CRTFilter_module.cc.

std::string filt::LArG4CRTFilter::fLArG4ModuleName
private

Definition at line 48 of file LArG4CRTFilter_module.cc.

std::vector<unsigned int> filt::LArG4CRTFilter::fLeftCRTAuxDetIDs
private

Definition at line 35 of file LArG4CRTFilter_module.cc.

std::vector<double> filt::LArG4CRTFilter::fMaxMomentums
private

Definition at line 47 of file LArG4CRTFilter_module.cc.

std::vector<double> filt::LArG4CRTFilter::fMinMomentums
private

Definition at line 46 of file LArG4CRTFilter_module.cc.

std::vector<int> filt::LArG4CRTFilter::fPDGs
private

Definition at line 45 of file LArG4CRTFilter_module.cc.

std::vector<unsigned int> filt::LArG4CRTFilter::fRightCRTAuxDetIDs
private

Definition at line 36 of file LArG4CRTFilter_module.cc.

std::vector<unsigned int> filt::LArG4CRTFilter::fTopHighCRTAuxDetIDs
private

Definition at line 30 of file LArG4CRTFilter_module.cc.

std::vector<unsigned int> filt::LArG4CRTFilter::fTopLowCRTAuxDetIDs
private

Definition at line 31 of file LArG4CRTFilter_module.cc.

bool filt::LArG4CRTFilter::fUseBackCRTs
private

Definition at line 42 of file LArG4CRTFilter_module.cc.

bool filt::LArG4CRTFilter::fUseBottomCRTs
private

Definition at line 40 of file LArG4CRTFilter_module.cc.

bool filt::LArG4CRTFilter::fUseFrontCRTs
private

Definition at line 41 of file LArG4CRTFilter_module.cc.

bool filt::LArG4CRTFilter::fUseLeftCRTs
private

Definition at line 43 of file LArG4CRTFilter_module.cc.

bool filt::LArG4CRTFilter::fUseReadoutWindow
private

Definition at line 49 of file LArG4CRTFilter_module.cc.

bool filt::LArG4CRTFilter::fUseRightCRTs
private

Definition at line 44 of file LArG4CRTFilter_module.cc.

bool filt::LArG4CRTFilter::fUseTightReadoutWindow
private

Definition at line 50 of file LArG4CRTFilter_module.cc.

bool filt::LArG4CRTFilter::fUseTopHighCRTs
private

Definition at line 38 of file LArG4CRTFilter_module.cc.

bool filt::LArG4CRTFilter::fUseTopLowCRTs
private

Definition at line 39 of file LArG4CRTFilter_module.cc.

bool filt::LArG4CRTFilter::fUseTPC
private

Definition at line 51 of file LArG4CRTFilter_module.cc.

double filt::LArG4CRTFilter::readoutWindow
private

Definition at line 54 of file LArG4CRTFilter_module.cc.


The documentation for this class was generated from the following file: