ShowerAnaCheck_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \version 1
3 /// \brief
4 /// \author Dan Pershey - pershdawg@gmail.com
5 ////////////////////////////////////////////////////////////////////////
6 #include <cmath>
7 #include <iostream>
8 #include <string>
9 #include <vector>
10 
11 #include "TH1.h"
12 
13 #include "Geometry/Geometry.h"
15 #include "MCCheater/BackTracker.h"
16 #include "Simulation/FLSHit.h"
17 #include "Simulation/FLSHitList.h"
18 #include "Simulation/Particle.h"
20 #include "Utilities/func/MathUtil.h"
21 
23 // Framework includes
27 #include "fhiclcpp/ParameterSet.h"
35 
36 
37 namespace mcchk {
38  /// check the results from the Monte Carlo generator
40 
41  public:
42  explicit ShowerAnaCheck(fhicl::ParameterSet const& pset);
43  virtual ~ShowerAnaCheck();
44 
45  void analyze(const art::Event& evt);
46  void reconfigure(const fhicl::ParameterSet& pset);
47  void beginJob();
48  void EMShowerWidth(float& Event_distance, int MomsID, const sim::Particle* part, const std::vector<sim::FLSHit>& hits);
49  void EMMoliereRadius(float& Moliere_radius, int MomsID, const sim::Particle* part, const std::vector<sim::FLSHit>& hits);
50  float PointLineDistance(float x, float y, float z, float x1, float y1, float z1, float x2, float y2, float z2);
51 
52  private:
53 
55  std::vector<float> Half_Width_of_Detector;
56  std::vector<float> Half_Height_of_Detector;
57  std::vector<float> Begin_of_Detector_in_Z;
58  std::vector<float> End_of_Detector_in_Z;
59  std::vector<float> Energy_Percent;
60  std::vector<float> Min_Energy;
61 
62  double XYZ[3];
63  double XYZWorld[3];
64 
65 
67  TH1F* fVx;
68  TH1F* fVy;
69  TH1F* fVz;
70  TH1F* fTruthEnergy;
71  TH1F* fCosTheta;
74  TH1F* fEventWidth;
76 
77  };
78 }
79 
80 
81 ////////////////////////////////////////////////////////////////////////
82 namespace mcchk {
83 
85  : EDAnalyzer(pset)
86  {
87  this->reconfigure(pset);
88  }
89 
90  ////////////////////////////////////////////////////////////////////////////
91 
93 
94  ////////////////////////////////////////////////////////////////////////////
95 
97  {
98  fFLSHitListLabel = pset.get< std::string >("FLSHitListLabel");
99  Half_Width_of_Detector = pset.get<std::vector<float> >("sHalf_Width_of_Detector");
100  Half_Height_of_Detector = pset.get<std::vector<float> >("sHalf_Height_of_Detector");
101  Begin_of_Detector_in_Z = pset.get<std::vector<float> >("sBegin_of_Detector_in_Z");
102  End_of_Detector_in_Z = pset.get<std::vector<float> >("sEnd_of_Detector_in_Z");
103  Energy_Percent = pset.get<std::vector<float> >("sEnergy_Percent");
104  Min_Energy = pset.get<std::vector<float> >("sMin_Energy");
105  }
106 
108  {
110 
111  fShowerMotherPDG = tfs->make<TH1F>("fFirstPDG","Plot of PDG Code of Shower Mother; PDG Code; Entries",120,0,120);
112  fVx = tfs->make<TH1F>("fVx","X Coordinate of Shower Mother Vertex; X (cm); Entries",30,-150,150);
113  fVy = tfs->make<TH1F>("fVy","Y Coordinate of Shower Mother Vertex; Y (cm); Entries",40,-200,200);
114  fVz = tfs->make<TH1F>("fVz","Z Coordinate of Shower Mother Vertex; Z (cm); Entries",100,0,1000);
115  fTruthEnergy = tfs->make<TH1F>("fTruthEnergy","MC Truth Energy of Shower Mother; Energy (GeV); Entries",40,0,8);
116  fCosTheta = tfs->make<TH1F>("fCosTheta","Cosine of Angle Between the Shower Mother Momentum and Z Axis; radians;Entries",100,-1,1);
117  fNumberPerEvent = tfs->make<TH1F>("fNumberPerEvent","Number of Particles in Each Event; Number in Event; Entries",15,0,15);
118  fEnergyofShower = tfs->make<TH1F>("fEnergyofShower","Energy of All FLSHITS; Energy (GeV); Entries",100,0,12);
119  fEventWidth = tfs->make<TH1F>("fEventWidth","Average Distance to FLSHit per Event; Distance (cm); Entries",50,0,15);
120  fMoliereRadius = tfs->make<TH1F>("fMoliereRadius","Plot of the Moliere Radius of Each Event, 10-11 cm for NOvA; R_{Moliere}; Entries",25,0,20);
121 
122  }
123 
124  ///////////////////////////////////////////////////////////////////////////
125 
127  {
128 
129  mf::LogInfo("ShowerAnaCheck") << "Calculating Shower Parameters";
130 
131  // Get the FLS list information out of the event
133  evt.getByLabel(fFLSHitListLabel, hitlist);
134  if (hitlist->empty()){
135  std::cerr << "no flshit list for this event" << std::endl;
136  return;
137  }
138 
139  art::Ptr<sim::FLSHitList> hlist(hitlist, 0);
140  const std::vector<sim::FLSHit>& hits = hlist->fHits;
141 
142  int Number_of_showers_per_event = 0;
143 
145  const sim::ParticleNavigator& pnav = bt->ParticleNavigator();
146 
147  // loop over particles, we test if each makes the cuts, then fill histograms
148  for (sim::ParticleNavigator::const_iterator k = pnav.begin(); k != pnav.end(); ++k){
149 
150  const sim::Particle* part = (*k).second;
151 
152  // Yay! Cuts!
153  bool isElectronorPhoton = 0;
154  bool isVertexContained = 0;
155  bool isEnergyContained = 0;
156  bool isEnergeticEnough = 0;
157  float Energy_of_shower = 0;
158  for (unsigned int j=0; j<hits.size(); ++j) if(hits[j].GetTrackID() == part->TrackId()) Energy_of_shower += hits[j].GetEdep();
159  // run a cut on particle type +-11 = electron/positron, 22 = photon
160  if (abs(part->PdgCode()) == 11 || part->PdgCode() == 22)
161  isElectronorPhoton = 1;
162  // run a cut on the vertex, require it to be in the interior of the detector
163  if (abs(part->Vx()) < Half_Width_of_Detector[0] &&
164  abs(part->Vy()) < Half_Height_of_Detector[0] &&
165  part->Vz() < End_of_Detector_in_Z[0] &&
166  part->Vz() > Begin_of_Detector_in_Z[0])
167  isVertexContained = 1;
168  // run a cut on how much energy is deposited by flshits, a measure of containment
169  // the detector is 70% active, so if we have ~70% deposited, the shower is essentially contained
170  if (Energy_of_shower/part->E() > Energy_Percent[0])
171  isEnergyContained = 1;
172  // cut on energy, don't want particles without enough energy to develope a shower
173  if (part->E() > Min_Energy[0])
174  isEnergeticEnough = 1;
175 
176  // If the particle passes the cuts, run this so that we can fill histograms
177  if (isElectronorPhoton && isVertexContained && isEnergyContained && isEnergeticEnough){
178  ++ Number_of_showers_per_event;
179  float Event_distance = 0;
180  float Moliere_radius = 0;
181  for (unsigned int j=0; j<hits.size(); ++j)
182  if(hits[j].GetTrackID() == part->TrackId()) Energy_of_shower += hits[j].GetEdep();
183  // Fill hists with quick MCParticle variables for MC validation
184  fShowerMotherPDG->Fill(part->PdgCode());
185  fVx->Fill(part->Vx());
186  fVy->Fill(part->Vy());
187  fVz->Fill(part->Vz());
188  fTruthEnergy->Fill(part->E());
189  fCosTheta->Fill(part->Pz()/part->P());
190  fEnergyofShower->Fill(Energy_of_shower);
191  // Call the shower width function to calculate the average distance to FLS Hits
192  EMShowerWidth(Event_distance, part->TrackId(), part, hits);
193  // Fill histogram with average distance for event
194  fEventWidth->Fill(Event_distance);
195  // Call the Moliere Radius function to calculate Moliere Radius
196  EMMoliereRadius(Moliere_radius, part->TrackId(), part, hits);
197  // Fill histogram with Moliere radius for event
198  fMoliereRadius->Fill(Moliere_radius);
199  }
200 
201  } // end loop over particles
202 
203  // After running through all particles, we want to fill the number of showers in the event here
204  if (Number_of_showers_per_event > 0) fNumberPerEvent->Fill(Number_of_showers_per_event);
205 
206  } // end of main analyze function
207 
208 
209  ////////////////////////////////////////////////////////////////////////////
210 
211 
212  /// When this function is run, it finds the average distance, weighted by energy, to flshits over the shower
213  void ShowerAnaCheck::EMShowerWidth(float& Event_distance, int MomsID, const sim::Particle* part, const std::vector<sim::FLSHit>& hits){
214 
216 
217  float Energy_of_shower = 0;
218  float X_cent = 0;
219  float Y_cent = 0;
220  float Z_cent = 0;
221 
222  // this for loop calculates the centroid of the shower, weights by energy
223  for (unsigned int index1=0; index1<hits.size(); ++index1){
224 
225  if (hits[index1].GetTrackID() == MomsID){
226  // here XYZ is the average of where the hit enters and where it exits in the cell's frame
227  // XYZ is the global position
228  XYZ[0] = 0.5*(hits[index1].GetEntryX() + hits[index1].GetExitX());
229  XYZ[1] = 0.5*(hits[index1].GetEntryY() + hits[index1].GetExitY());
230  XYZ[2] = 0.5*(hits[index1].GetEntryZ() + hits[index1].GetExitZ());
231  geom->Plane(hits[index1].GetPlaneID())->Cell(hits[index1].GetCellID())->LocalToWorld(XYZ,XYZWorld);
232  // recursively add the coordinates of the hit
233  X_cent += hits[index1].GetEdep()*XYZWorld[0];
234  Y_cent += hits[index1].GetEdep()*XYZWorld[1];
235  Z_cent += hits[index1].GetEdep()*XYZWorld[2];
236  Energy_of_shower += hits[index1].GetEdep();
237  }
238 
239  }
240 
241  // normalize the centroid by the overall shower Energy
242  X_cent = X_cent/Energy_of_shower;
243  Y_cent = Y_cent/Energy_of_shower;
244  Z_cent = Z_cent/Energy_of_shower;
245 
246  // here we find the average distance to the fls hits, weighted by energy
247  for (unsigned int j=0; j<hits.size(); ++j)
248  {
249  if (hits[j].GetTrackID() == MomsID){
250  // here XYZ is the average of where the hit enters and where it exits in the cell's frame
251  // XYZ is the global position
252  XYZ[0] = 0.5*(hits[j].GetEntryX() + hits[j].GetExitX());
253  XYZ[1] = 0.5*(hits[j].GetEntryY() + hits[j].GetExitY());
254  XYZ[2] = 0.5*(hits[j].GetEntryZ() + hits[j].GetExitZ());
255  geom->Plane(hits[j].GetPlaneID())->Cell(hits[j].GetCellID())->LocalToWorld(XYZ,XYZWorld);
256  // Add the distance to each hit * energy in that hit / total energy = energy weighting scheme
257  Event_distance = Event_distance + hits[j].GetEdep()*PointLineDistance(XYZWorld[0],XYZWorld[1],XYZWorld[2],
258  part->Vx(),part->Vy(),part->Vz(),
259  X_cent,Y_cent,Z_cent)/Energy_of_shower;
260  }
261 
262  } // end of FLSHit hits loop
263 
264  } // end of EMShowerWidth()
265 
266 
267  ////////////////////////////////////////////////////////////////////////////
268 
269  /// This function calculates the Moliere Radius
270  void ShowerAnaCheck::EMMoliereRadius(float& Moliere_radius, int MomsID, const sim::Particle* part, const std::vector<sim::FLSHit>& hits){
271 
273 
274  float Energy_of_shower = 0;
275  float X_cent = 0;
276  float Y_cent = 0;
277  float Z_cent = 0;
278 
279  // Again, calculate the centroid of the showewr
280  for (unsigned int index1=0; index1<hits.size(); ++index1){
281  if (hits[index1].GetTrackID() == MomsID){
282  // here XYZ is the average of where the hit enters and where it exits in the cell's frame
283  // XYZ is the global position
284  XYZ[0] = 0.5*(hits[index1].GetEntryX() + hits[index1].GetExitX());
285  XYZ[1] = 0.5*(hits[index1].GetEntryY() + hits[index1].GetExitY());
286  XYZ[2] = 0.5*(hits[index1].GetEntryZ() + hits[index1].GetExitZ());
287  geom->Plane(hits[index1].GetPlaneID())->Cell(hits[index1].GetCellID())->LocalToWorld(XYZ,XYZWorld);
288  // recursively add the coordinates of the hit
289  X_cent += hits[index1].GetEdep()*XYZWorld[0];
290  Y_cent += hits[index1].GetEdep()*XYZWorld[1];
291  Z_cent += hits[index1].GetEdep()*XYZWorld[2];
292  Energy_of_shower += hits[index1].GetEdep();
293  }
294  }
295  // Normalize the centroid of the shower by the overall energy
296  X_cent = X_cent/Energy_of_shower;
297  Y_cent = Y_cent/Energy_of_shower;
298  Z_cent = Z_cent/Energy_of_shower;
299 
300  // Here we have a crude way to find the Moliere radius by counting the energy deposited in a cylinder
301  // around the shower axis. The first for loop slowly increases the radius, and the second goes over
302  // all hits and finds how much energy is deposited within the shell
303  for (float Moliere_parameter = 0; Moliere_parameter < 100; Moliere_parameter += .1){
304  float Energy_contained = 0;
305  for (unsigned int j = 0; j<hits.size(); ++j){
306  // here XYZ is the average of where the hit enters and where it exits in the cell's frame
307  // XYZ is the global position
308  XYZ[0] = 0.5*(hits[j].GetEntryX() + hits[j].GetExitX());
309  XYZ[1] = 0.5*(hits[j].GetEntryY() + hits[j].GetExitY());
310  XYZ[2] = 0.5*(hits[j].GetEntryZ() + hits[j].GetExitZ());
311  geom->Plane(hits[j].GetPlaneID())->Cell(hits[j].GetCellID())->LocalToWorld(XYZ,XYZWorld);
312  // if the hit is within a cylinder of radius moliere_parameter,
313  if (hits[j].GetTrackID() == MomsID && PointLineDistance(XYZWorld[0],XYZWorld[1],XYZWorld[2],
314  part->Vx(),part->Vy(),part->Vz(),X_cent,Y_cent,Z_cent) < Moliere_parameter)
315  Energy_contained += hits[j].GetEdep();
316  }
317  // Set Moliere_radius to the iterative constant above if <90% of energy is contained, else break
318  if (Energy_contained/Energy_of_shower <= .9) Moliere_radius = Moliere_parameter;
319  else break;
320  }
321 
322  }
323 
324  ///////////////////////////////////////////////////////////////////////////
325 
326  /// This function calculates the distance between any point and a line.
327  /// The point is the first three arguments (x,y,z), then the line is the
328  /// line passing through the other two points, (x1,y1,z1) and (x2,y2,z2)
329  float ShowerAnaCheck::PointLineDistance(float x, float y, float z, float x1, float y1, float z1, float x2, float y2, float z2){
330  float Lineone_point[3];
331  float Lineone_linetwo[3];
332  float Distance_one = 0;
333  float Distance_two = 0;
334  float Cosine = 0;
335  float Dot_product = 0;
336  Lineone_point[0]=x-x1;
337  Lineone_point[1]=y-y1;
338  Lineone_point[2]=z-z1;
339  Lineone_linetwo[0]=x2-x1;
340  Lineone_linetwo[1]=y2-y1;
341  Lineone_linetwo[2]=z2-z1;
342  Distance_one = util::pythag(Lineone_point[0], Lineone_point[1], Lineone_point[2]);
343  Distance_two = util::pythag(Lineone_linetwo[0], Lineone_linetwo[1], Lineone_linetwo[2]);
344  Dot_product = Lineone_point[0]*Lineone_linetwo[0]+Lineone_point[1]*Lineone_linetwo[1]+Lineone_point[2]*Lineone_linetwo[2];
345  Cosine = Dot_product/Distance_one/Distance_two;
346  return Distance_one*sqrt(1-util::sqr(Cosine));
347  }
348 
349 } // end mcchk
350 
351 ////////////////////////////////////////////////////////////////////////
352 namespace mcchk
353 {
355 }
double E(const int i=0) const
Definition: MCParticle.h:232
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
std::vector< sim::FLSHit > fHits
Definition: FLSHitList.h:21
Simple module to analyze MC cosmics distributions.
Float_t y1[n_points_granero]
Definition: compare.C:5
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
void LocalToWorld(const double *local, double *world) const
Definition: CellGeo.cxx:80
std::vector< float > Half_Height_of_Detector
Float_t x1[n_points_granero]
Definition: compare.C:5
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
T sqrt(T number)
Definition: d0nt_math.hpp:156
list_type::const_iterator const_iterator
std::vector< float > Half_Width_of_Detector
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
void reconfigure(const fhicl::ParameterSet &pset)
int TrackId() const
Definition: MCParticle.h:209
TString part[npart]
Definition: Style.C:32
void hits()
Definition: readHits.C:15
std::vector< float > End_of_Detector_in_Z
std::vector< float > Energy_Percent
double P(const int i=0) const
Definition: MCParticle.h:233
T get(std::string const &key) const
Definition: ParameterSet.h:231
std::vector< float > Min_Energy
const double j
Definition: BetheBloch.cxx:29
void EMShowerWidth(float &Event_distance, int MomsID, const sim::Particle *part, const std::vector< sim::FLSHit > &hits)
When this function is run, it finds the average distance, weighted by energy, to flshits over the sho...
void EMMoliereRadius(float &Moliere_radius, int MomsID, const sim::Particle *part, const std::vector< sim::FLSHit > &hits)
This function calculates the Moliere Radius.
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
ShowerAnaCheck(fhicl::ParameterSet const &pset)
z
Definition: test.py:28
check the results from the Monte Carlo generator
std::vector< float > Begin_of_Detector_in_Z
double Vx(const int i=0) const
Definition: MCParticle.h:220
::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
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
void geom(int which=0)
Definition: geom.C:163
double Pz(const int i=0) const
Definition: MCParticle.h:231
double Vz(const int i=0) const
Definition: MCParticle.h:222
float PointLineDistance(float x, float y, float z, float x1, float y1, float z1, float x2, float y2, float z2)
void analyze(const art::Event &evt)
Definition: fwd.h:28
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)