All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CRTSimAnalysis_module.cc
Go to the documentation of this file.
1 /**
2  * @file CRTSimAnalysis_module.cc
3  * @brief Access CRT data and reco products and compare to MCTruth info
4  * @author Chris Hilgenberg (Chris.Hilgenberg@colostate.edu)
5  *
6  * The last revision of this code was done in October 2018 with LArSoft v07_06_01.
7  */
8 
9 // LArSoft includes
10 #include "lardataobj/Simulation/SimChannel.h"
11 #include "lardataobj/Simulation/AuxDetSimChannel.h"
12 #include "lardataobj/RecoBase/Hit.h"
13 #include "lardataobj/RecoBase/Cluster.h"
14 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
15 #include "larcore/Geometry/Geometry.h"
16 #include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom()
17 #include "larcorealg/Geometry/GeometryCore.h"
18 #include "larcorealg/Geometry/AuxDetGeometryCore.h"
19 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
20 #include "nusimdata/SimulationBase/MCParticle.h"
21 #include "nusimdata/SimulationBase/MCTruth.h"
22 #include "larsim/Simulation/LArG4Parameters.h"
23 
24 // Framework includes
25 #include "art/Framework/Core/EDAnalyzer.h"
26 #include "art/Framework/Principal/Event.h"
27 #include "art/Framework/Principal/Handle.h"
28 #include "art/Framework/Services/Registry/ServiceHandle.h"
29 #include "art_root_io/TFileService.h"
30 #include "art/Framework/Core/ModuleMacros.h"
31 #include "canvas/Persistency/Common/FindManyP.h"
32 #include "canvas/Utilities/Exception.h"
33 
34 // Utility libraries
35 #include "messagefacility/MessageLogger/MessageLogger.h"
36 #include "fhiclcpp/ParameterSet.h"
37 #include "fhiclcpp/types/Table.h"
38 #include "fhiclcpp/types/Atom.h"
39 #include "cetlib/pow.h" // cet::sum_of_squares()
40 
41 // ROOT includes
42 #include "TH1.h"
43 #include "TH2.h"
44 #include "TTree.h"
45 #include "TLorentzVector.h"
46 #include "TVector3.h"
47 #include "TGeoManager.h"
48 #include "TMath.h"
49 #include "TROOT.h"
50 
51 // C++ includes
52 #include <map>
53 #include <vector>
54 #include <string>
55 #include <set>
56 #include <cmath>
57 #include <iostream>
58 #include <utility>
59 #include <array>
60 
61 // CRT data products
62 #include "sbnobj/ICARUS/CRT/CRTData.hh"
63 #include "sbnobj/Common/CRT/CRTHit.hh"
64 #include "sbnobj/Common/CRT/CRTTrack.hh"
65 #include "icaruscode/CRT/CRTUtils/CRTCommonUtils.h"
66 #include "icaruscode/CRT/CRTUtils/CRTBackTracker.h"
67 
68 using std::string;
69 using std::vector;
70 using std::map;
71 using std::set;
72 using std::pair;
73 
74 namespace {
75  int ProcessToICode(string const& p);
76 
77 } // local namespace
78 
79 
80 namespace icarus {
81 namespace crt {
82 
83  //-----------------------------------------------------------------------
84  //-----------------------------------------------------------------------
85  // class definition
86  /**
87  * @brief Example analyzer
88  *
89  * This class extracts information from the generated and reconstructed
90  * particles.
91  *
92  * It produces histograms for the simulated particles in the input file:
93  * - PDG ID (flavor) of all particles
94  * - momentum of the primary particles selected to have a specific PDG ID
95  * - length of the selected particle trajectory
96  *
97  * It also produces two ROOT trees.
98  *
99  * The first ROOT tree contains information on the simulated
100  * particles, including "dEdx", a binned histogram of collected
101  * charge as function of track range.
102  *
103  * Configuration parameters
104  * =========================
105  *
106  * - *SimulationLabel* (string, default: "largeant"): tag of the input data
107  * product with the detector simulation information (typically an instance
108  * of the LArG4 module)
109  *
110  */
111  class CRTSimAnalysis : public art::EDAnalyzer
112  {
113  public:
114 
115  using CRTHit = sbn::crt::CRTHit;
116 
117  // -------------------------------------------------------------------
118  // -------------------------------------------------------------------
119  // Standard constructor for an ART module with configuration validation;
120  // we don't need a special destructor here.
121 
122  /// Constructor: configures the module (see the Config structure above)
123  explicit CRTSimAnalysis(fhicl::ParameterSet const& p);
124 
125  virtual void beginJob() override;
126  virtual void beginRun(const art::Run& run) override;
127  virtual void analyze (const art::Event& event) override;
128 
129  private:
130 
131  // The parameters we'll read from the .fcl file.
132  art::InputTag fSimulationProducerLabel; ///< The name of the producer that tracked simulated particles through the detector
133  art::InputTag fCRTSimHitProducerLabel; ///< The name of the producer that created hits
134  art::InputTag fCRTTrueHitProducerLabel;
135  art::InputTag fCRTDetSimProducerLabel;
136  art::InputTag fCRTSimTrackProducerLabel;
137  vector<int> fPDGs; ///< PDG code of particle we'll focus on
138  vector<float> fMinMomenta;
139  vector<float> fMaxMomenta;
140 
141  // The n-tuples we'll create.
142  TTree* fCosmicDisplayNtuple; ///< for ROOT based event display
143  TTree* fGenNtuple;
144  TTree* fSimulationNtuple; ///< tuple with simulated data
145  TTree* fRegionsNtuple;
146  TTree* fDetSimNtuple;
147  TTree* fSimHitNtuple;
148  TTree* fTrueCRTHitNtuple;
149  TTree* fSimTrackNtuple;
150 
151  // The comment lines with the @ symbols define groups in doxygen.
152  /// @name The variables that will go into both n-tuples.
153  /// @{
154  int fEvent; ///< number of the event being processed
155  int fRun; ///< number of the run being processed
156  int fSubRun; ///< number of the sub-run being processed
157  /// @}
158 
159  /// @name The variables that will go into the CosmicDisplay n-tuple.
160  /// @{
161  static const int LAR_PROP_DELAY = 1.0/(30.0/1.38); //[ns/cm]
162  int fCDTrackID;
163  int fNCD;
164  int fCDpdg;
165  vector<vector<double>> fCDSlopes; ///< direction cosines
166  vector<vector<double>> fCDpe; ///< 4-momentum
167  vector<vector<double>> fCDxyzt; ///< 4-position
168  /// @}
169 
170  /// @name The variables that will go into the Gen n-tuple.
171  /// @{
172  int fNGen;
173  vector<int> fGenTrack;
174  vector<int> fGenPDG;
175  vector<vector<double>> fGenStartXYZT;
176  vector<vector<double>> fGenEndXYZT;
177  vector<vector<double>> fGenStartPE;
178  vector<vector<double>> fGenEndPE;
179  /// @}
180 
181  /// @name The variables that will go into the Simulation n-tuple.
182  /// @{
183  uint32_t fSimHits; ///< number of trajectory points for each MCParticle
184  float fTrackLength; ///< total track length for each MCParticle
185  int fSimPDG; ///< PDG ID of the particle being processed
186  int fSimProcess; ///< process that created the particle (e.g. brehmstralung)
187  int fSimEndProcess; ///< process the killed the particle (e.g. annihilation)
188  int fSimTrackID; ///< GEANT ID of the particle being processed
189  uint32_t fNAuxDet; ///< Number of scintillator strips hit
190 
191  vector<uint32_t> fAuxDetID; ///< Global CRT module ID
192  vector<uint32_t> fAuxDetSensitiveID; ///< Strip ID in module
193  vector<double> fADEDep; ///< Energy deposited in CRT strip (GeV)
194  vector<double> fADdEdx; ///< average dEdx for particle traversing CRT strip
195  vector<double> fADTrackLength; ///< Track length in CRT strip (cm)
196  vector<uint32_t> fAuxDetReg; ///< CRT region code
197  vector<uint32_t> fADMac; ///< Mac5 address of the CRT module
198  vector<uint32_t> fADType;
199  vector<vector<double>> fADEnterXYZT; ///< 4-position of entry into CRT strip
200  vector<vector<double>> fADExitXYZT; ///< 4-position of exit from CRT strip
201  vector<vector<double>> fADEnterPE; ///< 4-position of entry into CRT strip
202  vector<vector<double>> fADExitPE; ///< 4-position of exit from CRT strip
203 
204  int fParentPDG;
205  float fParentE;
206 
207  // Arrays for 4-vectors: (x,y,z,t) and (Px,Py,Pz,E).
208  float fStartXYZT[4];///< (x,y,z,t) of the true start of the particle
209  float fEndXYZT[4]; ///< (x,y,z,t) of the true end of the particle
210  float fStartPE[4]; ///< (Px,Py,Pz,E) at the true start of the particle
211  float fEndPE[4]; ///< (Px,Py,Pz,E) at the true end of the particle
212 
213  int fProgenitor; ///< G4 track ID of the primary particle that ultimately led to this one
214  int fMother; ///< G4 track ID of mother that directly produced this MCParticle
215  int fNDaught; ///< number of daughters belonging to this MCParticle
216 
217  //Regions tree vars (1 track per entry)
218  int fNReg;
219  int fRegFid;
220  int fRegActive;
221  int fRegInactive;
222  int fRegCRTs;
223  int fRegTrkID;
224  int fRegPDG;
225  vector<int> fRegRegions;
226  vector<double> fRegEDep;
227  vector<double> fRegdL;
228  vector<vector<double>> fRegEntryPE;
229  vector<vector<double>> fRegExitPE;
230  vector<vector<double>> fRegEntryXYZT;
231  vector<vector<double>> fRegExitXYZT;
232  vector<vector<double>> fRegEntrySlope;
233  vector<vector<double>> fRegExitSlope;
234  vector<int> fRegOpDetID;
235  vector<double> fRegDistToOpDet;
236  vector<vector<double>> fRegOpDetXYZT;
237 
238  //CRT data product vars
239  //int fNChan; ///< number of channels above threshold for this front-end board readout
240  int fEntry; ///< front-end board entry number (reset for each event)
241  int fFEBReg; ///< CRT region for this front-end board
242  int fMac5; ///< Mac5 address for this front-end board
243  int fDetSubSys;
244  int fT0;///< signal time w.r.t. global event time
245  int fT1;///< signal time w.r.t. PPS
246  int fADC[64];///< signal amplitude
247  int fMaxAdc;
248  int fMaxChan;
249  int fNAbove;
250  vector<int> fTrackID;///< track ID(s) of particle that produced the signal
251  vector<int> fDetPDG; /// signal inducing particle(s)' PDG code
252 
253  // sim CRT hit product vars
254  float fXHit; ///< reconstructed X position of CRT hit (cm)
255  float fYHit; ///< reconstructed Y position of CRT hit (cm)
256  float fZHit; ///< reconstructed Z position of CRT hit (cm)
257  float fXHitErr; ///< stat error of CRT hit reco X (cm)
258  float fYHitErr; ///< stat error of CRT hit reco Y (cm)
259  float fZHitErr; ///< stat error of CRT hit reco Z (cm)
260  float fT0Hit; ///< hit time w.r.t. global event time
261  float fT1Hit; ///< hit time w.r.t. PPS
262  int fHitReg; ///< region code of CRT hit
263  int fHitSubSys;
264  int fNHit; ///< number of CRT hits for this event
265  vector<int> fHitTrk;
266  vector<int> fHitPDG;
267  vector<int> fHitMod;
268  vector<int> fHitStrip;
269  int fNHitFeb;
270  vector<float> fHitPe;
271  float fHitTotPe;
272  float fHitPeRms;
273 
274  // truth CRT hit vars
275  float fTrueXHit; ///< reconstructed X position of CRT hit (cm)
276  float fTrueYHit; ///< reconstructed Y position of CRT hit (cm)
277  float fTrueZHit; ///< reconstructed Z position of CRT hit (cm)
278  float fTrueXHitErr; ///< stat error of CRT hit reco X (cm)
279  float fTrueYHitErr; ///< stat error of CRT hit reco Y (cm)
280  float fTrueZHitErr; ///< stat error of CRT hit reco Z (cm)
281  float fTrueT0Hit; ///< hit time w.r.t. global event time
282  float fTrueT1Hit; ///< hit time w.r.t. global event time
283  int fTrueHitReg; ///< region code of CRT hit
284  int fTrueHitSubSys;
285  uint32_t fTrueNHit; ///< number of CRT hits for this event
286  vector<int> fTrueHitTrk;
287  vector<int> fTrueHitPDG;
288  vector<int> fTrueHitStrip;
289  vector<int> fTrueHitMod;
290  int fTrueNHitFeb;
291  vector<float> fTrueHitPe;
292  float fTrueHitTotPe;
293  float fTrueHitPeRms;
294 
295  //CRTSimTrack vars
296  int fNSimTrack; ///< number of simulated CRT tracks for this event
297  float fSimTrackPE;
298  double fSimTrackT0;
299  float fSimTrackStart[3];
300  float fSimTrackEnd[3];
301  float fSimTrackL;
302  float fSimTrackTheta;
303  float fSimTrackPhi;
304  int fNHitSimTrack;
305  float fSimTrackHitStart[4];
306  float fSimTrackHitEnd[4];
307  int fSimTrackRegStart;
308  int fSimTrackRegEnd;
309 
310  /// @}
311 
312  // Other variables that will be shared between different methods.
313  geo::GeometryCore const* fGeometryService; ///< pointer to Geometry provider
314  int fTriggerOffset; ///< (units of ticks) time of expected neutrino event
315  CRTBackTracker bt;
316  CRTCommonUtils* fCrtutils;
317 
318  }; // class CRTSimAnalysis
319 
320 
321  //-----------------------------------------------------------------------
322  //-----------------------------------------------------------------------
323  // class implementation
324 
325  //-----------------------------------------------------------------------
326  // Constructor
327  //
328  // Note that config is a Table<Config>, and to access the Config
329  // value we need to use an operator: "config()". In the same way,
330  // each element in Config is an Atom<Type>, so to access the type we
331  // again use the call operator, e.g. "SimulationLabel()".
332  CRTSimAnalysis::CRTSimAnalysis(fhicl::ParameterSet const& p)
333  : EDAnalyzer{p}
334  , fSimulationProducerLabel(p.get<art::InputTag>("SimulationLabel","largeant"))
335  , fCRTSimHitProducerLabel(p.get<art::InputTag>("CRTSimHitLabel","crthit"))
336  , fCRTTrueHitProducerLabel(p.get<art::InputTag>("CRTTrueHitLabel","crttruehit"))
337  , fCRTDetSimProducerLabel(p.get<art::InputTag>("CRTDetSimLabel","crtdaq"))
338  , fCRTSimTrackProducerLabel(p.get<art::InputTag>("CRTSimTrackLabel","crttrack"))
339  , fPDGs(p.get<vector<int>>("PDGs"))
340  , fMinMomenta(p.get<vector<float>>("MinMomenta"))
341  , fMaxMomenta(p.get<vector<float>>("MaxMomenta"))
342  , bt(p.get<fhicl::ParameterSet>("CRTBackTrack"))
343  , fCrtutils(new CRTCommonUtils())
344  {
345  fGeometryService = lar::providerFrom<geo::Geometry>();
346  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
347  fTriggerOffset = sampling_rate(clockData);
348  }
349 
350  //-----------------------------------------------------------------------
351  void CRTSimAnalysis::beginJob()
352  {
353  std::cout << " starting analysis job" << std::endl;
354 
355  // Access ART's TFileService, which will handle creating and writing
356  // histograms and n-tuples for us.
357  art::ServiceHandle<art::TFileService> tfs;
358 
359  // Define our n-tuples
360  fCosmicDisplayNtuple = tfs->make<TTree>("DisplayTree", "track information for ROOT event display");
361  fGenNtuple = tfs->make<TTree>("GenTree", "truth information from the generator");
362  fSimulationNtuple = tfs->make<TTree>("SimTree", "MyCRTSimulation");
363  fRegionsNtuple = tfs->make<TTree>("RegTree", "Info about particles crossing boundaries");
364  fDetSimNtuple = tfs->make<TTree>("DetTree", "MyCRTDetSim");
365  fSimHitNtuple = tfs->make<TTree>("HitTree", "MyCRTSimHit");
366  fTrueCRTHitNtuple = tfs->make<TTree>("TrueCRTHitTree", "CRT hits from truth info");
367  fSimTrackNtuple = tfs->make<TTree>("SimTrackTree", "Simulated CRTTracks");
368 
369  // Define the branches of our event display n-tuple
370  fCosmicDisplayNtuple->Branch("event", &fEvent, "event/I");
371  fCosmicDisplayNtuple->Branch("run", &fRun, "run/I");
372  fCosmicDisplayNtuple->Branch("subRun", &fSubRun, "subRun/I");
373  fCosmicDisplayNtuple->Branch("trackID", &fCDTrackID, "trackID/I");
374  fCosmicDisplayNtuple->Branch("nSeg", &fNCD, "nSeg/I");
375  fCosmicDisplayNtuple->Branch("pdg", &fCDpdg, "pdg/I");
376  //fCosmicDisplayNtuple->Branch("regions", &fCDRegions);
377  fCosmicDisplayNtuple->Branch("slopes", &fCDSlopes);
378  fCosmicDisplayNtuple->Branch("pe", &fCDpe);
379  fCosmicDisplayNtuple->Branch("xyzt", &fCDxyzt);
380 
381  // Define the branches of our Gen n-tuple
382  fGenNtuple->Branch("event", &fEvent, "event/I");
383  fGenNtuple->Branch("run", &fRun, "run/I");
384  fGenNtuple->Branch("subRun", &fSubRun, "subRun/I");
385  fGenNtuple->Branch("nGen", &fNGen, "nGen/I");
386  fGenNtuple->Branch("trackID", &fGenTrack);
387  fGenNtuple->Branch("pdg", &fGenPDG);
388  fGenNtuple->Branch("startXYZT", &fGenStartXYZT);
389  fGenNtuple->Branch("endXYZT", &fGenEndXYZT);
390  fGenNtuple->Branch("startPE", &fGenStartPE);
391  fGenNtuple->Branch("endPE", &fGenEndPE);
392 
393  // Define the branches of our simulation n-tuple
394  fSimulationNtuple->Branch("event", &fEvent, "event/I");
395  fSimulationNtuple->Branch("run", &fRun, "run/I");
396  fSimulationNtuple->Branch("subRun", &fSubRun, "subrun/I");
397  fSimulationNtuple->Branch("nPoints" , &fSimHits, "nPoints/I");
398  fSimulationNtuple->Branch("trackID", &fSimTrackID, "trackID/I");
399  fSimulationNtuple->Branch("pdg", &fSimPDG, "pdg/I");
400  fSimulationNtuple->Branch("trackLength", &fTrackLength, "trackLenth/F");
401  fSimulationNtuple->Branch("process", &fSimProcess, "process/I");
402  fSimulationNtuple->Branch("endProcess", &fSimEndProcess, "endProcess/I");
403  fSimulationNtuple->Branch("parentPDG", &fParentPDG, "parentPDG/I");
404  fSimulationNtuple->Branch("parentE", &fParentE, "parentE/F");
405  fSimulationNtuple->Branch("progenitor", &fProgenitor, "progenitor/I");
406 
407  // CRT hits
408  fSimulationNtuple->Branch("auxDetSensitiveID", &fAuxDetSensitiveID);
409  fSimulationNtuple->Branch("auxDetID", &fAuxDetID);
410  fSimulationNtuple->Branch("auxDetEDep", &fADEDep);
411  fSimulationNtuple->Branch("auxDetdEdx", &fADdEdx);
412  fSimulationNtuple->Branch("auxDetTrackLength", &fADTrackLength);
413  fSimulationNtuple->Branch("auxDetEnterXYZT", &fADEnterXYZT);
414  fSimulationNtuple->Branch("auxDetExitXYZT", &fADExitXYZT);
415  fSimulationNtuple->Branch("auxDetEnterPE", &fADEnterPE);
416  fSimulationNtuple->Branch("auxDetExitPE", &fADExitPE);
417  fSimulationNtuple->Branch("auxDetRegion", &fAuxDetReg);
418  fSimulationNtuple->Branch("mac5", &fADMac);
419  fSimulationNtuple->Branch("adType", &fADType);
420 
421  fSimulationNtuple->Branch("startXYZT", fStartXYZT, "startXYZT[4]/F");
422  fSimulationNtuple->Branch("endXYZT", fEndXYZT, "endXYZT[4]/F");
423  fSimulationNtuple->Branch("startPE", fStartPE, "startPE[4]/F");
424  fSimulationNtuple->Branch("endPE", fEndPE, "endPE[4]/F");
425  fSimulationNtuple->Branch("nChan", &fNAuxDet, "nChan/I");
426  fSimulationNtuple->Branch("mother", &fMother, "mother/I");
427  fSimulationNtuple->Branch("nDaught", &fNDaught, "nDaught/I");
428 
429  //regions tree
430  fRegionsNtuple->Branch("event", &fEvent, "event/I");
431  fRegionsNtuple->Branch("run", &fRun, "run/I");
432  fRegionsNtuple->Branch("subRun", &fSubRun, "subRun/I");
433  fRegionsNtuple->Branch("nReg", &fNReg, "nReg/I");
434  fRegionsNtuple->Branch("fiducial", &fRegFid, "fiducial/I");
435  fRegionsNtuple->Branch("active", &fRegActive, "active/I");
436  fRegionsNtuple->Branch("inactive", &fRegInactive, "inactive/I");
437  fRegionsNtuple->Branch("crts", &fRegCRTs, "crts/I");
438  fRegionsNtuple->Branch("regions", &fRegRegions);
439  fRegionsNtuple->Branch("pdg", &fRegPDG, "pdg/I");
440  fRegionsNtuple->Branch("trackID", &fRegTrkID, "trackID/I");
441  fRegionsNtuple->Branch("eDep", &fRegEDep);
442  fRegionsNtuple->Branch("dL", &fRegdL);
443  fRegionsNtuple->Branch("opDetID", &fRegOpDetID);
444  fRegionsNtuple->Branch("distToOpDet", &fRegDistToOpDet);
445  fRegionsNtuple->Branch("opDetXYZT", &fRegOpDetXYZT);
446  fRegionsNtuple->Branch("entryPE", &fRegEntryPE);
447  fRegionsNtuple->Branch("exitPE", &fRegExitPE);
448  fRegionsNtuple->Branch("entryXYZT", &fRegEntryXYZT);
449  fRegionsNtuple->Branch("exitXYZT", &fRegExitXYZT);
450  fRegionsNtuple->Branch("entrySlope", &fRegEntrySlope);
451  fRegionsNtuple->Branch("exitSlope", &fRegExitSlope);
452 
453  // Define the branches of our DetSim n-tuple
454  fDetSimNtuple->Branch("event", &fEvent, "event/I");
455  fDetSimNtuple->Branch("run", &fRun, "run/I");
456  fDetSimNtuple->Branch("subRun", &fSubRun, "subRun/I");
457  fDetSimNtuple->Branch("nAbove", &fNAbove, "nAbove/I");
458  fDetSimNtuple->Branch("t0", &fT0, "t0/I");
459  fDetSimNtuple->Branch("t1", &fT1, "t1/I");
460  fDetSimNtuple->Branch("adc", fADC, "adc[64]/I");
461  fDetSimNtuple->Branch("maxAdc", &fMaxAdc, "maxAdc/I");
462  fDetSimNtuple->Branch("maxChan", &fMaxChan, "maxChan/I");
463  fDetSimNtuple->Branch("trackID", &fTrackID);
464  fDetSimNtuple->Branch("detPDG", &fDetPDG);
465  fDetSimNtuple->Branch("entry", &fEntry, "entry/I");
466  fDetSimNtuple->Branch("mac5", &fMac5, "mac5/I");
467  fDetSimNtuple->Branch("region", &fFEBReg, "region/I");
468  fDetSimNtuple->Branch("subSys", &fDetSubSys, "subSys/I");
469 
470  // Define the branches of our SimHit n-tuple
471  fSimHitNtuple->Branch("event", &fEvent, "event/I");
472  fSimHitNtuple->Branch("run", &fRun, "run/I");
473  fSimHitNtuple->Branch("subRun", &fSubRun, "subRun/I");
474  fSimHitNtuple->Branch("nHit", &fNHit, "nHit/I");
475  fSimHitNtuple->Branch("x", &fXHit, "x/F");
476  fSimHitNtuple->Branch("y", &fYHit, "y/F");
477  fSimHitNtuple->Branch("z", &fZHit, "z/F");
478  fSimHitNtuple->Branch("xErr", &fXHitErr, "xErr/F");
479  fSimHitNtuple->Branch("yErr", &fYHitErr, "yErr/F");
480  fSimHitNtuple->Branch("zErr", &fZHitErr, "zErr/F");
481  fSimHitNtuple->Branch("t0", &fT0Hit, "t0/F");
482  fSimHitNtuple->Branch("t1", &fT1Hit, "t1/F");
483  fSimHitNtuple->Branch("region", &fHitReg, "region/I");
484  fSimHitNtuple->Branch("subSys", &fHitSubSys, "subSys/I");
485  fSimHitNtuple->Branch("trackID", &fHitTrk);
486  fSimHitNtuple->Branch("pdg", &fHitPDG);
487  fSimHitNtuple->Branch("modID", &fHitMod);
488  fSimHitNtuple->Branch("stripID", &fHitStrip);
489  fSimHitNtuple->Branch("nFeb", &fNHitFeb, "nFeb/I");
490  fSimHitNtuple->Branch("hitPe", &fHitPe);
491  fSimHitNtuple->Branch("totPe", &fHitTotPe, "totPe/F");
492  fSimHitNtuple->Branch("rmsPe", &fHitPeRms, "rmsPe/F");
493 
494  // Define the branches of our SimTrueHit n-tuple
495  fTrueCRTHitNtuple->Branch("event", &fEvent, "event/I");
496  fTrueCRTHitNtuple->Branch("run", &fRun, "run/I");
497  fTrueCRTHitNtuple->Branch("subRun", &fSubRun, "subRun/I");
498  fTrueCRTHitNtuple->Branch("nHit", &fTrueNHit, "nHit/I");
499  fTrueCRTHitNtuple->Branch("x", &fTrueXHit, "x/F");
500  fTrueCRTHitNtuple->Branch("y", &fTrueYHit, "y/F");
501  fTrueCRTHitNtuple->Branch("z", &fTrueZHit, "z/F");
502  fTrueCRTHitNtuple->Branch("xErr", &fTrueXHitErr, "xErr/F");
503  fTrueCRTHitNtuple->Branch("yErr", &fTrueYHitErr, "yErr/F");
504  fTrueCRTHitNtuple->Branch("zErr", &fTrueZHitErr, "zErr/F");
505  fTrueCRTHitNtuple->Branch("t0", &fTrueT0Hit, "t0/F");
506  fTrueCRTHitNtuple->Branch("t1", &fTrueT1Hit, "t1/F");
507  fTrueCRTHitNtuple->Branch("region", &fTrueHitReg, "region/I");
508  fTrueCRTHitNtuple->Branch("subSys", &fTrueHitSubSys, "subSys/I");
509  fTrueCRTHitNtuple->Branch("trackID", &fTrueHitTrk);
510  fTrueCRTHitNtuple->Branch("pdg", &fTrueHitPDG);
511  fTrueCRTHitNtuple->Branch("modID", &fTrueHitMod);
512  fTrueCRTHitNtuple->Branch("stripID", &fTrueHitStrip);
513  fTrueCRTHitNtuple->Branch("nFeb", &fTrueNHitFeb, "nFeb/I");
514  fTrueCRTHitNtuple->Branch("hitPe", &fTrueHitPe);
515  fTrueCRTHitNtuple->Branch("totPe", &fTrueHitTotPe, "totPe/F");
516  fTrueCRTHitNtuple->Branch("rmsPe", &fTrueHitPeRms, "rmsPe/F");
517 
518  fSimTrackNtuple->Branch("ntrack", &fNSimTrack, "ntrack/I");
519  fSimTrackNtuple->Branch("pe", &fSimTrackPE, "pe/F");
520  fSimTrackNtuple->Branch("t", &fSimTrackT0, "t/D");
521  fSimTrackNtuple->Branch("start", fSimTrackStart, "start[3]/F");
522  fSimTrackNtuple->Branch("end", fSimTrackEnd, "end[3]/F");
523  fSimTrackNtuple->Branch("l", &fSimTrackL, "l/F");
524  fSimTrackNtuple->Branch("theta", &fSimTrackTheta, "theta/F");
525  fSimTrackNtuple->Branch("phi", &fSimTrackPhi, "phi/F");
526  fSimTrackNtuple->Branch("nhit", &fNHitSimTrack, "nhit/I");
527  fSimTrackNtuple->Branch("hitstart", fSimTrackHitStart, "hitstart[4]/F");
528  fSimTrackNtuple->Branch("hitend", fSimTrackHitEnd, "hitend[4]/F");
529  fSimTrackNtuple->Branch("regstart", &fSimTrackRegStart, "regstart/I");
530  fSimTrackNtuple->Branch("regend", &fSimTrackRegEnd, "regend/I");
531 
532 
533 }
534 
535  void CRTSimAnalysis::beginRun(const art::Run& /*run*/)
536  {
537  //art::ServiceHandle<sim::LArG4Parameters> larParameters;
538  //fElectronsToGeV = 1./larParameters->GeVToElectrons();
539  std::cout << "beginning run" << std::endl;
540  }
541 
542  //-----------------------------------------------------------------------
543  void CRTSimAnalysis::analyze(const art::Event& event)
544  {
545  MF_LOG_DEBUG("CRT") << "beginning analyis" << '\n';
546 
547  // Check that lists for momenta limits is same size as last of PDGs from FHiCL
548  if (fPDGs.size() != fMinMomenta.size() || fPDGs.size() != fMaxMomenta.size())
549  throw cet::exception("CRTSimAnalysis")
550  << " PDG/Momtenta values not set correctly in fhicl - lists have different sizes"
551  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
552 
553 
554  // Start by fetching some basic event information for our n-tuple.
555  fEvent = event.id().event();
556  fRun = event.run();
557  fSubRun = event.subRun();
558 
559  //CRTBackTracker matches CRTProducts to the true trackIDs
560  //bt.Initialize(event);
561 
562  // Define "handle" to Generator level MCTruth objects
563  art::Handle< vector<simb::MCTruth>> genHandle;
564 
565  // Define a "handle" to point to a vector of MCParticle objects.
566  art::Handle< vector<simb::MCParticle> > particleHandle;
567  map< int, const simb::MCParticle*> particleMap;
568 
569  if (!event.getByLabel("generator", genHandle)) {
570  std::cout << "could not get handle to gen objects!!!" << std::endl;
571  }
572 
573  // If there aren't any simb::MCParticle object art will
574  // display a "ProductNotFound" exception message and may skip
575  // all processing for the rest of this event or stop the execution.
576  if (!event.getByLabel(fSimulationProducerLabel, particleHandle))
577  {
578  // If we have no MCParticles at all in an event, then we're in
579  // big trouble. Throw an exception.
580  throw cet::exception("CRTSimAnalysis")
581  << " No simb::MCParticle objects in this event - "
582  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
583  }
584 
585  // Handle to AuxDetSimChannel (CRT module) objects generated by LArG4
586  art::Handle<vector<sim::AuxDetSimChannel> > auxDetSimChannelHandle;
587  if (!event.getByLabel(fSimulationProducerLabel, auxDetSimChannelHandle)) {
588  throw cet::exception("CRTSimAnalysis")
589  << " No sim::AuxDetSimChannel objects in this event - "
590  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
591  }
592 
593  if((*genHandle).size()>1)
594  throw cet::exception("CRTSimAnalysis") << "gen stage MCParticle vector has more than 1 entry!" << '\n';
595 
596  auto const& truth = (*genHandle)[0];
597  fNGen = truth.NParticles();
598  fGenTrack.clear();
599  fGenPDG.clear();
600  fGenStartXYZT.clear();
601  fGenEndXYZT.clear();
602  fGenStartPE.clear();
603  fGenEndPE.clear();
604 
605  for ( int i=0; i<fNGen; i++ )
606  {
607  auto const& part = truth.GetParticle(i); //simb::MCParticle
608 
609  fGenTrack.push_back(part.TrackId());
610  fGenPDG.push_back(part.PdgCode());
611 
612  const TLorentzVector startPos = part.Position(0);
613  const TLorentzVector endPos = part.EndPosition();
614  const TLorentzVector startMom = part.Momentum(0);
615  const TLorentzVector endMom = part.EndMomentum();
616 
617  fGenStartXYZT.push_back({startPos.X(),startPos.Y(),startPos.Z(),startPos.T()});
618  fGenEndXYZT.push_back({endPos.X(),endPos.Y(),endPos.Z(),endPos.T()});
619  fGenStartPE.push_back({startMom.Px(),startMom.Py(),startMom.Pz(),startMom.E()});
620  fGenEndPE.push_back({endMom.Px(),endMom.Py(),endMom.Pz(),endMom.E()});
621  }
622 
623  fGenNtuple->Fill();
624 
625  // The MCParticle objects are not necessarily in any particular
626  // order. Since we may have to search the list of particles, let's
627  // put them into a map. To save both space and time, the map
628  // will not contain a copy of the MCParticle, but a pointer to it.
629  for ( auto const& particle : (*particleHandle) )
630  {
631  // Add the address of the MCParticle to the map, with the
632  // track ID as the key.
633  particleMap.insert(std::make_pair(particle.TrackId(),&particle));
634  }
635 
636  std::cout << "event " << fEvent << " with " << particleMap.size() << " MCParticles" << std::endl;
637 
638  //get TPC objects
639  geo::CryostatGeo const& cryo0 = fGeometryService->Cryostat(0);
640  geo::CryostatGeo const& cryo1 = fGeometryService->Cryostat(1);
641 
642  geo::TPCGeo const& tpc00 = cryo0.TPC(0);
643  geo::TPCGeo const& tpc01 = cryo0.TPC(1);
644  geo::TPCGeo const& tpc10 = cryo1.TPC(0);
645  geo::TPCGeo const& tpc11 = cryo1.TPC(1);
646 
647  //loop over MCParticles
648  for ( auto const& particle : (*particleHandle) )
649  {
650  fSimPDG = particle.PdgCode();
651  vector<int>::iterator it = fPDGs.begin(); //iterator to list of interested PDGs from FHiCL
652  const TLorentzVector& momentumStart = particle.Momentum(0);//initial momentum
653  const double p = (momentumStart.Vect()).Mag();
654  size_t index = 0; //index in PDG list
655 
656  //if ( (particle.Process() != "primary" &&
657  //find matching PDG code given current MCParticle PDG
658  // use for momentum cuts
659  while (it!=fPDGs.end()) {
660  if (*it==fSimPDG) {
661  index = (size_t)(it - fPDGs.begin());
662  break;
663  }
664  it++;
665  }
666  //if PDG not included in interest list, skip to next
667  if (!(fPDGs.size()==1 && fPDGs[0]==0) && it == fPDGs.end()) continue;
668  //check momentum is within region of interest
669  if ( fMinMomenta[index] != 0 && p < fMinMomenta[index]) continue;
670  if ( fMaxMomenta[index] != 0 && p > fMaxMomenta[index]) continue;
671  //if ( abs(fSimPDG)==11 && particle.Process()!="compt" && particle.Process()!="conv" )
672  // continue;
673 
674  //count total number of muons present in event
675  //if (abs(fSimPDG)==13){
676  // fNmuTruth++;
677  // }
678  fSimTrackID = particle.TrackId();
679 
680  // the following bit attempts to establish the ancestry
681  // of each MCParticle of interest
682  // this is useful for matching gammas produced by muons
683  // for evaluation of removal algorithms
684  fMother = particle.Mother();
685  fNDaught = particle.NumberDaughters();
686  fSimProcess = ProcessToICode(particle.Process());
687  fSimEndProcess = ProcessToICode(particle.EndProcess());
688 
689  if(fMother!=0){ //if not primary
690  map<int,const simb::MCParticle*>::const_iterator it = particleMap.find(fMother);
691  if(it==particleMap.end()){
692  fParentPDG=INT_MAX;
693  fParentE = FLT_MAX;
694  fProgenitor = INT_MAX;
695  }//if mother not found in particle list
696  else{ //otherwise try to find primary muon
697  fParentPDG = it->second->PdgCode();
698  fParentE = it->second->E(0);
699  int tmpID=it->second->Mother();
700  size_t ctr=0;
701  map<int,const simb::MCParticle*>::iterator it2 = particleMap.begin();
702 
703  if(fParentPDG==13||fParentPDG==-13) fProgenitor = it->second->TrackId();
704  else while(it2!=particleMap.end()&&ctr<particleMap.size()){
705  it2=particleMap.find(tmpID);
706  if(it2!=particleMap.end()){
707  if(it2->second->PdgCode()==13||it2->second->PdgCode()==-13){
708  fProgenitor=tmpID;
709  break;
710  }
711  tmpID = it2->second->Mother();
712  }
713 
714  ctr++;
715  if(ctr==particleMap.size()) std::cout<<"manual break!"<<std::endl;
716  }
717  //std::cout<<"out of progenitor search loop" << std::endl;
718  }//else mother is found
719  }//if particle has mother (not a primary)
720  else{ //if current particle in loop is primary
721  fParentPDG=0;
722  fProgenitor=-10;
723  fParentE=-1.0;
724  }
725  //end of ancestry matching
726  //now get some other useful info about the trajectory
727 
728  // A particle has a trajectory, consisting of a set of
729  // 4-positions and 4-mommenta.
730  fSimHits = particle.NumberTrajectoryPoints();
731  if(fSimHits==1) continue;
732 
733  // For trajectories, as for vectors and arrays, the first
734  // point is #0, not #1.
735  const int last = fSimHits - 1;
736  const TLorentzVector& positionStart = particle.Position(0);
737  const TLorentzVector& positionEnd = particle.Position(last);
738  //const TLorentzVector& momentumStart = particle.Momentum(0);
739  const TLorentzVector& momentumEnd = particle.Momentum(last);
740 
741  // Fill arrays with the 4-values.
742  positionStart.GetXYZT( fStartXYZT );
743  positionEnd.GetXYZT( fEndXYZT );
744  momentumStart.GetXYZT( fStartPE );
745  momentumEnd.GetXYZT( fEndPE );
746  fTrackLength = ( positionEnd - positionStart ).Rho();
747 
748  fNCD = 0;
749  fCDpdg = fSimPDG;
750  fCDTrackID = fSimTrackID;
751  //fCDRegions.clear();
752  fCDSlopes.clear();
753  fCDpe.clear();
754  fCDxyzt.clear();
755 
756  fNReg = 0;
757  fRegFid = 0;
758  fRegActive = 0;
759  fRegInactive = 0;
760  fRegCRTs = 0;
761  fRegPDG = fSimPDG;;
762  fRegTrkID = fSimTrackID;
763  fRegRegions.clear();
764  fRegEDep.clear();
765  fRegDistToOpDet.clear();
766  fRegOpDetID.clear();
767  fRegEntryXYZT.clear();
768  fRegExitXYZT.clear();
769  fRegEntryPE.clear();
770  fRegExitPE.clear();
771  fRegOpDetXYZT.clear();
772  fRegEntrySlope.clear();
773  fRegExitSlope.clear();
774 
775  int oldreg = -1;
776 
777  //loop over trajectory points
778  for (unsigned int i=0; i<fSimHits; i++){
779  const TLorentzVector& pos = particle.Position(i); // 4-position in World coordinates
780  const TLorentzVector& posnext = particle.Position(i+1); // problem for last point???
781  const TLorentzVector& mom = particle.Momentum(i); // 4-momentum
782  const double point[3] = {pos.X(),pos.Y(),pos.Z()};
783  const double pointnext[3] = {posnext.X(),posnext.Y(),posnext.Z()};
784  double opDetPos[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
785  double entryPos[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
786  double entryT = -FLT_MAX;
787  bool active0 = false, active1 = false, activenext0 = false, activenext1 = false;
788 
789  // CosmicDisplay info
790  fCDxyzt.push_back({pos.X(),pos.Y(),pos.Z(),pos.T()});
791  fCDpe.push_back({mom.Px(),mom.Py(),mom.Pz(),mom.E()});
792  fCDSlopes.push_back({mom.Px()/mom.P(), mom.Py()/mom.P(), mom.Pz()/mom.P()});
793  fNCD++;
794 
795  // Regions info
796  // Check if trajectory points are in cryostats (active + inactve LAr )
797  if(cryo0.ContainsPosition(point)) {
798  active0 = tpc00.ContainsPosition(point);
799  active1 = tpc01.ContainsPosition(point);
800  activenext0 = tpc00.ContainsPosition(pointnext);
801  activenext1 = tpc01.ContainsPosition(pointnext);
802 
803  // if last point was not in this cryostat or is now entering AV
804  if ( (oldreg!=10&&!active0&&!active1) || (active0&&oldreg!=5) || (active1&&oldreg!=6)) {
805  fRegEntryXYZT.push_back({pos.X(),pos.Y(),pos.Z(),pos.T()});
806  fRegEntryPE.push_back({mom.Px(),mom.Py(),mom.Pz(),mom.E()});
807  fRegEntrySlope.push_back({mom.Px()/mom.P(), mom.Py()/mom.P(), mom.Pz()/mom.P()});
808  oldreg = 10;
809  if (active0) oldreg = 5;
810  if (active1) oldreg = 6;
811  }
812 
813  // if next point is outside of this volume or is last traj. point
814  if (!cryo0.ContainsPosition(pointnext) || (oldreg==10&&(activenext0||activenext1))
815  || i==fSimHits-1
816  || (active0 && !activenext0) || (active1&&!activenext1) ){
817 
818  if (!cryo0.ContainsPosition(pointnext))
819  oldreg=-1;
820 
821  fRegExitXYZT.push_back({pos.X(),pos.Y(),pos.Z(),pos.T()});
822  fRegExitPE.push_back({mom.Px(),mom.Py(),mom.Pz(),mom.E()});
823  fRegExitSlope.push_back({mom.Px()/mom.P(), mom.Py()/mom.P(), mom.Pz()/mom.P()});
824  if (active0) {
825  fRegRegions.push_back(5);
826  fRegActive++;
827  if(tpc00.InFiducialX(point[0],25,0) && tpc00.InFiducialY(point[1],25,25)
828  && tpc00.InFiducialZ(point[2],30,50))
829  fRegFid++;
830 
831  }
832  else if (active1) {
833  fRegRegions.push_back(6);
834  fRegActive++;
835  if(tpc01.InFiducialX(point[0],25,0) && tpc01.InFiducialY(point[1],25,25)
836  && tpc01.InFiducialZ(point[2],30,50))
837  fRegFid++;
838  }
839  else {
840  fRegRegions.push_back(10);
841  fRegInactive++;
842  }
843  if(fRegExitXYZT.size()!=fRegEntryXYZT.size())
844  std::cout << "entry/exit point size mismatch! " <<
845  fRegEntryXYZT.size() << " vs. " << fRegExitXYZT.size() << std::endl;
846  fRegdL.push_back(sqrt(pow(fRegExitXYZT[fNReg][0]-fRegEntryXYZT[fNReg][0],2)
847  +pow(fRegExitXYZT[fNReg][1]-fRegEntryXYZT[fNReg][1],2)
848  +pow(fRegExitXYZT[fNReg][2]-fRegEntryXYZT[fNReg][2],2)));
849  fRegEDep.push_back(fRegEntryPE[fNReg][3] - fRegExitPE[fNReg][3]);
850  for (int index=0; index<3; index++) entryPos[index] = fRegEntryXYZT[fNReg][index];
851  entryT = fRegEntryXYZT[fNReg][3];
852  fRegOpDetID.push_back(cryo0.GetClosestOpDet(entryPos));
853  geo::OpDetGeo const& opDet = cryo0.OpDet(fRegOpDetID[fNReg]);
854  opDet.GetCenter(opDetPos);
855  fRegDistToOpDet.push_back(sqrt(pow(opDetPos[0]-entryPos[0],2)
856  + pow(opDetPos[1]-entryPos[1],2)
857  + pow(opDetPos[2]-entryPos[2],2)));
858  fRegOpDetXYZT.push_back({});
859  for (int index=0; index<3; index++) fRegOpDetXYZT[fNReg].push_back(opDetPos[index]);
860  fRegOpDetXYZT[fNReg].push_back(entryT + fRegDistToOpDet[fNReg]*LAR_PROP_DELAY);
861  fNReg++;
862  }
863  } //if cryo0
864 
865  // if this point in the other cryostat
866  if(cryo1.ContainsPosition(point)) {
867  //check if this or next points are in active volumes
868  active0 = tpc10.ContainsPosition(point);
869  active1 = tpc11.ContainsPosition(point);
870  activenext0 = tpc10.ContainsPosition(pointnext);
871  activenext1 = tpc11.ContainsPosition(pointnext);
872 
873  // if last point was not in this cryostat or is now entering AV
874  if ( (oldreg!=12&&!active0&&!active1) || (active0&&oldreg!=7) || (active1&&oldreg!=8)) {
875  fRegEntryXYZT.push_back({pos.X(),pos.Y(),pos.Z(),pos.T()});
876  fRegEntryPE.push_back({mom.Px(),mom.Py(),mom.Pz(),mom.E()});
877  fRegEntrySlope.push_back({mom.Px()/mom.P(), mom.Py()/mom.P(), mom.Pz()/mom.P()});
878  oldreg = 12;
879  if (active0) oldreg = 7;
880  if (active1) oldreg = 8;
881  }
882 
883  if (!cryo1.ContainsPosition(pointnext) || (oldreg==12&&(activenext0||activenext1))
884  || i==fSimHits-1
885  || (active0 && !activenext0) || (active1&&!activenext1) ){
886 
887  if (!cryo1.ContainsPosition(pointnext))
888  oldreg=-1;
889 
890  fRegExitXYZT.push_back({pos.X(),pos.Y(),pos.Z(),pos.T()});
891  fRegExitPE.push_back({mom.Px(),mom.Py(),mom.Pz(),mom.E()});
892  fRegExitSlope.push_back({mom.Px()/mom.P(), mom.Py()/mom.P(), mom.Pz()/mom.P()});
893  if (active0) {
894  fRegRegions.push_back(7);
895  fRegActive++;
896  if(tpc10.InFiducialX(point[0],25,0) && tpc10.InFiducialY(point[1],25,25)
897  && tpc10.InFiducialZ(point[2],30,50))
898  fRegFid++;
899 
900  }
901  else if (active1) {
902  fRegRegions.push_back(8);
903  fRegActive++;
904  if(tpc11.InFiducialX(point[0],25,0) && tpc11.InFiducialY(point[1],25,25)
905  && tpc11.InFiducialZ(point[2],30,50))
906  fRegFid++;
907  }
908  else {
909  fRegRegions.push_back(12);
910  fRegInactive++;
911  }
912  fRegdL.push_back(sqrt(pow(fRegExitXYZT[fNReg][0]-fRegEntryXYZT[fNReg][0],2)
913  +pow(fRegExitXYZT[fNReg][1]-fRegEntryXYZT[fNReg][1],2)
914  +pow(fRegExitXYZT[fNReg][2]-fRegEntryXYZT[fNReg][2],2)));
915  fRegEDep.push_back(fRegEntryPE[fNReg][3] - fRegExitPE[fNReg][3]);
916  for (int index=0; index<3; index++) entryPos[index] = fRegEntryXYZT[fNReg][index];
917  entryT = fRegEntryXYZT[fNReg][3];
918  fRegOpDetID.push_back(cryo1.GetClosestOpDet(entryPos));
919  geo::OpDetGeo const& opDet = cryo1.OpDet(fRegOpDetID[fNReg]);
920  opDet.GetCenter(opDetPos);
921  fRegDistToOpDet.push_back(sqrt(pow(opDetPos[0]-entryPos[0],2)
922  + pow(opDetPos[1]-entryPos[1],2)
923  + pow(opDetPos[2]-entryPos[2],2)));
924  fRegOpDetXYZT.push_back({});
925  for (int index=0; index<3; index++) fRegOpDetXYZT[fNReg].push_back(opDetPos[index]);
926  fRegOpDetXYZT[fNReg].push_back(entryT + fRegDistToOpDet[fNReg]*LAR_PROP_DELAY);
927  fNReg++;
928  } // if exiting from volume
929  } //if cryo1
930 
931  }//for trajectory points
932 
933  fCosmicDisplayNtuple->Fill();
934 
935  //map module IDs to strip IDs hit by muons
936  //map< uint16_t,set<uint8_t>* > muHitMapC; //hits in C modules only
937  //map< uint16_t,set<uint8_t>* > muHitMapM; //hits in M modules only
938  //map< uint16_t,set<uint8_t>* > muHitMapD; //hits in D modules only
939  //map< uint16_t,set<uint8_t>* > muHitMap; //all hits
940 
941  map< int, vector<double> > regCRTEnter, regCRTExit, regCRTEnterPE, regCRTExitPE;
942 
943  //reinitialize ADChannel vars
944  fNAuxDet = 0;
945  fADType.clear();
946  fAuxDetReg.clear();
947  fADTrackLength.clear();
948  fADEDep.clear();
949  fADdEdx.clear();
950  fAuxDetID.clear();
951  fAuxDetSensitiveID.clear();
952  fADEnterXYZT.clear();
953  fADExitXYZT.clear();
954  fADEnterPE.clear();
955  fADExitPE.clear();
956  fADMac.clear();
957 
958  // To look at the energy deposited by this particle's track,
959  // we loop over the AuxDetSimChannel objects in the event.
960  // Note all volumes are included, not just ones with energy deps
961  for ( auto const& channel : (*auxDetSimChannelHandle) )
962  {
963  // Get vector of hits in this AuxDet channel
964  auto const& auxDetIDEs = channel.AuxDetIDEs();
965 
966  // For every hit in this channel:
967  for ( auto const& ide : auxDetIDEs )
968  {
969  // Check if the track that deposited the
970  // energy matches the track of the MCParticle.
971  if ( ide.trackID != fSimTrackID ) continue;
972  if ( ide.energyDeposited * 1.0e6 < 50 ) continue;
973 
974  size_t adid = channel.AuxDetID();
975  //size_t adsid = channel.AuxDetSensitiveID();
976  fADType.push_back(fCrtutils->GetAuxDetTypeCode(adid));
977  fAuxDetReg.push_back(
978  fCrtutils->AuxDetRegionNameToNum(fCrtutils->GetAuxDetRegion(adid)));
979 
980  // What is the distance from the hit (centroid of the entry
981  // and exit points) to the readout end?
982  /*double trueX = (ide.entryX + ide.exitX) / 2.0;
983  double trueY = (ide.entryY + ide.exitY) / 2.0;
984  double trueZ = (ide.entryZ + ide.exitZ) / 2.0;
985  //double trueT = (ide.entryT + ide.exitT) / 2.0;
986  double pos[3] = {trueX,trueY,trueZ};
987  size_t chanad, chanads;
988  int32_t adchan = fGeometryService->PositionToAuxDetChannel(pos,chanad,chanads);
989  std::cout << "retrieved AuxDetChannel from position (" << trueX << ", " << trueY
990  << ", " << trueZ << ")" << '\n'
991  << " true AuxDetID, AuxDetSensitiveID: " << adid << ", " << adsid << '\n'
992  << " from map-> chan, adid, adsid: " << adchan << ", " << chanad
993  << ", " << chanads << std::endl;*/
994 
995  //calculate track length in strip
996  double dx = ide.entryX-ide.exitX;
997  double dy = ide.entryY-ide.exitY;
998  double dz = ide.entryZ-ide.exitZ;
999  double adlength = sqrt(dx*dx+dy*dy+dz*dz);
1000  if ( adlength < 0.0001) continue;
1001 
1002  fADTrackLength.push_back(adlength);
1003  fADEDep.push_back(ide.energyDeposited);
1004  fADdEdx.push_back(ide.energyDeposited/fADTrackLength.back());
1005  fAuxDetID.push_back(channel.AuxDetID());
1006  fAuxDetSensitiveID.push_back(channel.AuxDetSensitiveID());
1007  fADEnterXYZT.push_back({ide.entryX,ide.entryY,ide.entryZ,ide.entryT});
1008  fADExitXYZT.push_back({ide.exitX,ide.exitY,ide.exitZ,ide.exitT});
1009  //we dont have entry mom. or total E info on adsc, so very rough approx. E~P
1010  double pmag = sqrt(pow(ide.exitMomentumX,2)+pow(ide.exitMomentumY,2)+pow(ide.exitMomentumZ,2));
1011  fADExitPE.push_back({ide.exitMomentumX,ide.exitMomentumY,ide.exitMomentumZ,pmag});
1012  fADEnterPE.push_back({fADExitPE[fNAuxDet][0]+pmag/3,fADExitPE[fNAuxDet][1]+pmag/3,
1013  fADExitPE[fNAuxDet][2]+pmag/3,pmag+ide.energyDeposited});
1014  fADMac.push_back(fCrtutils->ADToMac(channel.AuxDetID()).first);
1015 
1016  if (regCRTEnter.find(fAuxDetReg[fNAuxDet])!=regCRTEnter.end()) {
1017  if (regCRTEnter[fAuxDetReg[fNAuxDet]][3] > ide.entryT) {
1018  regCRTEnter[fAuxDetReg[fNAuxDet]] = fADEnterXYZT[fNAuxDet];
1019  regCRTEnterPE[fAuxDetReg[fNAuxDet]] = fADEnterPE[fNAuxDet];
1020  }
1021  }
1022  else {
1023  regCRTEnter[fAuxDetReg[fNAuxDet]] = fADEnterXYZT[fNAuxDet];
1024  regCRTEnterPE[fAuxDetReg[fNAuxDet]] = fADEnterPE[fNAuxDet];
1025  }
1026 
1027  if (regCRTExit.find(fAuxDetReg[fNAuxDet])!=regCRTExit.end()) {
1028  if (regCRTExit[fAuxDetReg[fNAuxDet]][3] < ide.exitT){
1029  regCRTExit[fAuxDetReg[fNAuxDet]] = fADExitXYZT[fNAuxDet];
1030  regCRTExitPE[fAuxDetReg[fNAuxDet]] = fADExitPE[fNAuxDet];
1031  }
1032  }
1033  else {
1034  regCRTExit[fAuxDetReg[fNAuxDet]] = fADExitXYZT[fNAuxDet];
1035  regCRTExitPE[fAuxDetReg[fNAuxDet]] = fADExitPE[fNAuxDet];
1036  }
1037 
1038  fNAuxDet++;
1039 
1040  } // For each IDE (strip hit by muon)
1041  } // For each SimChannel (module)
1042 
1043  // write values to tree for this event and particle
1044  fSimulationNtuple->Fill();
1045 
1046  for( auto it=regCRTEnter.begin(); it!=regCRTEnter.end(); it++) {
1047  //std::cout << "found CRT region " << it->first << std::endl;
1048  fRegRegions.push_back(it->first);
1049  fRegEntryXYZT.push_back({(it->second)[0],(it->second)[1],(it->second)[2],(it->second)[3]});
1050  fRegExitXYZT.push_back({regCRTExit[it->first][0],regCRTExit[it->first][1],regCRTExit[it->first][2],
1051  regCRTExit[it->first][3]});
1052  fRegdL.push_back(sqrt(pow(fRegExitXYZT[fNReg][0]-fRegEntryXYZT[fNReg][0],2)
1053  +pow(fRegExitXYZT[fNReg][1]-fRegEntryXYZT[fNReg][1],2)
1054  +pow(fRegExitXYZT[fNReg][2]-fRegEntryXYZT[fNReg][2],2)));
1055  fRegEntryPE.push_back({regCRTEnterPE[it->first][0],regCRTEnterPE[it->first][1],
1056  regCRTEnterPE[it->first][2], regCRTEnterPE[it->first][3]});
1057  fRegExitPE.push_back({regCRTExitPE[it->first][0],regCRTExitPE[it->first][1],
1058  regCRTExitPE[it->first][2], regCRTExitPE[it->first][3]});
1059  fRegEDep.push_back(fRegEntryPE[fNReg][3] - fRegExitPE[fNReg][3]);
1060  fRegDistToOpDet.push_back(-1);
1061  fRegOpDetID.push_back(-1);
1062  fRegOpDetXYZT.push_back({-1,-1,-1,-1});
1063  fRegEntrySlope.push_back({-1,-1,-1});
1064  fRegExitSlope.push_back({-1,-1,-1});
1065 
1066  fNReg++; fRegCRTs++;
1067  }
1068 
1069  //sort region tree entries by entry time
1070  int flag = 1; // set flag to 1 to start first pass
1071  int tempInt;
1072  double tempDoub;
1073 
1074  for(int i = 1; (i <= (int)fNReg) && flag; i++)
1075  {
1076  flag = 0;
1077  for (int j=0; j < ((int)fNReg -1); j++)
1078  {
1079  if (fRegEntryXYZT[j+1][3] < fRegEntryXYZT[j][3])
1080  {
1081  tempInt = fRegRegions[j]; // swap regions
1082  fRegRegions[j] = fRegRegions[j+1];
1083  fRegRegions[j+1] = tempInt;
1084  tempDoub = fRegEDep[j]; // swap EDep
1085  fRegEDep[j] = fRegEDep[j+1];
1086  fRegEDep[j+1] = tempDoub;
1087  tempDoub = fRegdL[j]; // swap dL
1088  fRegdL[j] = fRegdL[j+1];
1089  fRegdL[j+1] = tempDoub;
1090  tempDoub = fRegDistToOpDet[j]; //swap distToOpDet
1091  fRegDistToOpDet[j] = fRegDistToOpDet[j+1];
1092  fRegDistToOpDet[j+1] = tempDoub;
1093  tempInt = fRegOpDetID[j]; //swap opDetID
1094  fRegOpDetID[j] = fRegOpDetID[j+1];
1095  fRegOpDetID[j+1] = tempInt;
1096  for (int k=0; k<4; k++) {
1097  tempDoub = fRegEntryPE[j][k]; //swap entryPE
1098  fRegEntryPE[j][k] = fRegEntryPE[j+1][k];
1099  fRegEntryPE[j+1][k] = tempDoub;
1100  tempDoub = fRegExitPE[j][k]; //swap exitPE
1101  fRegExitPE[j][k] = fRegExitPE[j+1][k];
1102  fRegExitPE[j+1][k] = tempDoub;
1103  tempDoub = fRegEntryXYZT[j][k]; //swap entryXYZT
1104  fRegEntryXYZT[j][k] = fRegEntryXYZT[j+1][k];
1105  fRegEntryXYZT[j+1][k] = tempDoub;
1106  tempDoub = fRegExitXYZT[j][k]; //swap exitXYZT
1107  fRegExitXYZT[j][k] = fRegExitXYZT[j+1][k];
1108  fRegExitXYZT[j+1][k] = tempDoub;
1109  tempDoub = fRegOpDetXYZT[j][k]; //swap opDetXYZT
1110  fRegOpDetXYZT[j][k] = fRegOpDetXYZT[j+1][k];
1111  fRegOpDetXYZT[j+1][k] = tempDoub;
1112  if(k<3) {
1113  tempDoub = fRegEntrySlope[j][k];
1114  fRegEntrySlope[j][k] = fRegEntrySlope[j+1][k];
1115  fRegEntrySlope[j+1][k] = tempDoub;
1116  tempDoub = fRegExitSlope[j][k];
1117  fRegExitSlope[j][k] = fRegExitSlope[j+1][k];
1118  fRegExitSlope[j+1][k] = tempDoub;
1119  }
1120  }
1121  flag = 1; // indicates that a swap occurred.
1122  }
1123  }
1124  }
1125 
1126  fRegionsNtuple->Fill();
1127 
1128  } // loop over all particles in the event.
1129 
1130  //CRTData
1131  art::Handle<vector<CRTData>> crtDetSimHandle;
1132  if (event.getByLabel(fCRTDetSimProducerLabel, crtDetSimHandle)) {
1133 
1134  vector<art::Ptr<CRTData>> febdata;
1135  art::fill_ptr_vector(febdata,crtDetSimHandle);
1136  mf::LogPrint("CRTSimAnalysis") << "found " << febdata.size() << " CRTData" << '\n';
1137 
1138  for ( auto const& data : febdata ) {
1139  fMac5 = data->fMac5;
1140  fEntry = data->fEntry;
1141  fFEBReg = fCrtutils->AuxDetRegionNameToNum(fCrtutils->MacToRegion(fMac5));
1142  fDetSubSys = fCrtutils->MacToTypeCode(fMac5);
1143  fT0 = data->fTs0;
1144  fT1 = data->fTs1;
1145  fMaxAdc = 0;
1146  fMaxChan = -1;
1147  fNAbove = 0;
1148  fTrackID.clear();
1149  fDetPDG.clear();
1150  bool passfilter = false;
1151  if(fPDGs.size()==1 && fPDGs[0]==0)
1152  passfilter = true;
1153  //set max channel by subsystem (32 (c,m) or 64 (d))
1154  //int chmax = 0;
1155  //if(fDetSubSys!=2) chmax=32;
1156  //else chmax=64;
1157 
1158  //loop over all FEB channels
1159  for(int ch=0; ch<64; ch++) {
1160  fADC[ch] = data->fAdc[ch];
1161  if(fADC[ch] > 300) fNAbove++;
1162  if(fADC[ch]>fMaxAdc) {
1163  fMaxAdc = fADC[ch];
1164  fMaxChan = ch;
1165  }
1166  //loop over all track IDs
1167  }//loop over channels
1168  //loop over all track IDs
1169  for( int trk : bt.AllTrueIds(event,*data)) {
1170  fTrackID.push_back(trk);
1171  if (particleMap.find(trk) != particleMap.end() ){
1172  fDetPDG.push_back(particleMap[trk]->PdgCode());
1173  if(!passfilter)
1174  for(int const& pdg: fPDGs) {
1175  if(fDetPDG.back()==pdg)
1176  passfilter = true;
1177  }
1178  }
1179  else
1180  fDetPDG.push_back(INT_MAX);
1181  }//for trackIDs
1182 
1183  if(passfilter)
1184  fDetSimNtuple->Fill();
1185 
1186  } //for CRT FEB events
1187 
1188  }//if crtdetsim products present
1189 
1190  else
1191  mf::LogWarning("CRTSimAnalysis") << "CRTData products not found! (expected if gen/G4 step)" << '\n';
1192 
1193 
1194  //simulted CRTHits
1195  art::Handle<vector<CRTHit>> crtSimHitHandle;
1196  fNHit = 0;
1197  if (event.getByLabel(fCRTSimHitProducerLabel, crtSimHitHandle)) {
1198 
1199  vector<art::Ptr<CRTHit>> crtSimHits;
1200  art::fill_ptr_vector(crtSimHits,crtSimHitHandle);
1201 
1202  mf::LogPrint("CRTSimAnalysis") << "found " << crtSimHits.size() << " sim CRTHits" << '\n';
1203  for ( auto const& hit : crtSimHits )
1204  {
1205  fNHit++;
1206  fXHit = hit->x_pos;
1207  fYHit = hit->y_pos;
1208  fZHit = hit->z_pos;
1209  fXHitErr = hit->x_err;
1210  fYHitErr = hit->y_err;
1211  fZHitErr = hit->z_err;
1212  fT0Hit = hit->ts0_ns;
1213  fT1Hit = hit->ts1_ns;
1214  fNHitFeb = hit->feb_id.size();
1215  fHitTotPe = hit->peshit;
1216  fHitPeRms = 0.;
1217  fHitPe.clear();
1218  fHitTrk.clear();
1219  fHitPDG.clear();
1220  fHitMod.clear();
1221  fHitStrip.clear();
1222 
1223  uint8_t mactmp = hit->feb_id[0];
1224  fHitReg = fCrtutils->AuxDetRegionNameToNum(fCrtutils->MacToRegion(mactmp));
1225  fHitSubSys = fCrtutils->MacToTypeCode(mactmp);
1226 
1227  bool passfilter=false;
1228  if(fPDGs.size()==1 && fPDGs[0]==0)
1229  passfilter = true;
1230 
1231  for(auto const& mactopes : hit->pesmap){
1232  //std::cout << "SimHit with mac5 " << (int)mactopes.first << std::endl;
1233  for(auto const& chanpe : mactopes.second) {
1234  //std::cout << " chan: " << chanpe.first << std::endl;
1235  fHitMod.push_back(fCrtutils->MacToAuxDetID(mactopes.first, chanpe.first));
1236  fHitStrip.push_back(fCrtutils->ChannelToAuxDetSensitiveID(mactopes.first,chanpe.first));
1237  fHitPe.push_back(chanpe.second);
1238  fHitPeRms+=pow(chanpe.second-fHitTotPe,2);
1239  }
1240  }
1241  fHitPeRms = sqrt(fHitPeRms/(fHitPe.size()-1));
1242 
1243  //loop over all track IDs
1244  for( int trk : bt.AllTrueIds(event,*hit)) {
1245  fHitTrk.push_back(trk);
1246  if ( particleMap.find(trk) != particleMap.end()) {
1247  fHitPDG.push_back(particleMap[trk]->PdgCode());
1248  if(!passfilter)
1249  for(int const& pdg: fPDGs) {
1250  if(fHitPDG.back()==pdg)
1251  passfilter = true;
1252  }
1253  }
1254  else
1255  fHitPDG.push_back(INT_MAX);
1256  }
1257 
1258  if(passfilter)
1259  fSimHitNtuple->Fill();
1260  }//for CRTHits
1261  }//if sim CRTHits
1262 
1263  else
1264  mf::LogWarning("CRTSimAnalysis")
1265  << "CRTHit products not found! (expected if gen/G4/detsim step)" << '\n';
1266 
1267  //true CRTHits
1268  art::Handle<vector<CRTHit>> crtTrueHitHandle;
1269  fTrueNHit = 0;
1270  if (event.getByLabel(fCRTTrueHitProducerLabel, crtTrueHitHandle)) {
1271 
1272  vector<art::Ptr<CRTHit>> crtTrueHits;
1273  art::fill_ptr_vector(crtTrueHits,crtTrueHitHandle);
1274 
1275  mf::LogPrint("CRTSimAnalysis") << "found " << crtTrueHits.size() << " true CRTHits" << '\n';
1276  for ( auto const& hit : crtTrueHits )
1277  {
1278  fTrueNHit++;
1279  fTrueXHit = hit->x_pos;
1280  fTrueYHit = hit->y_pos;
1281  fTrueZHit = hit->z_pos;
1282  fTrueXHitErr = hit->x_err;
1283  fTrueYHitErr = hit->y_err;
1284  fTrueZHitErr = hit->z_err;
1285  fTrueT0Hit = hit->ts0_ns;
1286  fTrueT1Hit = hit->ts1_ns;
1287  fTrueNHitFeb = hit->feb_id.size();
1288  fTrueHitTotPe = hit->peshit;
1289  fTrueHitPeRms = 0.;
1290  fTrueHitPe.clear();
1291  fTrueHitTrk.clear();
1292  fTrueHitPDG.clear();
1293  fTrueHitMod.clear();
1294  fTrueHitStrip.clear();
1295 
1296  uint8_t mactmp = hit->feb_id[0];
1297  fTrueHitReg = fCrtutils->AuxDetRegionNameToNum(fCrtutils->MacToRegion(mactmp));
1298  fTrueHitSubSys = fCrtutils->MacToTypeCode(mactmp);
1299 
1300  bool passfilter=false;
1301  if(fPDGs.size()==1 && fPDGs[0]==0)
1302  passfilter = true;
1303 
1304  for(auto const& mactopes : hit->pesmap){
1305  for(auto const& chanpe : mactopes.second) {
1306  fTrueHitMod.push_back(fCrtutils->MacToAuxDetID(mactopes.first, chanpe.first));
1307  fTrueHitStrip.push_back(fCrtutils->ChannelToAuxDetSensitiveID(mactopes.first,chanpe.first));
1308  fTrueHitPe.push_back(chanpe.second);
1309  fTrueHitPeRms+=pow(chanpe.second-fHitTotPe,2);
1310  }
1311  }
1312  fTrueHitPeRms = sqrt(fTrueHitPeRms/(fTrueHitPe.size()-1));
1313 
1314  //loop over all track IDs
1315  for( int trk : bt.AllTrueIds(event,*hit)) {
1316  fTrueHitTrk.push_back(trk);
1317  if ( particleMap.find(trk) != particleMap.end()){
1318  fTrueHitPDG.push_back(particleMap[trk]->PdgCode());
1319  if(!passfilter)
1320  for(int const& pdg: fPDGs) {
1321  if(fTrueHitPDG.back()==pdg)
1322  passfilter = true;
1323  }
1324  }
1325  else
1326  fTrueHitPDG.push_back(INT_MAX);
1327  }
1328 
1329  if(passfilter)
1330  fTrueCRTHitNtuple->Fill();
1331  }//for CRTHits
1332  }//if true CRTHits
1333 
1334  else
1335  mf::LogWarning("CRTSimAnalysis") << "true CRTHit products not found" << '\n';
1336 
1337  //CRTSimTrack
1338  art::Handle<vector<sbn::crt::CRTTrack>> crtSimTrackHandle;
1339  fTrueNHit = 0;
1340  if (event.getByLabel(fCRTSimTrackProducerLabel, crtSimTrackHandle)) {
1341 
1342  vector<art::Ptr<sbn::crt::CRTTrack>> crtTracks;
1343  art::fill_ptr_vector(crtTracks,crtSimTrackHandle);
1344  fNSimTrack=crtTracks.size();
1345 
1346  art::FindManyP<CRTHit> findManyHits(
1347  crtSimTrackHandle, event, fCRTSimTrackProducerLabel);
1348 
1349  mf::LogPrint("CRTSimAnalysis") << "found " << crtTracks.size() << " sim CRTTracks" << '\n';
1350  for ( size_t itrk=0; itrk<crtTracks.size(); itrk++ )
1351  {
1352  auto const& trk = crtTracks[itrk];
1353  fSimTrackPE = trk->peshit;
1354  fSimTrackT0 = (double)trk->ts0_ns;
1355  fSimTrackStart[0] = trk->x1_pos;
1356  fSimTrackStart[1] = trk->y1_pos;
1357  fSimTrackStart[2] = trk->z1_pos;
1358  fSimTrackEnd[0] = trk->x2_pos;
1359  fSimTrackEnd[1] = trk->y2_pos;
1360  fSimTrackEnd[2] = trk->z2_pos;
1361  fSimTrackL = trk->length;
1362  fSimTrackTheta = trk->thetaxy;
1363  fSimTrackPhi = trk->phizy;
1364 
1365  auto const& trkhits = findManyHits.at(itrk);
1366  fNHitSimTrack = (int)trkhits.size();
1367  uint32_t tstart=UINT32_MAX;
1368  uint32_t tend=0;
1369  size_t ihit_start=0, ihit_end=0;
1370  for(size_t ihit=0; ihit<trkhits.size(); ihit++){
1371  if(trkhits[ihit]->ts0_ns<tstart) {
1372  ihit_start = ihit;
1373  tstart = trkhits[ihit]->ts0_ns;
1374  }
1375  if(trkhits[ihit]->ts0_ns>tend){
1376  ihit_end = ihit;
1377  tend = trkhits[ihit]->ts0_ns;
1378  }
1379  }//for track hits
1380  fSimTrackHitStart[0] = trkhits[ihit_start]->x_pos;
1381  fSimTrackHitStart[1] = trkhits[ihit_start]->y_pos;
1382  fSimTrackHitStart[2] = trkhits[ihit_start]->z_pos;
1383  fSimTrackHitStart[3] = tstart-1.6e6;
1384  fSimTrackHitEnd[0] = trkhits[ihit_end]->x_pos;
1385  fSimTrackHitEnd[1] = trkhits[ihit_end]->y_pos;
1386  fSimTrackHitEnd[2] = trkhits[ihit_end]->z_pos;
1387  fSimTrackHitEnd[3] = tend-1.6e6;
1388  //fSimTrackRegStart = fCrtutils->AuxDetRegionNameToNum(trkhits[ihit_start]->tagger);
1389  //fSimTrackRegEnd = fCrtutils->AuxDetRegionNameToNum(trkhits[ihit_end]->tagger);
1390  fSimTrackRegStart = trk->plane1;
1391  fSimTrackRegEnd = trk->plane2;
1392 
1393  fSimTrackNtuple->Fill();
1394  }
1395  }//if tracks found
1396  else
1397  mf::LogWarning("CRTSimAnalysis") << "no CRTTrack products not found" << '\n';
1398 
1399  } // CRTSimAnalysis::analyze()
1400 
1401  DEFINE_ART_MODULE(CRTSimAnalysis)
1402 
1403 } // namespace crt
1404 } // namespace icarus
1405 
1406 
1407 // Back to our local namespace.
1408 namespace {
1409 
1410  int ProcessToICode(string const& p)
1411  {
1412  int icode = -1;
1413 
1414  if(p=="primary") icode=0;
1415 
1416  //gamma (PDG 22)
1417  if(p=="eBrem") icode=1;
1418  if(p=="muBrems") icode=2;
1419  if(p=="annihil") icode=3;
1420  if(p=="muMinusCaptureAtRest") icode=4;
1421  if(p=="nCapture") icode=5;
1422  if(p=="hBertiniCaptureAtRest") icode=6;
1423  if(p=="photonNuclear") icode=7;
1424  if(p=="muonNuclear") icode=8;
1425  if(p=="neutronInelastic") icode=9;
1426  if(p=="protonInelastic") icode=10;
1427  if(p=="pi-Inelastic") icode=11;
1428  if(p=="Decay") icode=12;
1429 
1430  //e (PDG +/- 11)
1431  if (p=="muIoni") icode=13;
1432  if (p=="eIoni") icode=14;
1433  if (p=="hIoni") icode=15;
1434  if (p=="compt") icode=16;
1435  if (p=="conv") icode=17;
1436  if (p=="muPairProd") icode=18;
1437  if (p=="phot") icode=19;
1438 
1439  if (p=="hadElastic") icode=20;
1440 
1441  if (p=="LArVoxelReadoutScoringProcess") icode=30;
1442  if (p=="CoupledTransportation") icode=31;
1443 
1444  return icode;
1445 
1446  }
1447 
1448 }//local namespace
LArG4Parameters.h
Store parameters for running LArG4.
geo::CryostatGeo::GetClosestOpDet
unsigned int GetClosestOpDet(geo::Point_t const &point) const
Definition: CryostatGeo.cxx:239
icarus::crt::CRTSimAnalysis::fRegEntryXYZT
vector< vector< double > > fRegEntryXYZT
Definition: CRTSimAnalysis_module.cc:230
icarus::crt::CRTSimAnalysis::fTrueXHitErr
float fTrueXHitErr
stat error of CRT hit reco X (cm)
Definition: CRTSimAnalysis_module.cc:278
icarus::crt::CRTSimAnalysis::fSimTrackID
int fSimTrackID
GEANT ID of the particle being processed.
Definition: CRTSimAnalysis_module.cc:188
ServiceUtil.h
Utilities related to art service access.
pdg
var pdg
Definition: selectors.fcl:14
icarus::crt::CRTSimAnalysis::fSimulationProducerLabel
art::InputTag fSimulationProducerLabel
The name of the producer that tracked simulated particles through the detector.
Definition: CRTSimAnalysis_module.cc:132
icarus::crt::CRTSimAnalysis::fEntry
int fEntry
front-end board entry number (reset for each event)
Definition: CRTSimAnalysis_module.cc:240
icarus::crt::CRTSimAnalysis::fXHit
float fXHit
signal inducing particle(s)&#39; PDG code
Definition: CRTSimAnalysis_module.cc:254
flag
then echo unknown compiler flag
Definition: build-dunepdsprce.sh:66
icarus::crt::CRTSimAnalysis::fZHitErr
float fZHitErr
stat error of CRT hit reco Z (cm)
Definition: CRTSimAnalysis_module.cc:259
icarus::crt::CRTSimAnalysis::fAuxDetID
vector< uint32_t > fAuxDetID
Global CRT module ID.
Definition: CRTSimAnalysis_module.cc:191
icarus::crt::CRTSimAnalysis::fT1Hit
float fT1Hit
hit time w.r.t. PPS
Definition: CRTSimAnalysis_module.cc:261
icarus::crt::CRTSimAnalysis::fNSimTrack
int fNSimTrack
number of simulated CRT tracks for this event
Definition: CRTSimAnalysis_module.cc:296
Hit.h
Declaration of signal hit object.
icarus::crt::CRTSimAnalysis::fCDxyzt
vector< vector< double > > fCDxyzt
Definition: CRTSimAnalysis_module.cc:167
icarus::crt::CRTSimAnalysis::fADdEdx
vector< double > fADdEdx
average dEdx for particle traversing CRT strip
Definition: CRTSimAnalysis_module.cc:194
p
pdgs p
Definition: selectors.fcl:22
icarus::crt::CRTSimAnalysis::fGenStartXYZT
vector< vector< double > > fGenStartXYZT
Definition: CRTSimAnalysis_module.cc:175
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
icarus::crt::CRTSimAnalysis::fYHit
float fYHit
reconstructed Y position of CRT hit (cm)
Definition: CRTSimAnalysis_module.cc:255
icarus::crt::CRTSimAnalysis::fSimulationNtuple
TTree * fSimulationNtuple
tuple with simulated data
Definition: CRTSimAnalysis_module.cc:144
geo::OpDetGeo::GetCenter
void GetCenter(double *xyz, double localz=0.0) const
Definition: OpDetGeo.cxx:40
icarus::crt::CRTSimAnalysis::fGeometryService
geo::GeometryCore const * fGeometryService
pointer to Geometry provider
Definition: CRTSimAnalysis_module.cc:313
hit
process_name hit
Definition: cheaterreco.fcl:51
geo::CryostatGeo
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
icarus::crt::CRTSimAnalysis::fNHit
int fNHit
number of CRT hits for this event
Definition: CRTSimAnalysis_module.cc:264
icarus::crt::CRTSimAnalysis::CRTSimAnalysis
CRTSimAnalysis(fhicl::ParameterSet const &p)
Constructor: configures the module (see the Config structure above)
Definition: CRTSimAnalysis_module.cc:332
icarus::crt::CRTSimAnalysis::fNAuxDet
uint32_t fNAuxDet
Number of scintillator strips hit.
Definition: CRTSimAnalysis_module.cc:189
icarus::crt::CRTSimAnalysis::analyze
virtual void analyze(const art::Event &event) override
Definition: CRTSimAnalysis_module.cc:543
icarus::crt::CRTSimAnalysis::fEndPE
float fEndPE[4]
(Px,Py,Pz,E) at the true end of the particle
Definition: CRTSimAnalysis_module.cc:211
icarus::crt::CRTSimAnalysis::fHitReg
int fHitReg
region code of CRT hit
Definition: CRTSimAnalysis_module.cc:262
icarus::crt::CRTSimAnalysis::fT0
int fT0
signal time w.r.t. global event time
Definition: CRTSimAnalysis_module.cc:244
icarus::crt::CRTSimAnalysis::beginRun
virtual void beginRun(const art::Run &run) override
Definition: CRTSimAnalysis_module.cc:535
AuxDetGeometryCore.h
Access the description of auxiliary detector geometry.
icarus::crt::CRTSimAnalysis::fTrueHitReg
int fTrueHitReg
region code of CRT hit
Definition: CRTSimAnalysis_module.cc:283
icarus::crt::CRTSimAnalysis::fRegExitSlope
vector< vector< double > > fRegExitSlope
Definition: CRTSimAnalysis_module.cc:233
lar::dump::vector
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
GeometryCore.h
Access the description of detector geometry.
icarus::crt::CRTSimAnalysis::fEvent
int fEvent
number of the event being processed
Definition: CRTSimAnalysis_module.cc:154
icarus::crt::CRTSimAnalysis::fMother
int fMother
G4 track ID of mother that directly produced this MCParticle.
Definition: CRTSimAnalysis_module.cc:214
icarus::crt::CRTSimAnalysis::fCosmicDisplayNtuple
TTree * fCosmicDisplayNtuple
for ROOT based event display
Definition: CRTSimAnalysis_module.cc:142
icarus::crt::CRTSimAnalysis::fRegExitPE
vector< vector< double > > fRegExitPE
Definition: CRTSimAnalysis_module.cc:229
AuxDetSimChannel.h
object containing MC truth information necessary for making RawDigits and doing back tracking ...
icarus::crt::CRTSimAnalysis::fNDaught
int fNDaught
number of daughters belonging to this MCParticle
Definition: CRTSimAnalysis_module.cc:215
icarus::crt::CRTSimAnalysis::fTrueZHitErr
float fTrueZHitErr
stat error of CRT hit reco Z (cm)
Definition: CRTSimAnalysis_module.cc:280
icarus::crt::CRTSimAnalysis::fSimProcess
int fSimProcess
process that created the particle (e.g. brehmstralung)
Definition: CRTSimAnalysis_module.cc:186
icarus::crt::CRTSimAnalysis::fTrueYHitErr
float fTrueYHitErr
stat error of CRT hit reco Y (cm)
Definition: CRTSimAnalysis_module.cc:279
icarus::crt::CRTSimAnalysis::fSimHits
uint32_t fSimHits
number of trajectory points for each MCParticle
Definition: CRTSimAnalysis_module.cc:183
geo::CryostatGeo::OpDet
const OpDetGeo & OpDet(unsigned int iopdet) const
Return the iopdet&#39;th optical detector in the cryostat.
Definition: CryostatGeo.cxx:107
icarus::crt::CRTSimAnalysis::fProgenitor
int fProgenitor
G4 track ID of the primary particle that ultimately led to this one.
Definition: CRTSimAnalysis_module.cc:213
icarus::crt::CRTSimAnalysis::fRun
int fRun
number of the run being processed
Definition: CRTSimAnalysis_module.cc:155
icarus::crt::CRTSimAnalysis::fTrueYHit
float fTrueYHit
reconstructed Y position of CRT hit (cm)
Definition: CRTSimAnalysis_module.cc:276
icarus::crt::CRTSimAnalysis::fTriggerOffset
int fTriggerOffset
(units of ticks) time of expected neutrino event
Definition: CRTSimAnalysis_module.cc:314
icarus::crt::CRTSimAnalysis::fMac5
int fMac5
Mac5 address for this front-end board.
Definition: CRTSimAnalysis_module.cc:242
icarus::crt::CRTSimAnalysis::fADEnterXYZT
vector< vector< double > > fADEnterXYZT
4-position of entry into CRT strip
Definition: CRTSimAnalysis_module.cc:199
icarus::crt::CRTSimAnalysis::fPDGs
vector< int > fPDGs
PDG code of particle we&#39;ll focus on.
Definition: CRTSimAnalysis_module.cc:137
icarus::crt::CRTSimAnalysis::beginJob
virtual void beginJob() override
Definition: CRTSimAnalysis_module.cc:351
icarus::crt::CRTSimAnalysis::fADMac
vector< uint32_t > fADMac
Mac5 address of the CRT module.
Definition: CRTSimAnalysis_module.cc:197
icarus::crt::CRTSimAnalysis::fADTrackLength
vector< double > fADTrackLength
Track length in CRT strip (cm)
Definition: CRTSimAnalysis_module.cc:195
geo::GeometryCore::Cryostat
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
Definition: GeometryCore.cxx:321
icarus::crt::CRTSimAnalysis::fTrackLength
float fTrackLength
total track length for each MCParticle
Definition: CRTSimAnalysis_module.cc:184
icarus::crt::CRTSimAnalysis::fStartPE
float fStartPE[4]
(Px,Py,Pz,E) at the true start of the particle
Definition: CRTSimAnalysis_module.cc:210
icarus::crt::CRTSimAnalysis::fADEDep
vector< double > fADEDep
Energy deposited in CRT strip (GeV)
Definition: CRTSimAnalysis_module.cc:193
icarus::crt::CRTSimAnalysis::fGenEndXYZT
vector< vector< double > > fGenEndXYZT
Definition: CRTSimAnalysis_module.cc:176
icarus::crt::CRTSimAnalysis::fCDpe
vector< vector< double > > fCDpe
4-momentum
Definition: CRTSimAnalysis_module.cc:166
icarus::crt::CRTSimAnalysis::fSimPDG
int fSimPDG
PDG ID of the particle being processed.
Definition: CRTSimAnalysis_module.cc:185
icarus::crt::CRTSimAnalysis::fZHit
float fZHit
reconstructed Z position of CRT hit (cm)
Definition: CRTSimAnalysis_module.cc:256
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1476
Cluster.h
Declaration of cluster object.
geo_types.h
Definition of data types for geometry description.
icarus::crt::CRTSimAnalysis::fCDSlopes
vector< vector< double > > fCDSlopes
direction cosines
Definition: CRTSimAnalysis_module.cc:165
icarus::crt::CRTSimAnalysis::fADC
int fADC[64]
signal amplitude
Definition: CRTSimAnalysis_module.cc:246
icarus::crt::CRTSimAnalysis::fRegEntryPE
vector< vector< double > > fRegEntryPE
Definition: CRTSimAnalysis_module.cc:228
icarus::crt::CRTSimAnalysis::fTrueNHit
uint32_t fTrueNHit
number of CRT hits for this event
Definition: CRTSimAnalysis_module.cc:285
icarus::crt::CRTSimAnalysis::fADEnterPE
vector< vector< double > > fADEnterPE
4-position of entry into CRT strip
Definition: CRTSimAnalysis_module.cc:201
geo::BoxBoundedGeo::InFiducialY
bool InFiducialY(double y, double neg_margin, double pos_margin) const
Returns whether TPC fiducial volume contains world y coordinate.
Definition: BoxBoundedGeo.h:272
icarus::crt::CRTSimAnalysis::fT1
int fT1
signal time w.r.t. PPS
Definition: CRTSimAnalysis_module.cc:245
icarus::crt::CRTSimAnalysis::fTrueXHit
float fTrueXHit
reconstructed X position of CRT hit (cm)
Definition: CRTSimAnalysis_module.cc:275
icarus::crt::CRTSimAnalysis::fFEBReg
int fFEBReg
CRT region for this front-end board.
Definition: CRTSimAnalysis_module.cc:241
icarus::crt::CRTSimAnalysis::fAuxDetReg
vector< uint32_t > fAuxDetReg
CRT region code.
Definition: CRTSimAnalysis_module.cc:196
icarus::crt::CRTSimAnalysis::fEndXYZT
float fEndXYZT[4]
(x,y,z,t) of the true end of the particle
Definition: CRTSimAnalysis_module.cc:209
icarus::crt::CRTSimAnalysis::fXHitErr
float fXHitErr
stat error of CRT hit reco X (cm)
Definition: CRTSimAnalysis_module.cc:257
icarus::crt::CRTSimAnalysis::fTrackID
vector< int > fTrackID
track ID(s) of particle that produced the signal
Definition: CRTSimAnalysis_module.cc:250
geo::CryostatGeo::TPC
const TPCGeo & TPC(unsigned int itpc) const
Return the itpc&#39;th TPC in the cryostat.
Definition: CryostatGeo.cxx:93
icarus::crt::CRTSimAnalysis::fYHitErr
float fYHitErr
stat error of CRT hit reco Y (cm)
Definition: CRTSimAnalysis_module.cc:258
icarus::crt::CRTSimAnalysis::fADExitXYZT
vector< vector< double > > fADExitXYZT
4-position of exit from CRT strip
Definition: CRTSimAnalysis_module.cc:200
icarus::crt::CRTSimAnalysis::fRegEntrySlope
vector< vector< double > > fRegEntrySlope
Definition: CRTSimAnalysis_module.cc:232
icarus::crt::CRTSimAnalysis::fSimEndProcess
int fSimEndProcess
process the killed the particle (e.g. annihilation)
Definition: CRTSimAnalysis_module.cc:187
SimChannel.h
object containing MC truth information necessary for making RawDigits and doing back tracking ...
icarus::crt::CRTBackTracker::AllTrueIds
std::vector< int > AllTrueIds(const art::Event &event, const CRTData &data)
Definition: icaruscode/icaruscode/CRT/CRTUtils/CRTBackTracker.cc:208
icarus::crt::CRTSimAnalysis::fTrueT1Hit
float fTrueT1Hit
hit time w.r.t. global event time
Definition: CRTSimAnalysis_module.cc:282
icarus::crt::CRTSimAnalysis::fStartXYZT
float fStartXYZT[4]
(x,y,z,t) of the true start of the particle
Definition: CRTSimAnalysis_module.cc:208
icarus::crt::CRTSimAnalysis::fCRTSimHitProducerLabel
art::InputTag fCRTSimHitProducerLabel
The name of the producer that created hits.
Definition: CRTSimAnalysis_module.cc:133
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
icarus::crt::CRTSimAnalysis::fT0Hit
float fT0Hit
hit time w.r.t. global event time
Definition: CRTSimAnalysis_module.cc:260
geo::BoxBoundedGeo::InFiducialX
bool InFiducialX(double x, double neg_margin, double pos_margin) const
Returns whether TPC fiducial volume contains world x coordinate.
Definition: BoxBoundedGeo.h:241
k
pdgs k
Definition: selectors.fcl:22
icarus::crt::CRTSimAnalysis::fTrueZHit
float fTrueZHit
reconstructed Z position of CRT hit (cm)
Definition: CRTSimAnalysis_module.cc:277
geo::BoxBoundedGeo::InFiducialZ
bool InFiducialZ(double z, double neg_margin, double pos_margin) const
Returns whether TPC fiducial volume contains world z coordinate.
Definition: BoxBoundedGeo.h:303
detinfo::sampling_rate
double sampling_rate(DetectorClocksData const &data)
Returns the period of the TPC readout electronics clock.
Definition: DetectorClocksData.h:643
icarus::crt::CRTSimAnalysis::fADExitPE
vector< vector< double > > fADExitPE
4-position of exit from CRT strip
Definition: CRTSimAnalysis_module.cc:202
icarus::crt::CRTSimAnalysis::fTrueT0Hit
float fTrueT0Hit
hit time w.r.t. global event time
Definition: CRTSimAnalysis_module.cc:281
crt
process_name crt
Definition: run_genericCRT.fcl:82
icarus::crt::CRTSimAnalysis::fGenStartPE
vector< vector< double > > fGenStartPE
Definition: CRTSimAnalysis_module.cc:177
icarus::crt::CRTSimAnalysis::fRegExitXYZT
vector< vector< double > > fRegExitXYZT
Definition: CRTSimAnalysis_module.cc:231
geo::BoxBoundedGeo::ContainsPosition
bool ContainsPosition(geo::Point_t const &point, double wiggle=1.0) const
Returns whether this volume contains the specified point.
Definition: BoxBoundedGeo.h:205
Geometry.h
art framework interface to geometry description
cout
BEGIN_PROLOG could also be cout
Definition: messageservice.fcl:10
icarus::crt::CRTSimAnalysis::fRegOpDetXYZT
vector< vector< double > > fRegOpDetXYZT
Definition: CRTSimAnalysis_module.cc:236
icarus::crt::CRTSimAnalysis::fAuxDetSensitiveID
vector< uint32_t > fAuxDetSensitiveID
Strip ID in module.
Definition: CRTSimAnalysis_module.cc:192
icarus::crt::CRTSimAnalysis::fGenEndPE
vector< vector< double > > fGenEndPE
Definition: CRTSimAnalysis_module.cc:178
Generated by   doxygen 1.8.5