NumuEnergy_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: NumuEnergy
3 // Module Type: producer
4 // File: NumuEnergy_module.cc
5 //
6 // Generated at Fri Sep 21 08:53:53 2012 by Susan Lein using artmod
7 // from art v1_00_11.
8 ////////////////////////////////////////////////////////////////////////
9 
10 #include <string>
11 #include "TMath.h"
12 
13 // Framework includes
19 
20 
21 // NOvA includes
22 #include "NumuEnergy/NumuEAlg.h"
23 #include "NumuEnergy/NumuE.h"
24 #include "Geometry/Geometry.h"
26 #include "RecoBase/Cluster.h"
27 #include "RecoBase/Track.h"
28 #include "RecoBase/Energy.h"
29 #include "Utilities/AssociationUtil.h"
30 #include "Utilities/func/MathUtil.h"
31 #include "NovaDAQConventions/DAQConventions.h"
32 #include "MCCheater/BackTracker.h"
33 
35 #include "SummaryData/SpillData.h"
36 
37 
38 namespace numue {
39 
40  class NumuEnergy : public art::EDProducer {
41  public:
42  explicit NumuEnergy(fhicl::ParameterSet const & pset);
43  virtual ~NumuEnergy();
44 
45  void produce(art::Event & e);
46  void beginRun(art::Run& run);
47 
48  private:
49 
50  fhicl::ParameterSet fTrackCleanUpAlgPSet; ///< Parameter Set to configure the TrackCleanUpAlg object
51  murem::TrackCleanUpAlg* fTrackCleanUpAlg; ///< Track Clean Up Algorithm Object
52  std::string fHitModuleLabel; ///< Label for module making hits
53  std::string fSlicerModuleLabel; ///< Label for module doing slicing
54  std::string fTrackModuleLabel; ///< Label for module doing tracking
55  std::string fCosRejModuleLabel; ///< Label for module doing cosmic rejection
56  std::string fGeneratorLabel; ///< Generator label for the spillPot
57  std::string fNuMIBeamLabel; ///< NuMI Beam label for the spillPot
58  fhicl::ParameterSet fNumuEAlgPSet; ///< Parameter Set to configure the NumuEAlg object
59  NumuEAlg* fNumuEAlg; ///< Numu Energy Algorithm Object
60  bool isRHC;
62  };
63 
64  //----------------------------------------------------------------------
66  fTrackCleanUpAlgPSet(pset.get< fhicl::ParameterSet >("TrackCleanUpAlgPSet") ),
68  fHitModuleLabel(pset.get< std::string >("HitModuleLabel")),
69  fSlicerModuleLabel(pset.get< std::string >("SlicerModuleLabel")),
70  fTrackModuleLabel(pset.get< std::string >("TrackModuleLabel")),
71  fCosRejModuleLabel(pset.get< std::string >("CosRejModuleLabel")),
72  fGeneratorLabel (pset.get<std::string>("GeneratorLabel")),
73  fNuMIBeamLabel (pset.get<std::string>("NuMIBeamLabel")),
74  fNumuEAlgPSet(pset.get< fhicl::ParameterSet >("NumuEAlgPSet")),
75  fNumuEAlg(0)
76  {
77  produces< std::vector<NumuE> >();
78  produces< art::Assns<NumuE, rb::Cluster> >();
79 
80  produces< std::vector<rb::Energy> >();
81  produces< art::Assns<rb::Energy, rb::Track> >();
82  fLoadHornCurrent = false;
83  }
84 
85 
86  //----------------------------------------------------------------------
88  {
89  if (fNumuEAlg) delete fNumuEAlg;
91  }
92 
93  //----------------------------------------------------------------------
95  {
96  // Have to recreate these because behaviour can depend on run number
97  if(fNumuEAlg) delete fNumuEAlg;
101  }//End of Begin Run
102 
103  //----------------------------------------------------------------------
105  {
107  novadaq::cnv::DetId detID = geom->DetId();
108 
109  //Get all the hits from the event
111  e.getByLabel(fHitModuleLabel, allHits);
112 
114 
115  //Load the slicer list from the event
117  e.getByLabel(fSlicerModuleLabel,slicecol);
118 
119  if(slicecol->empty()){
120  mf::LogWarning ("No Slices")<<"No Slices in the input file";
121  return;
122  }
123 
124  //Getting tracks associated with slices
125  art::FindManyP<rb::Track> sliceToTracks(slicecol, e, fTrackModuleLabel);
126 
127  //Declaring containers for things to be stored in event
128  std::unique_ptr< std::vector<NumuE> > energyCol(new std::vector<NumuE>);
129  std::unique_ptr< art::Assns<NumuE, rb::Cluster> > assoc(new art::Assns<NumuE, rb::Cluster>);
130 
131  std::unique_ptr< std::vector<rb::Energy> > trackECol(new std::vector<rb::Energy>);
132  std::unique_ptr< art::Assns<rb::Energy, rb::Track> > trackAssoc(new art::Assns<rb::Energy, rb::Track>);
133 
135 
136  if (!fLoadHornCurrent){
138  if (!e.isRealData()) {
139  e.getByLabel(fGeneratorLabel, spillPot);
140  }
141  else {
142  e.getByLabel(fNuMIBeamLabel, spillPot);
143  }
144 
145  if (spillPot.failedToGet()) {
146  mf::LogError("NumuEAlg") << "Spill Data not found, "
147  "aborting without horn current information";
148  abort();
149  }
150  isRHC = spillPot->isRHC;
151  mf::LogInfo("NumuEAlg") << "Setting Horn Current isRHC to: "<<isRHC;
152  fLoadHornCurrent = true;
153  }
154 
155  for(size_t s = 0; s < slicecol->size(); ++s){
156  const art::Ptr<rb::Cluster> slice(slicecol, s);
157 
158  if(slice->IsNoise()) continue;
159 
160  art::PtrVector<rb::CellHit> sliceHits = slice->AllCells();
161 
162  // Get all of the tracks in the slice
163  const std::vector< art::Ptr<rb::Track> > sliceTracks = sliceToTracks.at(s);
164  if(sliceTracks.empty()) continue;
165 
166  NumuE sliceEnergy;
167  rb::Energy angleqeE;
168  double errAngleQEE = -1.0;
169  NumuE truthEnergy;
170  std::vector<rb::Energy> trackEnergies;
171 
173  bool ismc = bt->HaveTruthInfo();
174 
175 
176 
177  trackEnergies = fNumuEAlg->MuonEnergies(sliceTracks,
178  e.run(),
179  ismc,
180  isRHC);
181  sliceEnergy = fNumuEAlg->Energy(sliceTracks,sliceHits,slice,e,fTrackCleanUpAlg,isRHC);
182  angleqeE = fNumuEAlg->QEFormulaEnergy(sliceTracks,sliceHits,errAngleQEE,e,isRHC);
183  sliceEnergy.SetAngleQEE(angleqeE.E());
184  sliceEnergy.SetAngleQEError(errAngleQEE);
185  truthEnergy = fNumuEAlg->MCTruthEnergyVariables(sliceTracks,slice,*allHits,e);
186  sliceEnergy.SetMCTrueMuonE(truthEnergy.MCTrueMuonE());
187  sliceEnergy.SetMCTrueMuonCatcherE(truthEnergy.MCTrueMuonCatcherE());
188  sliceEnergy.SetMCGoodTrueMuon(truthEnergy.MCGoodTrueMuon());
189 
190  if(nCosRej.isValid() && nCosRej.at(s).size() > 0) {
191  sliceEnergy.SetUCCCE(fNumuEAlg->GetUCE(nCosRej,s,slice->TotalGeV(),slice->NCell()));
192  if (detID == novadaq::cnv::kFARDET){
193  sliceEnergy.SetUCMuonESingle(fNumuEAlg->GetUCMuonESingle(sliceTracks,sliceHits,e,nCosRej,s));
194  sliceEnergy.SetUCMuonENonSingle(fNumuEAlg->GetUCMuonENonSingle(sliceTracks,sliceHits,e,nCosRej,s));
195  }
196  else{
197  sliceEnergy.SetUCMuonESingle(-5.0);
198  sliceEnergy.SetUCMuonENonSingle(-5.0);
199  }// Only for FarDet so far
200  }// End CosRej object conditional
201 
202  if((sliceEnergy.E() > 0) || (sliceEnergy.CalCCE() > 0)){
203  energyCol->push_back(sliceEnergy);
204  util::CreateAssn(*this,e,*(energyCol.get()),slice,*(assoc.get()));
205  }//End of loop over non-zero energy
206  for(size_t c = 0; c < trackEnergies.size(); ++c){
207  trackECol->push_back(trackEnergies[c]);
208  util::CreateAssn(*this,e,*(trackECol.get()),sliceTracks[c],*(trackAssoc.get()));
209  }//End of loop over track energies for the slice
210 
211  }//End of loop over slices
212 
213  //Put objects and associations in file
214  e.put(std::move(energyCol));
215  e.put(std::move(assoc));
216  e.put(std::move(trackECol));
217  e.put(std::move(trackAssoc));
218 
219  return;
220  }//End of Produce
221 
222  //----------------------------------------------------------------------
224 
225 }// End of namespace
bool isRHC
is the beam in antineutrino mode, aka RHC
Definition: SpillData.h:28
fhicl::ParameterSet fNumuEAlgPSet
Parameter Set to configure the NumuEAlg object.
double GetUCMuonENonSingle(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Event const &e, art::FindManyP< cosrej::CosRejObj > CRO, size_t slicenum)
Get values from cosrej object and use them to get TMVA trained value of uncontained muon energy for F...
Definition: NumuEAlg.cxx:1404
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.
std::string fTrackModuleLabel
Label for module doing tracking.
std::string fGeneratorLabel
Generator label for the spillPot.
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
NumuEAlg * fNumuEAlg
Numu Energy Algorithm Object.
Energy estimators for CC events.
Definition: FillEnergies.h:7
float E() const
Definition: Energy.cxx:27
void beginRun(art::Run &run)
void SetAngleQEE(float angleqe)
Definition: NumuE.cxx:102
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void SetMCGoodTrueMuon(bool mcgoodtruemuon)
Definition: NumuE.cxx:229
NumuE Energy(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Ptr< rb::Cluster > sliceCluster, art::Event const &e, murem::TrackCleanUpAlg *trkCleanUpAlg, bool isRHC)
Returns various energy estimations for the current detector. See comments at top of function in ...
Definition: NumuEAlg.cxx:639
void SetMCTrueMuonE(float mctruemuone)
Definition: NumuE.cxx:219
murem::TrackCleanUpAlg * fTrackCleanUpAlg
Track Clean Up Algorithm Object.
void SetAngleQEError(float angleqeerror)
Definition: NumuE.cxx:107
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
Definition: Run.h:31
std::string fSlicerModuleLabel
Label for module doing slicing.
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
const XML_Char * s
Definition: expat.h:262
void SetMCTrueMuonCatcherE(float mctruemuoncatchere)
Definition: NumuE.cxx:224
Far Detector at Ash River, MN.
NumuEnergy(fhicl::ParameterSet const &pset)
void SetUCCCE(float uccc)
Definition: NumuE.cxx:132
double GetUCMuonESingle(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, art::Event const &e, art::FindManyP< cosrej::CosRejObj > CRO, size_t slicenum)
Get values from cosrej object and use them to get TMVA trained value of uncontained muon energy for F...
Definition: NumuEAlg.cxx:1365
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
An algorithm to determine the energy of a slice.
Definition: NumuEAlg.h:28
std::string fCosRejModuleLabel
Label for module doing cosmic rejection.
rb::Energy QEFormulaEnergy(std::vector< art::Ptr< rb::Track > > const sliceTracks, art::PtrVector< rb::CellHit > sliceHits, double &error, art::Event const &e, bool isRHC)
Returns QE energy estimation using formula. See comments at top of function in .cxx for full explanat...
Definition: NumuEAlg.cxx:324
NumuE MCTruthEnergyVariables(const std::vector< art::Ptr< rb::Track >> &sliceTracks, const art::Ptr< rb::Cluster > &slice, const std::vector< rb::CellHit > &allHits, const art::Event &e)
Returns various useful truth energy values, good for fitting.
Definition: NumuEAlg.cxx:655
fhicl::ParameterSet fTrackCleanUpAlgPSet
Parameter Set to configure the TrackCleanUpAlg object.
std::string fNuMIBeamLabel
NuMI Beam label for the spillPot.
Definition: run.py:1
void produce(art::Event &e)
double GetUCE(art::FindManyP< cosrej::CosRejObj > CRO, size_t slicenum, float totalgev, int nhitslice)
Get values from cosrej object and use them to get TMVA trained value of uncontained energy for FD...
Definition: NumuEAlg.cxx:1350
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetUCMuonENonSingle(float ucmuonenonsingle)
Definition: NumuE.cxx:143
A container for energy information.
Definition: Energy.h:20
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
float CalCCE() const
Definition: NumuE.cxx:239
float MCTrueMuonE() const
Definition: NumuE.cxx:374
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
void SetUCMuonESingle(float ucmuonesingle)
Definition: NumuE.cxx:137
bool MCGoodTrueMuon() const
Definition: NumuE.cxx:384
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
float MCTrueMuonCatcherE() const
Definition: NumuE.cxx:379
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
std::vector< rb::Energy > MuonEnergies(std::vector< art::Ptr< rb::Track > > const sliceTracks, int run, bool ismc, bool isRHC) const
Function that returns energy for each track given the assumption that it is a muon. It uses the functions above, like MuonEFromTrackLength and NDMuonEFromTrackLength.
Definition: NumuEAlg.cxx:1137
std::string fHitModuleLabel
Label for module making hits.