MuonRemove_module.cc
Go to the documentation of this file.
1 
2 ///////////////////////////////////////////////////////////////////////
3 /// \file MuonRemove.cxx
4 /// \brief Remove raw digits corresponding to muon hits.
5 /// \author Kanika Sachdev (ksachdev@physics.umn.edu)
6 ////////////////////////////////////////////////////////////////////////
7 
8 
9 
10 // C/C++ includes
11 #include <cmath>
12 #include <iostream>
13 #include <map>
14 #include <vector>
15 #include <algorithm>
16 #include <list>
17 // ROOT includes
18 #include "TMath.h"
19 #include "TVector3.h"
20 #include "TDatabasePDG.h"
21 
22 // Framework includes
23 #include "fhiclcpp/ParameterSet.h"
33 #include "CVN/func/Result.h"
35 #include "RecoBase/FilterList.h"
36 
37 // NOvASoft includes
38 #include "NovaDAQConventions/DAQConventions.h"
39 #include "DAQDataFormats/NanoSliceVersionConvention.h"
41 #include "Geometry/Geometry.h"
44 #include "MCCheater/BackTracker.h"
45 #include "Simulation/Particle.h"
46 #include "Simulation/FLSHitList.h"
47 #include "Simulation/FLSHit.h"
49 #include "RawData/RawDigit.h"
50 #include "RecoBase/CellHit.h"
51 #include "RecoBase/WeightedHit.h"
52 #include "RecoBase/Track.h"
53 #include "RecoBase/Cluster.h"
54 #include "RecoBase/Vertex.h"
55 #include "RecoBase/FitSum.h"
56 #include "Calibrator/Calibrator.h"
57 #include "Utilities/AssociationUtil.h"
58 #include "ReMId/ReMId.h"
59 #include "NumuEnergy/NumuEAlg.h"
60 #include "NumuEnergy/NumuE.h"
61 
62 namespace murem
63 {
64  class MuonRemove : public art::EDProducer
65  {
66  public:
67 
68  explicit MuonRemove(fhicl::ParameterSet const& pset);
69  ~MuonRemove();
70 
71  void beginRun(art::Run& run);
72 
73  void produce(art::Event& evt);
74 
75  void reconfigure(const fhicl::ParameterSet& pset);
76 
77  ///\brief Returns a vector of Cellhits attributed to the muon by checking with truth.
78  /// For implementing muon removal by truth
79  std::vector< rb::WeightedHit > RemoveByTruth(const art::Ptr<rb::Cluster>& slice);
80 
81  ///\brief A function that returns the longest track amongst a vector of tracks
82  rb::Track LongestTrack( std::vector< art::Ptr<rb::Track> >& tracks, art::Event& evt, double& muonEnergy);
83 
84  ///\brief Returns the highest efficiency (by truth) muon track for removal;
85  rb::Track RemoveByEfficiency( std::vector< art::Ptr<rb::Track> >& tracks,
86  const rb::Cluster &slice, art::Event& evt, double& muonEnergy);
87 
88  ///\brief Returns the track with the highest PID value returned by the ReMId package
89  rb::Track RemoveByReMId( std::vector< art::Ptr<rb::Track> >& tracks,
90  art::Event& evt, double& muonEnergy);
91  ///\brief Returns the BPF muon track with the highest Prong CVN value
92  rb::Track RemoveByProngCVN( std::vector< art::Ptr<rb::Track> >& tracks, const art::Ptr<rb::Cluster>& slice,
93  art::Event& evt, double& muonEnergy);
94 
95  ///\brief Function for preselecting the slices for the further mrcc process
96  bool SliceSelection ( const art::Event& evt, const art::Ptr<rb::Cluster>& slice);
97 
98  ///\brief Saves the information about the removed muon in a simb::MCParticle object
99  /// The energy of the muon is obtained from NumuEAlg.
100  simb::MCParticle MuonInfo(rb::Track muonTrack, double muonEnergy);
101 
102  ///\brief Saves muon removed FLSHitLists, PhotonSignals, Particle list and
103  /// information about the removed muon in a simb::MCParticle.
104  void FillTruthInfo( const art::Event& evt,
105  std::unique_ptr< std::vector<sim::FLSHitList> >& mrccFLSHits,
106  std::unique_ptr< std::vector<sim::PhotonSignal> >& mrccPhotSig,
107  std::unique_ptr< std::vector<sim::Particle> >& mrccPartList,
108  std::unique_ptr< std::vector<simb::MCParticle> >& muonInfo) ;
109 
110  private:
111 
112  /// An object of TrackCleanUpAlg
114 
115  /// ParameterSet of TrackCleanUpAlg
117 
118  /// An object of NumuEAlg
120 
121  /// ParameterSet of NumuEAlg
123 
124  /// A FHiCL flag to turn of track clean up. Default if true.
126 
127  /// ADC thresholds
128 
133 
134  // Configurable module parameters
135  /// Label of module that produced verteces
137  /// Label of module that produced prongs
139  /// Label of module that produced Prong CVN PIDs
141  /// Label of module that produced BPF tracks
143  /// Label of module that produced slice CVN PID
145 
146  /// Label of module which will be used for track loading
148 
149  /// Label of module that produced raw digits
151  /// Label of module that produced tracks
153  /// Label of module that produced slice clusters
155  /// Label of module that produced muon pid object
157  /// Label of module that produced MC truth information
159  /// Label of module that produced Numu PID
161  /// Label of module that produced Energy object (in case of Kalman tracks)
163  /// Label of module that produced Energy object (in case of BPF tracks)
165 
166  /// Numu PID cut
167  double fNumuPIDCut;
168  /// DCS Distance: multipoint readout look back
170 
171  /// Handle to Geometry
173  };
174 
175  /// \brief Takes in a vector of raw digits
176  /// and sorts it by plane.
177  void SortDigitsByPlaneAndCell( std::list< rawdata::RawDigit >& digitcol );
178 
179  /// \brief Compare two rawdigits by plane. Return true if
180  /// plane of ra is greater then the plane of rb.
183 
184  void SortByPlaneAndCell( std::vector< rb::WeightedHit >& whits );
185 
186  /// \brief Compare two rawdigits by plane. Return true if
187  /// plane of ra is greater then the plane of rb.
189  rb::WeightedHit whitb);
190 
191 }
192 
193 
194 namespace murem
195 {
196 
197  /**********************************************************************/
198  ////////////////////////////////////////////////////////////////////////
199 
201  : fTrkCleanUpAlg(0)
202  , fTrkCleanUpPSet(pset.get< fhicl::ParameterSet >("TrkCleanUpPSet") )
203  , fNumuEAlg(0)
204  , fNumuEAlgPSet(pset.get< fhicl::ParameterSet >("NumuEAlgPSet") )
205  {
206  this->reconfigure(pset);
207 
208  // Produces a vector of muon removed raw digit
209  produces< std::vector <rawdata::RawDigit> >();
210  produces< std::vector< simb::MCParticle> >();
211  produces< art::Assns<simb::MCParticle, rb::Cluster> >();
212 
213  // If removing muon by truth, also want to store the mrcc
214  // FLSHit list to make BackTracking possible
215  if(fMuonID == "removeByTruth"){
216  produces< std::vector <sim::FLSHitList> >();
217  produces< std::vector <sim::PhotonSignal> >();
218  produces< std::vector <sim::Particle> >();
219  }
220 
221  mf::LogInfo("MuonRemove") << " MuonRemove::MuonRemove()\n";
222  }
223 
224  /**********************************************************************/
225  ////////////////////////////////////////////////////////////////////////
226 
228  {
229  if(fTrkCleanUpAlg) delete fTrkCleanUpAlg;
230  if(fNumuEAlg) delete fNumuEAlg;
231  }
232 
233  /**********************************************************************/
234  ////////////////////////////////////////////////////////////////////////
235 
237  {
238  fMCTruthModuleLabel = pset.get< std::string >("MCTruthModuleLabel");
239  fClusterModuleLabel = pset.get< std::string >("ClusterModuleLabel");
240  fKalmanTrackModuleLabel = pset.get< std::string >("KalmanTrackModuleLabel");
241  fRawDataLabel = pset.get< std::string >("RawDataLabel", "daq");
242  fMuonID = pset.get< std::string >("MuonID");
243  fNumuReMIdLabel = pset.get< std::string >("NumuReMIdLabel");
244  fEnergyLabel = pset.get< std::string >("EnergyLabel");
245  fEnergyLabelBPF = pset.get< std::string >("EnergyLabelBPF");
246  fElasticInput = pset.get< std::string >("ElasticInput");
247  fFuzzy3DInput = pset.get< std::string >("Fuzzy3DInput");
248  fCVNParticleLabel = pset.get< std::string >("CVNParticleLabel");
249  fBPFTrackModuleLabel = pset.get< std::string >("BPFTrackModuleLabel");
250  fWhichTrackModule = pset.get< std::string >("WhichTrackModule");
251  fCVNLabel = pset.get< std::string >("CVNLabel");
252 
253  fTrackCleanUp = pset.get< bool >("TrackCleanUp");
254  fADCThreshold_Far = pset.get< float >("ADCThreshold_Far");
255  fADCThreshold_NDOS = pset.get< float >("ADCThreshold_NDOS");
256  fADCThreshold_Near = pset.get< float >("ADCThreshold_Near");
257  fNumuPIDCut = pset.get< double >("NumuPIDCut");
258  fDCSDistance = pset.get< int >("DCSDistance");
259  }
260 
261  /**********************************************************************/
262  ////////////////////////////////////////////////////////////////////////
263 
265  {
266 
269 
270 
271  novadaq::cnv::DetId detID = fgeom->DetId();
272  if (detID == novadaq::cnv::kNEARDET)
274  if (detID == novadaq::cnv::kFARDET)
276  if (detID == novadaq::cnv::kNDOS)
278 
279  return;
280  }
281 
282  /**********************************************************************/
283  ////////////////////////////////////////////////////////////////////////
284 
286  {
287 
288 
289  std::unique_ptr< std::vector<rawdata::RawDigit> >
290  muonRemovedDigits(new std::vector<rawdata::RawDigit>);
291  std::unique_ptr< std::vector<simb::MCParticle > >
292  muonInfo(new std::vector< simb::MCParticle>);
293  std::unique_ptr< art::Assns<simb::MCParticle, rb::Cluster> >
295 
296  std::unique_ptr< std::vector<sim::Particle > >
297  mrccPartList(new std::vector< sim::Particle>);
298  std::unique_ptr< std::vector<sim::FLSHitList> >
299  mrccFLSHits(new std::vector<sim::FLSHitList>);
300  std::unique_ptr< std::vector<sim::PhotonSignal> >
301  mrccPhotSig(new std::vector<sim::PhotonSignal>);
302 
303  std::vector< rawdata::RawDigit> newDigits;
304  std::vector< rb::WeightedHit> allTrackHits;
306 
308  evt.getByLabel(fRawDataLabel, digitcol);
309  int nRaw = digitcol->size();
310  std::list< rawdata::RawDigit >digs;
311  for( int iRaw = 0; iRaw < nRaw; ++iRaw ){
312  art::Ptr< rawdata::RawDigit > iDig( digitcol, iRaw);
313  digs.push_back( *(iDig));
314  }// end loop over raw digits
315 
317  std::list< rawdata::RawDigit >::iterator digIter, digIterBegin, digIterEnd;
318 
319  // Get the slices in the event
321  evt.getByLabel(fClusterModuleLabel, slicecol);
323  for(unsigned int i = 0; i < slicecol->size(); ++i){
324  bool filt = true;
325  art::Ptr< rb::Cluster> sli(slicecol, i);
326  // Use only slices that passed preselection procedure and rockpresel filter
327  if(rb::IsFiltered(evt, slicecol, i, "rockpresel")) filt = false;
328  bool sel = SliceSelection(evt, sli);
329  if (sel && filt) slices.push_back(sli);
330  }
331 
332  rb::Track muonTrack;
333 
334  /// Specify lable which will be used for getting all tracks
335  std::string CallSpecificModule = "";
336  if (fWhichTrackModule == "kalman") CallSpecificModule = fKalmanTrackModuleLabel;
337  if (fWhichTrackModule == "bpf") CallSpecificModule = fBPFTrackModuleLabel;
338 
339  art::FindManyP<rb::Track> fmt(slices, evt, CallSpecificModule);
340 
341  int nSlices = slices.size();
342 
343  // Find the muon track in each slice. Longest track is assumed
344  // to be from the Muon for now.
345  for( int iSli = 0; iSli < nSlices; iSli++ ){
346 
347  // Find all the tracks in kth slice
348  std::vector< art::Ptr<rb::Track> > alltracks = fmt.at(iSli);
349  std::vector< art::Ptr<rb::Track> > tracks;
350  /// In case of BPF tracks there a re three tracks for each physical track
351  /// with different fit assumptions, loop tracks, check the fit sum and choose the muon ones
352  if (fWhichTrackModule == "bpf") {
353  art::FindManyP<rb::FitSum> fmBPFFitSums(alltracks, evt, fBPFTrackModuleLabel);
354  for (unsigned int trackId = 0; trackId < alltracks.size(); ++trackId){
355  if (fmBPFFitSums.isValid()){
356  std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(trackId);
357  if ( abs(bpfFitSums[0]->PDG()) == 13 ) tracks.push_back(alltracks[trackId]);
358  }
359  }
360  }
361  /// No special procedures for Kalman tracks
362  if (fWhichTrackModule == "kalman") tracks = alltracks;
363 
364  std::vector< rb::WeightedHit> trackHits;
365 
366  // Select the muon track based on the configuration in fhicl file.
367  //Muon energy estimation to be extract from Energy objects.
368  double MuEnergy = -5;
369 
370  if(fMuonID == "removeByLength")
371  muonTrack = LongestTrack( tracks, evt, MuEnergy );
372 
373  if(fMuonID == "removeByEfficiency")
374  muonTrack = RemoveByEfficiency( tracks, *slices[iSli], evt, MuEnergy );
375 
377  if(fMuonID == "removeByTruth"){
378  if( bt->HaveTruthInfo() )
379  trackHits = RemoveByTruth( slices[iSli]);
380  else
381  std::cerr<<" No MC truth info in this file.\n"
382  <<" Choose a different method for muon removal.\n";
383  }
384 
385  /// ReMId works only for Kalman tracks
386  if ( fMuonID == "removeByReMId" && fWhichTrackModule == "kalman" )
387  muonTrack = RemoveByReMId( tracks, evt, MuEnergy );
388 
389  // Works only for BPF tracks
390  if ( fMuonID == "removeByProngCVN" && fWhichTrackModule == "bpf" ){
391  muonTrack = RemoveByProngCVN( tracks, slices[iSli], evt, MuEnergy );
392  }
393 
394  if(fMuonID != "removeByTruth" && muonTrack.NCell() > 0){
395  muonInfo->push_back( MuonInfo(muonTrack, MuEnergy) );
396  util::CreateAssn(*this, evt, *muonInfo,
397  slices[iSli], *assns);
398  }
399 
400  if( muonTrack.NCell() != 0 && fMuonID != "removeByTruth" )
401  trackHits = fTrkCleanUpAlg->CleanUpTrack(muonTrack, *slices[iSli]);
402 
403  allTrackHits.insert( allTrackHits.end(), trackHits.begin(), trackHits.end() );
404 
405  }// End loop over slices
406 
407  SortByPlaneAndCell( allTrackHits );
408 
409  std::vector< rb::WeightedHit >::iterator trkIterStart = allTrackHits.begin(), trkIter;
410  for( auto digIter = digs.begin(); digIter != digs.end(); ++digIter){
411 
412  rawdata::RawDigit dig = *digIter;
413  rb::CellHit cHit(cal->MakeCellHit( &dig ));
414  int digPlane = cHit.Plane();
415  int digCell = cHit.Cell();
416  int digTDC = cHit.TDC();
417 
418  bool matched = false;
419 
420  for( trkIter = trkIterStart; trkIter != allTrackHits.end(); ++trkIter ){
421 
422  art::Ptr<rb::CellHit> hit = trkIter->hit;
423 
424  if( digPlane < hit->Plane() ){
425  trkIterStart = trkIter;
426  break;
427  }
428  else if( digPlane == hit->Plane() &&
429  digCell < hit->Cell() ){
430  trkIterStart = trkIter;
431  break;
432  }
433 
434  if(digPlane == hit->Plane() &&
435  digCell == hit->Cell() &&
436  digTDC == hit->TDC() ){
437 
438  float muonWeight = trkIter->weight;
439  float newADC = (1-muonWeight)*cHit.ADC();
440 
441  if( muonWeight == 1 || newADC < fADCThreshold ){
442  trkIterStart = trkIter;
443  matched = true;
444  break;
445  }
446  else if(muonWeight > 0 && muonWeight < 1){
447  // Need the following steps to be able to handle
448  // data/mc from before the multi-point readout era
449 
450  if( dig.IsMC() && dig.Version() == 0 && dig.NADC() == 3){
451  dig.SetADC(0, newADC );
452  dig.SetADC(1, (1-muonWeight)*cHit.ADC(1));
453  dig.SetADC(2, (1-muonWeight)*cHit.ADC(2));
454  }
455  else if( dig.Version() == 0 )
456  dig.SetADC(0, newADC );
457  else{
458 
460  unsigned int nPretrig = conv.getNPretriggeredSamples( dig.Version());
461  int16_t baseline = dig.ADC(nPretrig - fDCSDistance );
462  int nadc = dig.NADC();
463  for( int iadc = 0; iadc < nadc; ++iadc){
464  int16_t newadc = ((dig.ADC(iadc) - baseline )*( 1-muonWeight) ) + baseline;
465  dig.SetADC( iadc, newadc);
466  }
467  }
468  matched = true;
469  muonRemovedDigits->push_back( dig );
470  trkIterStart = trkIter;
471  break;
472  }// end if 0 < muonweight < 1 OC
473  }
474  }// end loop over all track hits
475  if( !matched)
476  muonRemovedDigits->push_back( *digIter );
477  }// end loop over all digits
478 
479 
480  std::cout<<"raw digits before and after removal: "
481  <<nRaw<<" "<<muonRemovedDigits->size()<<std::endl;
482 
483 
484  if( fMuonID == "removeByTruth"){
485  FillTruthInfo(evt, mrccFLSHits, mrccPhotSig, mrccPartList, muonInfo);
486 
487 
488  evt.put(std::move(mrccPartList));
489  evt.put(std::move(mrccPhotSig));
490  evt.put(std::move(mrccFLSHits));
491 
492  }
493 
494  evt.put(std::move(muonRemovedDigits));
495  evt.put(std::move(muonInfo));
496  evt.put(std::move(assns));
497 
498  return;
499 
500  }// end of producer
501 
502 
503 
504  /**********************************************************************/
505  ////////////////////////////////////////////////////////////////////////
506  /// Preslection procedure
508  {
509 
510  /// Get slice CVN PID
511  /// NB it can cause the bias in the further CVNe/presel ratio
512  /*art::FindManyP<cvn::Result> fmCVN ({slice}, evt, fCVNLabel);
513 
514  cvn::Result cvnResult;
515  float numuid = -1;
516  float nueid = -1;
517  float nutauid = -1;
518  float ncid = -1;
519  if (fmCVN.isValid()){
520  std::vector<art::Ptr<cvn::Result>> prods = fmCVN.at(0);
521  if (!prods.empty()){
522  cvnResult = *prods[0];
523  numuid = cvnResult.fOutput[0];
524  nueid = cvnResult.fOutput[1];
525  nutauid = -5;
526  ncid = cvnResult.fOutput[2];
527  }
528  }*/
529 
530  /// Loose quality
531  /// Vertex exists
532  bool vtxValid = false;
533  /// Number of prongs > 0
534  bool NProngValid = false;
535  /// at least one track with E > 0
536  double muonEmax = -1;
537 
538  art::FindManyP<rb::Vertex> fmElastic ({slice}, evt, fElasticInput);
539 
540  if (fmElastic.isValid()){
541  std::vector<art::Ptr<rb::Vertex>> elastics = fmElastic.at(0);
542  if(elastics.size() > 0) {
543  vtxValid = true; /// vertex exists
544  /// N Prongs > 0
545  art::FindManyP<rb::Prong> fmFuzzyProng3D ( elastics, evt, fFuzzy3DInput);
546  std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(0);
547  if (prongs.size() > 0) NProngValid = true; /// N prongs > 0 is valid
548  /// Check parent muon energy > 0
549  /// Check all tracks in this slice
550  std::string CallSpecificModule = "";
551  if (fWhichTrackModule == "kalman") CallSpecificModule = fKalmanTrackModuleLabel;
552  if (fWhichTrackModule == "bpf") CallSpecificModule = fBPFTrackModuleLabel;
553  /// Get all tracks
554  art::FindManyP<rb::Track> fmt({slice}, evt, CallSpecificModule);
555  std::vector< art::Ptr<rb::Track> > tracks = fmt.at(0);
556  for (unsigned int i = 0; i<tracks.size(); i++){
557  /// Use them as they are for Kalman tracks
558  if (fWhichTrackModule == "kalman") {
559  art::FindOneP<rb::Energy> muE({tracks[i]}, evt, fEnergyLabel);
560  double muonE = muE.at(0)->E();
561  if (muonE > muonEmax) muonEmax = muonE;
562  }
563  if (fWhichTrackModule == "bpf") {
564  /// Check BPF track's fit sum pdg (use only muons)
565  art::FindManyP<rb::FitSum> fmBPFFitSums({tracks[i]}, evt, fBPFTrackModuleLabel);
566  if (fmBPFFitSums.isValid()) {
567  std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(0);
568  if (abs(bpfFitSums[0]->PDG()) == 13) {
569  double muonE = bpfFitSums[0]->FourMom().E();
570  if (muonE > muonEmax) muonEmax = muonE;
571  }
572  }
573  }
574  } // end of loop over the tracks
575  } // end of elastics.size() > 0
576  } /// end of vtx is valid
577  /// resulting decision - all these things should exist and CVNm > 0.5 (?)
578  bool ret = vtxValid && NProngValid && (muonEmax > 0); //&& (numuid > 0.5);
579 
580  return ret;
581 
582  }
583 
584  /**********************************************************************/
585  ////////////////////////////////////////////////////////////////////////
586 
587  /// Returns the longest track amongst the given vector of tracks.
588 
589  rb::Track MuonRemove::LongestTrack( std::vector< art::Ptr<rb::Track> >& tracks, art::Event& evt, double& muonEnergy)
590  {
591  rb::Track longestTrack;
592  double maxLength = 0;
593  double length;
594  for(unsigned int i = 0; i < tracks.size(); i++){
595  length = tracks[i]->TotalLength();
596  if( length > maxLength){
597  maxLength = length;
598  longestTrack = *tracks[i];
599  std::string SpecificEnergyLabel = "";
600  /// check the muon energy
601  if (fWhichTrackModule == "kalman") {
602  SpecificEnergyLabel = fEnergyLabel;
603  art::FindOneP<rb::Energy> muE({tracks[i]}, evt, SpecificEnergyLabel);
604  muonEnergy = muE.at(0)->E();
605  }
606  if (fWhichTrackModule == "bpf") {
607  art::FindManyP<rb::FitSum> fmBPFFitSums({tracks[i]}, evt, fBPFTrackModuleLabel);
608  if (fmBPFFitSums.isValid()) {
609  std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(0);
610  muonEnergy = bpfFitSums[0]->FourMom().E();
611  }
612  }
613  } // end of if
614  } // end of loop
615 
616  return longestTrack;
617  }
618 
619 
620  /**********************************************************************/
621  ////////////////////////////////////////////////////////////////////////
622 
623  /// Finds muon hits by using truth
624  /// Loops through all the hits in a given slice and checks if
625  /// the photon signal was produced by a muon primary. If all
626  /// the light in a hit is not produced by the muon alone, then
627  /// only a fraction of the energy in the hit is attributed to the muon.
628 
629  std::vector<rb::WeightedHit> MuonRemove::RemoveByTruth( const art::Ptr<rb::Cluster>& slice)
630  {
631 
632  std::vector<rb::WeightedHit> muonHits;
634 
635  const sim::ParticleNavigator & pnav = bt->ParticleNavigator();
636  int nPrimaries = pnav.NumberOfPrimaries();
637  std::vector< int> primMuTrackIDs;
638 
639  for(int p = 0; p < nPrimaries; ++p){
640  const sim::Particle* prim = pnav.Primary(p);
641  if( abs(prim->PdgCode()) == 13 ){
642 
643  int muMotherID = prim->TrackId();
644  primMuTrackIDs.push_back( muMotherID);
645 
646  for( int i = 0; i < (int)pnav.size(); i++ ){
647 
648  const sim::Particle* part = pnav.Particle(i);
649  if(part->Mother() == muMotherID){
650  primMuTrackIDs.push_back( part->TrackId() );
651  }
652  }//end loop over particles
653  }
654  }// end loop over primaries
655 
656  for(int s = 0; s< (int)slice->NCell(); s++){
657 
658  if( !(bt->IsNoise( slice->Cell(s) ) ) ){
659 
660  const std::vector< sim::PhotonSignal> pSignal = bt->HitToPhotonSignal( slice->Cell(s));
661 
662  int muonPhotons = 0;
663  int totPhotons = 0;
664 
665  for( int p = 0; p < (int)pSignal.size(); p++){
666  totPhotons += pSignal[p].NPhoton();
667  int partID = pSignal[p].TrackId();
668 
669  // Find the photon signal produced by the muon
670  for(int i = 0; i < (int)primMuTrackIDs.size(); i++ ){
671  if( partID == primMuTrackIDs[i] ){
672  muonPhotons += pSignal[p].NPhoton();
673  }
674  }
675  }
676 
677  if(muonPhotons > 0){
678  if( muonPhotons == totPhotons ){
679  rb::WeightedHit newHit(slice->Cell(s), 1);
680  muonHits.push_back(newHit);
681  }
682  else {
683  rb::WeightedHit newHit(slice->Cell(s), (float)muonPhotons/totPhotons);
684  muonHits.push_back(newHit);
685  }
686  }
687  }// endif is noise
688  }// End loop over slice hits.
689 
690  return(muonHits);
691  }// End of RemoveByTruth
692 
693 
694  /**********************************************************************/
695  ////////////////////////////////////////////////////////////////////////
696 
697  /// Finds the highest efficiency muon track amongst a
698  /// vector of tracks. The efficiency is calculated based
699  /// on all hits in the slice that the track belongs to,
700  /// not on all hits in the event. BackTracker is used
701  /// to calculate efficiencies. Since there may be multiple
702  /// numuCC interactions in a slice, the process is twofold-
703  /// First find the highest efficiency track for every primary
704  /// muon, then select the track with highest effciency.
705 
706 
708  const rb::Cluster &slice, art::Event& evt, double& muonEnergy)
709  {
710 
711  rb::Track muonTrack;
713 
714  art::PtrVector< rb::CellHit > allHits = slice.AllCells();
715 
716  const sim::ParticleNavigator & pnav = bt->ParticleNavigator();
717  int nPrimaries = pnav.NumberOfPrimaries();
718  std::vector< std::set <int> > primMuTrackIDs;
719 
720  for(int n = 0; n < nPrimaries; ++n){
721  std::set <int> thisMuIDs;
722  const sim::Particle* prim = pnav.Primary(n);
723 
724  if( abs( prim->PdgCode() ) == 13 ){
725 
726  int muMotherID = prim->TrackId();
727  thisMuIDs.insert(muMotherID);
728 
729  // Find the track IDs of primary muons and its daughters
730  for(int p = 0; p < (int)pnav.size(); ++p){
731  const sim::Particle* part = pnav.Particle(p);
732  //Check if the primary mother is a Muon
733  if (part->Mother() == muMotherID )
734  thisMuIDs.insert( part->TrackId() );
735 
736  }// end loop over particles
737  }
738  primMuTrackIDs.push_back(thisMuIDs);
739  }// end loop over primaries
740 
741  double maxEff = -1;
742  double eff;
743 
744  for( int p = 0; p< (int)primMuTrackIDs.size(); p++){
745 
746  // The following loop finds the highest efficiency track for every primary muon.
747 
748  for( int i = 0; i< (int)tracks.size(); i++){
749 
750  art::PtrVector< rb::CellHit > trackHits = tracks[i]->AllCells();
751 
752  eff = bt->HitCollectionEfficiency( primMuTrackIDs[p], trackHits,
753  allHits, geo::kXorY);
754 
755  if( eff > maxEff){
756  maxEff = eff;
757  muonTrack = *tracks[i];
758  art::FindOneP<rb::Energy> muE({tracks[i]}, evt, fEnergyLabel);
759  muonEnergy = muE.at(0)->E();
760  }
761  }// end loop over tracks
762  }// end loop over primary muons
763 
764  return muonTrack;
765  } // End of RemoveByEfficiency
766 
767  /**********************************************************************/
768  ////////////////////////////////////////////////////////////////////////
769 
770 
771  /// A function to save the removed muon information in a simb::MCTrajectory
772  /// object. Note that only the first trajectory poin is filled put.
773 
774  simb::MCParticle MuonRemove::MuonInfo(rb::Track muonTrack, double muonEnergy)
775  {
776 
777  art::PtrVector< rb::CellHit> muonHits = muonTrack.AllCells();
778  rb::SortByPlane(muonHits);
779 
780  static TDatabasePDG* pdgt = TDatabasePDG::Instance();
781  TParticlePDG* pdgpm = pdgt->GetParticle(13);
782  assert(pdgpm);
783  const double muonMass = pdgpm->Mass();
784 
785 
786  // trackID, pdgCode, process, motherTrackID, status 1 = stableFinalState
787 
788  double muonMom = TMath::Sqrt((muonEnergy*muonEnergy) - (muonMass*muonMass));
789 
790  TLorentzVector startPosition( muonTrack.Start().X(), muonTrack.Start().Y(),
791  muonTrack.Start().Z(), muonTrack.MinTNS() );
792 
793  TLorentzVector startMomentum( muonTrack.Dir().X() * muonMom,
794  muonTrack.Dir().Y() * muonMom ,
795  muonTrack.Dir().Z() * muonMom, muonEnergy );
796 
797  TLorentzVector stopPosition( muonTrack.Stop().X(),
798  muonTrack.Stop().Y(),
799  muonTrack.Stop().Z(),
800  muonTrack.Cell( muonTrack.NCell() -1 )->TNS() );
801  TLorentzVector stopMomentum( 0,0,0,0 );
802 
803  simb::MCParticle muonPart( -1000, 13, "primary", -1000, muonMass, 1);
804  muonPart.AddTrajectoryPoint( startPosition, startMomentum );
805  muonPart.AddTrajectoryPoint( stopPosition, stopMomentum );
806  return muonPart;
807  }
808 
809 
810  /**********************************************************************/
811  ////////////////////////////////////////////////////////////////////////
812 
813 
814  /// A function to write muon removed truth information to the file.
815  /// It fills the unique_ptrs for FLSHitList, PhotonSignal, Particles
816  /// Also saves the removed muon's information in a vector of simb::MCTrajectory
817 
818 
820  std::unique_ptr< std::vector<sim::FLSHitList> >& mrccFLSHits,
821  std::unique_ptr< std::vector<sim::PhotonSignal> >& mrccPhotSig,
822  std::unique_ptr< std::vector<sim::Particle> >& mrccPartList,
823  std::unique_ptr< std::vector<simb::MCParticle> >& muonInfo)
824  {
825 
827 
829  evt.getByLabel(fMCTruthModuleLabel, MCTruthList);
830  if (MCTruthList->empty())
831  return;
832 
833  const sim::ParticleNavigator & pnav = bt->ParticleNavigator();
834  std::vector< int> primMuTrackIDs;
835 
836  // Get the TrackId's of Primary Muon and its daughters
837  for(int p = 0; p < (int)pnav.NumberOfPrimaries(); ++p){
838  const sim::Particle* prim = pnav.Primary(p);
839 
840  if( abs(prim->PdgCode()) == 13){
841 
842  muonInfo->push_back( *prim);
843  int muMotherID = prim->TrackId();
844  primMuTrackIDs.push_back(muMotherID);
845 
846  for( int i = 0; i < (int)pnav.size(); i++){
847  const sim::Particle* part = pnav.Particle(i);
848  if( part->Mother() == muMotherID ){
849  primMuTrackIDs.push_back( part->TrackId() );
850  }
851  }// end loop over particles
852  }
853  }// end loop over primaries
854 
855  // Make FLSHitLists for FLSHits not originating from Primary muon and its daughters.
856 
858  evt.getByLabel("geantgen", flslistHandle);
859 
860 
861  for( size_t l = 0; l < flslistHandle->size(); ++l){
862  sim::FLSHitList newFLSList;
863 
864  const std::vector<sim::FLSHit> &hits = flslistHandle->at(l).fHits;
865 
866  for( size_t f = 0; f < hits.size(); ++f){
867  bool hitIsFromMuon = false;
868 
869  for( int id = 0; id< (int)primMuTrackIDs.size(); id++){
870  if(hits[f].GetTrackID() == primMuTrackIDs[id])
871  hitIsFromMuon = true;
872  }
873  if(!hitIsFromMuon)
874  newFLSList.fHits.push_back(hits[f]);
875  }// end loop over flshits in this flshitlist
876 
877 
878  mrccFLSHits->push_back(newFLSList);
879  }// end loop over flshit lists
880 
881 
882  // Make vector of photon signals not originating from primary muon and its daughters.
883 
885  evt.getByLabel("photrans", phots);
886 
887  for(size_t i = 0; i < phots->size(); ++i){
888  const sim::PhotonSignal& phot = (*phots)[i];
889  bool sigIsFromMuon = false;
890 
891  for(int id = 0; id<(int)primMuTrackIDs.size(); id++){
892  if( phot.TrackId() == primMuTrackIDs[id])
893  sigIsFromMuon = true;
894  }
895 
896  if(!sigIsFromMuon)
897  mrccPhotSig->push_back(phot);
898  } // end of loop over Photon Signals
899 
900  // Make particle lists for all things besides muon and its daughters.
901  // Also, the muonInfo vector get s filled at this point.
902 
903  for (unsigned int s = 0; s < MCTruthList->size(); s++){
904 
905  art::Ptr<simb::MCTruth> mc (MCTruthList, s);
906  const art::Ptr<simb::MCTruth> neutrinoInteraction = mc;
907 
908  std::vector<const sim::Particle*> particleList = bt->MCTruthToParticles(neutrinoInteraction);
909 
910  for (unsigned int i = 0; i < particleList.size(); ++i){
911  bool partIsFromMuon = false;
912  for(int id = 0; id<(int)primMuTrackIDs.size(); id++){
913  if(particleList[i]->TrackId() == primMuTrackIDs[id])
914  partIsFromMuon = true;
915  }
916  if(!partIsFromMuon)
917  mrccPartList->push_back(*particleList[i]);
918  }
919  }// end loop over MCTruth
920 
921  }// End of FillTruthInfo
922 
923  /**********************************************************************/
924  ////////////////////////////////////////////////////////////////////////
925 
927  art::Event &evt, double& muonEnergy)
928  {
929  float maxPidVal = -1;
930  rb::Track mostMuLikeTrack;
932  if(fo.isValid()){
933  for(int i = 0; i < (int)tracks.size(); i++){
934  // Find the ReMId object associated with this track
935  art::Ptr<remid::ReMId> pid = fo.at(i);
936  float pidVal = pid->Value();
937  if( pidVal > maxPidVal ){
938  maxPidVal = pidVal;
939  mostMuLikeTrack = *tracks[i];
940  art::FindOneP<rb::Energy> muE({tracks[i]}, evt, fEnergyLabel);
941  muonEnergy = muE.at(0)->E();
942  }
943  }//end loop over tracks
944  }
945  else{
946  std::cerr<<" There are no remid objects in this file.\n"
947  <<" Choose another method of muon removal.\n";
948  }
949 
950  if( maxPidVal > fNumuPIDCut )
951  return mostMuLikeTrack;
952  else{
953  rb::Track voidTrack;
954  return voidTrack;
955  }
956  }// End of ReMId
957 
958  /**********************************************************************/
959  ////////////////////////////////////////////////////////////////////////
960 
962  art::Event &evt, double& muonEnergy)
963  {
964  /// Prong CVN works for prongs only, so need to find the most muon like prong first
965  /// and then use respective track for the mrcc purposes
966  float maxPidVal = -1;
967  art::Ptr<rb::Prong> mostMuLikeProng;
968  rb::Track mostMuLikeTrack;
969  /// Get the slice's vrtx
970  art::FindManyP<rb::Vertex> fmElastic ({slice}, evt, fElasticInput);
971  std::vector<art::Ptr<rb::Vertex>> elastics = fmElastic.at(0);
972  if( elastics.size() > 0 ){
973  /// Get all prongs for this vtx
974  art::FindManyP<rb::Prong> fmFuzzyProng3D ( elastics, evt, fFuzzy3DInput);
975  std::vector< art::Ptr<rb::Prong> > prongs = fmFuzzyProng3D.at(0);
976  /// Get the collection of CVNs for these prongs
977  art::FindManyP<rb::PID> fmCVNParticleResult(prongs, evt, fCVNParticleLabel);
978  for ( unsigned int iPng = 0; iPng < prongs.size(); ++iPng ) {
979  art::Ptr<rb::Prong> prong = prongs[iPng];
980  std::vector<art::Ptr<rb::PID>> cvnparts = fmCVNParticleResult.at(iPng);
981  for ( unsigned int j=0; j<cvnparts.size(); ++j ){
982  if ( abs(cvnparts[j]->Pdg()) == 13 ){
983  double thisPid = cvnparts[j]->Value();
984  if ( thisPid > maxPidVal) {
985  maxPidVal = thisPid;
986  mostMuLikeProng = prongs[iPng];
987  }
988  }
989  }
990  /// Prong CVN is empty for the prongs with L > 500 cm, assume these are definitely muons
991  if ( cvnparts.size() == 0 && fmCVNParticleResult.isValid() && prong->TotalLength() > 500){
992  maxPidVal = 1; ///?
993  mostMuLikeProng = prongs[iPng];
994  }
995  } // end of loop over the prongs
996  }
997  /// Done with searhing the most muon like prong
998  /// Find respective tracks (again x3 more tracks in BPF since there are three fit assumptions)
999  art::FindManyP<rb::Track> fmMuLikeTracks ( {mostMuLikeProng}, evt, fBPFTrackModuleLabel);
1000  std::vector< art::Ptr<rb::Track> > MuLikeTracks = fmMuLikeTracks.at(0);
1001  art::FindManyP<rb::FitSum> fmBPFFitSums ( MuLikeTracks, evt, fBPFTrackModuleLabel );
1002  /// Loop over all BPF tracks, look for the muon one
1003  for (unsigned int trackId = 0; trackId < MuLikeTracks.size(); ++trackId){
1004  if (fmBPFFitSums.isValid()){
1005  std::vector<art::Ptr<rb::FitSum>> bpfFitSums = fmBPFFitSums.at(trackId);
1006  if ( abs(bpfFitSums[0]->PDG()) == 13 ) {
1007  muonEnergy = bpfFitSums[0]->FourMom().E();
1008  mostMuLikeTrack=*MuLikeTracks[trackId];
1009  }
1010  }
1011  }
1012  if ( maxPidVal > fNumuPIDCut )
1013  return mostMuLikeTrack;
1014  else{
1015  rb::Track voidTrack;
1016  return voidTrack;
1017  }
1018 
1019  } // End remove by prong CVN
1020 
1021  /**********************************************************************/
1022  ////////////////////////////////////////////////////////////////////////
1023 
1024  // Compares two raw digits by plane. Same as CompareByPlane
1025  // function for CellHits in CellHit.h/.cxx. Helper to
1026  // SortDigitsByPlane.
1029  {
1031  if( cmap->GetPlane( &ra ) < cmap->GetPlane( &rb ) )
1032  return true;
1033  if( cmap->GetPlane( &ra ) > cmap->GetPlane( &rb ) )
1034  return false;
1035  return cmap->GetCell( &ra ) < cmap->GetCell( &rb );
1036 
1037  }// End of CompareByPlane
1038 
1039  /**********************************************************************/
1040  ////////////////////////////////////////////////////////////////////////
1041 
1042  // Sorts a vector of raw digits in order of increasing plane
1043  // number. Similar to SortByPlane function for CellHits.
1044  void SortDigitsByPlaneAndCell( std::list< rawdata::RawDigit >& digs)
1045  {
1046  digs.sort( CompareDigitsByPlaneAndCell );
1047  }// End of SortDigitsByPlane
1048 
1049 
1050  /**********************************************************************/
1051  ////////////////////////////////////////////////////////////////////////
1052 
1054  rb::WeightedHit whitb)
1055  {
1056  if( whita.hit->Plane() < whitb.hit->Plane() )
1057  return true;
1058  if( whita.hit->Plane() > whitb.hit->Plane() )
1059  return false;
1060  return whita.hit->Cell() < whitb.hit->Cell();
1061  }// End of CompareByPlane
1062 
1063  /**********************************************************************/
1064  ////////////////////////////////////////////////////////////////////////
1065 
1066  void SortByPlaneAndCell( std::vector< rb::WeightedHit >& whits)
1067  {
1068  std::sort( whits.begin(), whits.end(), CompareByPlaneAndCell );
1069  }// End of SortDigitsByPlane
1070 
1071 
1073 
1074 } // End namespace murem
1075 
1076 /**********************************************************************/
1077 ////////////////////////////////////////////////////////////////////////
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
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.
int PdgCode() const
Definition: MCParticle.h:211
bool fTrackCleanUp
A FHiCL flag to turn of track clean up. Default if true.
std::vector< sim::FLSHit > fHits
Definition: FLSHitList.h:21
std::vector< rb::WeightedHit > RemoveByTruth(const art::Ptr< rb::Cluster > &slice)
Returns a vector of Cellhits attributed to the muon by checking with truth. For implementing muon rem...
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void produce(art::Event &evt)
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:256
int32_t TDC() const
The time of the last baseline sample.
Definition: RawDigit.h:94
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
std::string fFuzzy3DInput
Label of module that produced prongs.
int TrackId() const
Definition: PhotonSignal.h:33
X or Y views.
Definition: PlaneGeo.h:30
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
size_type size() const
unsigned short Plane() const
Definition: CellHit.h:39
int Mother() const
Definition: MCParticle.h:212
const char * p
Definition: xmltok.h:285
std::string fCVNLabel
Label of module that produced slice CVN PID.
A collection of associated CellHits.
Definition: Cluster.h:47
rb::Track RemoveByProngCVN(std::vector< art::Ptr< rb::Track > > &tracks, const art::Ptr< rb::Cluster > &slice, art::Event &evt, double &muonEnergy)
Returns the BPF muon track with the highest Prong CVN value.
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
std::string fKalmanTrackModuleLabel
Label of module that produced tracks.
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
std::string fEnergyLabelBPF
Label of module that produced Energy object (in case of BPF tracks)
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
std::string fCVNParticleLabel
Label of module that produced Prong CVN PIDs.
double fNumuPIDCut
Numu PID cut.
rb::Track LongestTrack(std::vector< art::Ptr< rb::Track > > &tracks, art::Event &evt, double &muonEnergy)
A function that returns the longest track amongst a vector of tracks.
rb::CellHit MakeCellHit(const rawdata::RawDigit *rawdigit)
double Value() const
Definition: PID.h:22
void FillTruthInfo(const art::Event &evt, std::unique_ptr< std::vector< sim::FLSHitList > > &mrccFLSHits, std::unique_ptr< std::vector< sim::PhotonSignal > > &mrccPhotSig, std::unique_ptr< std::vector< sim::Particle > > &mrccPartList, std::unique_ptr< std::vector< simb::MCParticle > > &muonInfo)
Saves muon removed FLSHitLists, PhotonSignals, Particle list and information about the removed muon i...
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
Definition: Run.h:31
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
int TrackId() const
Definition: MCParticle.h:209
bool SliceSelection(const art::Event &evt, const art::Ptr< rb::Cluster > &slice)
Function for preselecting the slices for the further mrcc process.
fhicl::ParameterSet fTrkCleanUpPSet
ParameterSet of TrackCleanUpAlg.
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
const XML_Char * s
Definition: expat.h:262
double MinTNS() const
Definition: Cluster.cxx:482
void SetADC(uint32_t i, int16_t iADC)
Definition: RawDigit.cxx:108
TString part[npart]
Definition: Style.C:32
virtual double TotalLength() const
Distance along prong to reach last cell hit.
Definition: Prong.cxx:186
Far Detector at Ash River, MN.
void hits()
Definition: readHits.C:15
int fDCSDistance
DCS Distance: multipoint readout look back.
std::string fClusterModuleLabel
Label of module that produced slice clusters.
length
Definition: demo0.py:21
Result for CVN.
const sim::Particle * Primary(const int) const
void beginRun(art::Run &run)
uint32_t getNPretriggeredSamples(const version_t ver) const
Get number of pretriggered samples.
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
bool CompareDigitsByPlaneAndCell(rawdata::RawDigit ra, rawdata::RawDigit rb)
Compare two rawdigits by plane. Return true if plane of ra is greater then the plane of rb...
std::string fBPFTrackModuleLabel
Label of module that produced BPF tracks.
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
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
An algorithm to determine the energy of a slice.
Definition: NumuEAlg.h:28
Near Detector in the NuMI cavern.
bool IsMC() const
Definition: RawDigit.h:108
unsigned short GetPlane(const rawdata::RawDigit *dig)
Definition: CMap.cxx:285
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
float fADCThreshold_Far
ADC thresholds.
const double j
Definition: BetheBloch.cxx:29
void SortDigitsByPlaneAndCell(std::list< rawdata::RawDigit > &digitcol)
Takes in a vector of raw digits and sorts it by plane.
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
uint8_t Version() const
Definition: RawDigit.h:87
rb::Track RemoveByEfficiency(std::vector< art::Ptr< rb::Track > > &tracks, const rb::Cluster &slice, art::Event &evt, double &muonEnergy)
Returns the highest efficiency (by truth) muon track for removal;.
Definition: run.py:1
size_type size() const
Definition: PtrVector.h:308
OStream cout
Definition: OStream.cxx:6
std::vector< rb::WeightedHit > CleanUpTrack(const rb::Track &muonTrack, const rb::Cluster &slice, double *mipRangeHigh=NULL, double *mipRangeLow=NULL, double *mipValue=NULL, double *vertexRegionDeDxCutOff=NULL)
A function to disentangle the muon track from the hadronic energy close to the vertex. Returns a vector of rb::WeightedHits where the hit is the cellhits that were on the given track and weight is the fraction of the energy in the hit that came from the track. mipRangeLow, mipRangeHigh, mipValue and vertexRegionDeDxCutOff are pointers that point to the value of mip range limits, peak values andthe dedx cut off used to find the vertex region. If these pointers are not provided, the default values are used. The default values are good for muons, not for other particles.
bool CompareByPlaneAndCell(rb::WeightedHit whita, rb::WeightedHit whitb)
Compare two rawdigits by plane. Return true if plane of ra is greater then the plane of rb...
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const std::vector< sim::PhotonSignal > HitToPhotonSignal(const art::Ptr< rb::CellHit > &hit) const
Returns the PhotonSignals contributing the signal in the specified hit.
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
std::string fEnergyLabel
Label of module that produced Energy object (in case of Kalman tracks)
numue::NumuEAlg * fNumuEAlg
An object of NumuEAlg.
rb::Track RemoveByReMId(std::vector< art::Ptr< rb::Track > > &tracks, art::Event &evt, double &muonEnergy)
Returns the track with the highest PID value returned by the ReMId package.
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
uint32_t NADC() const
Definition: RawDigit.h:81
A vector of FLSHit from single neutrino interaction.
Definition: FLSHitList.h:13
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
MuonRemove(fhicl::ParameterSet const &pset)
fhicl::ParameterSet fNumuEAlgPSet
ParameterSet of NumuEAlg.
art::Ptr< rb::CellHit > hit
Definition: WeightedHit.h:23
Definition: event.h:1
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
std::string fMCTruthModuleLabel
Label of module that produced MC truth information.
Simple hit+weight pair, returned from rb::Cluster::WeightedHits.
Definition: WeightedHit.h:18
std::string fMuonID
Label of module that produced muon pid object.
assert(nhit_max >=nhit_nbins)
simb::MCParticle MuonInfo(rb::Track muonTrack, double muonEnergy)
Saves the information about the removed muon in a simb::MCParticle object The energy of the muon is o...
std::string fElasticInput
Label of module that produced verteces.
unsigned short GetCell(const rawdata::RawDigit *dig)
Definition: CMap.cxx:327
particleList
Definition: demo.py:41
std::string fNumuReMIdLabel
Label of module that produced Numu PID.
void reconfigure(const fhicl::ParameterSet &pset)
void SortByPlaneAndCell(std::vector< rb::WeightedHit > &whits)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
TString filt[ntarg]
Definition: Style.C:28
TrackCleanUpAlg * fTrkCleanUpAlg
An object of TrackCleanUpAlg.
Encapsulate the geometry of one entire detector (near, far, ndos)
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
Definition: CellHit.cxx:124
std::string fWhichTrackModule
Label of module which will be used for track loading.
std::string fRawDataLabel
Label of module that produced raw digits.
art::ServiceHandle< geo::Geometry > fgeom
Handle to Geometry.
mapped_type Particle(const size_type) const
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const