CosmicVeto_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 // \file CosmicVeto_module.cc
3 // \brief Preselection to filter cosmic events for reconstruction/PID
4 // for everything
5 // \author Gavin Davies, Kirk Bays
6 ///////////////////////////////////////////////////////////////////////////
7 #include "TMath.h"
8 #include "TVector3.h"
9 
10 //Novasoft includes
11 #include "CosRej/CosRejFxs.h"
12 #include "Geometry/Geometry.h"
13 #include "Geometry/LiveGeometry.h"
14 #include "MCCheater/BackTracker.h"
15 #include "NovaDAQConventions/DAQConventions.h"
16 #include "RecoBase/Cluster.h"
17 #include "RecoBase/FilterList.h"
18 #include "RecoBase/Track.h"
19 #include "Preselection/Veto.h"
20 #include "Utilities/AssociationUtil.h"
23 // Framework includes
30 #include "fhiclcpp/ParameterSet.h"
32 
33 namespace presel
34 {
36  {
37  public:
38  explicit CosmicVeto(fhicl::ParameterSet const &pset);
39  virtual ~CosmicVeto();
40  void beginRun(art::Run& run);
41  virtual void produce(art::Event& evt);
42 
43  protected:
44 
47  bool fRunFilter;
48 
56  bool fDoMinPCut;
58 
59  unsigned int fNCellLow;
60  unsigned int fNCellHigh;
61  unsigned int fFwdCut;
62  unsigned int fBakCut;
63  unsigned int fMinPZ;
64  double fRatio;
66 
68  private:
69 
71 
72  };
73 
74 }//End namespace presel
75 
76 namespace presel
77 {
78  //----------------------------------------------------------------------
80  fSlicerLabel (pset.get< std::string >("SlicerLabel")),
81  fTrackType (pset.get< std::string >("TrackType") ),
82  fRunFilter (pset.get< bool >("RunFilter")),
84  fDoNCellLowCut (0),
85  fDoNCellHighCut (0),
86  fDoContainment (0),
87  fDoOutsideCut (0),
89  fUseBetterCut (0),
90  fDoMinPCut (0),
91  fKeepNuMIInTime (0),
92  fNCellLow (0),
93  fNCellHigh (0),
94  fFwdCut (0),
95  fBakCut (0),
96  fMinPZ (0),
97  fRatio (0),
98  fAngVarCutValue (0),
99  fPSetNDOS (pset.get<fhicl::ParameterSet>("ndos")),
100  fPSetND (pset.get<fhicl::ParameterSet>("nd")),
101  fPSetFD (pset.get<fhicl::ParameterSet>("fd")),
102  fPSetTB (pset.get<fhicl::ParameterSet>("tb"))
103 
104  {
105  produces< rb::FilterList<rb::Cluster> >();
106  produces< std::vector<Veto> >();
107  produces< art::Assns<Veto, rb::Cluster> >();
108  }
109 
110  //----------------------------------------------------------------------
112  {
113  }
114 
115  //----------------------------------------------------------------------
117  {
119 
120  fhicl::ParameterSet pset;
121 
122  switch(geom->DetId()){
123  case novadaq::cnv::kNDOS:
124  pset = fPSetNDOS;
125  break;
127  pset = fPSetND;
128  break;
130  pset = fPSetFD;
131  break;
133  pset = fPSetTB;
134  break;
135  default:
136  assert(0 && "Unknown detector");
137  }
138  fDoTrackAngleCut = pset.get< bool >("DoTrackAngleCut");
139  fDoNCellLowCut = pset.get< bool >("DoNCellLowCut");
140  fDoNCellHighCut = pset.get< bool >("DoNCellHighCut");
141  fDoContainment = pset.get< bool >("DoContainment");
142  fDoOutsideCut = pset.get< bool >("DoOutsideCut");
143  fDoSlicePlanesCut= pset.get< bool >("DoSlicePlanesCut");
144 
145  fUseBetterCut = pset.get< bool >("UseBetterCut");
146  fDoMinPCut = pset.get< bool >("DoMinPCut");
147  fKeepNuMIInTime = pset.get< bool >("KeepNuMIInTime");
148  fNCellLow = pset.get< unsigned int >("NCellLow");
149  fNCellHigh = pset.get< unsigned int >("NCellHigh");
150  fFwdCut = pset.get< unsigned int >("FwdCut");
151  fBakCut = pset.get< unsigned int >("BakCut");
152  fMinPZ = pset.get< unsigned int >("MinPZ");
153  fRatio = pset.get< double >("Ratio");
154  fAngVarCutValue = pset.get< double >("AngVarCutValue");
155  std::cout << "CosmicVeto settings: angle cut: "; // add one line of config to log file
156  if(fUseBetterCut) std::cout << "v2, throughgoing cut: v1, nhit cut: ";
157  if(!fUseBetterCut) std::cout << "v1, throughgoing cut: v1, nhit cut: ";
158  if(!fDoNCellLowCut) std::cout << "none" << std::endl;
159  if(fDoNCellLowCut) std::cout << fNCellLow << std::endl;
160  //if(fKeepNuMIInTime)
161  if(fKeepNuMIInTime)
162  std::cout<<"Not applying cuts to events in time with NuMI\n";
163  }
164 
165  //----------------------------------------------------------------------
167  {
168  std::unique_ptr< rb::FilterList<rb::Cluster> > filtcol(new rb::FilterList<rb::Cluster>);
169 
170  std::unique_ptr< std::vector<presel::Veto> > vetocol(new std::vector<presel::Veto>);
171 
172  std::unique_ptr< art::Assns<presel::Veto, rb::Cluster> > vetoinfo(new art::Assns<presel::Veto, rb::Cluster>);
173 
174  //Grab event's slices
176  evt.getByLabel(fSlicerLabel, slices);
177 
178  // Object holding tracks associated with slice
179  art::FindManyP<rb::Track> fmpTrack(slices, evt, fTrackType);
180 
181  const int sliceMax = slices->size();
182 
183  for(int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
184  art::Ptr<rb::Cluster> slicePtr(slices,sliceIdx); // Used to write association
185 
186  const rb::Cluster& slice = (*slices)[sliceIdx];
187 
188  Veto sliceinfo;
189  if(slice.IsNoise()) continue;
190 
191  float fwddist = 100.0, bakdist = 100.0;
192  double ratio = 1, diry = 0, angle = 1;
193  bool passanglenew=true, passangleold=true, passthru=true;
194 
195  const double nCells = slice.NCell(); //Number of cells in slice
196  sliceinfo.SetNCell(nCells);
197  // Selected unless we fail any cuts below
198  bool sel = true;
199 
200  if(fmpTrack.isValid()){
201  std::vector< art::Ptr<rb::Track> > tracks = fmpTrack.at(sliceIdx);
202  if(tracks.size()>0) {
203  const art::Ptr<rb::Track> track = tracks[0];
204 
206 
207  // get all the important track information
208  TVector3 start = track->Start();
209  TVector3 end = track->Stop();
210  TVector3 dirfor = track->Dir().Unit();
211  TVector3 dirbak = -dirfor;
212 
213  fwddist = livegeom->ProjectedLiveDistanceToEdge(end,dirfor);
214  bakdist = livegeom->ProjectedLiveDistanceToEdge(start,dirbak);
215  sliceinfo.SetFwdDist(fwddist);
216  sliceinfo.SetBakDist(bakdist);
217  diry = track->Dir().Y();
218  sliceinfo.SetDirY(diry);
219  angle = fCosRejFxs->getAngle(track->Dir());
220  sliceinfo.SetAngle(angle);
221  ratio = (track->NCell()*1.0) / (1.0*nCells);
222  sliceinfo.SetRatio(ratio);
223  double anglediff = (angle*angle) - (diry*diry);
224  double anglevar = fabs(angle) * (diry+1.0);
225  sliceinfo.SetAngleVarOld(anglediff);
226  sliceinfo.SetAngleVar(anglevar);
227 
228  // would this pass the old, flawed, through-going check?
229  if(fwddist < fFwdCut &&
230  bakdist < fBakCut)
231  passthru = false;
232  sliceinfo.SetPassThruOld(passthru);
233 
234  // how about the fixed cut we actually use?
235  if(!passthru &&
236  ratio > 0.8)
237  sliceinfo.SetPassThru(false);
238 
239  if(fDoContainment){
240  // apply cut if asked for (fixed only)
241  if(!passthru && ratio > 0.8) sel = false;
242  }
243 
244  // would this pass the old angle cut used for first analysis?
245  if(fmpTrack.isValid()) {
246  if(ratio > 0.88) {
247  if(slice.NCell() > 50 && slice.NCell() <= 135 &&
248  anglediff < ((1.1/85)*slice.NCell()-1.65)) passangleold = false;
249  if(slice.NCell() > 135 && slice.NCell() <= 250 &&
250  anglediff < ((.6/115)*slice.NCell()-.6)) passangleold = false;
251  if(slice.NCell() > 250 && slice.NCell() <= 500 &&
252  anglediff < ((.1/250)*slice.NCell()+0.6)) passangleold = false;
253  if(slice.NCell() > 500 && anglediff < 0.8) passangleold = false;
254  }
255  sliceinfo.SetPassAngleFirstAna(passangleold);
256 
257  // how about the latest version of the angle cut?
258  if( !fDoOutsideCut){
259  if(ratio > 0.75 &&
260  slice.MaxY() > 400 &&
261  anglevar < fAngVarCutValue)
262  passanglenew = false;
263  }
264  else{
265  // For nue cosmic veto, we want the angle cut applied
266  // to not only the events that are close to the dete-
267  // -ctor top, but close to any detector edge.
268  // see doc 14654
269 
270  // For Z containment cut, use LiveGeometry otherwise
271  // cosmic rejection efficiency drops for fewer numb-
272  // -ers of diblocks.
273  double zDistToFront = livegeom->DistToFront(slice.MinXYZ());
274  double zDistToBack = livegeom->DistToBack(slice.MaxXYZ());
275 
276  bool fromOutside = (slice.MaxY() > 400 ||
277  slice.MinY() < -700 ||
278  slice.MaxX() > 650 ||
279  slice.MinX() < -650 ||
280  zDistToFront < 100 ||
281  zDistToBack < 150 );
282 
283  if(ratio > 0.75 &&
284  anglevar < fAngVarCutValue &&
285  fromOutside)
286  passanglenew = false;
287  }
288  sliceinfo.SetPassAngle(passanglenew);
289  }
290 
291  // do the track angle cut if asked for:
292 
293  if (fDoTrackAngleCut){
294  if(!fUseBetterCut) {
295  // old, first analysis tuning
296  if(!passangleold) sel = false;
297  }
298  if(fUseBetterCut) {
299  // better tuning, much better cosmic removal, much simpler
300  if(!passanglenew) sel = false;
301  }
302  }
303 
304  } // end loop over tracks (should only be one)
305  } // end if track object is valid
306 
307  // do the NCell cuts if asked for:
308  if(fDoNCellLowCut){
309  if(nCells < fNCellLow) sel = false;
310  }
311  if(fDoNCellHighCut){
312  if(nCells > fNCellHigh) sel = false;
313  }
314  if( fDoSlicePlanesCut ){
315  // see doc 14654
316  if( slice.ExtentPlane(geo::kXorY) > 125 ||
317  slice.ExtentPlane(geo::kXorY) < 5)
318  sel = false;
319  }
320  // do ND (simulated half det, bad for full ND)
321  if(fDoMinPCut){
322  TVector3 minp = slice.MinXYZ();
323  if(minp.Z() < fMinPZ && ratio > fRatio) sel = false;
324  }
325 
326  // record whether or not we should keep this slice
327  // before checking the timing cuts. Timing cuts are
328  // tracked separately through the PassNuMICut flag.
329  sliceinfo.SetKeep(sel);
330 
331  if( util::IsInBeamWindow(evt.run(), slice.MeanTNS()) ){
332 
333  sliceinfo.SetPassNuMICut(true);
334 
335  if(fKeepNuMIInTime){
336  // if this slice is in the numi beam window,
337  // keep it.
338  sel = true;
339  }
340  }
341 
342  if(!sel && fRunFilter)
343  filtcol->Add(slices, sliceIdx);
344 
345  vetocol->push_back(sliceinfo);
346  util::CreateAssn(*this,evt,*vetocol,slicePtr,*vetoinfo);
347 
348  } // end for sliceIdx
349 
350  evt.put(std::move(filtcol));
351  evt.put(std::move(vetocol));
352  evt.put(std::move(vetoinfo));
353 
354  }
355  //............................................................................................
356 }//end namespace presel
357 
358 namespace presel
359 {
361 }
fhicl::ParameterSet fPSetFD
Preselection Object.
Definition: FillPIDs.h:20
A simple list of products that have been marked "filtered out".
Definition: FilterList.h:74
cosrej::CosRejFxs * fCosRejFxs
back track the reconstruction to the simulation
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.
Double_t angle
Definition: plot.C:86
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void SetBakDist(float)
Definition: Veto.cxx:140
void SetDirY(double)
Definition: Veto.cxx:150
X or Y views.
Definition: PlaneGeo.h:30
fhicl::ParameterSet fPSetTB
static bool IsInBeamWindow(const int run, const double time)
TH1 * ratio(TH1 *h1, TH1 *h2)
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
void SetAngleVarOld(double)
Definition: Veto.cxx:165
double DistToFront(TVector3 vertex)
double DistToBack(TVector3 vertex)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
void SetKeep(bool)
Definition: Veto.cxx:195
double MaxX() const
Definition: Cluster.h:217
void SetRatio(double)
Definition: Veto.cxx:145
Definition: Run.h:31
fhicl::ParameterSet fPSetND
CosmicVeto(fhicl::ParameterSet const &pset)
TVector3 MaxXYZ() const
Definition: Cluster.cxx:492
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
A class detailing the cuts made on a particular slice.
Definition: Veto.h:12
double ProjectedLiveDistanceToEdge(TVector3 vertex, TVector3 dir)
Far Detector at Ash River, MN.
void SetAngle(double)
Definition: Veto.cxx:155
Prototype Near Detector on the surface at FNAL.
T get(std::string const &key) const
Definition: ParameterSet.h:231
float getAngle(TVector3 trackdir)
Definition: CosRejFxs.cxx:33
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
int evt
Near Detector in the NuMI cavern.
void SetPassThru(bool)
Definition: Veto.cxx:170
void SetNCell(unsigned int)
Definition: Veto.cxx:130
double MinY() const
Definition: Cluster.h:205
void SetPassAngle(bool)
Definition: Veto.cxx:180
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
TVector3 MinXYZ() const
Definition: Cluster.cxx:446
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
void geom(int which=0)
Definition: geom.C:163
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
virtual void produce(art::Event &evt)
assert(nhit_max >=nhit_nbins)
void SetAngleVar(double)
Definition: Veto.cxx:160
void SetPassNuMICut(bool)
Definition: Veto.cxx:190
void SetPassAngleFirstAna(bool)
Definition: Veto.cxx:185
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
void SetPassThruOld(bool)
Definition: Veto.cxx:175
fhicl::ParameterSet fPSetNDOS
RunNumber_t run() const
Definition: Event.h:77
double MinX() const
Definition: Cluster.h:204
Encapsulate the geometry of one entire detector (near, far, ndos)
void SetFwdDist(float)
Definition: Veto.cxx:135
double MaxY() const
Definition: Cluster.h:218
enum BeamMode string
void beginRun(art::Run &run)