EMPi0Energy_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: EMPi0Energy
3 // Module Type: analyzer
4 // File: EMPi0Energy_module.cc
5 // Author: Mike Wallbank (m.wallbank@sheffield.ac.uk), August 2015
6 //
7 // Analyser module to produce useful information for characterising
8 // pi0 showers.
9 ////////////////////////////////////////////////////////////////////////
10 
11 // Framework includes:
14 #include "fhiclcpp/ParameterSet.h"
19 #include "art_root_io/TFileService.h"
20 #include "art_root_io/TFileDirectory.h"
23 
24 // LArSoft includes
37 #include "nug4/ParticleNavigation/ParticleList.h"
39 
40 // ROOT & C++ includes
41 #include <string>
42 #include <vector>
43 #include <map>
44 #include "TTree.h"
45 #include "TBranch.h"
46 #include "TLeaf.h"
47 #include "TVector3.h"
48 
49 const int kMaxClusters = 100;
50 
51 namespace emshower {
52  class EMPi0Energy;
53 }
54 
56 public:
57 
58  EMPi0Energy(fhicl::ParameterSet const& pset);
59  void analyze(art::Event const& evt);
60  void reconfigure(fhicl::ParameterSet const& p);
61  void reset();
62  double ConvertChargeToEnergy(double charge, int plane);
63  double FindDepositedEnergy(int trackID);
64  int FindTrackID(detinfo::DetectorClocksData const& clockData,
65  art::Ptr<recob::Hit> const& hit);
66  int FindTrueTrack(detinfo::DetectorClocksData const& clockData,
67  std::vector<art::Ptr<recob::Hit> > const& clusterHits);
68  double FindVertexDetectorDistance(const simb::MCParticle* particle);
69 
70 private:
71 
72  // Variables for the tree
73  TTree* fTree;
74  double trueEnergyPi0;
84  int nclusters;
90 
97 
98  // For converting charge to energy
102 
103 };
104 
106  this->reconfigure(pset);
107  fTree = tfs->make<TTree>("EMPi0Energy","EMPi0Energy");
108  fTree->Branch("TrueEnergyPi0", &trueEnergyPi0);
109  fTree->Branch("TrueEnergyHighEPhoton", &trueEnergyHighEPhoton);
110  fTree->Branch("TrueEnergyLowEPhoton", &trueEnergyLowEPhoton);
111  fTree->Branch("TrueTrackIDPi0", &trueTrackIDPi0);
112  fTree->Branch("TrueTrackIDHighEPhoton", &trueTrackIDHighEPhoton);
113  fTree->Branch("TrueTrackIDLowEPhoton", &trueTrackIDLowEPhoton);
114  fTree->Branch("DepositHighEPhoton", &depositHighEPhoton);
115  fTree->Branch("DepositLowEPhoton", &depositLowEPhoton);
116  fTree->Branch("VertexDetectorDistHighEPhoton",&vertexDetectorDistHighEPhoton);
117  fTree->Branch("VertexDetectorDistLowEPhoton", &vertexDetectorDistLowEPhoton);
118  fTree->Branch("NClusters", &nclusters);
119  fTree->Branch("Cluster_Plane", cluster_plane, "cluster_plane[NClusters]/I");
120  fTree->Branch("Cluster_Size", cluster_size, "cluster_size[NClusters]/I");
121  fTree->Branch("Cluster_TrueTrackID", cluster_truetrackid,"cluster_truetrackid[NClusters]/I");
122  fTree->Branch("Cluster_Charge", cluster_charge, "cluster_charge[NClusters]/D");
123  fTree->Branch("Cluster_Energy", cluster_energy, "cluster_energy[NClusters]/D");
124 
125  Uintercept = -1519.33; Ugradient = 148867;
126  Vintercept = -1234.91; Vgradient = 149458;
127  Zintercept = -1089.73; Zgradient = 145372;
128 
129 }
130 
132  fHitsModuleLabel = pset.get<std::string>("HitsModuleLabel");
133  fClusterModuleLabel = pset.get<std::string>("ClusterModuleLabel");
134 }
135 
137 
138  /// Analyse function to save information for calibrating shower energies
139  /// This is written assuming single pi0 per event
140 
141  this->reset();
142 
143  // Get the hits out of the event record
144  std::vector<art::Ptr<recob::Hit> > hits;
145  auto hitHandle = evt.getHandle<std::vector<recob::Hit> >(fHitsModuleLabel);
146  if (hitHandle)
147  art::fill_ptr_vector(hits, hitHandle);
148 
149  // Get the clusters out of the event record
150  std::vector<art::Ptr<recob::Cluster> > clusters;
151  auto clusterHandle = evt.getHandle<std::vector<recob::Cluster> >(fClusterModuleLabel);
152  if (clusterHandle)
153  art::fill_ptr_vector(clusters, clusterHandle);
154 
155  art::FindManyP<recob::Hit> fmh(clusterHandle, evt, fClusterModuleLabel);
156 
157  // Look at the clusters
158  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
159  for (unsigned int clus = 0; clus < clusters.size(); ++clus) {
160 
161  art::Ptr<recob::Cluster> cluster = clusters.at(clus);
162  std::vector<art::Ptr<recob::Hit> > hits = fmh.at(clus);
163 
164  // Find the charge deposited by hits in this cluster in this plane
165  double charge = 0;
166  for (std::vector<art::Ptr<recob::Hit> >::iterator hit = hits.begin(); hit != hits.end(); ++hit)
167  charge += ((*hit)->Integral() * TMath::Exp((500 * (*hit)->PeakTime())/3e6));
168 
169  cluster_plane [clus] = cluster->Plane().Plane;
170  cluster_size [clus] = hits.size();
171  cluster_truetrackid[clus] = this->FindTrueTrack(clockData, hits);
172  cluster_charge [clus] = charge;
173  cluster_energy [clus] = this->ConvertChargeToEnergy(charge, cluster->Plane().Plane);
174 
175  }
176 
177  // Event level information
178 
179  nclusters = clusters.size();
180 
181  // Get the pi0 and the decay photons
182  const sim::ParticleList& trueParticles = particleinventory->ParticleList();
183  const simb::MCParticle* truePi0 = trueParticles.Primary(0);
184  if (truePi0->NumberDaughters() != 2) return;
185  const simb::MCParticle* truePhoton1 = particleinventory->TrackIdToParticle_P(truePi0->Daughter(0));
186  const simb::MCParticle* truePhoton2 = particleinventory->TrackIdToParticle_P(truePi0->Daughter(1));
187 
188  // Make sure photon 1 energy > photon 2 energy
189  if (truePhoton1->Momentum().E() < truePhoton2->Momentum().E()) {
190  const simb::MCParticle* tmp = truePhoton2;
191  truePhoton2 = truePhoton1;
192  truePhoton1 = tmp;
193  }
194 
195  trueEnergyPi0 = truePi0->Momentum().E();
196  trueTrackIDPi0 = truePi0->TrackId();
197 
198  trueEnergyHighEPhoton = truePhoton1->Momentum().E();
199  trueEnergyLowEPhoton = truePhoton2->Momentum().E();
200  trueTrackIDHighEPhoton = truePhoton1->TrackId();
201  trueTrackIDLowEPhoton = truePhoton2->TrackId();
202 
203  // Find the energy deposited on each plane in the TPC
206 
207  // Find the distance between the particle vertex and the edge of the detector
210 
211  fTree->Fill();
212 
213  return;
214 
215 }
216 
217 double emshower::EMPi0Energy::ConvertChargeToEnergy(double charge, int plane) {
218 
219  /// Converts charge to energy
220 
221  double energy = 0;
222 
223  switch (plane) {
224  case 0:
225  energy = (double)(charge - Uintercept)/(double)Ugradient;
226  break;
227  case 1:
228  energy = (double)(charge - Vintercept)/(double)Vgradient;
229  break;
230  case 2:
231  energy = (double)(charge - Zintercept)/(double)Zgradient;
232  break;
233  }
234 
235  return energy;
236 
237 }
238 
240 
241  std::vector<sim::IDE> ides;
242  for( auto ide_P : backtracker->TrackIdToSimIDEs_Ps(trackID) ){
243  ides.push_back(*ide_P);
244  }
245 
246  double deposit = 0;
247 
248  for (std::vector<sim::IDE>::iterator ideIt = ides.begin(); ideIt != ides.end(); ++ideIt)
249  deposit += ideIt->energy;
250 
251  // Put energies in GeV units
252  deposit /= 1000;
253 
254  return deposit;
255 
256 }
257 
259  art::Ptr<recob::Hit> const& hit) {
260 
261  /// Find the true track ID this hit is associated with
262 
263  double particleEnergy = 0;
264  int likelyTrackID = 0;
265  std::vector<sim::TrackIDE> trackIDs = backtracker->HitToTrackIDEs(clockData, hit);
266  for (unsigned int idIt = 0; idIt < trackIDs.size(); ++idIt) {
267  if (trackIDs.at(idIt).energy > particleEnergy) {
268  particleEnergy = trackIDs.at(idIt).energy;
269  likelyTrackID = TMath::Abs(trackIDs.at(idIt).trackID);
270  }
271  }
272  return likelyTrackID;
273 }
274 
276  std::vector<art::Ptr<recob::Hit> > const& clusterHits) {
277 
278  /// Find the true track which is most associated to this cluster
279 
280  std::map<int,double> trackMap;
281  for (std::vector<art::Ptr<recob::Hit> >::const_iterator clusHitIt = clusterHits.begin(); clusHitIt != clusterHits.end(); ++clusHitIt) {
282  art::Ptr<recob::Hit> hit = *clusHitIt;
283  int trackID = FindTrackID(clockData, hit);
284  trackMap[trackID] += hit->Integral();
285  }
286  //return std::max_element(trackMap.begin(), trackMap.end(), [](const std::pair<int,double>& p1, const std::pair<int,double>& p2) {return p1.second < p2.second;} )->first;
287  double highestCharge = 0;
288  int clusterTrack = 0;
289  for (std::map<int,double>::iterator trackIt = trackMap.begin(); trackIt != trackMap.end(); ++trackIt)
290  if (trackIt->second > highestCharge) {
291  highestCharge = trackIt->second;
292  clusterTrack = trackIt->first;
293  }
294  return clusterTrack;
295 }
296 
298 
299  /// Finds the rough distance from the end of a particle track to the edge of the detector, along the direction of the particle
300 
301  TVector3 end = TVector3(particle->EndX(),particle->EndY(),particle->EndZ());
302  TVector3 direction = TVector3(particle->Px(),particle->Py(),particle->Pz()).Unit();
303 
304  int distanceStep = 1, steps = 0;
305  TVector3 pos;
306  bool inTPC = true;
307 
308  while (inTPC) {
309  pos = end + ( (steps*distanceStep) * direction );
310  double currentPos[3]; currentPos[0] = pos.X(); currentPos[1] = pos.Y(); currentPos[2] = pos.Z();
311  if (!geom->FindTPCAtPosition(currentPos).isValid)
312  inTPC = false;
313  ++steps;
314  }
315 
316  double distance = (end - pos).Mag();
317 
318  return distance;
319 
320 }
321 
323  trueEnergyPi0 = 0;
326  trueTrackIDPi0 = 0;
329  depositHighEPhoton = 0;
330  depositLowEPhoton = 0;
333  nclusters = 0;
334  for (int cluster = 0; cluster < kMaxClusters; ++cluster) {
335  cluster_plane [cluster] = 0;
336  cluster_size [cluster] = 0;
338  cluster_charge [cluster] = 0;
339  cluster_energy [cluster] = 0;
340  }
341 }
342 
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
intermediate_table::iterator iterator
double cluster_charge[kMaxClusters]
double Py(const int i=0) const
Definition: MCParticle.h:231
void cluster(In first, In last, Out result, Pred *pred)
Definition: NNClusters.h:41
double EndZ() const
Definition: MCParticle.h:228
Encapsulate the construction of a single cyostat.
std::vector< sim::TrackIDE > HitToTrackIDEs(detinfo::DetectorClocksData const &clockData, recob::Hit const &hit) const
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
art::ServiceHandle< art::TFileService > tfs
std::string string
Definition: nybbler.cc:12
art::ServiceHandle< cheat::BackTrackerService > backtracker
const simb::MCParticle * TrackIdToParticle_P(int id) const
bool isValid
Whether this ID points to a valid element.
Definition: geo_types.h:211
double Px(const int i=0) const
Definition: MCParticle.h:230
struct vector vector
std::vector< const sim::IDE * > TrackIdToSimIDEs_Ps(int const &id) const
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
intermediate_table::const_iterator const_iterator
geo::PlaneID Plane() const
Returns the plane ID this cluster lies on.
Definition: Cluster.h:744
Cluster finding and building.
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
Particle class.
double EndY() const
Definition: MCParticle.h:227
int NumberDaughters() const
Definition: MCParticle.h:217
art framework interface to geometry description
EMPi0Energy(fhicl::ParameterSet const &pset)
int TrackId() const
Definition: MCParticle.h:210
int Daughter(const int i) const
Definition: MCParticle.cxx:112
double FindVertexDetectorDistance(const simb::MCParticle *particle)
art::ServiceHandle< geo::Geometry > geom
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
art::ServiceHandle< cheat::ParticleInventoryService > particleinventory
double ConvertChargeToEnergy(double charge, int plane)
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
double FindDepositedEnergy(int trackID)
int cluster_truetrackid[kMaxClusters]
T get(std::string const &key) const
Definition: ParameterSet.h:271
p
Definition: test.py:223
const int kMaxClusters
string tmp
Definition: languages.py:63
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
double cluster_energy[kMaxClusters]
Detector simulation of raw signals on wires.
const sim::ParticleList & ParticleList() const
int cluster_size[kMaxClusters]
void analyze(art::Event const &evt)
Declaration of signal hit object.
Encapsulate the construction of a single detector plane.
int FindTrackID(detinfo::DetectorClocksData const &clockData, art::Ptr< recob::Hit > const &hit)
Contains all timing reference information for the detector.
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
int cluster_plane[kMaxClusters]
double Pz(const int i=0) const
Definition: MCParticle.h:232
int FindTrueTrack(detinfo::DetectorClocksData const &clockData, std::vector< art::Ptr< recob::Hit > > const &clusterHits)
Tools and modules for checking out the basics of the Monte Carlo.
TCEvent evt
Definition: DataStructs.cxx:7
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:297
double EndX() const
Definition: MCParticle.h:226
void reconfigure(fhicl::ParameterSet const &p)
Encapsulate the construction of a single detector plane.