MRE_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \file MRE.cxx
3 /// \brief Creates an electron in place on the muon removed by MuonRemove job
4 ///
5 /// \version $Id: MRE.cxx,v 1.3 2012-10-01 23:24:12 bckhouse Exp $
6 /// \author ksachdev@physics.umn.edu
7 ////////////////////////////////////////////////////////////////////////
8 
9 
10 
11 // C/C++ includes
12 #include <cmath>
13 #include <iostream>
14 #include <map>
15 #include <vector>
16 
17 // ROOT includes
18 #include "TMath.h"
19 #include "TVector3.h"
20 #include "TLorentzVector.h"
21 #include "TDatabasePDG.h"
22 
23 // Framework includes
25 #include "fhiclcpp/ParameterSet.h"
29 #include "cetlib/search_path.h"
30 #include "Utilities/AssociationUtil.h"
33 
34 // NOvASoft includes
35 #include "g4nova/G4Alg.h"
38 #include "Simulation/FLSHitList.h"
39 #include "Simulation/Particle.h"
40 #include "MCCheater/BackTracker.h"
41 
42 
43 namespace art
44 {
45  class Event;
46  class ParameterSet;
47 }
48 
49 namespace murem
50 {
51  class MRE : public art::EDProducer
52  {
53  public:
54  explicit MRE(fhicl::ParameterSet const& pset);
55  ~MRE();
56 
57  void produce(art::Event& evt);
58 
59  void reconfigure(const fhicl::ParameterSet& pset);
60 
61  void beginRun(art::Run& run);
62 
63  private:
64 
65  g4n::G4Alg* fG4Alg; ///< G4Helper object
66  fhicl::ParameterSet fG4AlgPSet; ///< parameter set to configure the G4Alg object
67 
70 
71  };
72 }
73 
74 
75 ////////////////////////////////////////////////////////////////////////
76 
77 
78 namespace murem
79 {
80 
81  /*********************************************************************************************/
82  ////////////////////////////////////////////////////////////////////////
83 
84  MRE::MRE(fhicl::ParameterSet const& pset)
85  : fG4Alg(0)
86  , fG4AlgPSet(pset.get< fhicl::ParameterSet >("G4AlgPSet") )
87  {
88  this->reconfigure(pset);
89 
90  produces< std::vector<simb::MCTruth> >();
91 
92  // association between mre electron truth
93  // and removed muon mcparticle
94  produces< art::Assns<simb::MCTruth, simb::MCParticle> >();
95 
96  // association between mre electron truth
97  // and electron sim particles
98  produces< art::Assns< simb::MCTruth, sim::Particle> >();
99 
100  produces< std::vector<sim::FLSHitList> >();
101  produces< std::vector<sim::Particle> >();
102 
103 
104  mf::LogInfo("MRE") << " MRE::MRE()\n";
105  }
106 
107  /*********************************************************************************************/
108  ////////////////////////////////////////////////////////////////////////
109 
111  {
112  if(fG4Alg) delete fG4Alg;
113  }
114 
115  /*********************************************************************************************/
116  ////////////////////////////////////////////////////////////////////////
117 
119  {
120  fMrccInfoLabel = pset.get< std::string >("MrccInfoLabel");
121  fMCTruthLabel = pset.get< std::string >("MCTruthLabel");
122  }
123 
124  /*********************************************************************************************/
125  ////////////////////////////////////////////////////////////////////////
127  {
128 
130  return;
131  }
132 
134  {
135 
136  std::unique_ptr<std::vector<simb::MCTruth> >
137  mctcol (new std::vector<simb::MCTruth> );
138  std::unique_ptr<std::vector<sim::FLSHitList> >
139  flshlcol(new std::vector<sim::FLSHitList> );
140  std::unique_ptr<std::vector<sim::Particle> >
141  pcol (new std::vector<sim::Particle> );
142  std::unique_ptr< art::Assns<simb::MCTruth, simb::MCParticle> >
144  std::unique_ptr< art::Assns<simb::MCTruth, sim::Particle> >
146 
147 
149  evt.getByLabel(fMrccInfoLabel, muonDataCol);
150 
151  static TDatabasePDG* pdgt = TDatabasePDG::Instance();
152  TParticlePDG* pdgpm = pdgt->GetParticle(13);
153  assert(pdgpm);
154 
155  TParticlePDG* pdgpe = pdgt->GetParticle(11);
156  assert(pdgpe);
157  const double electronMass = pdgpe->Mass();
158 
159  // Set the trackID offset for the electron.
160  // If file is MC, electron trackID's will be offset by
161  // the highest pre-existing trackID in the event
162  // If data is real, trackID offset is just -1000.
163 
164  int trackidOffset = -1000;
165 
167  evt.getByLabel(fMCTruthLabel, pListHandle);
168 
169  if( pListHandle.failedToGet() )
170  trackidOffset = -1000;
171  else{
172  for( int p = 0; p < (int) pListHandle->size(); p++){
173  art::Ptr<sim::Particle> part(pListHandle, p);
174  int partId = part->TrackId();
175  if( partId > trackidOffset)
176  trackidOffset = partId+1;
177  }
178  }
179 
180  // Loop over all the removed muons
181  for( int mData = 0; mData < (int)muonDataCol->size(); mData++){
182 
183  art::Ptr<simb::MCParticle> muonData(muonDataCol, mData);
184 
185  simb::MCTruth truth;
187 
188  double muonE = muonData->E(0);
189  double electronE = muonE;
190  double electronMom = TMath::Sqrt((electronE*electronE) -
191  (electronMass*electronMass));
192  TVector3 muonDir = muonData->Momentum(0).Vect().Unit();
193 
194  // A nearly vertical track usually has unreasonable
195  // reconstructed length and energy. Do not introduce
196  // a Mr.E for it.
197  if(muonDir.Z() < 1e-6)
198  continue;
199  TLorentzVector electron4Momentum( electronMom*muonDir.X(),
200  electronMom*muonDir.Y(),
201  electronMom*muonDir.Z(),
202  electronE);
203  TLorentzVector origin(muonData->Vx(0), muonData->Vy(0),
204  muonData->Vz(0), muonData->T(0));
205 
206  int id = trackidOffset;
207  simb::MCParticle part(id, 11, "primary");
208  part.AddTrajectoryPoint(origin, electron4Momentum);
209  truth.Add(part);
210  mctcol->push_back(truth);
211 
212  util::CreateAssn(*this, evt, *mctcol, muonData, *assns);
213 
214  std::vector< sim::Particle > particles;
215  fG4Alg->RunGeant(&truth, *flshlcol, particles,
216  id+1);
217 
218  int first = pcol->size();
219  int last = pcol->size() + particles.size();
220 
221  for(int p = 0; p < (int) particles.size(); p++){
222  pcol->push_back(particles[p]);
223  int partId = particles[p].TrackId();
224  if( partId > trackidOffset)
225  trackidOffset = partId;
226  }
227 
228  util::CreateAssn(*this, evt, *mctcol,
229  *pcol, *assnspart, first, last);
230 
231 
232  }// end loop over muon data
233 
234 
235  evt.put(std::move(mctcol));
236  evt.put(std::move(flshlcol));
237  evt.put(std::move(pcol));
238  evt.put(std::move(assns));
239  evt.put(std::move(assnspart));
240  return;
241 
242  }// end of producer
243 
244 
246 
247 
248 }// end of namespace murem
249 /////////////////////////////////////////////////////////////////
double E(const int i=0) const
Definition: MCParticle.h:232
back track the reconstruction to the simulation
void RunGeant(std::vector< art::Handle< std::vector< simb::MCTruth > > > &mclists, std::vector< sim::FLSHitList > &flshitlist, std::vector< sim::Particle > &particlelist, std::vector< std::vector< std::pair< size_t, size_t > > > &pListLimits)
Definition: G4Alg.cxx:280
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.
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:256
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:80
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void produce(art::Event &evt)
Definition: MRE_module.cc:133
const char * p
Definition: xmltok.h:285
void beginRun(art::Run &run)
Definition: MRE_module.cc:126
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
void Add(simb::MCParticle &part)
Definition: MCTruth.h:79
std::string fMrccInfoLabel
Definition: MRE_module.cc:69
Definition: Run.h:31
fhicl::ParameterSet fG4AlgPSet
parameter set to configure the G4Alg object
Definition: MRE_module.cc:66
int TrackId() const
Definition: MCParticle.h:209
void reconfigure(const fhicl::ParameterSet &pset)
Definition: MRE_module.cc:118
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
TString part[npart]
Definition: Style.C:32
g4n::G4Alg * fG4Alg
G4Helper object.
Definition: MRE_module.cc:65
single particles thrown at the detector
Definition: MCTruth.h:26
T get(std::string const &key) const
Definition: ParameterSet.h:231
double T(const int i=0) const
Definition: MCParticle.h:223
Definition: run.py:1
double Vx(const int i=0) const
Definition: MCParticle.h:220
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
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
std::string fMCTruthLabel
Definition: MRE_module.cc:68
An algorithm to pass interaction information to Geant4 and create hits and particle lists...
Definition: G4Alg.h:32
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
double Vz(const int i=0) const
Definition: MCParticle.h:222
assert(nhit_max >=nhit_nbins)
Service to store calibration data products (CDP) in the SQLite3 metadatabase of a file...
Definition: FillParentInfo.h:8
Event generator information.
Definition: MCTruth.h:32
Float_t e
Definition: plot.C:35
Definition: fwd.h:28
double Vy(const int i=0) const
Definition: MCParticle.h:221
bool failedToGet() const
Definition: Handle.h:196