NeutrinoAna_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 /// \brief
3 /// \author messier@indiana.edu
4 ///////////////////////////////////////////////////////////////////////////////
5 
6 // C/C++ Includes
7 #include <cmath>
8 #include <map>
9 #include <vector>
10 
11 // ART includes
20 #include "fhiclcpp/ParameterSet.h"
23 
24 // NOvA includes
25 #include "Geometry/Geometry.h"
30 #include "MCReweight/FluxWeights.h"
31 
32 // ROOT includes
33 #include "TF1.h"
34 #include "TH1F.h"
35 #include "TH2.h"
36 #include "TH3.h"
37 #include "TSystem.h"
38 
39 namespace mcchk {
40  /// Check the results from the Monte Carlo generator
41  /// This module analyzes neutrinos and neutrino interactions
42  class NeutrinoAna : public art::EDAnalyzer {
43 
44  public:
45  explicit NeutrinoAna(fhicl::ParameterSet const& pset);
46  virtual ~NeutrinoAna();
47 
48  void analyze(art::Event const& evt);
49  void beginJob();
50  void beginRun(art::Run const& run);
51 
52  double DotProduct(double*, double*); ///< Calculate dot product between two vectors
53 
54  private:
55  /// Calculate true L, E, L/E from neutrino interaction information
57  const art::Ptr<simb::MCFlux>& mcflux);
58 
59  bool fDetOnly;
60 
61  std::string fMCTruthModuleLabel; ///< label for module producing mc truth information
62  std::string fFluxReweightModuleLabel; ///< label for module producing flux reweight information
63 
64  int fEnu_Nbins; ///< Number of bins for energy histograms
65  double fEnu_min; ///< Minimum energy for energy histograms
66  double fEnu_max; ///< Maximum energy for energy histgorams
67 
68  double fVtxBinSize; ///< Size of bins in cm for neutrino interaction vertex histograms
69  double fRockDepthX; ///< Amount of space in cm outside of detector along x to plot neutrino interaction points
70  double fRockDepthY; ///< Amount of space in cm outside of detector along y to plot neutrino interaction points
71  double fRockDepthZ; ///< Amount of space in cm outside of detector along z to plot neutrino interaction points
72 
73  /// Neutrino information
74 
75  TH1F* fEnu_nue_CC; ///< Neutrino energy for nue CC events
76  TH1F* fEnu_nuebar_CC; ///< Nuetrino energy for anti-nue CC events
77  TH1F* fEnu_numu_CC; ///< Neutrino energy for numu CC events
78  TH1F* fEnu_numubar_CC; ///< Neutrino energy for anti-numu CC events
79  TH1F* fEnu_nu_NC; ///< Neutrino energy for neutrino NC events
80  TH1F* fEnu_nubar_NC; ///< Neutrino energy for anti-nu NC events
81 
82  std::map<int, TH1F*> fEnu_nue_ByParent; ///< Neutrino energy by parent species, nu_e(bar)s
83  std::map<int, TH1F*> fEnu_numu_ByParent; ///< Neutrino energy by parent species, nu_mu(bar)s
84  std::map<int, TH1F*> fEnu_other_ByParent; ///< Neutrino energy by parent species, NC/nu_tau(bar)s
85  std::map<int, TH1F*> fEnu_ByMode; ///< Neutrino energy by interaction mode
86 
87  TH1F* fNeutrinoStartingPositionX; ///< Neutrino Production Point X Coordinate
88  TH1F* fNeutrinoStartingPositionY; ///< Neutrino Production Point Y Coordinate
89  TH1F* fNeutrinoStartingPositionZ; ///< Neutrino Production Point Z Coordinate
90  TH2F* fNeutrinoStartingPositionXY; ///< Neutrino Production Point X vs Y Coordinates
91  TH2F* fNeutrinoStartingPositionXZ; ///< Neutrino Production Point X vs Z Coordinates
92  TH2F* fNeutrinoStartingPositionYZ; ///< Neutrino Production Point Y vs Z Coordinates
93 
94  TH1F* fNeutrinoInteractionVertexX; ///< Neutrino Interaction Vertex X Coordinate
95  TH1F* fNeutrinoInteractionVertexY; ///< Neutrino Interaction Vertex Y Coordinate
96  TH1F* fNeutrinoInteractionVertexZ; ///< Neutrino Interaction Vertex Z Coordinate
97  TH2F* fNeutrinoInteractionVertexXY; ///< Neutrino Interaction Vertex X vs Y Coordinates
98  TH2F* fNeutrinoInteractionVertexXZ; ///< Neutrino Interaction Vertex X vs Z Coordinates
99  TH2F* fNeutrinoInteractionVertexYZ; ///< Neutrino Interaction Vertex Y vs Z Coordinates
100  //TH3F* fNeutrinoInteractionVertex; // Neutrino Interaction Vertex
101 
102  TH1F* fNuCosX; ///< X direction cosine of incoming neutrino
103  TH1F* fNuCosY; ///< Y direction cosine of incoming neutrino
104  TH1F* fNuCosZ; ///< Z direction cosine of incoming neutrino
105 
106  TH1F* fMode; ///< Interaction mode
107 
108  TH1F* fNeutrinoTrueL; ///< True L for neutrino
109  TH2F* fNeutrinoTrueETrueL; ///< True L and E of the neutrino
110  TH1F* fNeutrinoTrueLOverE; ///< True L/E of the neutrino
111 
112  //TH1F* fProtonMultiplicityAllNucl; ///< Proton multiplicity
113  std::map<int, TH1F*> fProtonMultiplicity; ///< Proton multiplicity by nucleus PDG
114 
115  TH1F* fW; ///< W (shower mass)
116  TH1F* fX; ///< Bjorken X
117  TH1F* fY; ///< Bjorken Y
118  TH1F* fQSqr; ///< Q^2
119 
120  TH1F* fYUnscaled; ///< Difference between neutrino and charged lepton energies
121 
122  std::vector<std::string> fModeBinIndices; // used internally to label bins in interaction mode histogram
123 
124  const std::map<simb::int_type_, std::string> fModeNames
125  {
126  { simb::kQE, "QE" },
127  { simb::kMEC, "MEC" },
128  { simb::kRes, "Res" },
129  { simb::kDIS, "DIS" },
130  { simb::kCoh, "Coh" },
131  { simb::kUnknownInteraction, "Other" },
132  };
133 
134  //PPFX weights:
135  int fPPFXWgt_Nbins; ///< Number of bins for PPFX weight histograms
136  double fPPFXWgt_min; ///< Minimum weight PPFX weight histograms
137  double fPPFXWgt_max; ///< Maximum weight PPFX weight histograms
138  TH1F* fPPFXWgt_cv; ///< Central value weights
139  TH1F* fPPFXWgt_univ; ///< Multi-universe weights
140  TH1F* fPPFXWgtEnu_numu; ///< Neutrino energy for numu events weighted by PPFX
141  TH1F* fPPFXWgtEnu_numubar; ///< Neutrino energy for anti-numu events weighted by PPFX
142  TH1F* fPPFXWgtEnu_nue; ///< Neutrino energy for nue events weighted by PPFX
143  TH1F* fPPFXWgtEnu_nuebar; ///< Neutrino energy for nuebar events weighted by PPFX
144  };
145 }
146 
147 ///////////////////////////////////////////////////////////////////////////////
148 namespace mcchk {
149  static int msg1cnt = 0;
150 
151  //......................................................................
153  : EDAnalyzer(pset),
154  fDetOnly(pset.get<bool> ("DetOnly")),
155  fMCTruthModuleLabel(pset.get<std::string>("MCTruthModuleLabel")),
156  fFluxReweightModuleLabel(pset.get<std::string>("FluxReweightModuleLabel")),
157  fEnu_Nbins (pset.get<int> ("Enu_Nbins")),
158  fEnu_min (pset.get<double>("Enu_Emin")),
159  fEnu_max (pset.get<double>("Enu_Emax")),
160  fVtxBinSize(pset.get<double>("VtxBinSize")),
161  fRockDepthX(pset.get<double>("RockDepthX")),
162  fRockDepthY(pset.get<double>("RockDepthY")),
163  fRockDepthZ(pset.get<double>("RockDepthZ")),
164  fNeutrinoInteractionVertexX(0), // Flag that these histograms are uninitialized
165  fPPFXWgt_Nbins (pset.get<int> ("PPFXWgt_Nbins")),
166  fPPFXWgt_min (pset.get<double>("PPFXWgt_Emin")),
167  fPPFXWgt_max (pset.get<double>("PPFXWgt_Emax"))
168  {
169  }
170 
171  //......................................................................
173  {
174  mf::LogInfo("NeutrinoAna") << "~NeutrinoAna" << std::endl;
175  }
176 
177  //......................................................................
179  {
181 
182  /// Helper variables for creating histograms
183  std::string energyName = "fEnu_";
184  std::string energyLabel = ";E_{#nu} (GeV);Events";
185 
186  // Create histograms for neutrino energy by flavor and current
187  std::string energyTitle = "Neutrino energy, by flavor/current" + energyLabel;
188  fEnu_nue_CC = tfs->make<TH1F>("fEnuByNu{group=Flavor+Current,cat=nue CC}", energyTitle.c_str(),
190  fEnu_nuebar_CC = tfs->make<TH1F>("fEnuByNu{group=Flavor+Current,cat=nuebar CC}", energyTitle.c_str(),
192  fEnu_numu_CC = tfs->make<TH1F>("fEnuByNu{group=Flavor+Current,cat=numu CC}", energyTitle.c_str(),
194  fEnu_numubar_CC = tfs->make<TH1F>("fEnuByNu{group=Flavor+Current,cat=numubar CC}", energyTitle.c_str(),
196  fEnu_nu_NC = tfs->make<TH1F>("fEnuByNu{group=Flavor+Current,cat=nu NC}", energyTitle.c_str(),
198  fEnu_nubar_NC = tfs->make<TH1F>("fEnuByNu{group=Flavor+Current,cat=nubar NC}", energyTitle.c_str(),
200 
201  // Helper maps for creating histograms of neutrino energy by parent
202  std::map<int, std::string> decayMode;
203  decayMode[1] = "Klong";
204  decayMode[5] = "Kcharged";
205  decayMode[11] = "Muon";
206  decayMode[13] = "Pion";
207 
208  std::map<int, std::string> decayLabel;
209  decayLabel[1] = "K_{L}";
210  decayLabel[5] = "K^{+/-}";
211  decayLabel[11] = "#mu^{+/-}";
212  decayLabel[13] = "#pi^{+/-}";
213 
214  // Create histograms of neutrino energy by parent
215  for(std::map<int, std::string>::iterator it = decayMode.begin(); it != decayMode.end(); ++it) {
216  std::string histTitle = "Neutrino energy, by flavor/current/parent particle;E_{#nu} (GeV)";
217 
218  std::string histName = energyName + "ByParent{group=Flavor+Current,cat=nu_e CC}{group=Parent,cat=" + it->second + "}";
219  fEnu_nue_ByParent[it->first] = tfs->make<TH1F>(histName.c_str(), histTitle.c_str(),
221 
222  histName = energyName + "ByParent{group=Flavor+Current,cat=nu_mu CC}{group=Parent,cat=" + it->second + "}";
223  fEnu_numu_ByParent[it->first] = tfs->make<TH1F>(histName.c_str(), histTitle.c_str(),
225 
226  histName = energyName + "ByParent{group=Flavor+Current,cat=other}{group=Parent,cat=" + it->second + "}";
227  fEnu_other_ByParent[it->first] = tfs->make<TH1F>(histName.c_str(), histTitle.c_str(),
229  }
230 
231  // Helper vector for creating histograms of neutrino energy by interaction mode
232  fModeBinIndices.push_back("QE");
233  fModeBinIndices.push_back("MEC");
234  fModeBinIndices.push_back("Res");
235  fModeBinIndices.push_back("DIS");
236  fModeBinIndices.push_back("Coh");
237  fModeBinIndices.push_back("Other");
238 
239  // Create histograms of neutrino energy by interaction mode
240  for(const auto & modePair : fModeNames) {
241  std::string histName = energyName + "ByInteraction{group=Interaction,cat=" + modePair.second + "}";
242  std::string histTitle = "Neutrino energy, by interaction type" + energyLabel;
243  fEnu_ByMode[modePair.first] = tfs->make<TH1F>(histName.c_str(), histTitle.c_str(),
245  }
246 
247  std::string hNameStart = "fNeutrinoStartingPosition";
248  std::string hTitleStart = "Neutrino Production ";
249  std::string hLabel = " (cm);Events";
250 
251  int nuPosNBins = 200;
252  double nuPosXY = 100.;
253  double nuPosZMin = -150.;
254  double nuPosZMax = 75000.;
255 
256  std::string histName = hNameStart + "X";
257  std::string histTitle = hTitleStart + "X Coordinate;x" + hLabel;
258  fNeutrinoStartingPositionX = tfs->make<TH1F>(histName.c_str(), histTitle.c_str(),
259  nuPosNBins, -nuPosXY, nuPosXY);
260 
261  histName = hNameStart + "Y";
262  histTitle = hTitleStart + "Y Coordinate;y" + hLabel;
263  fNeutrinoStartingPositionY = tfs->make<TH1F>(histName.c_str(), histTitle.c_str(),
264  nuPosNBins, -nuPosXY, nuPosXY);
265 
266  histName = hNameStart + "Z";
267  histTitle = hTitleStart + "Z Coordinate;z" + hLabel;
268  fNeutrinoStartingPositionZ = tfs->make<TH1F>(histName.c_str(), histTitle.c_str(),
269  nuPosNBins, nuPosZMin, nuPosZMax);
270 
271  histName = hNameStart + "XY";
272  histTitle = hTitleStart + "Coordinates, XY View;x (cm);y (cm)";
273  fNeutrinoStartingPositionXY = tfs->make<TH2F>(histName.c_str(), histTitle.c_str(),
274  nuPosNBins, -nuPosXY, nuPosXY,
275  nuPosNBins, -nuPosXY, nuPosXY);
276 
277  histName = hNameStart + "XZ";
278  histTitle = hTitleStart + "Coordinates, XZ View;z (cm);x (cm)";
279  fNeutrinoStartingPositionXZ = tfs->make<TH2F>(histName.c_str(), histTitle.c_str(),
280  nuPosNBins, nuPosZMin, nuPosZMax,
281  nuPosNBins, -nuPosXY, nuPosXY);
282 
283  histName = hNameStart + "YZ";
284  histTitle = hTitleStart + "Coordinates, YZ View;z (cm);y (cm)";
285  fNeutrinoStartingPositionYZ = tfs->make<TH2F>(histName.c_str(), histTitle.c_str(),
286  nuPosNBins, nuPosZMin, nuPosZMax,
287  nuPosNBins, -nuPosXY, nuPosXY);
288 
289  fNuCosX = tfs->make<TH1F>("fNuCosX", "Incoming #nu X Directron Cosine;dx/ds;Events", 200, -1., 1.);
290  fNuCosY = tfs->make<TH1F>("fNuCosY", "Incoming #nu Y Directron Cosine;dy/ds;Events", 200, -1., 1.);
291  fNuCosZ = tfs->make<TH1F>("fNuCosZ", "Incoming #nu Z Directron Cosine;dz/ds;Events", 200, -1., 1.);
292 
293  fMode = tfs->make<TH1F>("fMode", "Interaction mode;Interaction Mode;", fModeNames.size(), 0., fModeNames.size());
294  for(const auto & modePair : fModeNames) {
295  auto binIdx = std::distance(fModeBinIndices.begin(), std::find(fModeBinIndices.begin(), fModeBinIndices.end(), modePair.second));
296  fMode->GetXaxis()->SetBinLabel(binIdx+1, modePair.second.c_str()); // ROOT bin indices starts at 1, not 0
297  }
298  fMode->GetXaxis()->CenterLabels();
299 
300  fNeutrinoTrueL = tfs->make<TH1F>("fNeutrinoTrueL", "Neutrino L;L (m);Events",
301  201, -0.05, 2000.05);
302  fNeutrinoTrueETrueL = tfs->make<TH2F>("fNeutrinoTrueETrueL", "Neutrino Energy vs Length;E (GeV);L (m)",
303  1000, 0., 10., 10000, 100., 1100.);
304  fNeutrinoTrueLOverE = tfs->make<TH1F>("fNeutrinoTrueLOverE", "Neutrino L/E;L/E (m/GeV);Events",
305  2000, 0., 2000.);
306 
307  // Create histograms of proton multiplicity in map, using the PDG of nuclear target as the key
308  fProtonMultiplicity[1000010010] = tfs->make<TH1F>("fProtonMultiplicity{group=target,cat=H}", "Proton Multiplicity;Proton Multiplicity;", 15, 0., 15.);
309  fProtonMultiplicity[1000060120] = tfs->make<TH1F>("fProtonMultiplicity{group=target,cat=C12}", "Proton Multiplicity;Proton Multiplicity;", 15, 0., 15.);
310  fProtonMultiplicity[1000070140] = tfs->make<TH1F>("fProtonMultiplicity{group=target,cat=N14}", "Proton Multiplicity;Proton Multiplicity;", 15, 0., 15.);
311  fProtonMultiplicity[1000080160] = tfs->make<TH1F>("fProtonMultiplicity{group=target,cat=O16}", "Proton Multiplicity;Proton Multiplicity;", 15, 0., 15.);
312  fProtonMultiplicity[1000160320] = tfs->make<TH1F>("fProtonMultiplicity{group=target,cat=S32}", "Proton Multiplicity;Proton Multiplicity;", 15, 0., 15.);
313  fProtonMultiplicity[1000170350] = tfs->make<TH1F>("fProtonMultiplicity{group=target,cat=Cl35}", "Proton Multiplicity;Proton Multiplicity;", 15, 0., 15.);
314  fProtonMultiplicity[1000220480] = tfs->make<TH1F>("fProtonMultiplicity{group=target,cat=Ti48}", "Proton Multiplicity;Proton Multiplicity;", 15, 0., 15.);
315  fProtonMultiplicity[1000260560] = tfs->make<TH1F>("fProtonMultiplicity{group=target,cat=Fe56}", "Proton Multiplicity;Proton Multiplicity;", 15, 0., 15.);
316  fProtonMultiplicity[0] = tfs->make<TH1F>("fProtonMultiplicity{group=target,cat=Other}", "Proton Multiplicity;Proton Multiplicity;", 15, 0., 15.);
317 
318  fW = tfs->make<TH1F>("fW", "Invariant hadronic mass;W (GeV);Events", 150, 0.,3.);
319  fX = tfs->make<TH1F>("fX", "Bjorken X;x{Bj};Events", 100, 0., 1.);
320  fY = tfs->make<TH1F>("fY", "Bjorken Y;y_{Bj};Events", 100, 0., 1.);
321  fQSqr = tfs->make<TH1F>("fQSqr", "Q2;Q^{2} (GeV^{2});Events", 100, 0., 2.);
322 
323  fYUnscaled = tfs->make<TH1F>("nu", "Energy transfer;E_{#nu}-E_{lep} (GeV);Events", 500, 0., 10.);
324 
325  //PPFX weights:
326  std::string energyTitlePPFXWgt = "Neutrino energy by flavor weighted by PPFX" + energyLabel;
327  std::string ppfxwgtLabel = ";Weights;Events";
328  fPPFXWgt_cv = tfs->make<TH1F>("fPPFXWgt_cv", ppfxwgtLabel.c_str(),
330  fPPFXWgt_univ = tfs->make<TH1F>("fPPFXWgt_univ", ppfxwgtLabel.c_str(),
332  fPPFXWgtEnu_numu = tfs->make<TH1F>("fPPFXWgtEnu_numuEnuByNu{numu}", energyTitlePPFXWgt.c_str(),
334  fPPFXWgtEnu_numubar = tfs->make<TH1F>("fPPFXWgtEnu_numuEnuByNu{numubar}", energyTitlePPFXWgt.c_str(),
336  fPPFXWgtEnu_nue = tfs->make<TH1F>("fPPFXWgtEnu_numuEnuByNu{nue}", energyTitlePPFXWgt.c_str(),
338  fPPFXWgtEnu_nuebar = tfs->make<TH1F>("fPPFXWgtEnu_numuEnuByNu{nuebar}", energyTitlePPFXWgt.c_str(),
340 
341  }
342 
343  //......................................................................
345  {
346  if(fNeutrinoInteractionVertexX) { return; }
347 
350 
351  // First get detector coordinates
352  double detX = geo->DetHalfWidth(); // -detX is minimum, +detX is maximum value
353  double detY = geo->DetHalfHeight(); // Likewise, use +/- detY
354  double detZ = geo->DetLength(); // Here, mininmum value is 0, max value is detZ
355 
356  // Add on the space for rock events
357  double nuPosX = detX + fRockDepthX;
358  double nuPosY = detY + fRockDepthY;
359  double nuPosZMin = -fRockDepthZ;
360  double nuPosZMax = detZ + fRockDepthZ;
361 
362  // This section adds a little extra space to make bins sizes exact
363  int nbinsHelper = 0; // Helper to make bin sizes exact
364 
365  // First calculate number of complete bins from min to max values
366  // Add one more bin for the last partial bin
367  // For X/Y, half of those bins are above 0, half below
368  // The max value is then the bin size times then number of bins over 2
369  // For Z, calculate the adjusted number of bins in the same way
370  // Leave the min value unchanged
371  // Set the max value as the min value plus the number of bins times the bin size
372 
373  nbinsHelper = (int)(2.*nuPosX/fVtxBinSize);
374  ++nbinsHelper; // Add one extra bin to cover truncated area
375  nuPosX = 0.5*fVtxBinSize*(double)nbinsHelper;
376 
377  nbinsHelper = (int)(2.*nuPosY/fVtxBinSize);
378  ++nbinsHelper; // Add one extra bin to cover truncated area
379  nuPosY = 0.5*fVtxBinSize*(double)nbinsHelper;
380 
381  nbinsHelper = (int)((nuPosZMax - nuPosZMin)/fVtxBinSize);
382  ++nbinsHelper; // Add one extra bin to cover truncated area
383  nuPosZMax = nuPosZMin + (fVtxBinSize*(double)nbinsHelper);
384 
385  int nBinsX = 2.*nuPosX/fVtxBinSize;
386  int nBinsY = 2.*nuPosY/fVtxBinSize;
387  int nBinsZ = (nuPosZMax - nuPosZMin)/fVtxBinSize;
388 
389  std::string hNameStart = "fNeutrinoInteractionVertex";
390  std::string hTitleStart = "Neutrino Vertex ";
391  std::string hLabel = " (cm);Events";
392 
393  std::string histName = hNameStart + "X";
394  std::string histTitle = hTitleStart + "X Coordinate;x" + hLabel;
395  fNeutrinoInteractionVertexX = tfs->make<TH1F>(histName.c_str(), histTitle.c_str(),
396  nBinsX, -nuPosX, nuPosX);
397 
398  histName = hNameStart + "Y";
399  histTitle = hTitleStart + "Y Coordinate;y" + hLabel;
400  fNeutrinoInteractionVertexY = tfs->make<TH1F>(histName.c_str(), histTitle.c_str(),
401  nBinsY, -nuPosY, nuPosY);
402 
403  histName = hNameStart + "Z";
404  histTitle = hTitleStart + "Z Coordinate;z" + hLabel;
405  fNeutrinoInteractionVertexZ = tfs->make<TH1F>(histName.c_str(), histTitle.c_str(),
406  nBinsZ, nuPosZMin, nuPosZMax);
407 
408  histName = hNameStart + "XY";
409  histTitle = hTitleStart + "Coordinates, XY View;x (cm);y (cm)";
410  fNeutrinoInteractionVertexXY = tfs->make<TH2F>(histName.c_str(), histTitle.c_str(),
411  nBinsX, -nuPosX, nuPosX,
412  nBinsY, -nuPosY, nuPosY);
413 
414  histName = hNameStart + "XZ";
415  histTitle = hTitleStart + "Coordinates, XZ View;z (cm);x (cm)";
416  fNeutrinoInteractionVertexXZ = tfs->make<TH2F>(histName.c_str(), histTitle.c_str(),
417  nBinsZ, nuPosZMin, nuPosZMax,
418  nBinsX, -nuPosX, nuPosX);
419 
420  histName = hNameStart + "YZ";
421  histTitle = hTitleStart + "Coordinates, YZ View;z (cm);y (cm)";
422  fNeutrinoInteractionVertexYZ = tfs->make<TH2F>(histName.c_str(), histTitle.c_str(),
423  nBinsZ, nuPosZMin, nuPosZMax,
424  nBinsY, -nuPosY, nuPosY);
425 
426  /*histName = hNameStart;
427  histTitle = hTitleStart + "Coordinates;x (cm);y (cm);z (cm)";
428  hNeutrinoInteractionVertex = tfs->make<TH3F>(histName.c_str(), histTitle.c_str(),
429  nBinsX, -nuPosX, nuPosX,
430  nBinsY, -nuPosY, nuPosY,
431  nBinsZ, nuPosZMin, nuPosZMax);*/
432  }
433 
434  //......................................................................
436  {
437  // Get the MC generator information out of the event
438  bool isMCTruthListEmpty = false;
440  try {
441  evt.getByLabel(fMCTruthModuleLabel, mclist);
442  if(mclist->empty()) { isMCTruthListEmpty = true; }
443  }
444  catch(...) { isMCTruthListEmpty = true; }
445  if(isMCTruthListEmpty) {
446  mf::LogError("MCTruth Error") << "Error retrieving MCTruth list." << std::endl;
447  return;
448  }
449 
450  // Get the MCFlux information out of the event
451  bool isMCFluxListEmpty = false;
453  try {
454  evt.getByLabel(fMCTruthModuleLabel, fllist);
455  if(fllist->empty()) { isMCFluxListEmpty = true; }
456  }
457  catch(...) { isMCFluxListEmpty = true; }
458  if(isMCFluxListEmpty) {
459  if(msg1cnt < 5) {
460  mf::LogWarning("MCFlux Error") << "Error retrieving MCFlux list. Let it slide." << std::endl;
461  ++msg1cnt;
462  }
463  }
464 
465  if(mclist->size() != fllist->size()) {
466  mf::LogWarning("List sizes") << "MCTruth: " << mclist->size()
467  << ", MCFlux: " << fllist->size() << std::endl;
468  }
469 
470  //Get the flux weights from mc list:
472 
473  /// Loop over neutrino interactions
474  for(unsigned int i_intx = 0, n_intx = mclist->size(); i_intx < n_intx; ++i_intx) {
475  /// Pointers to current MCTruth and MCFlux
476  art::Ptr<simb::MCTruth> mc (mclist, i_intx);
477  art::Ptr<simb::MCFlux> mcflux(fllist, i_intx);
478 
479  // If there is no MCNeutrino, skip the event
480  if(!mc->NeutrinoSet()) {
481  if(msg1cnt < 5) {
482  mf::LogError("MCNeutrino") << "There is no MCNeutrino in this event." << std::endl;
483  ++msg1cnt;
484  }
485 
486  return;
487  }
488 
489  // Calculating true L, L/E from neutrino interaction
490  // This will return false for non detector events if FCL says to skip them
491  if(!AnalyzeNeutrinoInteraction(mc, mcflux)) { continue; }
492 
493  // Link to the MCNeutrino class
494  // The class contains information not only about the incoming neutrino,
495  // but also about the products of the decay
496  const simb::MCNeutrino& mc_neutrino = mc->GetNeutrino();
497  const simb::MCParticle& nu = mc_neutrino.Nu(); // Get the actual neutrino
498  const simb::MCParticle& lep = mc_neutrino.Lepton(); // Get the outgoing lepton
499 
500  // Neutrino starting positions
501  fNeutrinoStartingPositionX ->Fill(mcflux->fvx);
502  fNeutrinoStartingPositionY ->Fill(mcflux->fvy);
503  fNeutrinoStartingPositionZ ->Fill(mcflux->fvz);
504  fNeutrinoStartingPositionXY->Fill(mcflux->fvx, mcflux->fvy);
505  fNeutrinoStartingPositionXZ->Fill(mcflux->fvz, mcflux->fvx);
506  fNeutrinoStartingPositionYZ->Fill(mcflux->fvz, mcflux->fvy);
507 
508  // Neutrino interaction vertex
509  fNeutrinoInteractionVertexX ->Fill(nu.Vx());
510  fNeutrinoInteractionVertexY ->Fill(nu.Vy());
511  fNeutrinoInteractionVertexZ ->Fill(nu.Vz());
512  fNeutrinoInteractionVertexXY->Fill(nu.Vx(), nu.Vy());
513  fNeutrinoInteractionVertexXZ->Fill(nu.Vz(), nu.Vx());
514  fNeutrinoInteractionVertexYZ->Fill(nu.Vz(), nu.Vy());
515  //fNeutrinoInteractionVertex->Fill(nu.Vx(), nu.Vy(), nu.Vz());
516 
517  // Store these locally to repetitive function calls
518  int nuPDG = nu.PdgCode();
519  double nuEnergy = nu.E();
520 
521  // Fill nu energy histogram by flavor and current
522  if(mc_neutrino.CCNC() == simb::kCC) {
523  if (nuPDG == +14) { fEnu_numu_CC ->Fill(nuEnergy); }
524  else if(nuPDG == -14) { fEnu_numubar_CC->Fill(nuEnergy); }
525  else if(nuPDG == +12) { fEnu_nue_CC ->Fill(nuEnergy); }
526  else if(nuPDG == -12) { fEnu_nuebar_CC ->Fill(nuEnergy); }
527  else {
528  mf::LogWarning("Unknown Flavor") << "Unknown neutrino flavor when trying to fill Enu CC histograms." << std::endl;
529  }
530  }
531  else {
532  if(nuPDG > 0) { fEnu_nu_NC ->Fill(nuEnergy); }
533  else { fEnu_nubar_NC->Fill(nuEnergy); }
534  }
535 
536  // Find the decay mode (corresponding to parent type)
537  // Then set the value to a corresponding map key
538  int decayMode = mcflux->fndecay;
539  if (decayMode >= 1 && decayMode <= 4) { decayMode = 1; }
540  else if(decayMode >= 5 && decayMode <= 10) { decayMode = 5; }
541  else if(decayMode >= 11 && decayMode <= 12) { decayMode = 11; }
542  else if(decayMode >= 13 && decayMode <= 14) { decayMode = 13; }
543  else { decayMode = 0; }
544 
545  // Fill nu energy by parent species histograms
546  if (decayMode != 0) {
547  if (abs(nuPDG) == 12)
548  fEnu_nue_ByParent[decayMode] ->Fill(nuEnergy);
549  else if (abs(nuPDG) == 14)
550  fEnu_numu_ByParent[decayMode]->Fill(nuEnergy);
551  else
552  fEnu_other_ByParent[decayMode]->Fill(nuEnergy);
553  }
554 
555  // Fill nu energy by interaction mode histograms
556  int nuMode = mc_neutrino.Mode();
557  if(fEnu_ByMode.find(nuMode) == fEnu_ByMode.end())
558  nuMode = simb::kUnknownInteraction;
559  fEnu_ByMode[nuMode]->Fill(nuEnergy);
560  int binIdx = std::distance(fModeBinIndices.begin(),
561  std::find(fModeBinIndices.begin(),
562  fModeBinIndices.end(),
563  fModeNames.at(static_cast<simb::int_type_>(nuMode))
564  )
565  );
566  fMode->Fill(binIdx);
567 
568  // Fill neutrino direction cosine histograms
569  if(fabs(nu.P()) > 0.) {
570  fNuCosX->Fill(nu.Px()/nu.P());
571  fNuCosY->Fill(nu.Py()/nu.P());
572  fNuCosZ->Fill(nu.Pz()/nu.P());
573  } // end of checking whether P is nonzero
574 
575  // Proton multiplicity -- Look for protons in the particle stack
576  int nprotons = 0; // number of protons
577 
578  // Loop over particles produced in the event
579  // Get number of protons produced
580  for(int i_particle = 0, n_particle = mc->NParticles(); i_particle < n_particle; ++i_particle) {
581  const simb::MCParticle& part = mc->GetParticle(i_particle);
582  if(part.PdgCode() == 2212) {
583  ++nprotons;
584  }
585  }
586 
587  // Fill histograms of proton multiplicity
588  int targetPDG = mc_neutrino.Target();
589  if(fProtonMultiplicity.find(targetPDG) == fProtonMultiplicity.end())
590  targetPDG = 0;
591  fProtonMultiplicity[targetPDG]->Fill(nprotons);
592 
593  fW->Fill(mc_neutrino.W());
594  fX->Fill(mc_neutrino.X());
595  fY->Fill(mc_neutrino.Y());
596  fQSqr->Fill(mc_neutrino.QSqr());
597  fYUnscaled->Fill(nuEnergy - lep.E());
598 
599  //PPFX:
600  std::vector<art::Ptr<fxwgt::FluxWeights> > vfluxwgts = fmFLXRW.at(i_intx);
601  if(vfluxwgts.size()>0){
602 
603  //Weights:
604  std::vector< float > vwgts_univ = vfluxwgts[0]->GetTotalMultiUniverses();
605  fPPFXWgt_cv->Fill(vfluxwgts[0]->GetTotalCentralValue());
606  for(unsigned int iuniv = 0;iuniv < vfluxwgts[0]->GetNumberOfUniverses();iuniv++){
607  fPPFXWgt_univ->Fill(vwgts_univ[iuniv]);
608  }
609  //Weighted spectra: nuPDG
610  if(nuPDG == 14)fPPFXWgtEnu_numu->Fill( nuEnergy,vfluxwgts[0]->GetTotalCentralValue());
611  if(nuPDG ==-14)fPPFXWgtEnu_numubar->Fill(nuEnergy,vfluxwgts[0]->GetTotalCentralValue());
612  if(nuPDG == 12)fPPFXWgtEnu_nue->Fill( nuEnergy,vfluxwgts[0]->GetTotalCentralValue());
613  if(nuPDG ==-12)fPPFXWgtEnu_nuebar->Fill( nuEnergy,vfluxwgts[0]->GetTotalCentralValue());
614  }
615  } // end loop over interactions
616 
617  }
618 
619  //......................................................................
621  const art::Ptr<simb::MCFlux>& mcflux)
622  {
624 
625  // Link to the MCNeutrino class
626  // The class contains information not only about the incoming neutrino,
627  // but also about the products of the decay
628  const simb::MCNeutrino& mc_neutrino = mc->GetNeutrino();
629  //const simb::MCParticle& incoming_nu = mc_neutrino.Nu(); // Incoming neutrino
630  const simb::MCParticle& outgoing_lepton = mc_neutrino.Lepton(); // Outgoing lepton
631 
632  double detHalfWidth = geo->DetHalfWidth();
633  double detHalfHeight = geo->DetHalfHeight();
634  double detLength = geo->DetLength();
635 
636  if(fDetOnly) {
637  const simb::MCParticle& nu = mc_neutrino.Nu();
638  if(std::abs(nu.Vx()) > detHalfWidth ||
639  std::abs(nu.Vy()) > detHalfHeight ||
640  nu.Vz() < 0. || nu.Vz() > detLength) {
641  return false;
642  }
643  }
644 
645  // Neutrino origin vertex (parent decay point) in beam coordinates
646  // The CLHEP::cm is required for the coordinate transformation
647  double neutrinoTrueStart_BeamCoords[3];
648  neutrinoTrueStart_BeamCoords[0] = mcflux->fvx * CLHEP::cm;
649  neutrinoTrueStart_BeamCoords[1] = mcflux->fvy * CLHEP::cm;
650  neutrinoTrueStart_BeamCoords[2] = mcflux->fvz * CLHEP::cm;
651 
652  // Neutrino interaction vertex in detector coordinates
653  double neutrinoTrueEnd_DetectorCoords[3];
654  neutrinoTrueEnd_DetectorCoords[0] = outgoing_lepton.Vx();
655  neutrinoTrueEnd_DetectorCoords[1] = outgoing_lepton.Vy();
656  neutrinoTrueEnd_DetectorCoords[2] = outgoing_lepton.Vz();
657 
658  // Neutrino origin vertex in detector coordinates
659  double neutrinoTrueStart_DetectorCoords[3];
660 
661  geo::CoordinateTransformation* coordinateTransformation = geo->getCoordinateTransformation();
662 
663  // Try to get start point in detector coordinates
664  // The output will have dimensions of CLHEP::cm
665  if(!coordinateTransformation->transformBeamToDetectorCoordinates(neutrinoTrueStart_BeamCoords,
666  neutrinoTrueStart_DetectorCoords)) {
667  return false;
668  }
669 
670  // Get neutrino L components
671  double vecNeutrinoL[3];
672  for(unsigned int i = 0, n = 3; i < n; ++i) {
673  // Calculate the distance traveled by the neutrino
674  vecNeutrinoL[i] = neutrinoTrueEnd_DetectorCoords[i] - neutrinoTrueStart_DetectorCoords[i]/CLHEP::cm;
675  vecNeutrinoL[i] /= 100.; // L is typically given in meters, not cm (especially for L/E)
676  }
677 
678  // Get neutrino L and E
679  double neutrinoL = sqrt(DotProduct(vecNeutrinoL, vecNeutrinoL));
680  double neutrinoE = mc_neutrino.Nu().E();
681 
682  // Fill histograms related to neutrino L, E
683  fNeutrinoTrueL->Fill(neutrinoL);
684  fNeutrinoTrueETrueL->Fill(neutrinoE, neutrinoL);
685  fNeutrinoTrueLOverE->Fill(neutrinoL/neutrinoE);
686 
687  return true;
688  }
689 
690  //......................................................................
691  /// Calculate dot product between two vectors
692  double NeutrinoAna::DotProduct(double* v1, double* v2)
693  {
694  double product = 0.0;
695 
696  for(int i = 0; i < 3; ++i) {
697  product += v1[i] * v2[i];
698  }
699 
700  return product;
701  }
702 } // end of namespace mcchk
703 
704 ///////////////////////////////////////////////////////////////////////////////
705 namespace mcchk {
707 }
TH1F * fNeutrinoStartingPositionZ
Neutrino Production Point Z Coordinate.
double E(const int i=0) const
Definition: MCParticle.h:232
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
double QSqr() const
Definition: MCNeutrino.h:157
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
double fVtxBinSize
Size of bins in cm for neutrino interaction vertex histograms.
double Py(const int i=0) const
Definition: MCParticle.h:230
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TH1F * fNeutrinoStartingPositionY
Neutrino Production Point Y Coordinate.
Simple module to analyze MC cosmics distributions.
TH1F * fPPFXWgt_cv
Central value weights.
TH1F * fNeutrinoInteractionVertexX
Neutrino Interaction Vertex X Coordinate.
TH2F * fNeutrinoInteractionVertexXZ
Neutrino Interaction Vertex X vs Z Coordinates.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
TH2F * fNeutrinoTrueETrueL
True L and E of the neutrino.
bool transformBeamToDetectorCoordinates(double *input_beam_coords, double *output_detector_coords) const
TH1F * fYUnscaled
Difference between neutrino and charged lepton energies.
TH1F * fNuCosZ
Z direction cosine of incoming neutrino.
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1F * fW
W (shower mass)
double Px(const int i=0) const
Definition: MCParticle.h:229
double fvx
Definition: MCFlux.h:52
TH1F * fNeutrinoInteractionVertexY
Neutrino Interaction Vertex Y Coordinate.
TH1F * fEnu_nue_CC
Neutrino information.
bool AnalyzeNeutrinoInteraction(const art::Ptr< simb::MCTruth > &mc, const art::Ptr< simb::MCFlux > &mcflux)
Calculate true L, E, L/E from neutrino interaction information.
double DetLength() const
std::vector< std::string > fModeBinIndices
void abs(TH1 *hist)
TH1F * fEnu_nu_NC
Neutrino energy for neutrino NC events.
TH1F * fNeutrinoInteractionVertexZ
Neutrino Interaction Vertex Z Coordinate.
TH1F * fNeutrinoStartingPositionX
Neutrino Production Point X Coordinate.
int NParticles() const
Definition: MCTruth.h:74
double fEnu_max
Maximum energy for energy histgorams.
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
DEFINE_ART_MODULE(TestTMapFile)
unsigned distance(const T &t1, const T &t2)
double fRockDepthX
Amount of space in cm outside of detector along x to plot neutrino interaction points.
float abs(float number)
Definition: d0nt_math.hpp:39
CoordinateTransformation * getCoordinateTransformation() const
Definition: GeometryBase.h:245
TH1F * fY
Bjorken Y.
TH1F * fPPFXWgtEnu_numubar
Neutrino energy for anti-numu events weighted by PPFX.
std::string fMCTruthModuleLabel
label for module producing mc truth information
object containing MC flux information
Definition: Run.h:31
TH1F * fPPFXWgtEnu_numu
Neutrino energy for numu events weighted by PPFX.
static constexpr double cm
Definition: SystemOfUnits.h:99
NeutrinoAna(fhicl::ParameterSet const &pset)
TString part[npart]
Definition: Style.C:32
const std::map< simb::int_type_, std::string > fModeNames
int fPPFXWgt_Nbins
Number of bins for PPFX weight histograms.
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
double W() const
Definition: MCNeutrino.h:154
double Y() const
Definition: MCNeutrino.h:156
int fndecay
Definition: MCFlux.h:50
TH1F * fMode
Interaction mode.
void analyze(art::Event const &evt)
std::map< int, TH1F * > fEnu_other_ByParent
Neutrino energy by parent species, NC/nu_tau(bar)s.
double fPPFXWgt_min
Minimum weight PPFX weight histograms.
int fEnu_Nbins
Number of bins for energy histograms.
double P(const int i=0) const
Definition: MCParticle.h:233
int evt
TH1F * fEnu_numu_CC
Neutrino energy for numu CC events.
double fPPFXWgt_max
Maximum weight PPFX weight histograms.
TH2F * fNeutrinoInteractionVertexYZ
Neutrino Interaction Vertex Y vs Z Coordinates.
double X() const
Definition: MCNeutrino.h:155
std::string fFluxReweightModuleLabel
label for module producing flux reweight information
double DotProduct(double *, double *)
Calculate dot product between two vectors.
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
void beginRun(art::Run const &run)
TH2F * fNeutrinoStartingPositionXY
Neutrino Production Point X vs Y Coordinates.
Definition: run.py:1
TH1F * fPPFXWgtEnu_nue
Neutrino energy for nue events weighted by PPFX.
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
double fEnu_min
Minimum energy for energy histograms.
TH1F * fNeutrinoTrueLOverE
True L/E of the neutrino.
TH1F * fNuCosY
Y direction cosine of incoming neutrino.
double fvy
Definition: MCFlux.h:53
std::map< int, TH1F * > fEnu_ByMode
Neutrino energy by interaction mode.
double Vx(const int i=0) const
Definition: MCParticle.h:220
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TH2F * fNeutrinoStartingPositionXZ
Neutrino Production Point X vs Z Coordinates.
T * make(ARGS...args) const
n.b.: this group is similar but not quite, entirely unlike GENIE ScatteringType convention ...
Definition: MCNeutrino.h:80
int Target() const
Definition: MCNeutrino.h:151
double DetHalfWidth() const
T product(std::vector< T > dims)
TH1F * fPPFXWgt_univ
Multi-universe weights.
std::map< int, TH1F * > fEnu_numu_ByParent
Neutrino energy by parent species, nu_mu(bar)s.
TH1F * fEnu_nuebar_CC
Nuetrino energy for anti-nue CC events.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
TH2F * fNeutrinoStartingPositionYZ
Neutrino Production Point Y vs Z Coordinates.
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
double fRockDepthY
Amount of space in cm outside of detector along y to plot neutrino interaction points.
double Pz(const int i=0) const
Definition: MCParticle.h:231
TH1F * fX
Bjorken X.
double Vz(const int i=0) const
Definition: MCParticle.h:222
std::map< int, TH1F * > fProtonMultiplicity
Proton multiplicity by nucleus PDG.
double fRockDepthZ
Amount of space in cm outside of detector along z to plot neutrino interaction points.
bool NeutrinoSet() const
Definition: MCTruth.h:77
double fvz
Definition: MCFlux.h:54
TH1F * fNeutrinoTrueL
True L for neutrino.
std::map< int, TH1F * > fEnu_nue_ByParent
Neutrino energy by parent species, nu_e(bar)s.
static int msg1cnt
Event generator information.
Definition: MCNeutrino.h:18
Helper for AttenCurve.
Definition: Path.h:10
int Mode() const
Definition: MCNeutrino.h:149
TH1F * fPPFXWgtEnu_nuebar
Neutrino energy for nuebar events weighted by PPFX.
TH1F * fEnu_numubar_CC
Neutrino energy for anti-numu CC events.
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)
TH2F * fNeutrinoInteractionVertexXY
Neutrino Interaction Vertex X vs Y Coordinates.
TH1F * fEnu_nubar_NC
Neutrino energy for anti-nu NC events.
TH1F * fNuCosX
X direction cosine of incoming neutrino.