CosmicAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Simple module to analyze MC cosmics distributions
3 /// \author messier@indiana.edu
4 ////////////////////////////////////////////////////////////////////////
5 
6 // C/C++ includes
7 #include <cmath>
8 #include <map>
9 #include <vector>
10 
11 // ART includes
21 #include "fhiclcpp/ParameterSet.h"
23 
24 // NOvA includes
25 #include "Geometry/Geometry.h"
26 #include "GeometryObjects/Geo.h"
28 #include "MCCheater/BackTracker.h"
29 #include "Simulation/FLSHit.h"
30 #include "Simulation/FLSHitList.h"
31 #include "Simulation/Particle.h"
33 #include "Utilities/func/MathUtil.h"
34 
35 // ROOT includes
36 #include "TH1F.h"
37 #include "TH2F.h"
38 
39 namespace mcchk {
40  /// A module to make some plots of generated cosmic-rays
41  class CosmicAna : public art::EDAnalyzer {
42 
43  public:
44  explicit CosmicAna(fhicl::ParameterSet const& pset);
45  virtual ~CosmicAna();
46 
47  void analyze(art::Event const& evt);
48  void beginRun(art::Run const& run);
49 
50  private:
51  double fVtxBinSize; ///< Size of bins in cm for neutrino interaction vertex histograms
52  double fRockDepthX; ///< Amount of space in cm outside of detector along x to plot interaction points
53  double fRockDepthY; ///< Amount of space in cm outside of detector along y to plot interaction points
54  double fRockDepthZ; ///< Amount of space in cm outside of detector along z to plot interaction points
55 
56  std::map<int, std::string> fParticleMap; ///< Map of particle PDG and name
57 
58  TH1F* fDmin; ///< Closest approach for particles leaving hits
59  TH1F* fDminNoHit; ///< Closest approach for particles leaving no hits
60 
61  std::map<int, TH2F*> fParticleAngles;
62  std::map<int, TH2F*> fParticleAnglesLo; ///< Particle rate vs angle, low momenta
63  std::map<int, TH2F*> fParticleAnglesMi; ///< Particle rate vs angle, middle momenta
64  std::map<int, TH2F*> fParticleAnglesHi; ///< Particle rate vs angle, high momenta
65  std::map<int, TH1F*> fParticleCosTheta; ///< Particle rate vs cos(Q)
66  std::map<int, TH1F*> fParticleEnergy; ///< Particle energy (GeV)
67  std::map<int, TH2F*> fParticleEnergyVsAngle; ///< Particle Energy vs angle
68  std::map<int, TH2F*> fParticleVertexXY; ///< Particle vertex x vs y view
69  std::map<int, TH1F*> fParticleVertexZ; ///< Particle vertex z view
70 
71  TH1F* fVerticalMuons; ///< Number of vertical muons
72 
73  /// For each particle from ParticleNavigator...
74  TH1F* fCRCosX; ///< Particle X Direction cosine
75  TH1F* fCRCosY; ///< Particle Y Direction cosine
76  TH1F* fCRCosZ; ///< Particle Z Direction cosine
77  TH1F* fCRVtxX;
78  TH1F* fCRVtxY;
79  TH1F* fCRVtxZ;
80  TH1F* fCREndX;
81  TH1F* fCREndY;
82  TH1F* fCREndZ;
83  TH1F* fCRAvgCellsX;
84  TH1F* fCRAvgCellsY;
85  TH1F* fCREnergy;
86  TH1F* fCRLowEnergy;
87  TH1F* fCRLength;
88 
89  TH1F* fSampleTime; ///< Total sample time, number of events times 500 usec
90  };
91 }
92 
93 ///////////////////////////////////////////////////////////////////////////////
94 namespace mcchk {
95  //......................................................................
97  : EDAnalyzer(pset),
98  fVtxBinSize(pset.get<double>("VtxBinSize")),
99  fRockDepthX(pset.get<double>("RockDepthX")),
100  fRockDepthY(pset.get<double>("RockDepthY")),
101  fRockDepthZ(pset.get<double>("RockDepthZ")),
102  fDmin(0) // Flag that histograms are uninitialized
103  {
104  }
105 
106  //......................................................................
108  {
109  }
110 
111  //......................................................................
113  {
114  // Only make histograms once
115  if(fDmin) { return; }
116 
119 
120  // Create bounds and bins for vertex plots
121  // Get detector coordinates
122  double detX = geo->DetHalfWidth(); // -detX is minimum, +detX is maximum value
123  double detY = geo->DetHalfHeight(); // Likewise, use +/- detY
124  double detZ = geo->DetLength(); // Here, mininmum value is 0, max value is detZ
125 
126  // Add on the space for rock events
127  double CRPosX = detX + fRockDepthX;
128  double CRPosY = detY + fRockDepthY;
129  double CRPosZMin = -fRockDepthZ;
130  double CRPosZMax = detZ + fRockDepthZ;
131 
132  // This section adds a little extra space to make bins sizes exact
133  int nbinsHelper = 0; // Helper to make bin sizes exact
134 
135  // First calculate number of complete bins from min to max values
136  // Add one more bin for the last partial bin
137  // For X/Y, half of those bins are above 0, half below
138  // The max value is then the bin size times then number of bins over 2
139  // For Z, calculate the adjusted number of bins in the same way
140  // Leave the min value unchanged
141  // Set the max value as the min value plus the number of bins times the bin size
142 
143  // Set the correct bounds
144  nbinsHelper = (int)(2.*CRPosX/fVtxBinSize);
145  ++nbinsHelper; // Add one extra bin to cover truncated area
146  CRPosX = 0.5*fVtxBinSize*(double)nbinsHelper;
147 
148  nbinsHelper = (int)(2.*CRPosY/fVtxBinSize);
149  ++nbinsHelper; // Add one extra bin to cover truncated area
150  CRPosY = 0.5*fVtxBinSize*(double)nbinsHelper;
151 
152  nbinsHelper = (int)((CRPosZMax - CRPosZMin)/fVtxBinSize);
153  ++nbinsHelper; // Add one extra bin to cover truncated area
154  CRPosZMax = CRPosZMin + (fVtxBinSize*(double)nbinsHelper);
155 
156  // Set the correct number of bins
157  int nBinsX = 2.*CRPosX/fVtxBinSize;
158  int nBinsY = 2.*CRPosY/fVtxBinSize;
159  int nBinsZ = (CRPosZMax - CRPosZMin)/fVtxBinSize;
160 
161  // Helpful for spacing/readability
162  // Also helpful in case histograms are made in loops over, say, axis labels
163  char HistoName[200];
164  char HistoTitle[200];
165  char TitleHelper[200];
166 
167  sprintf(TitleHelper, "Distance of Closest Approach to Detector Center for Particles That Left");
168 
169  sprintf(HistoName, "fDmin");
170  sprintf(HistoTitle, "%s Hits;d (cm);", TitleHelper);
171  fDmin = tfs->make<TH1F>(HistoName, HistoTitle, 100, 0., 500.E2);
172 
173  sprintf(HistoName, "fDminNoHit");
174  sprintf(HistoTitle, "%s No Hits;d (cm);", TitleHelper);
175  fDminNoHit = tfs->make<TH1F>(HistoName, HistoTitle, 100, 0., 500.E2);
176 
177  fParticleMap[11] = "Electron";
178  fParticleMap[13] = "Muon";
179  fParticleMap[22] = "Photon";
180  fParticleMap[2112] = "Neutron";
181  fParticleMap[2212] = "Proton";
182 
183  for(const auto& particlePair : fParticleMap) {
184  int pdg = particlePair.first;
185  std::string particleName = particlePair.second;
186 
187  sprintf(TitleHelper, "%s Angle Distribution", particleName.c_str());
188 
189  sprintf(HistoName, "f%sAngles", particleName.c_str());
190  sprintf(HistoTitle, "%s;#phi;cos#theta", TitleHelper);
191  fParticleAngles[pdg] = tfs->make<TH2F>(HistoName, HistoTitle,
192  36, -180., 180., 50, -1., 1.);
193 
194  sprintf(HistoName, "f%sAnglesLo", particleName.c_str());
195  sprintf(HistoTitle, "%s for Low Momentum %ss;#phi;cos#theta", TitleHelper, particleName.c_str());
196  fParticleAnglesLo[pdg] = tfs->make<TH2F>(HistoName, HistoTitle,
197  36, -180., 180., 50, -1., 1.);
198 
199  sprintf(HistoName, "f%sAnglesMi", particleName.c_str());
200  sprintf(HistoTitle, "%s for Medium Momentum %ss;#phi;cos#theta", TitleHelper, particleName.c_str());
201  fParticleAnglesMi[pdg] = tfs->make<TH2F>(HistoName, HistoTitle,
202  36, -180., 180., 50, -1., 1.);
203 
204  sprintf(HistoName, "f%sAnglesHi", particleName.c_str());
205  sprintf(HistoTitle, "%s for High Momentum %ss;#phi;cos#theta", TitleHelper, particleName.c_str());
206  fParticleAnglesHi[pdg] = tfs->make<TH2F>(HistoName, HistoTitle,
207  36, -180., 180., 50, -1., 1.);
208 
209  sprintf(HistoName, "f%sCosTheta", particleName.c_str());
210  sprintf(HistoTitle, "%s cos#theta;cos#theta;", TitleHelper);
211  fParticleCosTheta[pdg] = tfs->make<TH1F>(HistoName, HistoTitle, 50, -1., 1.);
212 
213  sprintf(HistoName, "f%sVertexXY", particleName.c_str());
214  sprintf(HistoTitle, "%s Vertex, X vs Y View;x (cm);y (cm)", particleName.c_str());
215  fParticleVertexXY[pdg] = tfs->make<TH2F>(HistoName, HistoTitle,
216  nBinsX, -CRPosX, CRPosX, nBinsY, -CRPosY, CRPosY);
217 
218  sprintf(HistoName, "f%sVertexZ", particleName.c_str());
219  sprintf(HistoTitle, "%s Vertex, Z View;z (cm);", particleName.c_str());
220  fParticleVertexZ[pdg] = tfs->make<TH1F>(HistoName, HistoTitle, nBinsZ, CRPosZMin, CRPosZMax);
221  }
222 
223  // Make plots with energy by hand, as particles have vastly different energy scales
224  sprintf(TitleHelper, "Energy;E (Gev);");
225 
226  sprintf(HistoName, "fElectronEnergy");
227  sprintf(HistoTitle, "Electron %s", TitleHelper);
228  fParticleEnergy[11] = tfs->make<TH1F>(HistoName, HistoTitle, 1000, 0., 1.);
229 
230  sprintf(HistoName, "fMuonEnergy");
231  sprintf(HistoTitle, "Muon %s", TitleHelper);
232  fParticleEnergy[13] = tfs->make<TH1F>(HistoName, HistoTitle, 10000, 0., 1000.);
233 
234  sprintf(HistoName, "fPhotonEnergy");
235  sprintf(HistoTitle, "Photon %s", TitleHelper);
236  fParticleEnergy[22] = tfs->make<TH1F>(HistoName, HistoTitle, 200, 0., 5.);
237 
238  sprintf(HistoName, "fNeutronEnergy");
239  sprintf(HistoTitle, "Neutron %s", TitleHelper);
240  fParticleEnergy[2112] = tfs->make<TH1F>(HistoName, HistoTitle, 100, 0., 10.);
241 
242  sprintf(HistoName, "fProtonEnergy");
243  sprintf(HistoTitle, "Proton %s", TitleHelper);
244  fParticleEnergy[2212] = tfs->make<TH1F>(HistoName, HistoTitle, 100, 0., 10.);
245 
246  sprintf(TitleHelper, "Energy vs cos#theta for ");
247 
248  sprintf(HistoName, "fElectronEnergyVsAngle");
249  sprintf(HistoTitle, "%s Electrons;cos#theta;Energy (GeV)", TitleHelper);
250  fParticleEnergyVsAngle[11] = tfs->make<TH2F>(HistoName, HistoTitle, 200, -1., 1., 100, 0., 0.1);
251 
252  sprintf(HistoName, "fMuonEnergyVsAngle");
253  sprintf(HistoTitle, "%s Muons;cos#theta;Energy (GeV)", TitleHelper);
254  fParticleEnergyVsAngle[13] = tfs->make<TH2F>(HistoName, HistoTitle, 200, -1., 1., 10000, 0., 1000.);
255 
256  sprintf(HistoName, "fPhotonEnergyVsAngle");
257  sprintf(HistoTitle, "%s Photons;cos#theta;Energy (GeV)", TitleHelper);
258  fParticleEnergyVsAngle[22] = tfs->make<TH2F>(HistoName, HistoTitle, 200, -1., 1., 100, 0., 0.1);
259 
260  sprintf(HistoName, "fNeutronEnergyVsAngle");
261  sprintf(HistoTitle, "%s Neutrons;cos#theta;Energy (GeV)", TitleHelper);
262  fParticleEnergyVsAngle[2112] = tfs->make<TH2F>(HistoName, HistoTitle, 200, -1., 1., 100, 0., 10.);
263 
264  sprintf(HistoName, "fProtonEnergyVsAngle");
265  sprintf(HistoTitle, "%s Protons;cos#theta;Energy (GeV)", TitleHelper);
266  fParticleEnergyVsAngle[2212] = tfs->make<TH2F>(HistoName, HistoTitle, 200, -1., 1., 100, 0., 10.);
267 
268  sprintf(HistoName, "fVerticalMuons");
269  sprintf(HistoTitle, "Number of Vertical Muons;");
270  fVerticalMuons = tfs->make<TH1F>(HistoName, HistoTitle, 1, -0.5, 0.5);
271 
272  sprintf(HistoName, "fCRCosX");
273  sprintf(HistoTitle, "Particle X Direction Cosine;dx/ds;Tracks");
274  fCRCosX = tfs->make<TH1F>(HistoName, HistoTitle, 200, -1., 1.);
275 
276  sprintf(HistoName, "fCRCosY");
277  sprintf(HistoTitle, "Particle Y Direction Cosine;dy/ds;Tracks");
278  fCRCosY = tfs->make<TH1F>(HistoName, HistoTitle, 200, -1., 1.);
279 
280  sprintf(HistoName, "fCRCosZ");
281  sprintf(HistoTitle, "Particle Z Direction Cosine;dz/ds;Tracks");
282  fCRCosZ = tfs->make<TH1F>(HistoName, HistoTitle, 200, -1., 1.);
283 
284  sprintf(TitleHelper, "Particle Entry Location, ");
285 
286  sprintf(HistoName, "fCRVtxX");
287  sprintf(HistoTitle, "%s X;x (cm);", TitleHelper);
288  fCRVtxX = tfs->make<TH1F>(HistoName, HistoTitle, nBinsX, -CRPosX, CRPosX);
289 
290  sprintf(HistoName, "fCRVtxY");
291  sprintf(HistoTitle, "%s Y;y (cm);", TitleHelper);
292  fCRVtxY = tfs->make<TH1F>(HistoName, HistoTitle, nBinsY, -CRPosY, CRPosY);
293 
294  sprintf(HistoName, "fCRVtxZ");
295  sprintf(HistoTitle, "%s Z;z (cm);", TitleHelper);
296  fCRVtxZ = tfs->make<TH1F>(HistoName, HistoTitle, nBinsZ, CRPosZMin, CRPosZMax);
297 
298  sprintf(TitleHelper, "Particle Exit Location, ");
299 
300  sprintf(HistoName, "fCREndX");
301  sprintf(HistoTitle, "%s X;x (cm);", TitleHelper);
302  fCREndX = tfs->make<TH1F>(HistoName, HistoTitle, nBinsX, -CRPosX, CRPosX);
303 
304  sprintf(HistoName, "fCREndY");
305  sprintf(HistoTitle, "%s Y;y (cm);", TitleHelper);
306  fCREndY = tfs->make<TH1F>(HistoName, HistoTitle, nBinsY, -CRPosY, CRPosY);
307 
308  sprintf(HistoName, "fCREndZ");
309  sprintf(HistoTitle, "%s Z;z (cm);", TitleHelper);
310  fCREndZ = tfs->make<TH1F>(HistoName, HistoTitle, nBinsZ, CRPosZMin, CRPosZMax);
311 
312  sprintf(TitleHelper, "Average Number of Cells Hit in ");
313 
314  sprintf(HistoName, "fCRAvgCellsX");
315  sprintf(HistoTitle, "%s X View;Cells;", TitleHelper);
316  fCRAvgCellsX = tfs->make<TH1F>(HistoName, HistoTitle, 100, 0., 100.);
317 
318  sprintf(HistoName, "fCRAvgCellsY");
319  sprintf(HistoTitle, "%s Y View;Cells;", TitleHelper);
320  fCRAvgCellsY = tfs->make<TH1F>(HistoName, HistoTitle, 100, 0., 100.);
321 
322  sprintf(HistoName, "fCREnergy");
323  sprintf(HistoTitle, "Energy of Particle;Energy (GeV);");
324  fCREnergy = tfs->make<TH1F>(HistoName, HistoTitle, 1000, 0., 1000.);
325 
326  sprintf(HistoName, "fCRLowEnergy");
327  sprintf(HistoTitle, "Energy of Low Energy Particles;Energy (GeV);");
328  fCRLowEnergy = tfs->make<TH1F>(HistoName, HistoTitle, 200, 0., 2.);
329 
330  sprintf(HistoName, "fCRLength");
331  sprintf(HistoTitle, "Length of Particle in Detector;Length (cm);");
332  fCRLength = tfs->make<TH1F>(HistoName, HistoTitle, 500, 0., 5000.);
333 
334  sprintf(HistoName, "fSampleTime");
335  sprintf(HistoTitle, "Total Sample Time;Time (s);");
336  fSampleTime = tfs->make<TH1F>(HistoName, HistoTitle, 1, -0.5, 0.5);
337  }
338 
339  //......................................................................
341  {
342  fSampleTime->Fill(0., 500.E-6); // Add 500 usec
343 
345 
347  const sim::ParticleNavigator &pnav = bt->ParticleNavigator();
348 
349  // Flags for "interesting" particles
350  double electronE = -1.0;
351  double photonE = -1.0;
352 
353  // Loop over particles in the event
354  for(sim::ParticleNavigator::const_iterator it = pnav.begin(); it != pnav.end(); ++it) {
355  const sim::Particle* particle = (*it).second;
356 
357  // Get a list of FLSHits corresponding to the current particle
358  const std::vector<sim::FLSHit>& flshits = bt->ParticleToFLSHit(particle->TrackId());
359  const int nflshits = flshits.size();
360 
361  const TLorentzVector& fourPosition = particle->Position();
362  const TLorentzVector& fourMomentum = particle->Momentum();
363 
364  // Inputs to ClosestApproach function
365  // Closest approach to this point
366  double detCenter[3] = {0., 0., 0.5*geo->DetLength()};
367  // Point on the line--calculate closest approach of this line to point above
368  double particleVtx[3] = {fourPosition.X(), fourPosition.Y(), fourPosition.Z()};
369  // Slope of line
370  double particleP[3] = {fourMomentum.Px(), fourMomentum.Py(), fourMomentum.Pz()};
371  double closestApp[3]; // Return vector
372  double distClosestApp = geo::ClosestApproach(detCenter, particleVtx, particleP, closestApp);
373 
374  if(nflshits == 0) {
375  fDminNoHit->Fill(distClosestApp);
376  }
377  else {
378  fDmin->Fill(distClosestApp);
379  if(distClosestApp > 10000.) {
380  mf::LogVerbatim("Dmin > 10000") << "Run " << evt.run() << ", Event " << evt.id().event() << std::endl
381  << "Particle PDG " << particle->PdgCode() << std::endl
382  << "Vertex " << fourPosition.X() << " " << fourPosition.Y() << " " << fourPosition.Z() << std::endl
383  << "Momentum" << fourMomentum.Px() << " " << fourMomentum.Py() << " " << fourMomentum.Pz() << std::endl
384  << "Distance of closest approach to detector center " << distClosestApp << std::endl;
385  }
386  } // end of (else) conditionals if size of FLSHit list from Particle is 0
387 
388  if(nflshits > 0) {
389  int pdg = abs(particle->PdgCode());
390 
391  double energy = fourMomentum.E();
392  double cosTheta = -fourMomentum.Py()/fourMomentum.P();
393  double phi = atan2(fourMomentum.Pz(), fourMomentum.Px());
394  phi *= 180./M_PI;
395 
396  // Add a vertical muon if angle is less than 10 degrees
397  if(pdg == 13 && cosTheta > 0.984808) {
398  fVerticalMuons->Fill(0.);
399  }
400 
401  // Proceed only if the current particle is one we care about
402  if(fParticleMap.find(pdg) != fParticleMap.end()) {
403  fParticleCosTheta[pdg]->Fill(cosTheta);
404  fParticleAngles[pdg] ->Fill(phi, cosTheta);
405 
406  // Fill angle plots based on momentum
407  if (energy < 1.) { fParticleAnglesLo[pdg]->Fill(phi, cosTheta); }
408  else if(energy < 10.) { fParticleAnglesMi[pdg]->Fill(phi, cosTheta); }
409  else { fParticleAnglesHi[pdg]->Fill(phi, cosTheta); }
410 
411  // Fill particle energy information
412  fParticleEnergy[pdg] ->Fill(energy);
413  fParticleEnergyVsAngle[pdg]->Fill(cosTheta, energy);
414 
415  // Fill particle vertex information
416  fParticleVertexXY[pdg]->Fill(fourPosition.X(), fourPosition.Y());
417  fParticleVertexZ[pdg] ->Fill(fourPosition.Z());
418 
419  // Create matrix (vector of vectors) of number of hits in each cell
420  std::vector<std::vector<double> > hitsByPlaneCell(geo->GetPlanesByView().size());
421  unsigned int numCells = geo->Plane(0)->Ncells();
422  if(geo->Plane(1)->Ncells() > numCells) { numCells = geo->Plane(1)->Ncells(); }
423 
424  // Resize each vector in the larger plane vector to have the number of cells, initialized to 0
425  for(unsigned int i_plane = 0, n_plane = hitsByPlaneCell.size(); i_plane < n_plane; ++i_plane) {
426  hitsByPlaneCell[i_plane].resize(numCells, 0.);
427  }
428 
429  // Vector matching a plane to a view
430  std::vector<geo::View_t> planeView(geo->GetPlanesByView().size());
431 
432  // Find earliest and latest plane/cell, entry/exit point
433  double startT = flshits[0].GetEntryT();
434  double endT = flshits[0].GetExitT();
435  int startPlane = flshits[0].GetPlaneID();
436  int endPlane = startPlane;
437  int startCell = flshits[0].GetCellID();
438  int endCell = startCell;
439 
440  double entryPos[3] = {0.};
441  double exitPos[3] = {0.};
442 
443  double localStartPos[3] = {flshits[0].GetEntryX(),
444  flshits[0].GetEntryY(),
445  flshits[0].GetEntryZ()};
446  double localEndPos[3] = {flshits[0].GetExitX(),
447  flshits[0].GetExitY(),
448  flshits[0].GetExitZ()};
449 
450  // Loop over FLSHits for this particle
451  for(const auto& hit : flshits) {
452  if (hit.GetEntryT() < startT) { // Update quantities if this is an earlier hit
453  startPlane = hit.GetPlaneID();
454  startCell = hit.GetCellID();
455  localStartPos[0] = hit.GetEntryX();
456  localStartPos[1] = hit.GetEntryY();
457  localStartPos[2] = hit.GetEntryZ();
458  }
459  else if(hit.GetExitT() > endT) { // Update quantities if this is a later hit
460  endPlane = hit.GetPlaneID();
461  endCell = hit.GetCellID();
462  localEndPos[0] = hit.GetExitX();
463  localEndPos[1] = hit.GetExitY();
464  localEndPos[2] = hit.GetExitZ();
465  }
466 
467  hitsByPlaneCell[hit.GetPlaneID()][hit.GetCellID()] += 1.; // Add a hit to the proper cell
468  planeView[hit.GetPlaneID()] = geo->Plane(hit.GetPlaneID())->View();
469  } // end of loop over FLSHits
470 
471  // Convert entry/exit locations to world coordinates
472  geo->Plane(startPlane)->Cell(startCell)->LocalToWorld(localStartPos, entryPos);
473  geo->Plane(endPlane) ->Cell(endCell) ->LocalToWorld(localEndPos, exitPos);
474 
475  double planesX = 0.;
476  double cellsX = 0.;
477  double planesY = 0.;
478  double cellsY = 0.;
479 
480  // Loop over planes, cells
481  for(unsigned int i_plane = 0, n_plane = hitsByPlaneCell.size(); i_plane < n_plane; ++i_plane) {
482  // Count the number of hit cells in this plane
483  double hitCells = 0.;
484  for(unsigned int i_cell = 0, n_cell = hitsByPlaneCell[i_plane].size(); i_cell < n_cell; ++i_cell) {
485  if(hitsByPlaneCell[i_plane][i_cell] > 0.) { hitCells += 1.; }
486  }
487 
488  if(planeView[i_plane] == geo::kX && hitCells > 0.) {
489  planesX += 1.;
490  cellsX += hitCells;
491  }
492  if(planeView[i_plane] == geo::kY && hitCells > 0.) {
493  planesY += 1.;
494  cellsY += hitCells;
495  }
496  } // end loop over hit planes/cells
497 
498  // Fill final plots
499  fCRCosX->Fill(fourMomentum.Px()/fourMomentum.P());
500  fCRCosY->Fill(fourMomentum.Py()/fourMomentum.P());
501  fCRCosZ->Fill(fourMomentum.Pz()/fourMomentum.P());
502 
503  // Cosmic rays always come from above
504  if(exitPos[1] > entryPos[1]) {
505  fCREndX->Fill(entryPos[0]);
506  fCREndY->Fill(entryPos[1]);
507  fCREndZ->Fill(entryPos[2]);
508  fCRVtxX->Fill(exitPos[0]);
509  fCRVtxY->Fill(exitPos[1]);
510  fCRVtxZ->Fill(exitPos[2]);
511  }
512  else { // Normal condition
513  fCRVtxX->Fill(entryPos[0]);
514  fCRVtxY->Fill(entryPos[1]);
515  fCRVtxZ->Fill(entryPos[2]);
516  fCREndX->Fill(exitPos[0]);
517  fCREndY->Fill(exitPos[1]);
518  fCREndZ->Fill(exitPos[2]);
519  }
520 
521  if(planesX > 0.) { fCRAvgCellsX->Fill(cellsX/planesX); }
522  else { fCRAvgCellsX->Fill(0.); }
523 
524  if(planesY > 0.) { fCRAvgCellsY->Fill(cellsY/planesY); }
525  else { fCRAvgCellsY->Fill(0.); }
526 
527  fCREnergy->Fill(energy); // Fill particle energy
528 
529  // Fill low energy particles for higher resolution
530  if(energy < 2.) {
531  fCRLowEnergy->Fill(energy);
532  }
533 
534  fCRLength->Fill(util::pythag(exitPos[0]-entryPos[0],
535  exitPos[1]-entryPos[1],
536  exitPos[2]-entryPos[2]));
537 
538  // Check for interesting particles
539  if(energy > 0.5 && energy < 5.) {
540  if(pdg == 11) { electronE = energy; }
541  if(pdg == 22) { photonE = energy; }
542  }
543  } // end of condition that current particle is in the particle map
544  } // end of condition that size of list of FLSHits is greater than 0
545  } // end of loop over particles in ParticleNavigator
546 
547  if(electronE > 0. || photonE > 0.) {
548  mf::LogVerbatim("CosmicAna") << "Interesting particle:" << std::endl
549  << "Run" << evt.getRun().run() << ", "
550  << "Event " << evt.id().event() << std::endl
551  << "Electron energy " << electronE << std::endl
552  << "Photon energy " << photonE << std::endl;
553  }
554  }
555 } // end of namespace mcchk
556 
557 ///////////////////////////////////////////////////////////////////////////////
558 namespace mcchk {
560 }
std::map< int, TH1F * > fParticleEnergy
Particle energy (GeV)
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
set< int >::iterator it
Simple module to analyze MC cosmics distributions.
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
void LocalToWorld(const double *local, double *world) const
Definition: CellGeo.cxx:80
std::map< int, TH2F * > fParticleAnglesHi
Particle rate vs angle, high momenta.
std::map< int, TH2F * > fParticleAngles
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
TH1F * fCRCosY
Particle Y Direction cosine.
list_type::const_iterator const_iterator
Vertical planes which measure X.
Definition: PlaneGeo.h:28
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
double fRockDepthZ
Amount of space in cm outside of detector along z to plot interaction points.
Run const & getRun() const
double DetLength() const
void abs(TH1 *hist)
TH1F * fDminNoHit
Closest approach for particles leaving no hits.
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
RunNumber_t run() const
Definition: Run.h:47
#define M_PI
Definition: SbMath.h:34
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
std::map< int, TH2F * > fParticleEnergyVsAngle
Particle Energy vs angle.
Definition: Run.h:31
int TrackId() const
Definition: MCParticle.h:209
double fRockDepthX
Amount of space in cm outside of detector along x to plot interaction points.
void analyze(art::Event const &evt)
void beginRun(art::Run const &run)
Float_t E
Definition: plot.C:20
A module to make some plots of generated cosmic-rays.
double fRockDepthY
Amount of space in cm outside of detector along y to plot interaction points.
double energy
Definition: plottest35.C:25
Collect Geo headers and supply basic geometry functions.
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
TH1F * fCRCosZ
Particle Z Direction cosine.
Definition: View.py:1
Definition: run.py:1
std::map< int, TH1F * > fParticleCosTheta
Particle rate vs cos(Q)
CosmicAna(fhicl::ParameterSet const &pset)
EventNumber_t event() const
Definition: EventID.h:116
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
double DetHalfWidth() const
double ClosestApproach(const double point[], const double intercept[], const double slopes[], double closest[])
Find the distance of closest approach between point and line.
Definition: Geo.cxx:222
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
Definition: event.h:1
TH1F * fDmin
Closest approach for particles leaving hits.
std::map< int, TH2F * > fParticleAnglesMi
Particle rate vs angle, middle momenta.
const std::set< unsigned int > & GetPlanesByView(View_t v=kXorY) const
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
std::map< int, std::string > fParticleMap
Map of particle PDG and name.
std::map< int, TH2F * > fParticleAnglesLo
Particle rate vs angle, low momenta.
std::map< int, TH1F * > fParticleVertexZ
Particle vertex z view.
TH1F * fSampleTime
Total sample time, number of events times 500 usec.
TH1F * fCRCosX
For each particle from ParticleNavigator...
RunNumber_t run() const
Definition: Event.h:77
Helper for AttenCurve.
Definition: Path.h:10
TH1F * fVerticalMuons
Number of vertical muons.
std::map< int, TH2F * > fParticleVertexXY
Particle vertex x vs y view.
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
double fVtxBinSize
Size of bins in cm for neutrino interaction vertex histograms.
T atan2(T number)
Definition: d0nt_math.hpp:72