NuePresel_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 /// \file NuePresel.cxx
3 /// \brief Preselection to filter events for nue analysis, pre-PID
4 /// \author Ruth Toner
5 ///////////////////////////////////////////////////////////////////////////
6 //Novasoft includes:
8 #include "Geometry/Geometry.h"
9 #include "NovaDAQConventions/DAQConventions.h"
10 #include "Preselection/PreselObj.h"
11 #include "RecoBase/CellHit.h"
12 #include "RecoBase/RecoHit.h"
13 #include "RecoBase/Cluster.h"
14 #include "RecoBase/Prong.h"
15 #include "RecoBase/Track.h"
16 #include "RecoBase/Shower.h"
17 #include "RecoBase/FilterList.h"
18 #include "RecoBase/PID.h"
21 #include "Utilities/AssociationUtil.h"
22 
23 // Framework includes
30 #include "fhiclcpp/ParameterSet.h"
31 
32 
33 namespace presel
34 {
35  class NuePresel: public art::EDProducer
36  {
37  public:
38 
39  explicit NuePresel(const fhicl::ParameterSet& pset);
40  virtual ~NuePresel();
41  void beginRun(art::Run& run);
42  virtual void produce(art::Event& evt);
43 
44  protected:
45 
56  double fNCellLow;
57  double fNCellHigh;
58  double fLengthLow;
59  double fLengthHigh;
60  double fPercMipLow;
61  double fPercMipHigh;
62  double fPlanesMax;
63  double fCPlanesLow;
64  double fCPlanesHigh;
65  double fRecoELow;
66  double fRecoEHigh;
67  bool fInFidCut;
68 
70  };
71 
72 } // end namespace presel
73 
74 namespace presel
75 {
76  //----------------------------------------------------------------------
78  fSlicerLabel (pset.get<std::string>("SlicerLabel")),
79  fRecoType (pset.get<std::string>("RecoType" )),
80  fAssn3D (pset.get<std::string>("Assn3D" )),
81  fAssn2D (pset.get<std::string>("Assn2D" )),
82  fDoLengthCut (0),
83  fDoEnergyCut (0),
84  fDoPercMipCut (0),
85  fDoNCellCut (0),
86  fDoCPlanesCut (0),
87  fDoPlanesCut (0),
88  fNCellLow (0),
89  fNCellHigh (0),
90  fLengthLow (0),
91  fLengthHigh (0),
92  fPercMipLow (0),
93  fPercMipHigh (0),
94  fPlanesMax (0),
95  fCPlanesLow (0),
96  fCPlanesHigh (0),
97  fRecoELow (0),
98  fRecoEHigh (0),
99  fInFidCut (0),
100  fPSetNDOS (pset.get<fhicl::ParameterSet>("ndos")),
101  fPSetND (pset.get<fhicl::ParameterSet>("nd")),
102  fPSetFD (pset.get<fhicl::ParameterSet>("fd"))
103  {
104  produces< rb::FilterList<rb::Cluster> >();
105  produces< std::vector<presel::PreselObj> >();
106  produces< std::vector<rb::PID> >();
107  produces< art::Assns<rb::PID, rb::Cluster> >();
108  produces< art::Assns<presel::PreselObj, rb::Cluster> >();
109  }
110 
111  //----------------------------------------------------------------------
113  {
114  }
115 
117  {
119 
120  fhicl::ParameterSet pset;
121 
122  switch(geom->DetId()){
123  case novadaq::cnv::kNDOS:
124  break;
126  pset = fPSetND;
127  break;
129  pset = fPSetFD;
130  break;
131  default:
132  assert(0 && "Unknown detector");
133  }
134 
135  fDoLengthCut = pset.get< bool >("DoLengthCut" );
136  fDoEnergyCut = pset.get< bool >("DoEnergyCut" );
137  fDoPercMipCut = pset.get< bool >("DoPercMipCut" );
138  fDoNCellCut = pset.get< bool >("DoNCellCut" );
139  fDoCPlanesCut = pset.get< bool >("DoCPlanesCut" );
140  fDoPlanesCut = pset.get< bool >("DoPlanesCut" );
141  fNCellLow = pset.get< signed int >("NCellLow" );
142  fNCellHigh = pset.get< signed int >("NCellHigh" );
143  fLengthLow = pset.get< signed int >("LengthLow" );
144  fLengthHigh = pset.get< signed int >("LengthHigh" );
145  fPercMipLow = pset.get< double >("PercMipLow" );
146  fPercMipHigh = pset.get< double >("PercMipHigh" );
147  fPlanesMax = pset.get< signed int >("PlanesMax" );
148  fCPlanesLow = pset.get< signed int >("CPlanesLow" );
149  fCPlanesHigh = pset.get< signed int >("CPlanesHigh" );
150  fRecoELow = pset.get< signed int >("RecoELow" );
151  fRecoEHigh = pset.get< signed int >("RecoEHigh" );
152  fInFidCut = pset.get< bool >("InFidCut" );
153  }
154 
155  //----------------------------------------------------------------------
157  {
158  std::unique_ptr< rb::FilterList<rb::Cluster> > filtcol(new rb::FilterList<rb::Cluster>);
159  std::unique_ptr< std::vector<presel::PreselObj> > nueprecol(new std::vector<presel::PreselObj>);
160  std::unique_ptr< std::vector<rb::PID> > pidcol(new std::vector<rb::PID>);
161  std::unique_ptr< art::Assns<rb::PID, rb::Cluster> > assns(new art::Assns<rb::PID, rb::Cluster>);
162  std::unique_ptr< art::Assns<presel::PreselObj, rb::Cluster> > nuepreassns(new art::Assns<presel::PreselObj, rb::Cluster>);
163 
164  //Set up Calibrator and Geometry services
167 
168  //Grab event's slices
170  evt.getByLabel(fSlicerLabel, slices);
171 
172  // Object holding showers associated with slice
173  art::FindManyP<rb::Shower> fmpShower3D(slices, evt, art::InputTag(fRecoType,fAssn3D));
174  art::FindManyP<rb::Shower> fmpShower2D(slices, evt, art::InputTag(fRecoType,fAssn2D));
175 
176  // Object holding prongs associated with slice
177  art::FindManyP<rb::Prong> fmpProng3D(slices, evt, art::InputTag(fRecoType,fAssn3D));
178  art::FindManyP<rb::Prong> fmpProng2D(slices, evt, art::InputTag(fRecoType,fAssn2D));
179 
180  // Object holding tracks associated with slice
181  art::FindManyP<rb::Track> fmpTrack3D(slices, evt, art::InputTag(fRecoType,fAssn3D));
182  art::FindManyP<rb::Track> fmpTrack2D(slices, evt, art::InputTag(fRecoType,fAssn2D));
183 
184  const int sliceMax = slices->size();
185 
186  for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
187  const rb::Cluster& slice = (*slices)[sliceIdx];
188 
189  if(slice.IsNoise()) continue;
190 
191  presel::PreselObj nuepreobj; // nue preselection object initialisation
192 
193  double longestObject = 0;
194  if(fDoLengthCut){
195  //see if input is tracks
196  if(fmpTrack3D.isValid()){
197  std::vector< art::Ptr<rb::Track> > tracks = fmpTrack3D.at(sliceIdx);
198  //if there is a defined 2D label, use that as well
199  if(fmpTrack2D.isValid() && fAssn2D != ""){
200  std::vector< art::Ptr<rb::Track> > tracks2D = fmpTrack2D.at(sliceIdx);
201  for (unsigned int ip=0; ip<tracks2D.size(); ++ip) tracks.push_back(tracks2D[ip]);
202  }
203 
204  //Find the longest track in the slice:
205  for(unsigned int n = 0; n < tracks.size(); ++n){
206  const double L = tracks[n]->TotalLength();
207  if(L > longestObject) longestObject = L;
208  }
209  }
210  //else fall back to prongs
211  else if(fmpProng3D.isValid()){
212  std::vector< art::Ptr<rb::Prong> > prongs = fmpProng3D.at(sliceIdx);
213  //if there is a defined 2D label, use that as well
214  if(fmpProng2D.isValid() && fAssn2D != ""){
215  std::vector< art::Ptr<rb::Prong> > prongs2D = fmpProng2D.at(sliceIdx);
216  for (unsigned int ip=0; ip<prongs2D.size(); ++ip) prongs.push_back(prongs2D[ip]);
217  }
218 
219  //Find the longest prong in the slice:
220  for(unsigned int n = 0; n < prongs.size(); ++n){
221  const double L = prongs[n]->TotalLength();
222  if(L > longestObject) longestObject = L;
223  }
224  }
225  //else fall back to showers
226  else if(fmpShower3D.isValid()){
227  std::vector< art::Ptr<rb::Shower> > showers = fmpShower3D.at(sliceIdx);
228  //if there is a defined 2D label, use that as well
229  if(fmpShower2D.isValid() && fAssn2D != ""){
230  std::vector< art::Ptr<rb::Shower> > showers2D = fmpShower2D.at(sliceIdx);
231  for (unsigned int ip=0; ip<showers2D.size(); ++ip) showers.push_back(showers2D[ip]);
232  }
233 
234  //Find the longest shower in the slice:
235  for(unsigned int n = 0; n < showers.size(); ++n){
236  const double L = showers[n]->TotalLength();
237  if(L > longestObject) longestObject = L;
238  }
239  }
240  }
241 
242  const double nCells = slice.NCell(); //Number of cells in slice
243  const int lengthPlanes = slice.ExtentPlane(); //Number of planes in the slice
244 
245  //Find whether the slice is in the fiducial volume (based on center of mass):
246  const TVector3 mean = slice.MeanXYZ();
247  const bool inFid = geom->isInsideFiducialVolume(mean.X(), mean.Y(), mean.Z());
248 
249  const double wx = (slice.NXCell() > 0) ? slice.W(slice.XCell(0).get()) : 0;
250  const double wy = (slice.NYCell() > 0) ? slice.W(slice.YCell(0).get()) : 0;
251 
252  //Holds contiguous planes info:
253  int planeNum=-1;
254  std::vector<int> planelist;
255 
256  //Hit counters:
257  double nTotHit = 0.0;
258  double nMipHit = 0.0;
259  double fMipFrac = -1;
260 
261  for(unsigned int hitIdx = 0; hitIdx < slice.NCell(); ++hitIdx){
262  // This is slow because we have to keep recalculating the mean
263  // positions. The manual version below is much better.
264  // const rb::RecoHit rhit = slice.RecoHit(hitIdx);
265 
266  const art::Ptr<rb::CellHit>& chit = slice.Cell(hitIdx);
267 
268  //Get calbrated reco hit:
269  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, chit->View() == geo::kX ? wx : wy));
270 
271  //Make sure things have been calibrated:
272  if(!rhit.IsCalibrated())
273  continue;
274 
275 
276  //Don't do this if we're not using contiguous planes cut
277  if (fDoCPlanesCut){
278 
279  //Make a list of the planes which have hits in the event:
280  planeNum = chit->Plane();
281  if (std::find (planelist.begin(),planelist.end(),planeNum) == planelist.end()) {
282  planelist.push_back(planeNum);
283  }
284 
285  }
286 
287  //Count up total hits:
288  nTotHit++;
289 
290  //Count up MIP hits:
291  if ( (rhit.MIP() > 0.5) && (rhit.MIP() < 1.5) ) nMipHit++;
292  }
293 
294  //Calculate fraction of MIP hits:
295  if (nTotHit>0) fMipFrac = nMipHit/nTotHit;
296 
297  // Selected unless we fail any cuts below
298  bool sel = true;
299 
300  //Do the Fiducial Cut:
301  if(!inFid && fInFidCut) {
302  sel = false;
303  }
304 
305  //Do the NCell cut:
306  if ( fDoNCellCut ){
307  if (nCells <= fNCellLow || nCells > fNCellHigh) sel = false;
308  }
309 
310  //Do the Longest reco object cut:
311  if ( fDoLengthCut ){
312  if (longestObject <= fLengthLow || longestObject > fLengthHigh) sel = false;
313  }
314 
315  //Do the MIP Fraction cut:
316  if ( fDoPercMipCut){
317  if ( fMipFrac <= fPercMipLow || fMipFrac > fPercMipHigh ) sel = false;
318  }
319 
320  //Do the Plane Length cut:
321  if ( fDoPlanesCut ){
322  if ( lengthPlanes > fPlanesMax ) sel = false;
323  }
324 
325  //Do the Energy Cut:
326  if (fDoEnergyCut){
327  const double calE = slice.CalorimetricEnergy();
328  if (calE <= fRecoELow || calE > fRecoEHigh) sel = false;
329  }
330 
331  //----Find the maximum number of contiguous planes, if we're doing this
332  //----cut and if it's still necessary:
333  if (fDoCPlanesCut && sel){
334 
335  int nContig=1;
336  int maxContig=-1;
337 
338  //Contiguous plane counter:
339  if (planelist.size()>0) std::sort(planelist.begin(),planelist.end());
340  for (unsigned int iplane=0; iplane<planelist.size(); iplane++){
341 
342  if (iplane!=0){
343  if (planelist[iplane] == (planelist[iplane-1]+1)){
344  nContig++;
345  if (nContig>maxContig) maxContig=nContig;
346  } else {
347  nContig=0;
348  }
349  }
350  }
351 
352  //Do the Contiguous Planes Cut
353  if (maxContig <= fCPlanesLow || maxContig > fCPlanesHigh) sel = false;
354  }
355 
356  if(!sel) filtcol->Add(slices, sliceIdx);
357 
358  // We're a nue preselection
359  pidcol->push_back(rb::PID(12, sel));
360  nuepreobj.SetPassPresel(sel);
361 
362  nueprecol->push_back(nuepreobj);
363 
364  util::CreateAssn(*this, evt, *pidcol, art::Ptr<rb::Cluster>(slices, sliceIdx), *assns);
365  util::CreateAssn(*this, evt, *nueprecol, art::Ptr<rb::Cluster>(slices, sliceIdx), *nuepreassns);
366 
367  } // end for sliceIdx
368 
369  evt.put(std::move(pidcol));
370  evt.put(std::move(nueprecol));
371  evt.put(std::move(assns));
372  evt.put(std::move(nuepreassns));
373  evt.put(std::move(filtcol));
374  }
375 
376 }//end namespace presel
377 
std::string fSlicerLabel
Preselection Object.
Definition: FillPIDs.h:20
A simple list of products that have been marked "filtered out".
Definition: FilterList.h:74
A pid value and corresponding pdg code.
Definition: PID.h:13
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
pdg code and pid value
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
A collection of associated CellHits.
Definition: Cluster.h:47
TString ip
Definition: loadincs.C:5
TVector3 MeanXYZ(rb::AveragingScheme=kDefaultScheme) const
Definition: Cluster.cxx:538
DEFINE_ART_MODULE(TestTMapFile)
Definition: Run.h:31
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
void beginRun(art::Run &run)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
Far Detector at Ash River, MN.
void SetPassPresel(bool pass)
holds pass/fail output of presel
Definition: PreselObj.cxx:24
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
static constexpr double L
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
NuePresel(const fhicl::ParameterSet &pset)
Prototype Near Detector on the surface at FNAL.
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
int evt
Near Detector in the NuMI cavern.
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
Definition: run.py:1
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
virtual void produce(art::Event &evt)
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
float MIP() const
Definition: RecoHit.cxx:58
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
fhicl::ParameterSet fPSetND
T const * get() const
Definition: Ptr.h:321
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
fhicl::ParameterSet fPSetNDOS
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
fhicl::ParameterSet fPSetFD
void geom(int which=0)
Definition: geom.C:163
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
assert(nhit_max >=nhit_nbins)
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
Object collecting Preselection variables.
Encapsulate the geometry of one entire detector (near, far, ndos)
bool isInsideFiducialVolume(const double x_cm, const double y_cm, const double z_cm) const
Is the particle inside the detector Fiducial Volume?