SupernovaMCClusters_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \brief
3 // \author Justin Vasel <justin.vasel@gmail.com>
4 // \date April 2020
5 ////////////////////////////////////////////////////////////////////////
6 
7 // Framework includes
8 // #include "art/Framework/Core/EDAnalyzer.h"
15 #include "fhiclcpp/ParameterSet.h"
16 
17 // ROOT includes
18 #include "TFile.h"
19 #include "TTree.h"
20 
21 // NOvA includes
22 #include "MCCheater/BackTracker.h"
23 #include "RawData/RawTrigger.h"
24 #include "RecoBase/CellHit.h"
25 #include "Simulation/Particle.h"
26 
27 
28 namespace sn {
29  class SupernovaMCCluster;
30  enum ProcessID {
31  kNone=-1,
32  kUnknown=0,
37  };
38 }
39 
41 public:
42  explicit SupernovaMCCluster(fhicl::ParameterSet const & p);
43 
44  SupernovaMCCluster(SupernovaMCCluster const &) = delete;
48 
50  void GetDaughters(sim::ParticleNavigator pNav, const sim::Particle* p, std::vector<sim::Particle>& dVec);
53 
54  void produce(art::Event& e);
55  void beginJob() override;
56 
57 private:
60 
61  unsigned int fEvt; ///< Event number
63 
66 
67  std::map<int, std::pair<int, int>> mapTrackIDInteraction;
68  // std::map<int, int> mapTrackIDParticleType;
69  std::map<std::pair<int, int>, std::vector<art::Ptr<rb::CellHit>>> mapInteractionToHits;
70 
74  float fClusterMeV;
79 
81 
82 };
83 
84 
85 //...........................................................................
87  // EDAnalyzer(p),
88  fRawDataLabel(p.get<std::string>("RawDataLabel")),
89  fCellHitLabel(p.get<std::string>("CellHitLabel")),
91  fInteractionId(0),
92  fClusterNHits(0),
95  fClusterMeV(0),
96  fClusterNViews(0),
97  fClusterMeanX(0),
98  fClusterMeanY(0),
99  fClusterMeanZ(0)
100 {
102  produces<std::vector<rb::Cluster>>();
103 }
104 
105 
106 //...........................................................................
108 {
109 
110  if (process == "Primary") return ProcessID::kPrimary;
111  if (process == "hadElastic") return ProcessID::kHadElastic;
112  if (process == "nCapture") return ProcessID::kNCapture;
113  if (process == "neutronInelastic") return ProcessID::kNeutronInelastic;
114 
115  std::cout << "Unknown process: " << process << std::endl;
116 
117  return ProcessID::kUnknown;
118 }
119 
120 // ..........................................................................
121 void sn::SupernovaMCCluster::GetDaughters(sim::ParticleNavigator pNav, const sim::Particle* p, std::vector<sim::Particle>& dVec)
122 {
123  for (int dId=0; dId<p->NumberDaughters(); ++dId) {
125  if (dIter == pNav.end()) continue;
126 
127  const sim::Particle* d = dIter->second;
128  this->GetDaughters(pNav, d, dVec);
129  dVec.push_back(*d);
130  }
131 
132  return;
133 }
134 
135 
136 //...........................................................................
138 {
139  //
140  // Interaction 1
141  // | Particle 1 (nuebar)
142  // | Particle 2 (proton)
143  // |_ Particle 3 (neutron)
144  // | |_ Daughter 1
145  // | | |_ Hit 1
146  // | | |_ Hit 2
147  // | |_ Daughter 2
148  // |_Particle 4 (positron)
149  // |_ Hit 1
150  //
151  // Interaction 2
152  // etc. etc.
153 
154  int nPrimaries = pNav.NumberOfPrimaries();
155  for (int idxPrimary=0; idxPrimary<nPrimaries; ++idxPrimary) {
156  // Check that this is an IBD event. The first primary particle should be a
157  // nuebar (-12), the second a proton (1221), the third a positron (-11),
158  // and the fourth, a neutron (2112).
159  bool kIsInverseBetaDecay = pNav.Primary(idxPrimary+0)->PdgCode() == -12 &&
160  pNav.Primary(idxPrimary+1)->PdgCode() == 2212 &&
161  pNav.Primary(idxPrimary+2)->PdgCode() == -11 &&
162  pNav.Primary(idxPrimary+3)->PdgCode() == 2112;
163 
164  // bool kIsElasticScattering = pNav.Primary(idxPrimary+0)->PdgCode() == -12 &&
165  // pNav.Primary(idxPrimary+1)->PdgCode() == -11;
166 
167 
168  if (kIsInverseBetaDecay) {
169  ++fInteractionId;
170  } else {
171  continue;
172  }
173 
174  int interactionSubId = -1;
175 
176  std::vector<sim::Particle> daughters_ebar;
177  std::vector<sim::Particle> daughters_n;
178 
179  // Get the positron primary & any daugthers
180  interactionSubId = 1;
181  const sim::Particle* positron = pNav.Primary(idxPrimary+2);
182  this->mapTrackIDInteraction[positron->TrackId()] = std::make_pair(fInteractionId, interactionSubId);
183 
184  this->GetDaughters(pNav, positron, daughters_ebar);
185  // std::cout << "Found " << daughters_ebar.size() << " positron daughters, subId = " << interactionSubId << std::endl;
186  for (const sim::Particle daughter : daughters_ebar ) {
187  this->mapTrackIDInteraction[daughter.TrackId()] = std::make_pair(fInteractionId, interactionSubId);
188  }
189 
190  // Get the neutron primary & any daugthers
191  interactionSubId = 2;
192  const sim::Particle* neutron = pNav.Primary(idxPrimary+3);
193  this->mapTrackIDInteraction[neutron->TrackId()] = std::make_pair(fInteractionId, interactionSubId);
194 
195 
196  this->GetDaughters(pNav, neutron, daughters_n);
197  // std::cout << "Found " << daughters_n.size() << " neutron daughters, subId = " << interactionSubId << std::endl;
198  for (const sim::Particle daughter : daughters_n ) {
199  this->mapTrackIDInteraction[daughter.TrackId()] = std::make_pair(fInteractionId, interactionSubId);
200  }
201 
202  idxPrimary += 3;
203  }
204 }
205 
206 //...........................................................................
208 {
209  for (art::Ptr<rb::CellHit> hit : hits) {
210  if (!hit->IsMC()) continue;
211 
212  // Get the particle that contributed the most to the hit
213  std::vector<cheat::TrackIDE> trackIDEs = bt->HitToTrackIDE(*hit);
214  cheat::TrackIDE trackIDE = trackIDEs.at(0);
215  const sim::Particle* particle = bt->TrackIDToParticle(trackIDE.trackID);
216 
217  // Get the interaction associated with this particle
218  int interactionId = this->mapTrackIDInteraction[particle->TrackId()].first;
219  int interactionSubId = this->mapTrackIDInteraction[particle->TrackId()].second;
220  if (interactionId == 0 || interactionSubId == 0) continue;
221 
222  // Add the hit to the list of hits associated with this interaction
223  std::vector<art::Ptr<rb::CellHit>> vecHits = this->mapInteractionToHits[std::make_pair(interactionId, interactionSubId)];
224  vecHits.push_back(hit);
225  this->mapInteractionToHits[std::make_pair(interactionId, interactionSubId)] = vecHits;
226  }
227 }
228 
229 
230 //...........................................................................
232 {
233  std::unique_ptr<std::vector<rb::Cluster>> product_clusters(new std::vector<rb::Cluster>);
234 
235  fEvt = e.id().event();
236 
239 
240  /* Get the cell hits from event */
242  e.getByLabel(fCellHitLabel, hitsHdl);
243  std::vector<art::Ptr<rb::CellHit>> hits;
244  art::fill_ptr_vector(hits, hitsHdl);
245 
246  // Build Interaction-TrackID and Interaction-Hits maps
247  mapTrackIDInteraction.clear();
248  mapInteractionToHits.clear();
249 
250  this->BuildInteractionTrackIDMap(pNav);
251  this->BuildInteractionHitMap(bt, hits);
252 
253 
254  for (auto it = this->mapInteractionToHits.begin(); it != this->mapInteractionToHits.end(); ++it) {
255  if ((it->second).size() == 0) continue;
256 
257  fInteractionId = it->first.first;
258  fInteractionSubId = it->first.second;
259 
260  rb::Cluster* intCluster = new rb::Cluster();
261  for (art::Ptr<rb::CellHit> h : it->second) {
262  intCluster->Add(h);
263  }
264 
265  if (intCluster->NCell()>0) {
266  product_clusters->push_back(*intCluster);
267  fClusterNHits = intCluster->NCell();
268  fClusterMeanTFile = intCluster->MeanTNS() * 1e-9 + fTimeSinceFileStart;
269  fClusterMeanTEvent = intCluster->MeanTNS() * 1e-9;
270  fClusterMeV = intCluster->TotalGeV() * 1e3;
271  fClusterNViews = (intCluster->NXCell()>0 && intCluster->NYCell()>0) ? 2 : 1;
272  fClusterMeanX = intCluster->MeanX();
273  fClusterMeanY = intCluster->MeanY();
274  fClusterMeanZ = intCluster->MeanZ();
275  fTreeParticles->Fill();
276  }
277  }
278 
279  std::cout << "Number of interactions: " << this->mapTrackIDInteraction.size() << std::endl;
280  std::cout << "Number of interactions w/ hits: " << this->mapInteractionToHits.size() << std::endl;
281  std::cout << "Number of clusters: " << product_clusters->size() << std::endl;
282 
283  e.put(std::move(product_clusters));
284 
285  /* Increment time */
287  e.getByLabel(fRawDataLabel, hdlTrigger);
288  const rawdata::RawTrigger rawTrigger = hdlTrigger->at(0);
289  fTimeSinceFileStart += (float)rawTrigger.fTriggerRange_TriggerLength * 500e-9;
290 
291  return;
292 }
293 
294 
295 //...........................................................................
297 {
298  /* Build TTree */
300 
301  fTreeParticles = tfs->make<TTree>("Clusters", "Clusters");
302  fTreeParticles->Branch("Event", &fEvt);
303  fTreeParticles->Branch("InteractionId", &fInteractionId);
304  fTreeParticles->Branch("InteractionSubId", &fInteractionSubId);
305  fTreeParticles->Branch("cluster_hits", &fClusterNHits);
306  fTreeParticles->Branch("cluster_meanT_file", &fClusterMeanTFile);
307  fTreeParticles->Branch("cluster_meanT_event", &fClusterMeanTEvent);
308  fTreeParticles->Branch("cluster_sumMeV", &fClusterMeV);
309  fTreeParticles->Branch("cluster_nviews", &fClusterNViews);
310  fTreeParticles->Branch("cluster_meanX", &fClusterMeanX);
311  fTreeParticles->Branch("cluster_meanY", &fClusterMeanY);
312  fTreeParticles->Branch("cluster_meanZ", &fClusterMeanZ);
313 
314 
315  return;
316 }
317 
318 
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
set< int >::iterator it
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
const char * p
Definition: xmltok.h:285
std::map< int, std::pair< int, int > > mapTrackIDInteraction
list_type::const_iterator const_iterator
SupernovaMCCluster(fhicl::ParameterSet const &p)
A collection of associated CellHits.
Definition: Cluster.h:47
double MeanZ(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:233
void GetDaughters(sim::ParticleNavigator pNav, const sim::Particle *p, std::vector< sim::Particle > &dVec)
std::map< std::pair< int, int >, std::vector< art::Ptr< rb::CellHit > > > mapInteractionToHits
void BuildInteractionTrackIDMap(sim::ParticleNavigator pNav)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
DEFINE_ART_MODULE(TestTMapFile)
Track type is none of the below.
int NumberDaughters() const
Definition: MCParticle.h:216
int TrackId() const
Definition: MCParticle.h:209
int Daughter(const int i) const
Definition: MCParticle.cxx:112
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
double MeanX(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:231
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
void hits()
Definition: readHits.C:15
int trackID
Geant4 supplied trackID.
Definition: BackTracker.h:40
double MeanY(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.h:232
const sim::Particle * Primary(const int) const
Float_t d
Definition: plot.C:236
void BuildInteractionHitMap(art::ServiceHandle< cheat::BackTracker > bt, std::vector< art::Ptr< rb::CellHit >> hits)
Remove hits from hot and cold channels.
SupernovaMCCluster & operator=(SupernovaMCCluster const &)=delete
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
OStream cout
Definition: OStream.cxx:6
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
EventNumber_t event() const
Definition: EventID.h:116
int ProcessIDFromString(std::string process)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Definition: structs.h:12
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
uint32_t fTriggerRange_TriggerLength
Definition: RawTrigger.h:40
unsigned int fEvt
Event number.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
Float_t e
Definition: plot.C:35
iterator find(const key_type &key)
EventID id() const
Definition: Event.h:56