TrueEnergy_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief TrueEnergy module
3 /// A module to make sim::TrueEnergy's for each sim::Particle
4 /// in the module. This will ultimately be removed when we have
5 /// the sim::TrueEnergy's produced by G4Gen_module, for Prod4/5.
6 ///
7 /// \author karlwarb@iastate.edu
8 ///
9 ////////////////////////////////////////////////////////////////////////
10 
11 #include <string>
12 #include <vector>
13 
14 // NOvASoft includes
15 #include "MCCheater/BackTracker.h"
16 #include "Simulation/Particle.h"
17 #include "Simulation/TrueEnergy.h"
18 #include "Utilities/AssociationUtil.h"
19 
20 // Framework includes
22 #include "fhiclcpp/ParameterSet.h"
29 
30 namespace cheat {
31  class TrueEnergy : public art::EDProducer {
32  public:
33  explicit TrueEnergy(fhicl::ParameterSet const& pset);
34  virtual ~TrueEnergy();
35 
36  void produce(art::Event& evt);
37 
38  void reconfigure(fhicl::ParameterSet const& pset);
39 
40  private:
41 
42  std::string fGeantLabel; ///< label for module running G4 and making particles, etc
43  bool fHeavyDebug; ///< Boolean flag for whether you want heavy debugging info.
44  int fNumWarnings;///< The number of warnings I have sent to the screen.
45  };
46 }
47 
48 namespace cheat{
49 
50  //--------------------------------------------------------------------
52  {
53  this->reconfigure(pset);
54 
55  produces< std::vector<sim::TrueEnergy> >();
56  produces< art::Assns<sim::TrueEnergy, sim::Particle> >();
57 
58  // Set the number of warning I've sent to the screen to 0.
59  fNumWarnings = 0;
60  }
61 
62  //--------------------------------------------------------------------
64  {
65  }
66 
67  //--------------------------------------------------------------------
69  {
70  fGeantLabel = pset.get< std::string >("G4ModuleLabel");
71  fHeavyDebug = pset.get< bool >("HeavyDebug");
72  return;
73  }
74 
75  //--------------------------------------------------------------------
77  {
78  // --- Delcare my new vector, and its association.
79  std::unique_ptr<std::vector<sim::TrueEnergy> > trueens (new std::vector<sim::TrueEnergy>);
80  std::unique_ptr<art::Assns<sim::TrueEnergy, sim::Particle> > treassn (new art::Assns<sim::TrueEnergy, sim::Particle>);
81 
82  // --- Make a handle to BackTracker
84 
85  // --- Get my list of sim::Particle's out of the event
87  // --- Make a vector of the sim::Particles.
88  std::vector< art::Ptr<sim::Particle> > particles;
89  if ( evt.getByLabel(fGeantLabel, PartHand) ) {
90  art::fill_ptr_vector(particles, PartHand);
91  } else {
92  if (!evt.isRealData() && fNumWarnings < 10) {
93  std::cout << "TrueEnergy_module::I can't find any sim::Paritcle made by " << fGeantLabel << " in the event, fix the label!" << std::endl;
94  ++fNumWarnings;
95  }
96  evt.put(std::move(trueens));
97  evt.put(std::move(treassn));
98  return;
99  }
100 
101  for ( auto ThisPart : particles ) {
102 
103  if (fHeavyDebug)
104  std::cout << "\n\n**********************\nLooking at particle with TrackID " << ThisPart->TrackId() << ", it is a " << ThisPart->PdgCode() << ".\n" << std::endl;
105  // --- Get the information I need for this particular particle.
106  int This_Ent, This_Exit;
107  bool EntDetector = bt->InterceptsDetector( *ThisPart, This_Ent, This_Exit );
108  float This_EnterEn = 0;
109  float This_EscapEn = 0;
110  if ( EntDetector ) {
111  This_EnterEn = ThisPart->E( This_Ent ) - ThisPart->Mass();
112  if ( This_Exit != (int)(ThisPart->NumberTrajectoryPoints()-1) )
113  This_EscapEn = ThisPart->E( This_Exit ) - ThisPart->Mass();
114  }
115 
116  //std::cout << "\n ********** Worked out the deps from this particle, now work out daughter energies. **********\n" << std::endl;
117  // --- Now work out the deps including the particles daughers.
118  float TotalEnDep = 0.;
119  float TotalEscEn = bt->CalcTotalEscapingEnergy( *ThisPart, TotalEnDep, true );
120 
121  // --- Write out what I got from all of that....
122  if (fHeavyDebug) {
123  std::cout << "\n********** Finished, now what did I get? **********\n"
124  << "For TrackID " << ThisPart->TrackId() << ", PDG " << ThisPart->PdgCode() << ", Mother " << ThisPart->Mother() << ", had " << ThisPart->NumberDaughters() << " daughters.\n"
125  << "\t Using just this particle I got: EnteringEn " << This_EnterEn << ", EscapingEn " << This_EscapEn << ".\n"
126  << "\t Including the daughters I got : TotalDepEn was " << TotalEnDep << ", and the TotalEscEn was " << TotalEscEn << ".\n"
127  << std::endl;
128  }
129  // --- Now, make my association!!
130  // Make my sim::TrueEnergy
131  sim::TrueEnergy ThisEn( ThisPart->TrackId(), This_EnterEn, This_EscapEn, TotalEnDep, TotalEscEn );
132  (*trueens).push_back( ThisEn );
133  // Make the association.
134  util::CreateAssn(*this, evt, *trueens, ThisPart, *treassn);
135  } // Loop over sim::Particles
136 
137  evt.put(std::move(trueens));
138  evt.put(std::move(treassn));
139 
140  return;
141 
142  } // end produce
143 
144 } // end namespace
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.
bool InterceptsDetector(const sim::Particle &Par, int &TrajEnt, int &TrajEx)
A function which returns the trajectory point at which the particle first enters and leave the detect...
int fNumWarnings
The number of warnings I have sent to the screen.
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
std::string fGeantLabel
label for module running G4 and making particles, etc
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
void reconfigure(fhicl::ParameterSet const &pset)
OStream cout
Definition: OStream.cxx:6
TrueEnergy(fhicl::ParameterSet const &pset)
code to link reconstructed objects back to the MC truth information
Definition: FillTruth.h:17
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
void produce(art::Event &evt)
bool fHeavyDebug
Boolean flag for whether you want heavy debugging info.
float CalcTotalEscapingEnergy(const sim::Particle &Par, bool UseMCPart=true, std::string G4Label="trueens")
An extension of CalcEscapingEnergy which also takes into account the energy which the daughters of gi...
enum BeamMode string