MakeNueCosRej_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: MakeNueCosRej
3 // Module Type: producer
4 // \file MakeNueCosRej_module.cc
5 // \brief A module to write out the NueCosRej object
6 //
7 // Generated at Mon Jul 28 16:46:22 2014 by Kanika Sachdev using artmod
8 // from cetpkgsupport v1_06_02.
9 ////////////////////////////////////////////////////////////////////////
10 
11 // art includes
12 
18 #include "fhiclcpp/ParameterSet.h"
24 #include "cetlib/search_path.h"
25 #include "cetlib/filesystem.h"
26 
27 #include <memory>
28 #include "TMath.h"
29 #include "TVector3.h"
30 #include "TParticlePDG.h"
31 #include "TDatabasePDG.h"
32 #include "TLorentzVector.h"
33 #include "TMVA/Reader.h"
34 // novasoft includes
35 
36 #include "CosRej/NueCosRej.h"
37 #include "CVN/func/Result.h"
39 #include "ShowerLID/ShowerLID.h"
40 #include "NovaDAQConventions/DAQConventions.h"
41 #include "Geometry/LiveGeometry.h"
42 #include "Geometry/Geometry.h"
44 #include "Utilities/AssociationUtil.h"
45 #include "GeometryObjects/Geo.h"
46 
47 #include "RecoBase/FilterList.h"
48 #include "RecoBase/Prong.h"
49 #include "RecoBase/Cluster.h"
50 #include "RecoBase/Shower.h"
51 #include "RecoBase/Vertex.h"
52 
53 #include "SummaryData/SpillData.h"
54 
55 #include <iostream>
56 
57 
58 namespace cosrej
59 {
60  class MakeNueCosRej : public art::EDProducer {
61  public:
62 
63  explicit MakeNueCosRej(fhicl::ParameterSet const & p);
65 
66  void beginRun(art::Run &evt) override;
67  void reconfigure(const fhicl::ParameterSet& pset);
68  void produce(art::Event & e) override;
69 
70 
71  private:
73 
82 
87 
92 
93  TMVA::Reader* fReaderContain;
94  TMVA::Reader* fReaderContainXY;
95  TMVA::Reader* fReaderLight;
96 
97 
98  TMVA::Reader* fReaderFHCPeriBDT;
99  TMVA::Reader* fReaderRHCPeriBDT;
100  TMVA::Reader* fReaderFHCCoreBDT;
101  TMVA::Reader* fReaderRHCCoreBDT;
102 
103 
104  // There are a lot ov variables we need to send to TMVA::Factory
105  float TMVAvars[13];
106 
107  float PeriTMVAvars[16];
108  float CoreTMVAvars[17];
109 
110 
111  /// Particle ID alorithm (loglikelihoods and dE/dx)
113  /// FCL parameters for particle ID alorithm (loglikelihoods and dE/dx)
115 
116  /// Return transverse momentum fraction of the event. Calculation based on
117  /// reconstructed prongs.
118  double TransMomFraction(const std::vector
119  < art::Ptr<rb::Prong> >& prongs) const;
120 
121  /// Compute transverse momentum fraction of the event. Based on showers
122  /// One calculation assumes 0 mass in computing particle momentum,
123  /// the other assigns mass based on which particle hypothesis has
124  /// maximum summed longitudinal and transverse log-likelihood.
125  void TransMomFractionShower(std::vector< const slid::ShowerLID* > shwlids,
127  TVector3 & nuRecoP,
128  TVector3 & nuRecoPLIDMass);
129 
130  /// Return particle mass based on particle type in
131  /// \ref slid::DedxParticleType
132  double Mass(const int type);
133 
134 
135  //double ShowerlidWidth(std::vector< const slid::ShowerLID* > shwlids, art::FindOneP<rb::Shower> fos);
136  // double ShowerlidcalE(std::vector< const slid::ShowerLID* > shwlids, art::FindOneP<rb::Shower> fos);
137 
138  /// \brief Return variables measuring forward-backward
139  /// asymmetry of the prong or slice.
140  /// \var sparsenessAsymm(Slice) is the difference divided by
141  /// sum of the number of zero hit planes in the first
142  /// and last 8 planes of the prong/slice
143  /// \var hitsPerPlaneAsymm(Slice) is the difference divided
144  /// by sum of the number of hits per plane in the
145  /// first and last 5 planes of the prong/slice
146  void SparsenessAsymmetry(const rb::Cluster* cluster,
147  float & sparsenessAsymm,
148  float & hitsPerPlaneAsymm);
149 
150  void MuonParentByDist(const art::Ptr< rb::Prong>& leadingPng,
151  art::FindManyP<rb::Track> pngAssn,
152  const int& idx,
153  const int& nslc,
154  double& mintimebydist,
155  double& minanglebydist,
156  double& closestapproachbydist,
157  int& musliceidxbydist);
158 
159  void MuonParentByTime(const art::Ptr< rb::Prong>& leadingPng,
161  const int& idx,
162  const int& nslc,
163  double& mintimebytime,
164  double& minanglebydist,
165  double& closestapproachbytime,
166  int& musliceidxbytime);
167 
168 
169 
170 
171 
172  bool IsRHC(const art::Event &evt);
173 
174  };
175 
176 
177 
179  const art::Ptr<rb::Prong> b);
180 
181 }
182 
183 ////////////////////////////////////////////////////////////////////////
184 ////////////////////////////////////////////////////////////////////////
185 
187  :fSliceToken(consumes<std::vector<rb::Cluster>>(p.get<std::string>("SliceLabel"))),
188  fReaderContain(nullptr),
189  fReaderContainXY(nullptr),
190  fReaderLight(nullptr),
191  fReaderFHCPeriBDT(nullptr),
192  fReaderRHCPeriBDT(nullptr),
193  fReaderFHCCoreBDT(nullptr),
194  fReaderRHCCoreBDT(nullptr),
195  fParticleIDAlg(0),
196  fParticleIDAlgPSet(p.get< fhicl::ParameterSet >("ParticleIDAlgPSet"))
197 {
198  produces< std::vector<cosrej::NueCosRej> > ();
199  produces< art::Assns< cosrej::NueCosRej, rb::Cluster> > ();
200  this->reconfigure(p);
201 }
202 
203 ////////////////////////////////////////////////////////////////////////
204 ////////////////////////////////////////////////////////////////////////
205 
207 {
208  if (fReaderContain) delete fReaderContain;
210  if (fReaderLight) delete fReaderLight;
215 }
216 
217 ////////////////////////////////////////////////////////////////////////
218 ////////////////////////////////////////////////////////////////////////
219 
221 {
223 
224  // Make sure we can find all files before we set up MVA
225  if (!cet::file_exists(pidlibpath+fCosRejPIDFileContain)){
226  throw cet::exception("MakeNueCosRej")
227  << "Couldn't find file " << pidlibpath+fCosRejPIDFileContain<<std::endl;
228  }
229  if (!cet::file_exists(pidlibpath+fCosRejPIDFileContainXY)){
230  throw cet::exception("MakeNueCosRej") << "Couldn't find file " <<
231  pidlibpath+fCosRejPIDFileContainXY << std::endl;
232  }
233 
234  if (!cet::file_exists(pidlibpath+fCosRejPIDFileLight)){
235  throw cet::exception("MakeNueCosRej")
236  << "Couldn't find file " << pidlibpath+fCosRejPIDFileLight << std::endl;
237  }
238 
239  if (!cet::file_exists(pidlibpath+fCosRejPIDFileFHCPeri)){
240  throw cet::exception("MakeNueCosRej")
241  << "Couldn't find file " << pidlibpath+fCosRejPIDFileFHCPeri << std::endl;
242  }
243 
244  if (!cet::file_exists(pidlibpath+fCosRejPIDFileRHCPeri)){
245  throw cet::exception("MakeNueCosRej")
246  << "Couldn't find file " << pidlibpath+fCosRejPIDFileRHCPeri << std::endl;
247  }
248 
249 
250  if (!cet::file_exists(pidlibpath+fCosRejPIDFileFHCCore)){
251  throw cet::exception("MakeNueCosRej")
252  << "Couldn't find file " << pidlibpath+fCosRejPIDFileFHCCore << std::endl;
253  }
254 
255  if (!cet::file_exists(pidlibpath+fCosRejPIDFileRHCCore)){
256  throw cet::exception("MakeNueCosRej")
257  << "Couldn't find file " << pidlibpath+fCosRejPIDFileRHCCore << std::endl;
258  }
259 
260 
262  // Not all variables are in every tree in the same order
263  fReaderContain = new TMVA::Reader;
264  fReaderContain->AddVariable ("nhit", &TMVAvars[0]);
265  fReaderContain->AddVariable ("sparseassym",&TMVAvars[2]);
266  fReaderContain->AddVariable ("ptp", &TMVAvars[3]);
267  fReaderContain->AddVariable ("disttop", &TMVAvars[6]);
268  fReaderContain->AddVariable ("distnottop", &TMVAvars[7]);
269  fReaderContain->AddSpectator("woscdumb", &TMVAvars[12]);
270  fReaderContain->AddSpectator("issig", &TMVAvars[9]);
271  fReaderContain->AddSpectator("isbb", &TMVAvars[10]);
272  fReaderContain->AddSpectator("iscos", &TMVAvars[11]);
273  fReaderContain->BookMVA("BDT_Contain",pidlibpath+fCosRejPIDFileContain);
274 
275  fReaderContainXY = new TMVA::Reader;
276  fReaderContainXY->AddVariable ("nhit", &TMVAvars[0]);
277  fReaderContainXY->AddSpectator("CVNe", &TMVAvars[8]);
278  fReaderContainXY->AddSpectator("ptp", &TMVAvars[3]);
279  fReaderContainXY->AddVariable ("pxp", &TMVAvars[4]);
280  fReaderContainXY->AddVariable ("pyp", &TMVAvars[5]);
281  fReaderContainXY->AddVariable ("sparseSlice",&TMVAvars[1]);
282  fReaderContainXY->AddSpectator("sparseProng",&TMVAvars[2]);
283  fReaderContainXY->AddVariable ("disttop", &TMVAvars[6]);
284  fReaderContainXY->AddVariable ("distnottop", &TMVAvars[7]);
285  fReaderContainXY->AddSpectator("issig", &TMVAvars[9]);
286  fReaderContainXY->AddSpectator("isbb", &TMVAvars[10]);
287  fReaderContainXY->BookMVA("BDTContain2018",pidlibpath+fCosRejPIDFileContainXY);
288 
289  fReaderLight = new TMVA::Reader;
290  fReaderLight->AddVariable ("nhit", &TMVAvars[0]);
291  fReaderLight->AddVariable ("CVNe", &TMVAvars[8]);
292  fReaderLight->AddSpectator("ptp", &TMVAvars[3]);
293  fReaderLight->AddVariable ("pxp", &TMVAvars[4]);
294  fReaderLight->AddVariable ("pyp", &TMVAvars[5]);
295  fReaderLight->AddVariable ("sparseSlice",&TMVAvars[1]);
296  fReaderLight->AddSpectator("sparseProng",&TMVAvars[2]);
297  fReaderLight->AddVariable ("disttop", &TMVAvars[6]);
298  fReaderLight->AddVariable ("distnottop", &TMVAvars[7]);
299  fReaderLight->AddSpectator("issig", &TMVAvars[9]);
300  fReaderLight->AddSpectator("isbb", &TMVAvars[10]);
301  fReaderLight->BookMVA("BDTLight2018",pidlibpath+fCosRejPIDFileLight);
302 
303 
304 
305 
306  fReaderFHCPeriBDT = new TMVA::Reader;
307  fReaderFHCPeriBDT->AddVariable ("nhit", &PeriTMVAvars[0]);
308  fReaderFHCPeriBDT->AddVariable ("sparseassym", &PeriTMVAvars[1]);
309  fReaderFHCPeriBDT->AddVariable ("pxp", &PeriTMVAvars[3]);
310  fReaderFHCPeriBDT->AddVariable ("pyp", &PeriTMVAvars[4]);
311  fReaderFHCPeriBDT->AddVariable ("disttop", &PeriTMVAvars[5]);
312  fReaderFHCPeriBDT->AddVariable ("distnottop", &PeriTMVAvars[6]);
313  fReaderFHCPeriBDT->AddVariable ("vtxx", &PeriTMVAvars[7]);
314  fReaderFHCPeriBDT->AddVariable ("vtxy", &PeriTMVAvars[8]);
315  fReaderFHCPeriBDT->AddVariable ("vtxz", &PeriTMVAvars[9]);
316 
317  fReaderFHCPeriBDT->AddSpectator("cvnsse", &PeriTMVAvars[10]);
318  fReaderFHCPeriBDT->AddSpectator("woscdumb", &PeriTMVAvars[11]);
319  fReaderFHCPeriBDT->AddSpectator("isGENIE", &PeriTMVAvars[12]);
320  fReaderFHCPeriBDT->AddSpectator("pdg", &PeriTMVAvars[13]);
321  fReaderFHCPeriBDT->AddSpectator("pdgorig", &PeriTMVAvars[14]);
322  fReaderFHCPeriBDT->AddSpectator("isCC", &PeriTMVAvars[15]);
323  fReaderFHCPeriBDT->BookMVA("FHCPeriBDT",pidlibpath+fCosRejPIDFileFHCPeri);
324 
325 
326  fReaderRHCPeriBDT = new TMVA::Reader;
327  fReaderRHCPeriBDT->AddVariable ("nhit", &PeriTMVAvars[0]);
328  fReaderRHCPeriBDT->AddVariable ("sparseassym", &PeriTMVAvars[1]);
329  fReaderRHCPeriBDT->AddVariable ("ptp", &PeriTMVAvars[2]);
330  fReaderRHCPeriBDT->AddVariable ("pxp", &PeriTMVAvars[3]);
331  fReaderRHCPeriBDT->AddVariable ("pyp", &PeriTMVAvars[4]);
332  fReaderRHCPeriBDT->AddVariable ("disttop", &PeriTMVAvars[5]);
333  fReaderRHCPeriBDT->AddVariable ("distnottop", &PeriTMVAvars[6]);
334  fReaderRHCPeriBDT->AddVariable ("vtxx", &PeriTMVAvars[7]);
335  fReaderRHCPeriBDT->AddVariable ("vtxy", &PeriTMVAvars[8]);
336  fReaderRHCPeriBDT->AddVariable ("vtxz", &PeriTMVAvars[9]);
337 
338  fReaderRHCPeriBDT->AddSpectator("cvnsse", &PeriTMVAvars[10]);
339  fReaderRHCPeriBDT->AddSpectator("woscdumb", &PeriTMVAvars[11]);
340  fReaderRHCPeriBDT->AddSpectator("isGENIE", &PeriTMVAvars[12]);
341  fReaderRHCPeriBDT->AddSpectator("pdg", &PeriTMVAvars[13]);
342  fReaderRHCPeriBDT->AddSpectator("pdgorig", &PeriTMVAvars[14]);
343  fReaderRHCPeriBDT->AddSpectator("isCC", &PeriTMVAvars[15]);
344  fReaderRHCPeriBDT->BookMVA("RHCPeriBDT",pidlibpath+fCosRejPIDFileRHCPeri);
345 
346 
347  fReaderFHCCoreBDT = new TMVA::Reader;
348  fReaderFHCCoreBDT->AddVariable ("nhit", &CoreTMVAvars[0]);
349  fReaderFHCCoreBDT->AddVariable ("sparseassym", &CoreTMVAvars[1]);
350  fReaderFHCCoreBDT->AddVariable ("ptp", &CoreTMVAvars[2]);
351  fReaderFHCCoreBDT->AddVariable ("disttop", &CoreTMVAvars[3]);
352  fReaderFHCCoreBDT->AddVariable ("distback", &CoreTMVAvars[4]);
353  fReaderFHCCoreBDT->AddVariable ("distwest", &CoreTMVAvars[5]);
354  fReaderFHCCoreBDT->AddVariable ("distfront", &CoreTMVAvars[6]);
355  fReaderFHCCoreBDT->AddVariable ("disteast", &CoreTMVAvars[7]);
356  fReaderFHCCoreBDT->AddVariable ("distbottom", &CoreTMVAvars[8]);
357  fReaderFHCCoreBDT->AddVariable ("inelast", &CoreTMVAvars[9]);
358  fReaderFHCCoreBDT->AddVariable ("widthofshow", &CoreTMVAvars[10]);
359 
360  fReaderFHCCoreBDT->AddSpectator("cvnsse", &CoreTMVAvars[11]);
361  fReaderFHCCoreBDT->AddSpectator("woscdumb", &CoreTMVAvars[12]);
362  fReaderFHCCoreBDT->AddSpectator("isGENIE", &CoreTMVAvars[13]);
363  fReaderFHCCoreBDT->AddSpectator("pdg", &CoreTMVAvars[14]);
364  fReaderFHCCoreBDT->AddSpectator("pdgorig", &CoreTMVAvars[15]);
365  fReaderFHCCoreBDT->AddSpectator("isCC", &CoreTMVAvars[16]);
366  fReaderFHCCoreBDT->BookMVA("FHCCoreBDT",pidlibpath+fCosRejPIDFileFHCCore);
367 
368 
369  fReaderRHCCoreBDT = new TMVA::Reader;
370  fReaderRHCCoreBDT->AddVariable ("nhit", &CoreTMVAvars[0]);
371  fReaderRHCCoreBDT->AddVariable ("sparseassym", &CoreTMVAvars[1]);
372  fReaderRHCCoreBDT->AddVariable ("ptp", &CoreTMVAvars[2]);
373  fReaderRHCCoreBDT->AddVariable ("disttop", &CoreTMVAvars[3]);
374  fReaderRHCCoreBDT->AddVariable ("distback", &CoreTMVAvars[4]);
375  fReaderRHCCoreBDT->AddVariable ("distwest", &CoreTMVAvars[5]);
376  fReaderRHCCoreBDT->AddVariable ("distfront", &CoreTMVAvars[6]);
377  fReaderRHCCoreBDT->AddVariable ("disteast", &CoreTMVAvars[7]);
378  fReaderRHCCoreBDT->AddVariable ("distbottom", &CoreTMVAvars[8]);
379  fReaderRHCCoreBDT->AddVariable ("inelast", &CoreTMVAvars[9]);
380  fReaderRHCCoreBDT->AddVariable ("widthofshow", &CoreTMVAvars[10]);
381 
382  fReaderRHCCoreBDT->AddSpectator("cvnsse", &CoreTMVAvars[11]);
383  fReaderRHCCoreBDT->AddSpectator("woscdumb", &CoreTMVAvars[12]);
384  fReaderRHCCoreBDT->AddSpectator("isGENIE", &CoreTMVAvars[13]);
385  fReaderRHCCoreBDT->AddSpectator("pdg", &CoreTMVAvars[14]);
386  fReaderRHCCoreBDT->AddSpectator("pdgorig", &CoreTMVAvars[15]);
387  fReaderRHCCoreBDT->AddSpectator("isCC", &CoreTMVAvars[16]);
388  fReaderRHCCoreBDT->BookMVA("RHCCoreBDT",pidlibpath+fCosRejPIDFileRHCCore);
389 
390 }
391 
392 ////////////////////////////////////////////////////////////////////////
393 ////////////////////////////////////////////////////////////////////////
394 
396 {
397  fCVNLabel = p.get< std::string >("CVNLabel");
398  fVertexLabel = p.get< std::string >("VertexLabel");
399  fProngLabel = p.get< std::string >("ProngLabel");
400  fProngInstance = p.get< std::string >("ProngInstance");
401  fShowerLIDLabel = p.get< std::string >("ShowerLIDLabel");
402  fCosmicTrackLabel = p.get< std::string >("CosmicTrackLabel");
403 
404  fGeneratorLabel = p.get< std::string >("GeneratorLabel");
405  fNuMILabel = p.get< std::string >("NuMILabel");
406 
407 
408  fPIDLibPath = p.get< std::string >("PIDLibPath");
409  fCosRejPIDFileContain = p.get< std::string >("CosRejPIDFileContain");
410  fCosRejPIDFileContainXY= p.get< std::string >("CosRejPIDFileContainXY");
411  fCosRejPIDFileLight = p.get< std::string >("CosRejPIDFileLight");
412  fCosRejPIDFileFHCPeri = p.get< std::string >("CosRejPIDFileFHCPeri");
413  fCosRejPIDFileRHCPeri = p.get< std::string >("CosRejPIDFileRHCPeri");
414  fCosRejPIDFileFHCCore = p.get< std::string >("CosRejPIDFileFHCCore");
415  fCosRejPIDFileRHCCore = p.get< std::string >("CosRejPIDFileRHCCore");
416  fParticleIDAlgPSet = p.get< fhicl::ParameterSet >("ParticleIDAlgPSet");
417 }
418 
419 ////////////////////////////////////////////////////////////////////////
420 ////////////////////////////////////////////////////////////////////////
421 
423 {
424 
425  std::unique_ptr< std::vector< cosrej::NueCosRej > >
426  nueCosRejjies(new std::vector< cosrej::NueCosRej >);
427 
428  std::unique_ptr< art::Assns< cosrej::NueCosRej, rb::Cluster > >
430 
432 
433  // get the slices in the event
434 
436  evt.getByToken( fSliceToken, sliceHandle);
437  std::vector<art::Ptr<rb::Cluster> > slices;
438  art::fill_ptr_vector(slices, sliceHandle);
439 
440  // get things associated with slices
441 
443  prongAssn(sliceHandle, evt, art::InputTag(fProngLabel, fProngInstance));
444 
446  trackAssn(sliceHandle, evt, fCosmicTrackLabel);
447 
449  showerLidAssn(sliceHandle, evt, fShowerLIDLabel);
450 
452  vertexAssn(sliceHandle, evt, fVertexLabel);
453 
454  const int nslices = slices.size();
455  for( int isli = 0; isli< nslices; ++isli){
456  art::Ptr<rb::Cluster> thisSlice = slices.at(isli);
457  if( thisSlice->IsNoise())
458  continue;
459 
460  if(!prongAssn.isValid())
461  continue;
462 
463  cosrej::NueCosRej nuecosrej;
464 
465  std::vector<art::Ptr<rb::Prong>> prongs = prongAssn.at(isli);
466  if(prongs.empty()) continue;
467 
468  nuecosrej.SetProngTransMom( TransMomFraction(prongs) );
469 
470  double mintopStart = 9999;
471  double minbotStart = 9999;
472  double mineastStart = 9999;
473  double minwestStart = 9999;
474  double minfrontStart = 9999;
475  double minbackStart = 9999;
476 
477  double mintopStop = 9999;
478  double minbotStop = 9999;
479  double mineastStop = 9999;
480  double minwestStop = 9999;
481  double minfrontStop = 9999;
482  double minbackStop = 9999;
483 
484  const int npng = prongs.size();
485  for(int ipng = 0; ipng < npng; ++ipng){
486 
487  TVector3 start = prongs[ipng]->Start();
488  TVector3 stop = prongs[ipng]->Start() +
489  (prongs[ipng]->TotalLength() * prongs[ipng]->Dir() );
490 
491  mintopStart = std::min(mintopStart, livegeo->DistToTop(start));
492  mintopStop = std::min(mintopStop, livegeo->DistToTop(stop ));
493 
494  minbotStart = std::min(minbotStart, livegeo->DistToBottom(start));
495  minbotStop = std::min(minbotStop, livegeo->DistToBottom(stop ));
496 
497  mineastStart = std::min(mineastStart, livegeo->DistToEdgeXMinus(start));
498  mineastStop = std::min(mineastStop, livegeo->DistToEdgeXMinus(stop ));
499 
500  minwestStart = std::min(minwestStart, livegeo->DistToEdgeXPlus(start));
501  minwestStop = std::min(minwestStop, livegeo->DistToEdgeXPlus(stop ));
502 
503  minfrontStart = std::min(minfrontStart, livegeo->DistToFront(start));
504  minfrontStop = std::min(minfrontStop, livegeo->DistToFront(stop ));
505 
506  minbackStart = std::min(minbackStart, livegeo->DistToBack(start));
507  minbackStop = std::min(minbackStop, livegeo->DistToBack(stop ));
508  } // end for ipng
509 
510  nuecosrej.SetStartMinDistToTop( mintopStart );
511  nuecosrej.SetStartMinDistToBottom( minbotStart );
512  nuecosrej.SetStartMinDistToFront( minfrontStart );
513  nuecosrej.SetStartMinDistToBack( minbackStart );
514  nuecosrej.SetStartMinDistToEast( mineastStart );
515  nuecosrej.SetStartMinDistToWest( minwestStart );
516 
517  nuecosrej.SetStopMinDistToTop( mintopStop );
518  nuecosrej.SetStopMinDistToBottom( minbotStop );
519  nuecosrej.SetStopMinDistToFront( minfrontStop );
520  nuecosrej.SetStopMinDistToBack( minbackStop );
521  nuecosrej.SetStopMinDistToEast( mineastStop );
522  nuecosrej.SetStopMinDistToWest( minwestStop );
523 
524  // Find the most energetic prong
525  std::sort( prongs.begin(),
526  prongs.end(),
528 
529  art::Ptr< rb::Prong> leadingProng = prongs[0];
530 
531  if(trackAssn.isValid()){
532  int musliceidxbydist, musliceidxbytime;
533  double mintimebydist, minanglebydist, closestapproachbydist;
534  double mintimebytime, minanglebytime, closestapproachbytime;
535 
536  MuonParentByDist(leadingProng, trackAssn, isli, nslices,
537  mintimebydist, minanglebydist,
538  closestapproachbydist, musliceidxbydist);
539 
540  MuonParentByTime(leadingProng, trackAssn, isli, nslices,
541  mintimebytime, minanglebytime,
542  closestapproachbytime, musliceidxbytime);
543 
544  nuecosrej.SetMuSliceIdxByDist(musliceidxbydist);
545  nuecosrej.SetMuAngleDiffByDist(minanglebydist);
546  nuecosrej.SetMuTimeDiffByDist(mintimebydist);
547  nuecosrej.SetMuClosestApproachByDist(closestapproachbydist);
548 
549  nuecosrej.SetMuSliceIdxByTime(musliceidxbytime);
550  nuecosrej.SetMuAngleDiffByTime(minanglebytime);
551  nuecosrej.SetMuTimeDiffByTime(mintimebytime);
552  nuecosrej.SetMuClosestApproachByTime(closestapproachbytime);
553  }
554 
555  // Fill forward backward prong asymmetry variables
556  float sparsenessAsymm = -5, hitsperplaneAsymm = -5, sparsenessAsymmSlice = -5, hitsperplaneAsymmSlice = -5 ;
557 
558  SparsenessAsymmetry(leadingProng.get(),
559  sparsenessAsymm,
560  hitsperplaneAsymm);
561  nuecosrej.SetSparsenessAsymm(sparsenessAsymm);
562  nuecosrej.SetHitsPerPlaneAsymm(hitsperplaneAsymm);
563 
564  SparsenessAsymmetry(thisSlice.get(),
565  sparsenessAsymmSlice,
566  hitsperplaneAsymmSlice);
567  nuecosrej.SetSparsenessAsymmSlice(sparsenessAsymmSlice);
568  nuecosrej.SetHitsPerPlaneAsymmSlice(hitsperplaneAsymmSlice);
569 
570  TVector3 start = leadingProng->Start();
571  TVector3 stop = leadingProng->Start() +
572  (leadingProng->TotalLength() * leadingProng->Dir() );
573 
574  nuecosrej.SetStartDistToTop( livegeo->DistToTop(start) );
575  nuecosrej.SetStartDistToBottom( livegeo->DistToBottom(start) );
576  nuecosrej.SetStartDistToFront( livegeo->DistToFront(start) );
577  nuecosrej.SetStartDistToBack( livegeo->DistToBack(start) );
578  nuecosrej.SetStartDistToEast( livegeo->DistToEdgeXMinus(start) );
579  nuecosrej.SetStartDistToWest( livegeo->DistToEdgeXPlus(start) );
580 
581  nuecosrej.SetStopDistToTop( livegeo->DistToTop(stop) );
582  nuecosrej.SetStopDistToBottom( livegeo->DistToBottom(stop) );
583  nuecosrej.SetStopDistToFront( livegeo->DistToFront(stop) );
584  nuecosrej.SetStopDistToBack( livegeo->DistToBack(stop) );
585  nuecosrej.SetStopDistToEast( livegeo->DistToEdgeXMinus(stop) );
586  nuecosrej.SetStopDistToWest( livegeo->DistToEdgeXPlus(stop) );
587 
588  nuecosrej.SetProngMaxX( std::max(start.X(), stop.X()) );
589  nuecosrej.SetProngMaxY( std::max(start.Y(), stop.Y()) );
590  nuecosrej.SetProngMaxZ( std::max(start.Z(), stop.Z()) );
591  nuecosrej.SetProngMinX( std::min(start.X(), stop.X()) );
592  nuecosrej.SetProngMinY( std::min(start.Y(), stop.Y()) );
593  nuecosrej.SetProngMinZ( std::min(start.Z(), stop.Z()) );
594 
595  nuecosrej.SetHitsPerPlane((float)thisSlice->NCell()/
596  leadingProng->ExtentPlane());
597 
598  if(vertexAssn.isValid()){
599  const std::vector<art::Ptr<rb::Vertex>> vertCol =
600  vertexAssn.at(isli);
601  //currently only one vertex per slice
602  TVector3 vtx = vertCol[0]->GetXYZ();
603  // TODO: so far as I (CJB) can see fParticleIDAlg is never initialized,
604  // so presumably this code is never reached and can be removed.
606  start,
607  stop) );
608  }
609 
610  if(prongs.size() > 1){
611  // some cuts need the presence of two or more showers
612  art::Ptr<rb::Prong> png2 = prongs[1];
613 
614  nuecosrej.SetCosAngleToNextProng( leadingProng->Dir()
615  .Dot( png2->Dir() ) );
616  }
617  else{
618  nuecosrej.SetCosAngleToNextProng(-5);
619  }
620 
621 
622 
623  // Only execute when there are showers associated with the shwlids
624 
625  double widthofshower=0;
626  double shwlidCalE=0;
627 
628  if( showerLidAssn.isValid() ){
629  std::vector< const slid::ShowerLID* > shwlids = showerLidAssn.at(isli);
630  if( shwlids.empty() )
631  continue;
632 
633  sort(shwlids.begin(), shwlids.end(), slid::CompareByEnergy);
634  art::FindOneP<rb::Shower> fos(shwlids, evt, fShowerLIDLabel);
635 
636  widthofshower=shwlids[0]->Radius();
637 
638 
639  if(fos.at(0)){
640 
641  // widthofshower=ShowerlidWidth(shwlids, fos);
642 
643  shwlidCalE=fos.at(0)->CalorimetricEnergy();
644  TVector3 nuRecoP(0,0,0);
645  TVector3 nuRecoPLIDMass(0,0,0);
646  TransMomFractionShower( shwlids, fos, nuRecoP, nuRecoPLIDMass);
647  nuRecoP = nuRecoP.Unit(); // CosRej only cares about dirs of these vars
648  nuRecoPLIDMass = nuRecoPLIDMass.Unit();
650  const TVector3 beamDir = geom->NuMIBeamDirection();
651 
652  double nuPX = nuRecoP[0];
653  double nuPY = nuRecoP[1];
654  double nuPT = sqrt( 1 - util::sqr(beamDir.Dot(nuRecoP)) );
655  double nuPTLID = sqrt( 1 - util::sqr(beamDir.Dot(nuRecoPLIDMass)) );
656 
657  nuecosrej.SetPhotonShowerMomX(nuPX);
658  nuecosrej.SetPhotonShowerMomY(nuPY);
659  nuecosrej.SetParticleShowerTransMom( nuPT );
660  nuecosrej.SetPhotonShowerTransMom( nuPTLID );
661  }
662  }
663 
664 
665  // Fill the BDTs
666  // Need CVN, input to some BDTs
667  art::FindManyP<cvn::Result> fmCVN(slices, evt, fCVNLabel);
668  if (fmCVN.isValid() && fReaderContain && fReaderLight && fReaderContainXY){
669  std::vector<art::Ptr<cvn::Result> > cvns = fmCVN.at(isli);
670 
671  // Some containment vars already set up in nuecosrej
672  TMVAvars[0] = thisSlice->NCell();
673  TMVAvars[1] = sparsenessAsymmSlice; // Defined above
674  TMVAvars[2] = sparsenessAsymm; // Defined above
675  TMVAvars[3] = nuecosrej.PhotonShowerTransMom(); // Grab the massless one
676  TMVAvars[4] = nuecosrej.PhotonShowerMomX();
677  TMVAvars[5] = nuecosrej.PhotonShowerMomY();
678  TMVAvars[6] = std::min(nuecosrej.StartMinDistToTop(),
679  nuecosrej.StopMinDistToTop());
680  // Need algebra, nuecosrej dosn't store this var
681  double distnottop = nuecosrej.StartMinDistToBottom();
682  distnottop = std::min(distnottop,nuecosrej.StopMinDistToBottom());
683  distnottop = std::min(distnottop,nuecosrej.StartMinDistToFront());
684  distnottop = std::min(distnottop,nuecosrej.StopMinDistToFront());
685  distnottop = std::min(distnottop,nuecosrej.StartMinDistToBack());
686  distnottop = std::min(distnottop,nuecosrej.StopMinDistToBack());
687  distnottop = std::min(distnottop,nuecosrej.StartMinDistToEast());
688  distnottop = std::min(distnottop,nuecosrej.StopMinDistToEast());
689  distnottop = std::min(distnottop,nuecosrej.StartMinDistToWest());
690  distnottop = std::min(distnottop,nuecosrej.StopMinDistToWest());
691  TMVAvars[7] = distnottop;
692 
693  // CVN values
694  float cvne = 0;
695  if(!cvns.empty()){
696  cvne = cvns[0]->fOutput[1];
697  }
698  TMVAvars[8] = cvne;
699 
700  // The rest are true vars only. Ignore
701  TMVAvars[9] = 0;
702  TMVAvars[10]= 0;
703  TMVAvars[11]= 0;
704  TMVAvars[12]= 0;
705 
706  nuecosrej.SetCosPIDContain(fReaderContain->EvaluateMVA("BDT_Contain"));
707  nuecosrej.SetCosPIDContainXY(fReaderContainXY->EvaluateMVA("BDTContain2018"));
708  nuecosrej.SetCosPIDLight(fReaderLight ->EvaluateMVA("BDTLight2018"));
709 
710 
711 
712 
713  CoreTMVAvars[0] = thisSlice->NCell();
714  CoreTMVAvars[1] = sparsenessAsymm; // Defined above
715  CoreTMVAvars[2] = nuecosrej.PhotonShowerTransMom(); // Grab the massless one
716  CoreTMVAvars[3] = std::min(nuecosrej.StartMinDistToTop(),
717  nuecosrej.StopMinDistToTop());
718  CoreTMVAvars[4] = std::min(nuecosrej.StartMinDistToBack(),
719  nuecosrej.StopMinDistToBack());
720  CoreTMVAvars[5] = std::min(nuecosrej.StartMinDistToWest(),
721  nuecosrej.StopMinDistToWest());
722  CoreTMVAvars[6] = std::min(nuecosrej.StartMinDistToFront(),
723  nuecosrej.StopMinDistToFront());
724  CoreTMVAvars[7]= std::min(nuecosrej.StartMinDistToEast(),
725  nuecosrej.StopMinDistToEast());
726  CoreTMVAvars[8]= std::min(nuecosrej.StartMinDistToBottom(),
727  nuecosrej.StopMinDistToBottom());
728  double inel = (thisSlice->CalorimetricEnergy()-shwlidCalE)/thisSlice->CalorimetricEnergy();// inelasticity=(sliceCalE-shwlidCalE)/sliceCalE
729  CoreTMVAvars[9] = inel;
730 
731 
732  CoreTMVAvars[10]=widthofshower;
733  //spectator vars set zero
734  CoreTMVAvars[11]= 0;//cvnsse
735  CoreTMVAvars[12]= 0;//woscdumb
736  CoreTMVAvars[13]= 0;//isGENIE
737  CoreTMVAvars[14]= 0;//pdg
738  CoreTMVAvars[15]= 0;//pdgorig
739  CoreTMVAvars[16]= 0;//isCC
740 
741 
742 
743  if(!IsRHC(evt))
744  nuecosrej.SetCoreBDT(fReaderFHCCoreBDT ->EvaluateMVA("FHCCoreBDT"));
745  else if(IsRHC(evt))
746  nuecosrej.SetCoreBDT(fReaderRHCCoreBDT ->EvaluateMVA("RHCCoreBDT"));
747 
748 
749  TVector3 vtxXYZ(0,0,0);
750 
751  if(vertexAssn.isValid()){
752  const std::vector<art::Ptr<rb::Vertex>> vertCol = vertexAssn.at(isli);//currently only one vertex per slice
753  vtxXYZ = vertCol[0]->GetXYZ();
754  }
755 
756 
757 
758  PeriTMVAvars[0] = thisSlice->NCell();
759  PeriTMVAvars[1] = sparsenessAsymm; // Defined above
760  PeriTMVAvars[2] = nuecosrej.PhotonShowerTransMom(); // Grab the massless one
761  PeriTMVAvars[3] = nuecosrej.PhotonShowerMomX();
762  PeriTMVAvars[4] = nuecosrej.PhotonShowerMomY();
763  PeriTMVAvars[5] = std::min(nuecosrej.StartMinDistToTop(),
764  nuecosrej.StopMinDistToTop());
765  PeriTMVAvars[6] = distnottop;
766  PeriTMVAvars[7]= vtxXYZ[0];//vtxx
767  PeriTMVAvars[8]= vtxXYZ[1];//vtxy
768  PeriTMVAvars[9]= vtxXYZ[2];//vtxz
769 
770  //spectator vars set zero
771  PeriTMVAvars[10]= 0;//cvnsse
772  PeriTMVAvars[11]= 0;//woscdumb
773  PeriTMVAvars[12]= 0;//isGENIE
774  PeriTMVAvars[13]= 0;//pdg
775  PeriTMVAvars[14]= 0;//pdgorig
776  PeriTMVAvars[15]= 0;//isCC
777 
778  if(!IsRHC(evt))
779  nuecosrej.SetPeriBDT(fReaderFHCPeriBDT ->EvaluateMVA("FHCPeriBDT"));
780  else if(IsRHC(evt))
781  nuecosrej.SetPeriBDT(fReaderRHCPeriBDT ->EvaluateMVA("RHCPeriBDT"));
782 
783 
784  }
785  nueCosRejjies->push_back(nuecosrej);
786  util::CreateAssn( *this, evt, *nueCosRejjies, thisSlice ,*sliceAssn);
787 
788  }// end loop over slices
789 
790  evt.put( std::move(sliceAssn) );
791  evt.put( std::move(nueCosRejjies) );
792 
793 }// end of producer
794 
795 ////////////////////////////////////////////////////////////////////////
796 ////////////////////////////////////////////////////////////////////////
797 
798 double cosrej::MakeNueCosRej::TransMomFraction(const std::vector< art::Ptr<rb::Prong> >& prongs) const
799 {
800  if(prongs.empty()) return 0;
801 
802  TVector3 ret;
803  for(const art::Ptr<rb::Prong>& prong: prongs){
804  if(!prong->Is3D()) continue;
805  const double w = prong->CalorimetricEnergy();
806  ret += w*prong->Dir();
807  }
808 
809  ret = ret.Unit();
810 
812  const TVector3 beamDir = geom->NuMIBeamDirection();
813  const double beamProj = ret.Dot(beamDir);
814 
815  return (ret-beamProj*beamDir).Mag();
816 }
817 
818 ////////////////////////////////////////////////////////////////////////
819 ////////////////////////////////////////////////////////////////////////
820 
821 void cosrej::MakeNueCosRej::TransMomFractionShower(std::vector< const slid::ShowerLID* > shwlids,
823  TVector3 &nuRecoP,
824  TVector3 &nuRecoPLIDMass)
825 {
826 
827  int nsh = shwlids.size();
828 
829  TLorentzVector nuP, nuP0;
830 
831  for(int ish =0; ish < nsh; ++ish){
832  cet::maybe_ref<art::Ptr<rb::Shower> const> shwref(fos.at(ish));
833  art::Ptr<rb::Shower> shw = shwref.ref();
834  if(!shw)
835  continue;
836 
837  // figure out for which particle the LL is maximum
838  std::map<int, float> partLongLL = shwlids[ish]->fPartLongLL;
839  std::map<int, float> partTransLL = shwlids[ish]->fPartTransLL;
840 
841  int maxllType = 0;
842  int maxll = -1;
843 
844  for( auto const & itype : partLongLL){
845  if(itype.second + partTransLL[itype.first] > maxll){
846  maxllType = itype.first;
847  maxll = itype.second + partTransLL[itype.first];
848  }
849  }
850 
851 
852  TVector3 dir = shw->Dir();
853  double showerE = shwlids[ish]->ShowerEnergy();
854  double mass = Mass( maxllType );
855  // if the mass is non zero, the measured energy
856  // is the kineticenergy of the particle.
857  // so the momentum is p = (T^2 + 2Tm )^1/2
858  double momMagMass = sqrt( (showerE*showerE) + (2*showerE*mass) );
859 
860  nuRecoP += showerE*dir;
861  nuRecoPLIDMass += momMagMass*dir;
862 
863  }// end loop over showerlids
864 
865 }
866 
867 
868 ////////////////////////////////////////////////////////////////////////
869 ////////////////////////////////////////////////////////////////////////
870 
872 {
873  int pdg;
875  pdg = 11;
876  else if (type == slid::DedxParticleType::kPHOTON)
877  pdg = 22;
878  else if (type == slid::DedxParticleType::kMUON)
879  pdg = 13;
880  else if (type == slid::DedxParticleType::kPI0)
881  pdg = 111;
882  else if (type == slid::DedxParticleType::kPION)
883  pdg = 211;
884  else if (type == slid::DedxParticleType::kPROTON)
885  pdg = 2212;
886  else if (type == slid::DedxParticleType::kNEUTRON)
887  pdg = 2112;
888  else
889  return 0;
890 
891  static TDatabasePDG* pdgdb = TDatabasePDG::Instance();
892  return pdgdb->GetParticle(pdg)->Mass();
893 }
894 
895 
896 
897 
898 
899 ////////////////////////////////////////////////////////////////////////
900 
902  float & sparsenessAsymm,
903  float & hitsPerPlaneAsymm)
904 {
905 
907 
908  art::PtrVector<rb::CellHit> cells = cluster->AllCells();
909  unsigned int ncells = cells.size();
910 
911  if( ncells < 1 )
912  return;
913 
914  std::set<int> setplanes;
915  int minplane = cluster->MinPlane();
916  int maxplane = cluster->MaxPlane();
917 
918  // create a set of planes first; it is sorted
919  // and has unique entries
920  for(unsigned int i = 0; i < ncells; ++i)
921  setplanes.insert( cells[i]->Plane());
922 
923  if(setplanes.empty())
924  return;
925 
926  std::vector<int> planes;
927  planes.assign(setplanes.begin(), setplanes.end());
928  int last = planes.size()-1;
929 
930  float hitsperplanefront = 0;
931  float hitsperplaneback = 0;
932  float sparsenessfront = 0;
933  float sparsenessback = 0;
934 
935  sparsenessAsymm = 0;
936 
937  if( planes.size() > 8){
938  sparsenessfront = planes[7]-planes[0]-7;
939  sparsenessback = planes[last]-planes[last-7]-7;
940 
941  if(sparsenessback == 0 && sparsenessfront == 0)
942  sparsenessAsymm = 0;
943  else
944  sparsenessAsymm = (sparsenessfront - sparsenessback)/
945  (sparsenessfront + sparsenessback);
946  }
947 
948  for(unsigned int i = 0; i < ncells; ++i){
949  int pl = cells[i]->Plane();
950 
951  if(pl < minplane+5)
952  hitsperplanefront ++;
953 
954  if(pl > maxplane-5)
955  hitsperplaneback ++;
956 
957  }// end loop over cells
958 
959  hitsperplaneback = hitsperplaneback/5;
960  hitsperplanefront = hitsperplanefront/5;
961 
962  hitsPerPlaneAsymm = (hitsperplanefront - hitsperplaneback)/
963  (hitsperplanefront + hitsperplaneback);
964 
965  return;
966 }// End of SparsenessAsymmetry
967 
968 ////////////////////////////////////////////////////////////////////////
969 
970 
972  const art::Ptr<rb::Prong> b)
973 {
974  return a->CalorimetricEnergy() > b->CalorimetricEnergy();
975 }
976 
977 
979 {
981  if (!evt.isRealData())
982  {
983  evt.getByLabel(fGeneratorLabel, spillPot);
984  }else
985  {
986  evt.getByLabel(fNuMILabel, spillPot);
987  }
988  if (spillPot.failedToGet())
989  {
990  mf::LogError("NueCosRej::spillpot") <<
991  "Spill Data not found, aborting without horn current information";
992  abort();
993  }
994 
995  return spillPot->isRHC;
996 }
997 
998 
999 ////////////////////////////////////////////////////////////////////////
1000 
1002  art::FindManyP<rb::Track> trkAssn,
1003  const int& idx,
1004  const int& nsli,
1005  double& mintimebydist,
1006  double& minanglebydist,
1007  double& closestapproachbydist,
1008  int& musliceidxbydist)
1009 {
1010 
1012 
1013  closestapproachbydist = 1.0e10;
1014  minanglebydist = 1.0e10;
1015  mintimebydist = 1.0e10;
1016  musliceidxbydist = -999;
1017 
1018  TVector3 stop = leadingPng->Start() +
1019  leadingPng->TotalLength()*leadingPng->Dir();
1020 
1021 
1022  for(int jsli = 0; jsli < nsli; jsli++){
1023  if(jsli == idx)
1024  continue;
1025 
1026  std::vector<art::Ptr<rb::Track>> tracks = trkAssn.at(jsli);
1027  if(tracks.empty()) continue;
1028 
1029 
1030  art::Ptr< rb::Track> leadingtrk = tracks[0];
1031 
1032  double dist = geo::ClosestApproach( leadingPng->Start(),
1033  stop,
1034  leadingtrk->Start(),
1035  leadingtrk->Stop());
1036 
1037  if( dist < closestapproachbydist ){
1038  mintimebydist = std::abs(leadingPng->MeanTNS() - leadingtrk->MeanTNS());
1039  minanglebydist =
1040  leadingPng->Dir().Angle(leadingtrk->Dir())*TMath::RadToDeg();
1041  closestapproachbydist = dist;
1042  musliceidxbydist = jsli;
1043  }
1044  }// end loop over slices
1045  return;
1046 }// End of MuonParentByDist
1047 
1048 
1050  art::FindManyP<rb::Track> trkAssn,
1051  const int& idx,
1052  const int& nsli,
1053  double& mintimebytime,
1054  double& minanglebytime,
1055  double& closestapproachbytime,
1056  int& musliceidxbytime)
1057 {
1059 
1060  closestapproachbytime = 1.0e10;
1061  minanglebytime = 1.0e10;
1062  mintimebytime = 1.0e10;
1063  musliceidxbytime = -999;
1064 
1065  TVector3 stop = leadingPng->Start() +
1066  leadingPng->TotalLength()*leadingPng->Dir();
1067 
1068  for(int jsli = 0; jsli < nsli; jsli++){
1069  if(jsli == idx) continue;
1070 
1071  std::vector<art::Ptr<rb::Track>> tracks = trkAssn.at(jsli);
1072  if(tracks.empty()) continue;
1073 
1074  art::Ptr< rb::Track> leadingtrk = tracks[0];
1075 
1076  double time = std::abs(leadingPng->MeanTNS() - leadingtrk->MeanTNS());
1077 
1078  if( time < mintimebytime ){
1079  musliceidxbytime = jsli;
1080  mintimebytime = time;
1081  minanglebytime =
1082  leadingPng->Dir().Angle(leadingtrk->Dir())*TMath::RadToDeg();
1083  closestapproachbytime = geo::ClosestApproach( leadingPng->Start(),
1084  stop,
1085  leadingtrk->Start(),
1086  leadingtrk->Stop());
1087  }
1088  }// end loop over slices
1089  return;
1090 }// End of MuonParentByTime
1091 
1092 
1093 
1094 
1095 
1096 ////////////////////////////////////////////////////////////////////////
1097 
void SetMuClosestApproachByTime(double input)
Definition: NueCosRej.h:409
T max(const caf::Proxy< T > &a, T b)
void SetMuAngleDiffByDist(double input)
Definition: NueCosRej.h:394
TMVA::Reader * fReaderRHCPeriBDT
bool isRHC
is the beam in antineutrino mode, aka RHC
Definition: SpillData.h:28
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.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void SetProngMaxY(double input)
Definition: NueCosRej.h:361
double Mass(const int type)
void SetParticleShowerTransMom(double input)
Definition: NueCosRej.h:251
void MuonParentByTime(const art::Ptr< rb::Prong > &leadingPng, art::FindManyP< rb::Track > pngAssn, const int &idx, const int &nslc, double &mintimebytime, double &minanglebydist, double &closestapproachbytime, int &musliceidxbytime)
void SetStartMinDistToWest(double input)
Definition: NueCosRej.h:312
double StopMinDistToBottom() const
Minimum perpendicular distance of all prongs stop points to the bottom edge of the detector...
Definition: NueCosRej.h:124
void SetPhotonShowerMomX(double input)
Definition: NueCosRej.h:257
void SetStartDistToFront(double input)
Definition: NueCosRej.h:269
void SetStartDistToEast(double input)
Definition: NueCosRej.h:278
fhicl::ParameterSet fParticleIDAlgPSet
FCL parameters for particle ID alorithm (loglikelihoods and dE/dx)
const char * p
Definition: xmltok.h:285
double StopMinDistToWest() const
Minimum perpendicular distance of all prongs stop points to the west edge of the detector (west is po...
Definition: NueCosRej.h:136
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SetHitsPerPlaneAsymmSlice(double input)
Definition: NueCosRej.h:385
double StartMinDistToTop() const
Minimum perpendicular distance of all prongs start points.
Definition: NueCosRej.h:96
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
A collection of associated CellHits.
Definition: Cluster.h:47
void SetStartMinDistToTop(double input)
Definition: NueCosRej.h:300
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void beginRun(art::Run &evt) override
bool IsRHC(const art::Event &evt)
double DistToFront(TVector3 vertex)
double DistToBack(TVector3 vertex)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
bool isRealData() const
Definition: Event.h:83
double PhotonShowerTransMom() const
Transverse component of the event momentum. Particle masses in the momentum calculation are assumed t...
Definition: NueCosRej.h:36
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
void SetMuTimeDiffByTime(double input)
Definition: NueCosRej.h:403
TMVA::Reader * fReaderFHCPeriBDT
::xsd::cxx::tree::time< char, simple_type > time
Definition: Database.h:194
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
void SetMuSliceIdxByTime(int input)
Definition: NueCosRej.h:400
bool CompareByEnergy(const slid::ShowerLID *a, const slid::ShowerLID *b)
Definition: ShowerLID.cxx:51
float abs(float number)
Definition: d0nt_math.hpp:39
void SetStopMinDistToEast(double input)
Definition: NueCosRej.h:333
void MuonParentByDist(const art::Ptr< rb::Prong > &leadingPng, art::FindManyP< rb::Track > pngAssn, const int &idx, const int &nslc, double &mintimebydist, double &minanglebydist, double &closestapproachbydist, int &musliceidxbydist)
bool CompareByCalEnergy(const art::Ptr< rb::Prong > a, const art::Ptr< rb::Prong > b)
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
const art::ProductToken< std::vector< rb::Cluster > > fSliceToken
Definition: Run.h:31
void SetCosPIDContainXY(double input)
Definition: NueCosRej.h:339
::xsd::cxx::tree::type type
Definition: Database.h:110
void SetStopMinDistToFront(double input)
Definition: NueCosRej.h:324
double dist
Definition: runWimpSim.h:113
void SetCosPIDLight(double input)
Definition: NueCosRej.h:342
void SetPhotonShowerTransMom(double input)
Definition: NueCosRej.h:254
void reconfigure(const fhicl::ParameterSet &pset)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void SetStartDistToBack(double input)
Definition: NueCosRej.h:272
void SetProngTransMom(double input)
Definition: NueCosRej.h:248
void SetStopDistToFront(double input)
Definition: NueCosRej.h:287
double DistToTop(TVector3 vertex)
double StopMinDistToTop() const
Minimum perpendicular distance of all prongs stop points to the top edge of the detector.
Definition: NueCosRej.h:120
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Definition: Prong.cxx:186
double StartMinDistToEast() const
Minimum perpendicular distance of all prongs start points to the east edge of the detector (east is n...
Definition: NueCosRej.h:116
double StopMinDistToEast() const
Minimum perpendicular distance of all prongs stop points to the east edge of the detector (east is ne...
Definition: NueCosRej.h:140
double DistToBottom(TVector3 vertex)
Result for CVN.
double CalorimetricEnergy(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple estimate of neutrino energy.
Definition: Cluster.cxx:439
double StartMinDistToBottom() const
Minimum perpendicular distance of all prongs start points to the bottom edge of the detector...
Definition: NueCosRej.h:100
double StartMinDistToBack() const
Minimum perpendicular distance of all prongs start points to the back edge of the detector...
Definition: NueCosRej.h:108
void produce(art::Event &e) override
const double a
void SetMuAngleDiffByTime(double input)
Definition: NueCosRej.h:406
void SetProngMinX(double input)
Definition: NueCosRej.h:367
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
void SetProngMinY(double input)
Definition: NueCosRej.h:370
slid::ParticleIDAlg * fParticleIDAlg
Particle ID alorithm (loglikelihoods and dE/dx)
void SetProngDistToVtx(double input)
Definition: NueCosRej.h:355
void SparsenessAsymmetry(const rb::Cluster *cluster, float &sparsenessAsymm, float &hitsPerPlaneAsymm)
void SetStartDistToBottom(double input)
Definition: NueCosRej.h:266
Collect Geo headers and supply basic geometry functions.
void SetSparsenessAsymm(double input)
Definition: NueCosRej.h:376
void SetStopDistToEast(double input)
Definition: NueCosRej.h:296
void SetMuTimeDiffByDist(double input)
Definition: NueCosRej.h:391
double StartMinDistToWest() const
Minimum perpendicular distance of all prongs start points to the west edge of the detector (west is p...
Definition: NueCosRej.h:112
double StopMinDistToFront() const
Minimum perpendicular distance of all prongs stop points to the front edge of the detector...
Definition: NueCosRej.h:128
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
double PointDoca(TVector3 vtx, TVector3 start, TVector3 stop)
Get distance of closest approach between shower and vertex.
void SetStopDistToTop(double input)
Definition: NueCosRej.h:281
void SetProngMaxZ(double input)
Definition: NueCosRej.h:364
void SetHitsPerPlane(double input)
Definition: NueCosRej.h:245
double PhotonShowerMomY() const
The cosine between the reconstructed momentum and the y-axis. Particle masses in the momentum calcula...
Definition: NueCosRej.h:44
void SetStopMinDistToBottom(double input)
Definition: NueCosRej.h:321
void SetStartMinDistToEast(double input)
Definition: NueCosRej.h:315
void SetProngMaxX(double input)
Definition: NueCosRej.h:358
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
double StartMinDistToFront() const
Minimum perpendicular distance of all prongs start points to the front edge of the detector...
Definition: NueCosRej.h:104
void SetStopMinDistToWest(double input)
Definition: NueCosRej.h:330
void SetPhotonShowerMomY(double input)
Definition: NueCosRej.h:260
void SetStartMinDistToFront(double input)
Definition: NueCosRej.h:306
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
TVector3 NuMIBeamDirection() const
Direction of neutrinos from the NuMI beam (unit vector)
void SetCosPIDContain(double input)
Definition: NueCosRej.h:336
Cosmic Rejection PIDs for Numu analysis.
Definition: FillParentInfo.h:9
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetStartMinDistToBack(double input)
Definition: NueCosRej.h:309
void SetStopDistToWest(double input)
Definition: NueCosRej.h:293
Calculate dE/dx and log likelihoods for different particle hypotheses. This information will be of us...
double DistToEdgeXMinus(TVector3 vertex)
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
Definition: Geo.cxx:222
void SetStopMinDistToBack(double input)
Definition: NueCosRej.h:327
T const * get() const
Definition: Ptr.h:321
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void SetStopMinDistToTop(double input)
Definition: NueCosRej.h:318
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
void geom(int which=0)
Definition: geom.C:163
const hit & b
Definition: hits.cxx:21
double PhotonShowerMomX() const
The cosine between the reconstructed momentum and the x-axis. Particle masses in the momentum calcula...
Definition: NueCosRej.h:40
void SetStopDistToBottom(double input)
Definition: NueCosRej.h:284
int ncells
Definition: geom.C:124
void SetStopDistToBack(double input)
Definition: NueCosRej.h:290
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
void SetStartMinDistToBottom(double input)
Definition: NueCosRej.h:303
void SetCosAngleToNextProng(double input)
Definition: NueCosRej.h:352
bool file_exists(std::string const &qualified_filename)
void SetMuSliceIdxByDist(int input)
Definition: NueCosRej.h:388
void TransMomFractionShower(std::vector< const slid::ShowerLID * > shwlids, art::FindOneP< rb::Shower > fos, TVector3 &nuRecoP, TVector3 &nuRecoPLIDMass)
MakeNueCosRej(fhicl::ParameterSet const &p)
double StopMinDistToBack() const
Minimum perpendicular distance of all prongs stop points to the back edge of the detector.
Definition: NueCosRej.h:132
T min(const caf::Proxy< T > &a, T b)
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 fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
bool getByToken(ProductToken< PROD > const &token, Handle< PROD > &result) const
Definition: DataViewImpl.h:387
Float_t e
Definition: plot.C:35
void SetSparsenessAsymmSlice(double input)
Definition: NueCosRej.h:382
void SetProngMinZ(double input)
Definition: NueCosRej.h:373
Helper for AttenCurve.
Definition: Path.h:10
Float_t w
Definition: plot.C:20
ProductToken< T > consumes(InputTag const &)
Definition: fwd.h:28
void SetPeriBDT(double input)
Definition: NueCosRej.h:345
Encapsulate the geometry of one entire detector (near, far, ndos)
void SetMuClosestApproachByDist(double input)
Definition: NueCosRej.h:397
bool failedToGet() const
Definition: Handle.h:196
void SetCoreBDT(double input)
Definition: NueCosRej.h:348
double TransMomFraction(const std::vector< art::Ptr< rb::Prong > > &prongs) const
void SetStartDistToWest(double input)
Definition: NueCosRej.h:275
void SetStartDistToTop(double input)
Definition: NueCosRej.h:263
void SetHitsPerPlaneAsymm(double input)
Definition: NueCosRej.h:379
double DistToEdgeXPlus(TVector3 vertex)