CFilter_module.cc
Go to the documentation of this file.
1 //
2 // CFilter: Module to use M. Worcester's MuonCounter module as a filter
3 //
4 // M. Elnimr Dec. 2014
5 //
6 //
7 //
8 //
9 
17 #include "lardataobj/RawData/raw.h"
27 #include "TMath.h"
28 
31 
33 
34 
35 
36 namespace filt{
37 
38  class CFilter : public art::EDFilter {
39  public:
40  explicit CFilter(fhicl::ParameterSet const& pset);
41  virtual ~CFilter() { }
42  virtual bool filter(art::Event& e);
43  void reconfigure(fhicl::ParameterSet const& pset);
44 
45  private:
46 
47  // the parameters we'll read from the .fcl file
50  // std::string fCounterFile;
52  int fSelectedPDG; // PDG code for primary particle
53  int fTrigger; // as per Matthew which Trigger type
54  // log file
55  // ofstream outfile;
56 
57 
61 
62  };
63 
65  {
66  this->reconfigure(pset);
67 
68  }
69 
71  {
72 
73  fSimulationProducerLabel = p.get< std::string >("SimulationLabel");
74  // fCounterFile = p.get< std::string >("CounterFile");
75  fSelectedPDG = p.get< int >("PDGcode");
76  fCounterType=p.get< int >("CounterType");
77  fTrigger=p.get< int >("Trigger");
78  fnumTracks=1;
79 
80  // fnumTracks = pset.get<int>("numTracks");
81 
82  }
83 
85  {
86  //muon conter geometry
87  int counters_loaded=-1;
88  std::vector<std::vector<double> > countergeometry;
89  //load the muon counter positions from a text file
90 
91  char counterfile[]= "/afs/fnal.gov/files/home/room1/melnimr/lr_dev_7/srcs/dunetpc/dune/Geometry/muoncounters.txt";
92  counters_loaded=geo::MuonCounter35Alg::loadMuonCounterGeometry(counterfile,countergeometry);
93  bool keepFlag=0;
94 
95  // int icount=0; int trkind;
96  // define the handle as an MCParticle vector and fill it with events from the simulation
97  auto particleHandle = event.getHandle< std::vector<simb::MCParticle> >(fSimulationProducerLabel);
98 
99  // define a sorted map in which to put the particles
100  std::map< int, const simb::MCParticle* > particleMap;
101  // loop over all the particles, find the primary muon, and get the initial/final positions
102  for ( auto const& particle : (*particleHandle) )
103  {
104  // For the methods you can call to get particle information,
105  // see $NUTOOLS_INC/SimulationBase/MCParticle.h.
106  int fTrackID = particle.TrackId();
107  // Add the address of the MCParticle to the map, with the track ID as the key.
108  particleMap[fTrackID] = &particle;
109  // PDG code of every particle in the event.
110  // int fPDG = particle.PdgCode();
111  // locate the primary muon
112  // if ( particle.Process() == "primary" && fPDG == fSelectedPDG )
113  // {
114 
115  // A particle has a trajectory, consisting of a set of
116  // 4-positions and 4-mommenta.
117  size_t numberTrajectoryPoints = particle.NumberTrajectoryPoints();
118  // For trajectories, as for vectors and arrays, the
119  // first point is #0, not #1.
120  int last = numberTrajectoryPoints - 1;
121  const TLorentzVector& positionStart = particle.Position(0);
122  const TLorentzVector& positionEnd = particle.Position(last);
123  const TLorentzVector& momentumStart = particle.Momentum(0);
124  const TLorentzVector& momentumEnd = particle.Momentum(last);
125  // move the initial position (for muon counter studies)
126  double x_increment = 0.; // 349.666
127  double z_increment = 0.; // 231.065
128  TVector3 trackStart(positionStart.X()+x_increment,positionStart.Y(),positionStart.Z()+z_increment);
129  //}
130 
131  std::vector< std::vector<double> > hitcounters;
132  if(counters_loaded){
133  //unsigned int counters_hit=0;
134 
135  geo::MuonCounter35Alg::testTrackInAllCounters(fTrackID,trackStart,momentumStart.Vect(),
136  countergeometry,hitcounters);
137  }
138 
139 
140  // condition flags for each layer
141  bool Layer_1_2 = false;
142  bool Layer_3_4_5 = false;
143  bool Layer_E = false;
144  bool Layer_W = false;
145  bool Layer_N_U = false;
146  bool Layer_N_L = false;
147  bool Layer_S_U = false;
148  bool Layer_S_L = false;
149  int Trigger= 0;
150 
151  for(unsigned int ii=0;ii<hitcounters.size();ii++)
152  {
153  // for(unsigned int jj=0;jj<hitcounters[ii].size();jj++)
154  // {
155  std::cout<< "::Filter, counters triggered are: "<<hitcounters[ii][0] <<std::endl;
156  /* if (hitcounters[ii][0]==fCounterType)
157  {
158  keepFlag=1;
159  return keepFlag;
160  }*/
161  // check which layer is hit
162  if( 40 <= hitcounters[ii][0] && hitcounters[ii][0] <= 61){
163  Layer_1_2 = true;
164  }
165  if( hitcounters[ii][0] > 61){
166  Layer_3_4_5 = true;
167  }
168  if (14 <= hitcounters[ii][0] && hitcounters[ii][0] <=19){
169  Layer_N_U = true;
170  }
171  if ( 34 <= hitcounters[ii][0] && hitcounters[ii][0] <= 39){
172  Layer_S_L = true;
173  }
174  if (8 <= hitcounters[ii][0] && hitcounters[ii][0] <= 13){
175  Layer_N_L = true;
176  }
177  if (28 <= hitcounters[ii][0] && hitcounters[ii][0] <= 33){
178  Layer_S_U = true;
179  }
180  if (hitcounters[ii][0] <= 7){
181  Layer_E = true;
182  }
183  if (20 <= hitcounters[ii][0] && hitcounters[ii][0] <= 27){
184  Layer_W = true;
185  }
186  // }
187  }
188  // check for a satisfied trigger condition
189  if (Layer_1_2 && Layer_3_4_5){
190  Trigger = 1;
191  }
192  if (Layer_N_U && Layer_S_L){
193  Trigger = 2;
194  }
195  if (Layer_N_L && Layer_S_U){
196  Trigger = 3;
197  }
198  if (Layer_E && Layer_W){
199  Trigger = 4;
200  }
201  std::cout << "Trigger is ....." << Trigger << std::endl;
202  if(Trigger==fTrigger)
203  {
204  keepFlag=1;
205  return keepFlag;
206  }
207  }
208  return keepFlag;
209  }
210  // A macro required for a JobControl module.
212 
213 } // namespace filt
CFilter(fhicl::ParameterSet const &pset)
virtual ~CFilter()
std::string string
Definition: nybbler.cc:12
Declaration of signal hit object.
Definition of basic raw digits.
void reconfigure(fhicl::ParameterSet const &pset)
std::string fSimulationProducerLabel
Particle class.
static int loadMuonCounterGeometry(char *filename, std::vector< std::vector< double > > &geometry)
art framework interface to geometry description
virtual bool filter(art::Event &e)
const double e
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
Collect all the RawData header files together.
T get(std::string const &key) const
Definition: ParameterSet.h:271
std::string fHitsModuleLabel
Declaration of cluster object.
Provides recob::Track data product.
Encapsulate the geometry of a wire.
Interface to algorithm class for DUNE 35t muon counters.
p
Definition: test.py:223
std::string fRawDigitLabel
Utility object to perform functions of association.
Encapsulate the construction of a single detector plane.
EDFilter(fhicl::ParameterSet const &pset)
Definition: EDFilter.h:21
static int testTrackInAllCounters(int trackID, TVector3 trackpoint, TVector3 trackvector, std::vector< std::vector< double > > &geometry, std::vector< std::vector< double > > &hitcounters)
Declaration of basic channel signal object.
QTextStream & endl(QTextStream &s)
Event finding and building.
std::string fTrackModuleLabel