All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ThroughgoingmuonAnalyzer_module.cc
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////////////
2 // Class: ThroughgoingmuonAnalyzer
3 // Plugin Type: analyzer (art v3_02_06)
4 // File: ThroughgoingmuonAnalyzer_module.cc
5 //
6 // This is based on the Tracy's TrackHitEfficiency module.
7 //
8 // Generated at Thu Nov 14 15:21:00 2019 by Biswaranjan Behera using cetskelgen
9 // from cetlib version v3_07_02.
10 /////////////////////////////////////////////////////////////////////////////////
11 
12 #include "art/Framework/Core/EDAnalyzer.h"
13 #include "art/Framework/Core/ModuleMacros.h"
14 #include "art/Framework/Principal/Event.h"
15 #include "art/Framework/Principal/Handle.h"
16 #include "art/Framework/Principal/Run.h"
17 #include "art/Framework/Principal/SubRun.h"
18 #include "canvas/Utilities/InputTag.h"
19 #include "fhiclcpp/ParameterSet.h"
20 #include "messagefacility/MessageLogger/MessageLogger.h"
21 
22 #include "canvas/Persistency/Common/FindMany.h"
23 #include "fhiclcpp/ParameterSet.h"
24 #include "art/Utilities/ToolMacros.h"
25 #include "art/Framework/Services/Registry/ServiceHandle.h"
26 #include "art_root_io/TFileService.h"
27 #include "art/Framework/Core/ModuleMacros.h"
28 #include "art_root_io/TFileDirectory.h"
29 
30 
31 #include "lardata/DetectorInfoServices/DetectorClocksService.h"
32 #include "larcore/Geometry/Geometry.h"
33 #include "larcore/CoreUtils/ServiceUtil.h" // lar::providerFrom()
34 #include "larcorealg/Geometry/CryostatGeo.h"
35 #include "larevt/CalibrationDBI/Interface/ChannelStatusService.h"
36 #include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h"
37 #include "lardataobj/RawData/RawDigit.h"
38 #include "lardataobj/RecoBase/Wire.h"
39 #include "lardataobj/RecoBase/Hit.h"
40 #include "lardataobj/RecoBase/Track.h"
41 #include "lardataobj/Simulation/SimChannel.h"
42 #include "lardataobj/Simulation/SimEnergyDeposit.h"
43 #include "nusimdata/SimulationBase/MCParticle.h"
44 #include "larsim/Simulation/LArG4Parameters.h"
45 
46 // Eigen
47 #include <Eigen/Dense>
48 
49 //Root
50 #include "TH1.h"
51 #include "TH2.h"
52 #include "TProfile.h"
53 #include "TProfile2D.h"
54 #include "TTree.h"
55 
56 #include <cmath>
57 #include <algorithm>
58 
59 namespace thrugoingmuon {
60 
61 class ThroughgoingmuonAnalyzer : public art::EDAnalyzer
62 {
63 public:
64  explicit ThroughgoingmuonAnalyzer(fhicl::ParameterSet const& pset);
65  // The compiler-generated destructor is fine for non-base
66  // classes without bare pointers or other resource use.
67 
68  // Plugins should not be copied or assigned.
69  ThroughgoingmuonAnalyzer(ThroughgoingmuonAnalyzer const&) = delete;
70  ThroughgoingmuonAnalyzer(ThroughgoingmuonAnalyzer&&) = delete;
71  ThroughgoingmuonAnalyzer& operator=(ThroughgoingmuonAnalyzer const&) = delete;
72  ThroughgoingmuonAnalyzer& operator=(ThroughgoingmuonAnalyzer&&) = delete;
73 
74  // Required functions.
75  void analyze(art::Event const& evt) override;
76 
77 private:
78 
79 
80  // Declare member data here.
81  // Fcl parameters.
82  std::vector<art::InputTag> fRawDigitProducerLabelVec;
83  std::vector<art::InputTag> fWireProducerLabelVec;
84  std::vector<art::InputTag> fHitProducerLabelVec;
85  std::vector<art::InputTag> fTrackProducerLabelVec;
86  art::InputTag fMCParticleProducerLabel;
87  art::InputTag fSimChannelProducerLabel;
88  art::InputTag fSimEnergyProducerLabel;
89  art::InputTag fBadChannelProducerLabel;
90  bool fUseBadChannelDB;
91  bool fUseICARUSGeometry;
92  std::vector<int> fOffsetVec; ///< Allow offsets for each plane
93  std::vector<float> fSigmaVec; ///< Window size for matching to SimChannels
94  int fMinAllowedChanStatus; ///< Don't consider channels with lower status
95  float ftotalelctroncut; ///< Threshold for electron (10k)
96 
97  // double fElectronsToGeV; ///< conversion factor
98  // TTree variables
99  TTree* fTree;
100  int fEvent; ///< number of the event being processed
101  int fRun; ///< number of the run being processed
102  int fSubRun; ///< number of the sub-run being processed
103  /// @}
104 
105  /// @name The variables that will go into the simulation n-tuple.
106  /// @{
107  int fSimPDG; ///< PDG ID of the particle being processed
108  int fSimTrackID; ///< GEANT ID of the particle being processed
109  int fNtracks_reco;
110  int fNtracks_primary;
111 
112  //mutable std::vector<float> fHitSummedADCVec;
113  mutable std::vector<int> fTPCVec;
114  mutable std::vector<int> fCryoVec;
115  mutable std::vector<int> fPlaneVec;
116  mutable std::vector<int> fWireVec;
117 
118  mutable std::vector<float> fmcpartvx;
119  mutable std::vector<float> fmcpartvy;
120  mutable std::vector<float> fmcpartvz;
121 
122  mutable std::vector<float> fmcpartpx;
123  mutable std::vector<float> fmcpartpy;
124  mutable std::vector<float> fmcpartpz;
125  mutable std::vector<float> fmcparte;
126 
127  mutable std::vector<float> fTotalElectronsVec;
128  mutable std::vector<float> fTotalElecEnergyVec;
129  mutable std::vector<float> fMaxElectronsVec;
130  mutable std::vector<int> fStartTickVec;
131  mutable std::vector<int> fStopTickVec;
132  mutable std::vector<int> fMaxETickVec;
133  mutable std::vector<float> fPartDirX;
134  mutable std::vector<float> fPartDirY;
135  mutable std::vector<float> fPartDirZ;
136 
137  mutable std::vector<int> fNMatchedWires;
138  mutable std::vector<int> fNMatchedHits;
139 
140  mutable std::vector<float> fHitPeakTimeVec;
141  mutable std::vector<float> fHitPeakAmpVec;
142  mutable std::vector<float> fHitPeakRMSVec;
143  mutable std::vector<float> fHitBaselinevec;
144  mutable std::vector<float> fHitSummedADCVec;
145  mutable std::vector<float> fHitIntegralVec;
146  mutable std::vector<int> fHitStartTickVec;
147  mutable std::vector<int> fHitStopTickVec;
148  mutable std::vector<int> fHitMultiplicityVec;
149  mutable std::vector<int> fHitLocalIndexVec;
150  mutable std::vector<float> fHitGoodnessVec;
151  mutable std::vector<int> fNumDegreesVec;
152 
153  mutable std::vector<float> fmcstartx;
154  mutable std::vector<float> fmcstarty;
155  mutable std::vector<float> fmcstartz;
156 
157  mutable std::vector<float> fmcendx;
158  mutable std::vector<float> fmcendy;
159  mutable std::vector<float> fmcendz;
160 
161  mutable std::vector<float> fmcstartdirx;
162  mutable std::vector<float> fmcstartdiry;
163  mutable std::vector<float> fmcstartdirz;
164 
165  mutable std::vector<float> fmcenddirx;
166  mutable std::vector<float> fmcenddiry;
167  mutable std::vector<float> fmcenddirz;
168 
169  mutable std::vector<float> fmctstartx;
170  mutable std::vector<float> fmctstarty;
171  mutable std::vector<float> fmctstartz;
172 
173  mutable std::vector<float> fmctendx;
174  mutable std::vector<float> fmctendy;
175  mutable std::vector<float> fmctendz;
176 
177  mutable std::vector<float> fmctstartdirx;
178  mutable std::vector<float> fmctstartdiry;
179  mutable std::vector<float> fmctstartdirz;
180 
181  mutable std::vector<float> fmctenddirx;
182  mutable std::vector<float> fmctenddiry;
183  mutable std::vector<float> fmctenddirz;
184 
185  mutable std::vector<float> frecostartx;
186  mutable std::vector<float> frecostarty;
187  mutable std::vector<float> frecostartz;
188 
189  mutable std::vector<float> frecoendx;
190  mutable std::vector<float> frecoendy;
191  mutable std::vector<float> frecoendz;
192 
193  mutable std::vector<float> frecostartdirx;
194  mutable std::vector<float> frecostartdiry;
195  mutable std::vector<float> frecostartdirz;
196 
197  mutable std::vector<float> frecoenddirx;
198  mutable std::vector<float> frecoenddiry;
199  mutable std::vector<float> frecoenddirz;
200 
201  mutable std::vector<float> fLength;
202  mutable std::vector<float> fThetaXZ;
203  mutable std::vector<float> fThetaYZ;
204  mutable std::vector<float> fTheta;
205  mutable std::vector<float> fPhi;
206 
207  mutable std::vector<float> fmcThetaXZ;
208  mutable std::vector<float> fmcThetaYZ;
209  mutable std::vector<float> fmcTheta;
210  mutable std::vector<float> fmcPhi;
211 
212  mutable std::vector<float> fmctThetaXZ;
213  mutable std::vector<float> fmctThetaYZ;
214  mutable std::vector<float> fmctTheta;
215  mutable std::vector<float> fmctPhi;
216 
217  mutable std::vector<float> fNSimChannelHitsVec;
218  mutable std::vector<float> fNRecobHitVec;
219  mutable std::vector<float> fHitEfficiencyVec;
220  mutable std::vector<float> fNFakeHitVec;
221 
222  mutable std::vector<float> flocx;
223  mutable std::vector<float> flocy;
224  mutable std::vector<float> flocz;
225 
226  std::vector<TProfile*> fHitEffvselecVec;
227 
228  // Implementation of required member function here.
229  art::ServiceHandle< art::TFileService > tfs;
230 
231  // How to convert from number of electrons to GeV. The ultimate
232  // source of this conversion factor is
233  // ${LARCOREOBJ_INC}/larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h.
234  // But sim::LArG4Parameters might in principle ask a database for it.
235  //art::ServiceHandle<sim::LArG4Parameters const> larParameters;
236 
237 
238  // Useful services, keep copies for now (we can update during begin run periods)
239  const geo::GeometryCore* fGeometry; ///< pointer to Geometry service
240 
241  // Get geometry.
242  // art::ServiceHandle<geo::Geometry> geom;
243 };
244 
245 
246 ThroughgoingmuonAnalyzer::ThroughgoingmuonAnalyzer(fhicl::ParameterSet const& pset)
247  : EDAnalyzer{pset} // ,
248 // More initializers here.
249 {
250  std::cout << "In analyzer constructor" << std::endl;
251 
252  // Call appropriate consumes<>() for any products to be retrieved by this module.
253  fRawDigitProducerLabelVec = pset.get< std::vector<art::InputTag>>("RawDigitLabelVec", std::vector<art::InputTag>() = {"rawdigitfilter"});
254  fWireProducerLabelVec = pset.get< std::vector<art::InputTag>>("WireModuleLabelVec", std::vector<art::InputTag>() = {"decon1droi"});
255  fHitProducerLabelVec = pset.get< std::vector<art::InputTag>>("HitModuleLabelVec", std::vector<art::InputTag>() = {"gaushit"});
256  fTrackProducerLabelVec = pset.get< std::vector<art::InputTag>>("TrackModuleLabelVec",std::vector<art::InputTag>() = {"pandoraTrackGausCryo0"});
257  fMCParticleProducerLabel = pset.get< art::InputTag >("MCParticleLabel", "largeant");
258  fSimChannelProducerLabel = pset.get< art::InputTag >("SimChannelLabel", "largeant");
259  fSimEnergyProducerLabel = pset.get< art::InputTag >("SimEnergyLabel", "ionization");
260  fBadChannelProducerLabel = pset.get< art::InputTag >("BadChannelLabel", "simnfspl1:badchannels");
261  fUseBadChannelDB = pset.get< bool >("UseBadChannelDB", false);
262  fUseICARUSGeometry = pset.get< bool >("UseICARUSGeometry", false);
263  // fLocalDirName = pset.get<std::string >("LocalDirName", std::string("wow"));
264  fOffsetVec = pset.get<std::vector<int> >("OffsetVec", std::vector<int>()={0,0,0});
265  fSigmaVec = pset.get<std::vector<float> >("SigmaVec", std::vector<float>()={1.,1.,1.});
266  fMinAllowedChanStatus = pset.get< int >("MinAllowedChannelStatus");
267  ftotalelctroncut = pset.get<float >("totalelctroncut", 10000.);
268 
269  std::cout << "UseICARUSGeometry: " << fUseICARUSGeometry << std::endl;
270  std::cout << "Read the fhicl parameters, recovering services" << std::endl;
271 
272  // fHitProducerLabelVec = pset.get< std::vector<art::InputTag>>("HitModuleLabelVec", std::vector<art::InputTag>() = {"gauss"});
273 
274  fGeometry = lar::providerFrom<geo::Geometry>();
275 
276  art::TFileDirectory dir = tfs->mkdir("histos");
277  fHitEffvselecVec.resize(fGeometry->Nplanes());
278  for(size_t plane = 0; plane < fGeometry->Nplanes(); plane++)
279  {
280  fHitEffvselecVec.at(plane) = dir.make<TProfile> (("HitEffvselec" + std::to_string(plane)).c_str(), "Hit Efficiency; Total # e-", 200, 0., 200000., 0., 1.);
281  }
282 
283  std::cout << "Services recovered, setting up output tree" << std::endl;
284 
285  fTree = tfs->make<TTree>("Analysis", "Analysis tree");
286 
287  fTree->Branch("Event", &fEvent, "Event/I");
288  fTree->Branch("SubRun", &fSubRun, "SubRun/I");
289  fTree->Branch("Run", &fRun, "Run/I");
290  fTree->Branch("TrackID", &fSimTrackID, "TrackID/I");
291  fTree->Branch("PDG", &fSimPDG, "PDG/I");
292  fTree->Branch("ntracks_reco", &fNtracks_reco, "ntracks_reco/I");
293  fTree->Branch("ntracks_primary", &fNtracks_primary, "ntracks_primary/I");
294 
295  fTree->Branch("CryostataVec", "std::vector<int>", &fCryoVec);
296  fTree->Branch("TPCVec", "std::vector<int>", &fTPCVec);
297  fTree->Branch("PlaneVec", "std::vector<int>", &fPlaneVec);
298  fTree->Branch("WireVec", "std::vector<int>", &fWireVec);
299 
300  fTree->Branch("mcpartvx", "std::vector<float>", &fmcpartvx);
301  fTree->Branch("mcpartvy", "std::vector<float>", &fmcpartvy);
302  fTree->Branch("mcpartvz", "std::vector<float>", &fmcpartvz);
303 
304  fTree->Branch("mcpartpx", "std::vector<float>", &fmcpartpx);
305  fTree->Branch("mcpartpy", "std::vector<float>", &fmcpartpy);
306  fTree->Branch("mcpartpz", "std::vector<float>", &fmcpartpz);
307  fTree->Branch("mcparte", "std::vector<float>", &fmcparte);
308 
309  fTree->Branch("TotalElectronsVec", "std::vector<float>", &fTotalElectronsVec);
310  fTree->Branch("TotalElecEnergyVec","std::vector<float>", &fTotalElecEnergyVec);
311  fTree->Branch("MaxElectronsVec", "std::vector<float>", &fMaxElectronsVec);
312  fTree->Branch("StartTick", "std::vector<int>", &fStartTickVec);
313  fTree->Branch("StopTick", "std::vector<int>", &fStopTickVec);
314  fTree->Branch("MaxETick", "std::vector<int>", &fMaxETickVec);
315  fTree->Branch("PartDirX", "std::vector<float>", &fPartDirX);
316  fTree->Branch("PartDirY", "std::vector<float>", &fPartDirY);
317  fTree->Branch("PartDirZ", "std::vector<float>", &fPartDirZ);
318 
319  fTree->Branch("NMatchedWires", "std::vector<int>", &fNMatchedWires);
320  fTree->Branch("NMatchedHits", "std::vector<int>", &fNMatchedHits);
321 
322  fTree->Branch("HitPeakTimeVec", "std::vector<float>", &fHitPeakTimeVec);
323  fTree->Branch("HitPeakAmpVec", "std::vector<float>", &fHitPeakAmpVec);
324  fTree->Branch("HitPeakRMSVec", "std::vector<float>", &fHitPeakRMSVec);
325  fTree->Branch("HitBaselineVec", "std::vector<float>", &fHitBaselinevec);
326  fTree->Branch("HitSummedADCVec", "std::vector<float>", &fHitSummedADCVec);
327  fTree->Branch("HitIntegralVec", "std::vector<float>", &fHitIntegralVec);
328  fTree->Branch("HitStartTickVec", "std::vector<int>", &fHitStartTickVec);
329  fTree->Branch("HitStopTickVec", "std::vector<int>", &fHitStopTickVec);
330  fTree->Branch("HitMultiplicity", "std::vector<int>", &fHitMultiplicityVec);
331  fTree->Branch("HitLocalIndex", "std::vector<int>", &fHitLocalIndexVec);
332  fTree->Branch("HitGoodness", "std::vector<float>", &fHitGoodnessVec);
333  fTree->Branch("HitNumDegrees", "std::vector<int>", &fNumDegreesVec);
334 
335  fTree->Branch("NSimChannelHitsVec", "std::vector<float>", &fNSimChannelHitsVec);
336  fTree->Branch("NRecobHitVec", "std::vector<float>", &fNRecobHitVec);
337  fTree->Branch("HitEfficiencyVec", "std::vector<float>", &fHitEfficiencyVec);
338  fTree->Branch("NFakeHitVec", "std::vector<float>", &fNFakeHitVec);
339 
340  fTree->Branch("mcstartx", "std::vector<float>", &fmcstartx);
341  fTree->Branch("mcstarty", "std::vector<float>", &fmcstarty);
342  fTree->Branch("mcstartz", "std::vector<float>", &fmcstartz);
343  fTree->Branch("mcendx", "std::vector<float>", &fmcendx);
344  fTree->Branch("mcendy", "std::vector<float>", &fmcendy);
345  fTree->Branch("mcendz", "std::vector<float>", &fmcendz);
346 
347  fTree->Branch("mcstartdirx", "std::vector<float>", &fmcstartdirx);
348  fTree->Branch("mcstartdiry", "std::vector<float>", &fmcstartdiry);
349  fTree->Branch("mcstartdirz", "std::vector<float>", &fmcstartdirz);
350  fTree->Branch("mcenddirx", "std::vector<float>", &fmcenddirx);
351  fTree->Branch("mcenddiry", "std::vector<float>", &fmcenddiry);
352  fTree->Branch("mcenddirz", "std::vector<float>", &fmcenddirz);
353 
354  fTree->Branch("mctstartx", "std::vector<float>", &fmctstartx);
355  fTree->Branch("mctstarty", "std::vector<float>", &fmctstarty);
356  fTree->Branch("mctstartz", "std::vector<float>", &fmctstartz);
357  fTree->Branch("mctendx", "std::vector<float>", &fmctendx);
358  fTree->Branch("mctendy", "std::vector<float>", &fmctendy);
359  fTree->Branch("mctendz", "std::vector<float>", &fmctendz);
360 
361  fTree->Branch("mctstartdirx", "std::vector<float>", &fmctstartdirx);
362  fTree->Branch("mctstartdiry", "std::vector<float>", &fmctstartdiry);
363  fTree->Branch("mctstartdirz", "std::vector<float>", &fmctstartdirz);
364  fTree->Branch("mctenddirx", "std::vector<float>", &fmctenddirx);
365  fTree->Branch("mctenddiry", "std::vector<float>", &fmctenddiry);
366  fTree->Branch("mctenddirz", "std::vector<float>", &fmctenddirz);
367 
368  fTree->Branch("recostartx", "std::vector<float>", &frecostartx);
369  fTree->Branch("recostarty", "std::vector<float>", &frecostarty);
370  fTree->Branch("recostartz", "std::vector<float>", &frecostartz);
371  fTree->Branch("recoendx", "std::vector<float>", &frecoendx);
372  fTree->Branch("recoendy", "std::vector<float>", &frecoendy);
373  fTree->Branch("recoendz", "std::vector<float>", &frecoendz);
374 
375  fTree->Branch("recostartdirx", "std::vector<float>", &frecostartdirx);
376  fTree->Branch("recostartdiry", "std::vector<float>", &frecostartdiry);
377  fTree->Branch("recostartdirz", "std::vector<float>", &frecostartdirz);
378  fTree->Branch("recoenddirx", "std::vector<float>", &frecoenddirx);
379  fTree->Branch("recoenddiry", "std::vector<float>", &frecoenddiry);
380  fTree->Branch("recoenddirz", "std::vector<float>", &frecoenddirz);
381 
382  fTree->Branch("length", "std::vector<float>", &fLength);
383  fTree->Branch("thetaxz", "std::vector<float>", &fThetaXZ);
384  fTree->Branch("thetayz", "std::vector<float>", &fThetaYZ);
385  fTree->Branch("theta", "std::vector<float>", &fTheta);
386  fTree->Branch("phi", "std::vector<float>", &fPhi);
387  fTree->Branch("mctheta", "std::vector<float>", &fmcTheta);
388  fTree->Branch("mcphi", "std::vector<float>", &fmcPhi);
389  fTree->Branch("mcthetaxz", "std::vector<float>", &fmcThetaXZ);
390  fTree->Branch("mcthetayz", "std::vector<float>", &fmcThetaYZ);
391  fTree->Branch("mcttheta", "std::vector<float>", &fmctTheta);
392  fTree->Branch("mctphi", "std::vector<float>", &fmctPhi);
393  fTree->Branch("mctthetaxz", "std::vector<float>", &fmctThetaXZ);
394  fTree->Branch("mctthetayz", "std::vector<float>", &fmctThetaYZ);
395 
396  fTree->Branch("locx", "std::vector<float>", &flocx);
397  fTree->Branch("locy", "std::vector<float>", &flocy);
398  fTree->Branch("locz", "std::vector<float>", &flocz);
399 
400  std::cout << "Returning" << std::endl;
401  return;
402 }
403 
404 void ThroughgoingmuonAnalyzer::analyze(art::Event const& evt)
405 {
406  // Implementation of required member function here.
407  // Always clear the tuple
408  // clear();
409 
410  fNtracks_reco = 0;
411  fNtracks_primary = 0;
412 
413  fTPCVec.clear();
414  fCryoVec.clear();
415  fPlaneVec.clear();
416  fWireVec.clear();
417 
418  fmcpartvx.clear();
419  fmcpartvy.clear();
420  fmcpartvz.clear();
421 
422  fmcpartpx.clear();
423  fmcpartpy.clear();
424  fmcpartpz.clear();
425  fmcparte.clear();
426 
427  fTotalElectronsVec. clear();
428  fTotalElecEnergyVec.clear();
429  fMaxElectronsVec. clear();
430  fStartTickVec. clear();
431  fStopTickVec. clear();
432  fMaxETickVec. clear();
433  fPartDirX. clear();
434  fPartDirY. clear();
435  fPartDirZ. clear();
436 
437  fNMatchedWires. clear();
438  fNMatchedHits. clear();
439 
440  fHitPeakTimeVec. clear();
441  fHitPeakAmpVec. clear();
442  fHitPeakRMSVec. clear();
443  fHitBaselinevec. clear();
444  fHitSummedADCVec. clear();
445  fHitIntegralVec. clear();
446  fHitStartTickVec. clear();
447  fHitStopTickVec. clear();
448  fHitMultiplicityVec.clear();
449  fHitLocalIndexVec. clear();
450  fHitGoodnessVec. clear();
451  fNumDegreesVec. clear();
452  fmcstartx. clear();
453  fmcstarty. clear();
454  fmcstartz. clear();
455  fmcendx. clear();
456  fmcendy. clear();
457  fmcendz. clear();
458  fmcstartdirx. clear();
459  fmcstartdiry. clear();
460  fmcstartdirz. clear();
461  fmctstartx. clear();
462  fmctstarty. clear();
463  fmctstartz. clear();
464  fmctendx. clear();
465  fmctendy. clear();
466  fmctendz. clear();
467  fmctstartdirx. clear();
468  fmctstartdiry. clear();
469  fmctstartdirz. clear();
470 
471  flocx. clear();
472  flocy. clear();
473  flocz. clear();
474 
475  fmcenddirx.clear();
476  fmcenddiry.clear();
477  fmcenddirz.clear();
478  frecostartx.clear();
479  frecostarty.clear();
480  frecostartz.clear();
481  frecoendx.clear();
482  frecoendy.clear();
483  frecoendz.clear();
484  frecostartdirx.clear();
485  frecostartdiry.clear();
486  frecostartdirz.clear();
487  frecoenddirx. clear();
488  frecoenddiry. clear();
489  frecoenddirz. clear();
490  fLength. clear();
491  fThetaXZ.clear();
492  fThetaYZ.clear();
493  fTheta. clear();
494  fPhi. clear();
495  fmcThetaXZ.clear();
496  fmcThetaYZ.clear();
497  fmcTheta. clear();
498  fmcPhi. clear();
499  fmctThetaXZ.clear();
500  fmctThetaYZ.clear();
501  fmctTheta. clear();
502  fmctPhi. clear();
503 
504  fNSimChannelHitsVec.clear();
505  fNRecobHitVec. clear();
506  fHitEfficiencyVec. clear();
507  fNFakeHitVec. clear();
508 
509  fHitEffvselecVec. clear();
510 
511  //*/
512  // Start by fetching some basic event information for our n-tuple.
513  fEvent = evt.id().event();
514  fRun = evt.run();
515  fSubRun = evt.subRun();
516 
517  int ntrk_reco = 0;
518 
519  art::Handle< std::vector<sim::SimChannel>> simChannelHandle;
520  evt.getByLabel(fSimChannelProducerLabel, simChannelHandle);
521 
522 
523  art::Handle< std::vector<simb::MCParticle>> mcParticleHandle;
524  evt.getByLabel(fMCParticleProducerLabel, mcParticleHandle);
525 
526  art::Handle<std::vector<sim::SimEnergyDeposit>> simEnergyHandle;
527  evt.getByLabel(fSimEnergyProducerLabel, simEnergyHandle);
528 
529  //fElectronsToGeV = 1./larParameters->GeVToElectrons();
530 
531  // If there is no sim channel informaton then exit
532  if (!simChannelHandle.isValid() || simChannelHandle->empty() ||
533  !simEnergyHandle.isValid() || simEnergyHandle->empty() ||
534  !mcParticleHandle.isValid() ) return;
535 
536  /* std::cout << " Event: " << fEvent
537  << " Run: " << fRun
538  << " SubRun: "<< fSubRun<< std::endl;
539  */
540 
541  // There are several things going on here... for each channel we have particles (track id's) depositing energy in a range to ticks
542  // So... for each channel we want to build a structure that relates particles to tdc ranges and deposited energy (or electrons)
543  // Here is a complicated structure:
544  using TDCToIDEMap = std::map<unsigned short, sim::IDE>; // We need this one in order
545  using ChanToTDCToIDEMap = std::map<raw::ChannelID_t, TDCToIDEMap>;
546  using PartToChanToTDCToIDEMap = std::unordered_map<int, ChanToTDCToIDEMap>;
547 
548  PartToChanToTDCToIDEMap partToChanToTDCToIDEMap;
549 
550  // Build out the above data structure
551  for(const auto& simChannel : *simChannelHandle)
552  {
553  for(const auto& tdcide : simChannel.TDCIDEMap())
554  {
555  for(const auto& ide : tdcide.second) partToChanToTDCToIDEMap[ide.trackID][simChannel.Channel()][tdcide.first] = ide;
556  }
557  }
558 
559  // Here we make a map between track ID and associatied SimEnergyDeposit objects
560  // We'll need this for sorting out the track direction at each hit
561  using SimEnergyDepositVec = std::vector<const sim::SimEnergyDeposit*>;
562  using PartToSimEnergyMap = std::unordered_map<int, SimEnergyDepositVec>;
563 
564  PartToSimEnergyMap partToSimEnergyMap;
565 
566  for(const auto& simEnergy : *simEnergyHandle)
567  {
568  partToSimEnergyMap[simEnergy.TrackID()].push_back(&simEnergy);
569  }
570 
571  // Now we create a data structure to relate hits to their channel ID
572  using ChanToHitVecMap = std::unordered_map<raw::ChannelID_t,std::vector<const recob::Hit*>>;
573 
574  ChanToHitVecMap channelToHitVec;
575 
576  // And now fill it
577  for(const auto& hitLabel : fHitProducerLabelVec)
578  {
579  art::Handle<std::vector<recob::Hit>> hitHandle;
580  evt.getByLabel(hitLabel, hitHandle);
581 
582  for(const auto& hit : *hitHandle)
583  channelToHitVec[hit.Channel()].push_back(&hit);
584  }
585 
586  // We need a structure to relate hits to MCParticles and their trajectory ID
587  using MCParticleTrajIdxPair = std::pair<const simb::MCParticle*,size_t>;
588  using MCPartPointPairIdxVec = std::vector<MCParticleTrajIdxPair>;
589  using HitToMCPartToIdxMap = std::unordered_map<const recob::Hit*,MCPartPointPairIdxVec>;
590 
591  HitToMCPartToIdxMap hitToMcPartToIdxMap;
592 
593  // It is useful to create a mapping between trackID and MCParticle
594  using TrackIDToMCParticleMap = std::unordered_map<int, const simb::MCParticle*>;
595 
596  TrackIDToMCParticleMap trackIDToMCParticleMap;
597 
598  for(const auto& mcParticle : *mcParticleHandle)
599  trackIDToMCParticleMap[mcParticle.TrackId()] = &mcParticle;
600 
601  const lariov::ChannelStatusProvider& chanFilt = art::ServiceHandle<lariov::ChannelStatusService>()->GetProvider();
602 
603  // Look up the list of bad channels
604  art::Handle< std::vector<int>> badChannelHandle;
605  evt.getByLabel(fBadChannelProducerLabel, badChannelHandle);
606 
607  std::vector<int> nSimChannelHitVec = {0,0,0};
608  std::vector<int> nRecobHitVec = {0,0,0};
609  std::vector<int> nFakeHitVec = {0,0,0};
610  std::vector<int> nSimulatedWiresVec = {0,0,0};
611 
612  std::cout << "***************** EVENT " << fEvent << " ******************" << std::endl;
613  std::cout << "-- Looping over channels for hit efficiency, # MC Track IDs: " << partToChanToTDCToIDEMap.size() << std::endl;
614 
615  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
616 
617  // Initiate loop over MC Track IDs <--> SimChannel information by wire and TDC
618  for(const auto& partToChanInfo : partToChanToTDCToIDEMap)
619  {
620  TrackIDToMCParticleMap::const_iterator trackIDToMCPartItr = trackIDToMCParticleMap.find(partToChanInfo.first);
621 
622  if (trackIDToMCPartItr == trackIDToMCParticleMap.end()) continue;
623 
624  const simb::MCParticle* mcParticle = trackIDToMCPartItr->second;
625 
626  int trackPDGCode = mcParticle->PdgCode();
627  std::string processName = mcParticle->Process();
628 
629  // Looking for primary muons (e.g. CR Tracks)
630  // fTotalElectronsVec.push_back(totalElectrons);
631  if (fabs(trackPDGCode) != 13 || processName != "primary") continue;
632  // if (fabs(trackPDGCode) != 13) continue;
633  //if (fEvent != 44) continue;
634  std::cout << ">>>> Loop on MCParticle: " << mcParticle << ", pdg: " << trackPDGCode << ", process: " << processName << std::endl;
635 
636  // Let's recover the SimEnergyDeposit vector for this track
637 // PartToSimEnergyMap::iterator simEneDepItr = partToSimEnergyMap.find(partToChanInfo.first);
638 
639 // if (simEneDepItr == partToSimEnergyMap.end())
640 // {
641 // std::cout << "No match for SimEnergyDeposit" << std::endl;
642 // continue;
643 // }
644 
645  // We should sort this vector by time, this will facilitate searching for matches to SimChannels later
646 // SimEnergyDepositVec& simEnergyDepositVec = simEneDepItr->second;
647 
648 // std::sort(simEnergyDepositVec.begin(),simEnergyDepositVec.end(),[](const auto& left,const auto& right){return left->T() < right->T();});
649 
650 // std::cout << " -- Processing track id: " << partToChanInfo.first << ", with " << simEnergyDepositVec.size() << " SimEnergyDeposit objects" << ", # channels: " << partToChanInfo.second.size() << std::endl;
651 
652  fSimPDG = trackPDGCode;
653  fSimTrackID = mcParticle->TrackId();
654 
655  fmcpartvx.push_back(mcParticle->Vx());
656  fmcpartvy.push_back(mcParticle->Vy());
657  fmcpartvz.push_back(mcParticle->Vz());
658 
659  fmcpartpx.push_back(mcParticle->Px());
660  fmcpartpy.push_back(mcParticle->Py());
661  fmcpartpz.push_back(mcParticle->Pz());
662  fmcparte.push_back(mcParticle->E());
663  //std::cout << "Energy: " << mcParticle->E() << std::endl;
664  // Recover particle position and angle information
665  Eigen::Vector3f partStartPos(mcParticle->Vx(),mcParticle->Vy(),mcParticle->Vz());
666  Eigen::Vector3f partStartDir(mcParticle->Px(),mcParticle->Py(),mcParticle->Pz());
667 
668  partStartDir.normalize();
669 
670  Eigen::Vector2f partStartDirVecXZ(partStartDir[0],partStartDir[2]);
671 
672  partStartDirVecXZ.normalize();
673 
674  // Assuming the SimChannels contain position information (currently not true for WC produced SimChannels)
675  // then we want to keep a running position
676  std::vector<Eigen::Vector3f> lastPositionVec = {partStartPos,partStartPos,partStartPos};
677 
678  std::cout << " *** Looping over channels, have " << partToChanInfo.second.size() << " channels" << std::endl;
679 
680  for(const auto& chanToTDCToIDEMap : partToChanInfo.second)
681  {
682  // skip bad channels
683  if (fUseBadChannelDB)
684  {
685  // This is the "correct" way to check and remove bad channels...
686  if( chanFilt.Status(chanToTDCToIDEMap.first) < fMinAllowedChanStatus)
687  {
688  std::vector<geo::WireID> wids = fGeometry->ChannelToWire(chanToTDCToIDEMap.first);
689  std::cout << "*** skipping bad channel with status: " << chanFilt.Status(chanToTDCToIDEMap.first)
690  << " for channel: " << chanToTDCToIDEMap.first
691  << ", plane: " << wids[0].Plane
692  << ", wire: " << wids[0].Wire << std::endl;
693  continue;
694  }
695  }
696 
697  // Was a list made available?
698  // If so then we try that
699  if (badChannelHandle.isValid())
700  {
701  std::vector<int>::const_iterator badItr = std::find(badChannelHandle->begin(),badChannelHandle->end(),chanToTDCToIDEMap.first);
702 
703  if (badItr != badChannelHandle->end()) continue;
704  }
705 
706  TDCToIDEMap tdcToIDEMap = chanToTDCToIDEMap.second;
707  float totalElectrons(0.);
708  float totalEnergy(0.);
709  float maxElectrons(0.);
710  unsigned short maxElectronsTDC(0);
711  int nMatchedWires(0);
712  int nMatchedHits(0);
713 
714  // The below try-catch block may no longer be necessary
715  // Decode the channel and make sure we have a valid one
716  std::vector<geo::WireID> wids = fGeometry->ChannelToWire(chanToTDCToIDEMap.first);
717 
718  // Recover plane and wire in the plane
719  unsigned int plane = wids[0].Plane;
720  //unsigned int wire = wids[0].Wire;
721 
722  Eigen::Vector3f avePosition(0.,0.,0.);
723 
724  nSimulatedWiresVec[plane]++; // Loop insures channels are unique here
725 
726  for(const auto& ideVal : tdcToIDEMap)
727  {
728  totalElectrons += ideVal.second.numElectrons;
729  totalEnergy += ideVal.second.energy;
730 
731  if (maxElectrons < ideVal.second.numElectrons)
732  {
733  maxElectrons = ideVal.second.numElectrons;
734  maxElectronsTDC = ideVal.first;
735  }
736 
737  avePosition += Eigen::Vector3f(ideVal.second.x,ideVal.second.y,ideVal.second.z);
738  }
739 
740  // Get local track direction by using the average position of deposited charge as the current position
741  // and then subtracting the last position
742  avePosition /= float(tdcToIDEMap.size());
743 
744  Eigen::Vector3f partDirVec = avePosition - lastPositionVec[plane];
745 
746  partDirVec.normalize();
747 
748  lastPositionVec[plane] = avePosition;
749 
750  // std::cout << "sim wire: "<< nSimulatedWiresVec[plane]++ << ", sim ch hit : " <<nSimChannelHitVec[plane]++ << std::endl;
751  // Threshold for electrons
752  if (totalElectrons > ftotalelctroncut) nSimChannelHitVec[plane]++;
753 
754  // std::cout <<" total electron cut : \t" << ftotalelctroncut << std::endl;
755  // totalElectrons = std::min(totalElectrons, float(99900.));
756 
757  //nSimChannelHitVec[plane]++;
758  // std::cout << "after the check --- line 469 "<< nSimChannelHitVec[plane] << std::endl;
759  unsigned short startTDC = tdcToIDEMap.begin()->first;
760  unsigned short stopTDC = tdcToIDEMap.rbegin()->first;
761 
762  // Convert to ticks to get in same units as hits
763  unsigned short startTick = clockData.TPCTDC2Tick(startTDC) + fOffsetVec[plane];
764  unsigned short stopTick = clockData.TPCTDC2Tick(stopTDC) + fOffsetVec[plane];
765  unsigned short maxETick = clockData.TPCTDC2Tick(maxElectronsTDC) + fOffsetVec[plane];
766 
767  //fSimNumTDCVec[plane]->Fill(stopTick - startTick, 1.);
768  //fSimNumTDCVec[plane]->Fill(stopTick - startTick, 1.);
769 
770  // Let's do the match to the SimChannel now just to see that it is working
771 // Eigen::Vector3f simChanPos(ideVal.second.x,ideVal.second.y,ideVal.second.z);
772 //
773 // float minDist = std::numeric_limits<float>::max();
774 // const sim::SimEnergyDeposit* simEDepMatch(0);
775 //
776 // for(const auto& simEDep : simEnergyDepositVec)
777 // {
778 // Eigen::Vector3f simEnePos(simEDep->MidPointX(),simEDep->MidPointY(),simEDep->MidPointZ());
779 //
780 // float sepDist = (simEnePos - simChanPos).norm();
781 //
782 // if (sepDist < minDist)
783 // {
784 // minDist = sepDist;
785 // simEDepMatch = simEDep;
786 // }
787 // else if (sepDist > 1.5 * minDist) break;
788 // }
789 //
790 // if (!simEDepMatch || minDist > 0.15)
791 // {
792 // std::cout << "Bad match, simEDepMatch: " << simEDepMatch << ", minDist: " << minDist << std::endl;
793 // continue;
794 // }
795 //
796 // // Recalculate the angles of the track at this point
797 // Eigen::Vector3f newPartDirVec(simEDepMatch->EndX()-simEDepMatch->StartX(),
798 // simEDepMatch->EndY()-simEDepMatch->StartY(),
799 // simEDepMatch->EndZ()-simEDepMatch->StartZ());
800 //
801 // newPartDirVec.normalize();
802 //
803 // // Overwrite the angles
804 // partDirVec = newPartDirVec;
805 // projPairDirVec = Eigen::Vector2f(partDirVec[0],partDirVec[2]);
806 //
807 // projPairDirVec.normalize();
808 
809  // Set up to extract the "best" parameters in the event of more than one hit for this pulse train
810  float nElectronsTotalBest(0.);
811  float hitSummedADCBest(0.);
812  float hitIntegralBest(0.);
813  float hitPeakTimeBest(0.);
814  float hitPeakAmpBest(-100.);
815  float hitRMSBest(0.);
816  int hitMultiplicityBest(0);
817  int hitLocalIndexBest(0);
818  float hitGoodnessBest(0.);
819  int hitNumDegreesBest(0);
820  float hitBaselineBest(0.);
821  //float hitSnippetLenBest(0.);
822  unsigned short hitStopTickBest(0);
823  unsigned short hitStartTickBest(0);
824 
825  const recob::Hit* rejectedHit = 0;
826  const recob::Hit* bestHit = 0;
827  // The next mission is to recover the hits associated to this Wire
828  // The easiest way to do this is to simply look up all the hits on this channel and then match
829  ChanToHitVecMap::iterator hitIter = channelToHitVec.find(chanToTDCToIDEMap.first);
830 
831  if (hitIter != channelToHitVec.end())
832  {
833  // Loop through the hits for this channel and look for matches
834  // In the event of more than one hit associated to the sim channel range, keep only
835  // the best match (assuming the nearby hits are "extra")
836  // Note that assumption breaks down for long pulse trains but worry about that later
837  for(const auto& hit : hitIter->second)
838  {
839  //std::cout <<" ........... sigma vector: \t " << fSigmaVec[plane] << std::endl;
840  unsigned short hitStartTick = hit->PeakTime() - fSigmaVec[plane] * hit->RMS();
841  unsigned short hitStopTick = hit->PeakTime() + fSigmaVec[plane] * hit->RMS();
842  // If hit is out of range then skip, it is not related to this particle
843  if (hitStartTick > stopTick || hitStopTick < startTick)
844  {
845  nFakeHitVec[plane]++;
846  rejectedHit = hit;
847  continue;
848  }
849 
850  float hitHeight = hit->PeakAmplitude();
851 
852  // Use the hit with the largest pulse height as the "best"
853  if (hitHeight < hitPeakAmpBest) continue;
854 
855  hitPeakAmpBest = hitHeight;
856  bestHit = hit;
857  hitStartTickBest = hitStartTick;
858  hitStopTickBest = hitStopTick;
859  }
860  // Find a match?
861  if (bestHit)
862  {
863  nElectronsTotalBest = 0.;
864  hitPeakTimeBest = bestHit->PeakTime();
865  hitIntegralBest = bestHit->Integral();
866  hitSummedADCBest = bestHit->SummedADC();
867  hitRMSBest = bestHit->RMS();
868  hitMultiplicityBest = bestHit->Multiplicity();
869  hitLocalIndexBest = bestHit->LocalIndex();
870  hitGoodnessBest = bestHit->GoodnessOfFit();
871  hitNumDegreesBest = bestHit->DegreesOfFreedom();
872  //hitSnippetLenBest = bestHit->EndTick() - bestHit->StartTick();
873  hitBaselineBest = 0.; // To do...
874 
875  // Get the number of electrons
876  for(unsigned short tick = hitStartTickBest; tick <= hitStopTickBest; tick++)
877  {
878  unsigned short hitTDC = clockData.TPCTick2TDC(tick - fOffsetVec[plane]);
879 
880  TDCToIDEMap::iterator ideIterator = tdcToIDEMap.find(hitTDC);
881 
882  if (ideIterator != tdcToIDEMap.end()) nElectronsTotalBest += ideIterator->second.numElectrons;
883  }
884  // Ok, now we need to figure out which trajectory point this hit is associated to
885  // Use the associated IDE to get the x,y,z position for this hit
886  Eigen::Vector3f hitIDEPos(avePosition[0],avePosition[1],avePosition[2]);
887  size_t matchIdx(0);
888  float bestMatch(std::numeric_limits<float>::max());
889 
890  for (size_t idx = 0; idx < mcParticle->NumberTrajectoryPoints(); idx++)
891  {
892  const TLorentzVector& pos = mcParticle->Position(idx); // 4-position in World coordinates
893  Eigen::Vector3f trackPos(pos.X(),pos.Y(),pos.Z());
894  float missDist = (trackPos - hitIDEPos).norm();
895  if (missDist < bestMatch)
896  {
897  bestMatch = missDist;
898  matchIdx = idx;
899  }
900  }
901 
902 // std::cout << " --> Hit: " << bestHit << ", traj idx: " << matchIdx << " with dist: " << bestMatch << std::endl;
903 // if (matchIdx > 0)
904 // {
905 // const TLorentzVector& pos = mcParticle->Position(matchIdx); // 4-position in World coordinates
906 // const TLorentzVector& last = mcParticle->Position(matchIdx-1); // 4-position in World coordinates
907 // float dLastPos = (Eigen::Vector3f(pos.X(),pos.Y(),pos.Z()) - Eigen::Vector3f(last.X(),last.Y(),last.Z())).norm();
908 // std::cout << " --> distance to last trajectory point: " << dLastPos << std::endl;
909 // }
910  if (bestMatch < std::numeric_limits<float>::max())
911  {
912  hitToMcPartToIdxMap[bestHit].push_back(MCParticleTrajIdxPair(mcParticle,matchIdx));
913  //std::cout << "-hit " << bestHit << " associated to MC " << mcParticle << " at point " << matchIdx << std::endl;
914 
915  partDirVec = Eigen::Vector3f(mcParticle->Px(matchIdx),mcParticle->Py(matchIdx),mcParticle->Pz(matchIdx));
916 
917  partDirVec.normalize();
918  }
919  } // end of besthit
920 
921  if (nMatchedHits > 0)
922  {
923  //float chgRatioADC = hitSummedADCBest > 0. ? totalElectrons / hitSummedADCBest : 0.;
924  //float chgRatioInt = hitIntegralBest > 0. ? totalElectrons / hitIntegralBest : 0.;
925  if (totalElectrons > ftotalelctroncut) nRecobHitVec[plane]++;
926  // nRecobHitVec[plane]++;
927  }
928  else if (rejectedHit)
929  {
930  //unsigned short hitStartTick = rejectedHit->PeakTime() - fSigmaVec[plane] * rejectedHit->RMS();
931  //unsigned short hitStopTick = rejectedHit->PeakTime() + fSigmaVec[plane] * rejectedHit->RMS();
932  std::cout << " -Rejected: " << rejectedHit->WireID().Plane << ", ph: " << rejectedHit->PeakAmplitude() << std::endl;
933  }
934  else
935  {
936  mf::LogDebug("ThroughgoingmuonAnalyzer") << "==> No match, TPC/Plane/Wire: "
937  << "/" << wids[0].TPC
938  << "/" << wids[0].Plane
939  << "/" << wids[0].Wire
940  << ", # electrons: " << totalElectrons
941  << ", startTick: " << startTick
942  << ", stopTick: " << stopTick << std::endl;
943  } // end of matched hit
944  } // end of channel to hit map iterator
945 
946  float matchHit = std::min(nMatchedHits,1);
947  fHitEffvselecVec[plane]->Fill(totalElectrons, matchHit, 1.);
948 
949  fTPCVec.push_back(wids[0].TPC);
950  fCryoVec.push_back(wids[0].Cryostat);
951  fPlaneVec.push_back(wids[0].Plane);
952  fWireVec.push_back(wids[0].Wire);
953 
954  fTotalElectronsVec.push_back(totalElectrons);
955  fTotalElecEnergyVec.push_back(totalEnergy);
956  fMaxElectronsVec.push_back(maxElectrons);
957  fStartTickVec.push_back(startTick);
958  fStopTickVec.push_back(stopTick);
959  fMaxETickVec.push_back(maxETick);
960  fPartDirX.push_back(partDirVec[0]);
961  fPartDirY.push_back(partDirVec[1]);
962  fPartDirZ.push_back(partDirVec[2]);
963 
964  fNMatchedWires.push_back(nMatchedWires);
965  fNMatchedHits.push_back(nMatchedHits);
966 
967  fHitPeakTimeVec.push_back(hitPeakTimeBest);
968  fHitPeakAmpVec.push_back(hitPeakAmpBest);
969  fHitPeakRMSVec.push_back(hitRMSBest);
970  fHitBaselinevec.push_back(hitBaselineBest);
971  fHitSummedADCVec.push_back(hitSummedADCBest);
972  fHitIntegralVec.push_back(hitIntegralBest);
973  fHitStartTickVec.push_back(hitStartTickBest);
974  fHitStopTickVec.push_back(hitStopTickBest);
975  fHitMultiplicityVec.push_back(hitMultiplicityBest);
976  fHitLocalIndexVec.push_back(hitLocalIndexBest);
977  fHitGoodnessVec.push_back(hitGoodnessBest);
978  fNumDegreesVec.push_back(hitNumDegreesBest);
979  } //end of partcle to id map info
980  } // Done looping over track <--> MCParticle associations
981 
982  std::cout << "-- Now starting loop over tracks" << std::endl;
983  std::cout << "-- Hit/MCParticle map size: " << hitToMcPartToIdxMap.size() << std::endl;
984 
985  // Let's keep track of track to MCParticle and hits.
986  using MCPartToHitVecMap = std::unordered_map<const simb::MCParticle*,std::vector<const recob::Hit*>>;
987  using TrackToMCMap = std::unordered_map<const recob::Track*, MCPartToHitVecMap>;
988 
989  TrackToMCMap trackToMCMap;
990 
991  std::cout << "** Starting outer loop over all sets of tracks" << std::endl;
992 
993  // The game plan for this module is to look at recob::Tracks and objects associated to tracks
994  // The first thing we need to do is match tracks to MCParticles which we can do using structures
995  // created in the above loops over hits.
996  for(const auto& trackLabel : fTrackProducerLabelVec)
997  {
998  art::Handle<std::vector<recob::Track> > trackHandle;
999  //std::vector<art::Ptr<recob::Track> > tracklist;
1000  evt.getByLabel(trackLabel, trackHandle);
1001  if (!trackHandle.isValid()) return;
1002 
1003  //track information
1004  ntrk_reco += trackHandle->size();
1005  fNtracks_reco = ntrk_reco;
1006 
1007  // Recover the collection of associations between tracks and hits
1008  art::FindMany<recob::Hit> trackHitAssns(trackHandle, evt, trackLabel);
1009 
1010  // First mission will be to find the tracks associated to our primary MCParticle...
1011  for(size_t trackIdx = 0; trackIdx < trackHandle->size(); trackIdx++)
1012  {
1013  art::Ptr<recob::Track> ptrack(trackHandle,trackIdx);
1014 
1015  // Associated hits
1016  const std::vector<const recob::Hit*>& trackHitVec = trackHitAssns.at(ptrack.key());
1017 
1018  // Keep track of the number of hits per MCParticle for this track
1019  MCPartToHitVecMap mcPartCountMap;
1020 
1021  std::cout << "- Track idx " << trackIdx << ", ptr: " << ptrack.get() << " has " << trackHitVec.size() << " hits" << std::endl;
1022 
1023  // Loop through hits and use associations tables to relate to MCParticle
1024  for(const auto& hit : trackHitVec)
1025  {
1026  // Get list of MCParticles/traj point ids
1027  HitToMCPartToIdxMap::iterator hitToMcPartIdxItr = hitToMcPartToIdxMap.find(hit);
1028 
1029  if (hitToMcPartIdxItr != hitToMcPartToIdxMap.end())
1030  {
1031  for(const auto& mcPartIdxPair : hitToMcPartIdxItr->second)
1032  {
1033  mcPartCountMap[mcPartIdxPair.first].push_back(hit);
1034  }
1035  }
1036  else
1037  mcPartCountMap[nullptr].push_back(hit);
1038  }
1039 
1040  std::cout << " - associated to " << mcPartCountMap.size() << " MCParticles" << std::endl;
1041 
1042  for(const auto& mcPair : mcPartCountMap)
1043  std::cout << " - MCParticle: " << mcPair.first << ", # hits: " << mcPair.second.size() << std::endl;
1044 
1045  // Find the best match (simple majority logic)
1046  const simb::MCParticle* bestMCMatch(nullptr);
1047  size_t bestCount(0);
1048 
1049  for(const auto& pair : mcPartCountMap)
1050  {
1051  //std::cout << " - MCParticle: " << pair.first <<", pdg: " << pair.first->PdgCode() << ", # hits: " << pair.second.size() << std::endl;
1052  if (pair.first && pair.second.size() > bestCount)
1053  {
1054  bestMCMatch = pair.first;
1055  bestCount = pair.second.size();
1056  }
1057  }
1058 
1059  // std::cout << " ==> bestMCMatch: " << bestMCMatch << ", pdg:" << bestMCMatch->PdgCode()<< ", bestCount: " << bestCount << std::endl;
1060  //std::cout << " ==> bestMCMatch: " << bestMCMatch << ", bestCount: " << bestCount << std::endl;
1061 
1062  if (bestCount > 0)
1063  {
1064  MCPartToHitVecMap::iterator bestItr = mcPartCountMap.find(bestMCMatch);
1065 
1066  trackToMCMap[ptrack.get()][bestItr->first] = bestItr->second;
1067 
1068  }
1069  }
1070  }
1071 
1072  std::cout << "====> trackToMCMap size: " << trackToMCMap.size() << std::endl;
1073  // std::cout << "====> mcPartCountMap size: " << MCPartToHitVecMap.size() << std::endl;
1074 
1075  fNtracks_primary = trackToMCMap.size();
1076  std::cout << "==================================================================================> " << std::endl;
1077  if (fNtracks_primary != 1)
1078  std::cout << "**> tracks associated to primary = " << fNtracks_primary << ", event: " << evt.id().event() << std::endl;
1079 
1080 
1081  // Go through the tracks matched to the primary
1082  for(const auto& trackMCPair : trackToMCMap)
1083  {
1084  const recob::Track& track = *trackMCPair.first;
1085 
1086  // std::cout << track.ParticleId()<< std::endl;
1087  unsigned int ntraj = track.NumberTrajectoryPoints();
1088  //std::cout << "ntraj: \t"<< ntraj << std::endl;
1089 
1090  if(ntraj <= 0) continue;
1091 
1092 
1093  // const geo::TPCGeo* tpcGeo = fGeometry->PositionToTPCptr(track.Start());
1094 
1095  //if (!tpcGeo) continue;
1096 
1097  //std::cout << tpcGeo->ActiveBoundingBox().ContainsPosition(track.Start()) << std::endl;
1098  //std::cout << "track length : "<< track.Length() << std::endl;
1099 
1100 
1101  // const recob::tracking::Point_t& pos
1102  // TVector3 pos = track.Vertex<TVector3>();
1103  //TVector3 pos = track.Start<TVector3>();
1104  //TVector3 dir = track.StartDirection<TVector3>();
1105  //TVector3 end = track.End<TVector3>();
1106  // TVector3 enddir = track.EndDirection<TVector3>();
1107  // double pstart = track.VertexMomentum();
1108  //double pend = track.EndMomentum();
1109 
1110  // double length = track.Length();
1111  // double theta_xz = std::atan2(dir.X(), dir.Z());
1112  //double theta_yz = std::atan2(dir.Y(), dir.Z());
1113  // std::cout << end.X() << "\t" << end.Y() << "\t" << end.Z()<<std::endl;
1114  // TVector3 pos(0,2.43493,851.593);
1115  // if(!cryo0.ContainsPosition(pos) && !cryo1.ContainsPosition(pos)) continue;
1116  // std::cout << trackLabel << " no.of track:\t" <<trackHandle->size() << "\t"<< ntracks_reco<< std::endl;
1117  // This seems a redundant test for fit tracks?
1118  // geo::Point_t trackPoint(pos.X(),pos.Y(),pos.Z());
1119  // std::cout << track.X() << "\t" << track.Y() << "\t" << track.Z()<<std::endl;
1120 
1121  //std::cout << "Reco Track--- Before: \t" << pos.X() <<"\t"<< pos.Y() <<"\t"<< pos.Z() << std::endl;
1122  // geo::Point_t const trackPoint(track.Start());
1123  //geo::Point_t const trackEndPoint(track.End());
1124 
1125  //const geo::TPCGeo* tpcGeo = fGeometry->PositionToTPCptr(trackPoint);
1126 
1127  //if (!tpcGeo) continue;
1128  //if (!tpcGeo->ActiveBoundingBox().ContainsPosition(trackPoint)) continue; // out of active volume
1129  //if (!tpcGeo->ActiveBoundingBox().ContainsPosition(trackEndPoint)) continue; // out of active volume
1130 
1131  // std::cout << "Reco Track start point : \t" << trackPoint.X() <<"\t"<< trackPoint.Y() <<"\t"<< trackPoint.Z() << std::endl;
1132 
1133 
1134 
1135  // if (tpcGeo->ActiveBoundingBox().ContainsPosition(trackPoint)) trackEndPoint = trackPoint;
1136  //else break; // this means that once you are out of the active volume the first time, you are done; think if that's what you want,
1137  //since tracks scatter a lot in LAr
1138  // std::cout << "Reco Track start point : \t" << trackPoint.X() <<"\t"<< trackPoint.Y() <<"\t"<< trackPoint.Z() << std::endl;
1139  //std::cout << "Reco Track end point : \t" << trackEndPoint.X() <<"\t"<< trackEndPoint.Y() <<"\t"<< trackEndPoint.Z() << std::endl;
1140  // std::cout <<: \t" << end.X() <<"\t"<< end.Y() <<"\t"<< end.Z() << std::endl;
1141 
1142  //geo::Point_t trackstartPoint ;
1143  // geo::Point_t trackendPoint ;
1144  std::vector<double> posx, posy, posz;
1145  std::vector<double> dirx, diry, dirz;
1146  // const geo::CryostatID* cryo;
1147  // std::cout << geo::CryostatID(0)<< std::endl;
1148 
1149  // std::cout << "cryostat: " << cryostat << std::endl;
1150  const geo::CryostatID C0id { 0 };
1151  // ...
1152  //std::cout << "cryostat: " << C0id << std::endl;
1153  for (unsigned int i= 0; i < ntraj; i++)
1154  {
1155  //geo::Point_t trackstartPoint = track.LocationAtPoint(i);
1156  //geo::Point_t trackendPoint = track.LocationAtPoint(i);
1157  TVector3 mom = track.MomentumVectorAtPoint<TVector3>(i);
1158  geo::Point_t trackstartPoint = track.LocationAtPoint(i);
1159  //trackendPoint = track.LocationAtPoint(i);
1160  const geo::TPCGeo* tpcGeom = fGeometry->PositionToTPCptr(trackstartPoint);
1161 
1162  if (!tpcGeom) continue;
1163  if (tpcGeom->ID() != C0id) continue; // point not in cryostat 0
1164  if (!tpcGeom->ActiveBoundingBox().ContainsPosition(trackstartPoint)) continue;
1165  //if (tpcGeom->ActiveBoundingBox().ContainsPosition(trackstartPoint)) trackendPoint = trackstartPoint;
1166  //else break;
1167  posx.push_back( trackstartPoint.X() );
1168  posy.push_back( trackstartPoint.Y() );
1169  posz.push_back( trackstartPoint.Z() );
1170 
1171  dirx.push_back( mom.X() );
1172  diry.push_back( mom.Y() );
1173  dirz.push_back( mom.Z() );
1174 
1175  // if(i == 10) std::cout<< "track length inside loop: "<< track.Length() << std::endl;
1176  }
1177 
1178  //std::cout << "Reco Track start point : \t" << posx.front() <<"\t"<< posy.front() <<"\t"<< posz.front() << std::endl;
1179  //std::cout << "Reco Track end point : \t" << posx.back() <<"\t"<< posy.back() <<"\t"<< posz.back() << std::endl;
1180 
1181 
1182  // frecostartx.push_back( pos.X() );
1183  // frecostarty.push_back( pos.Y() );
1184  // frecostartz.push_back( pos.Z() );
1185 
1186  // frecoendx.push_back( end.X() );
1187  // frecoendy.push_back( end.Y() );
1188  // frecoendz.push_back( end.Z() );
1189  if (!posx.empty())
1190  {
1191  frecostartx.push_back( posx.front() );
1192  frecostarty.push_back( posy.front() );
1193  frecostartz.push_back( posz.front() );
1194 
1195  frecoendx.push_back( posx.back() );
1196  frecoendy.push_back( posy.back() );
1197  frecoendz.push_back( posz.back() );
1198 
1199  frecostartdirx.push_back( dirx.front() );
1200  frecostartdiry.push_back( diry.front() );
1201  frecostartdirz.push_back( dirz.front() );
1202 
1203  frecoenddirx.push_back( dirx.back() );
1204  frecoenddiry.push_back( diry.back() );
1205  frecoenddirz.push_back( dirz.back() );
1206 
1207  double theta_xz = std::atan2(dirx.front(), dirz.front());
1208  double theta_yz = std::atan2(diry.front(), dirz.front());
1209 
1210  fThetaXZ.push_back( theta_xz );
1211  fThetaYZ.push_back( theta_yz );
1212  }
1213  // frecostartdirx.push_back( dir.X() );
1214  // frecostartdiry.push_back( dir.Y() );
1215  // frecostartdirz.push_back( dir.Z() );
1216 
1217  // frecoenddirx.push_back( enddir.X() );
1218  //frecoenddiry.push_back( enddir.Y() );
1219  //frecoenddirz.push_back( enddir.Z() );
1220 
1221  // fLength. push_back( track.Length() );
1222  // fThetaXZ.push_back( theta_xz );
1223  //fThetaYZ.push_back( theta_yz );
1224 
1225  fLength. push_back( track.Length() );
1226  fTheta. push_back( track.Theta() );
1227  fPhi. push_back( track.Phi() );
1228 
1229  posx.clear(); posy.clear(); posz.clear();
1230  dirx.clear(); diry.clear(); dirz.clear();
1231  // }
1232  //std::cout << pos.X() << "\t" << pos.Y() << "\t" << pos.Z()<<std::endl;
1233  // std::cout << end.X() << "\t" << end.Y() << "\t" << end.Z()<<std::endl;
1234 
1235  //std::cout << "-- Now starting loop over trajectory points for track " << track.ID() << std::endl;
1236 
1237  // Now start going through the associated MCParticles, focus is on the primary for now
1238  for(const auto& mcPartHitVecPair : trackMCPair.second)
1239  {
1240  const simb::MCParticle* mcParticle = mcPartHitVecPair.first;
1241 
1242  // Skip the unassociated hits for now
1243  if (!mcParticle) continue;
1244 
1245 
1246  geo::Point_t vtxPoint(mcParticle->Vx(), mcParticle->Vy(), mcParticle->Vz());
1247 
1248  const geo::TPCGeo* tpcGeo = fGeometry->PositionToTPCptr(vtxPoint);
1249 
1250  if (!tpcGeo) continue;
1251  if (tpcGeo->ID() != C0id) continue; // point not in cryostat 0
1252  if (!tpcGeo->ActiveBoundingBox().ContainsPosition(vtxPoint)) continue; // out of active volume
1253 
1254  // Find the extreme trajectory points associated to hits for this track
1255  size_t minTrajIdx = std::numeric_limits<size_t>::max();
1256  size_t maxTrajIdx = 0;
1257 
1258  for(const auto& hit : mcPartHitVecPair.second)
1259  {
1260  // std::cout << " # Hit: "<< mcPartHitVecPair.second.size()<< "; Hit Index: " << &hit - &mcPartHitVecPair.second[0] << std::endl;
1261 
1262  geo::Point_t loc = track.LocationAtPoint(&hit - &mcPartHitVecPair.second[0]);
1263 
1264  // std::cout << "location at x: " << loc.X() << " ; location at y: " << loc.Y()<< " ; location at z: " << loc.Z() << std::endl;
1265 
1266  flocx.push_back(loc.X());
1267  flocy.push_back(loc.Y());
1268  flocz.push_back(loc.Z());
1269 
1270  // **THIS NEEDS TO BE CHECKED **
1271  // There should be only one entry for the trajectory index the way the code is set up to run now
1272  // (Only looking at the primary)
1273  size_t trajIdx = hitToMcPartToIdxMap[hit].front().second;
1274 
1275  minTrajIdx = std::min(minTrajIdx,trajIdx);
1276  maxTrajIdx = std::max(maxTrajIdx,trajIdx);
1277  }
1278 
1279  // std::cout << " # MCParticle trajectory points: " << mcParticle->NumberTrajectoryPoints() << ", max: " << maxTrajIdx << ", min: " << minTrajIdx << std::endl;
1280 
1281  if (mcParticle->Px(maxTrajIdx) == 0 && mcParticle->Py(maxTrajIdx) == 0 && mcParticle->Pz(maxTrajIdx) == 0) maxTrajIdx -= 1;
1282  //if (mcParticle->Px(minTrajIdx) == 0 && mcParticle->Py(minTrajIdx) == 0 && mcParticle->Pz(minTrajIdx) == 0){minTrajIdx = minTrajIdx - 1;}
1283 
1284  // Fill MC info for trajectory point limits
1285 
1286  TVector3 minPointDir(mcParticle->Px(minTrajIdx),mcParticle->Py(minTrajIdx),mcParticle->Pz(minTrajIdx));
1287  TVector3 maxPointDir(mcParticle->Px(maxTrajIdx),mcParticle->Py(maxTrajIdx),mcParticle->Pz(maxTrajIdx));
1288 
1289  // Px,Py,Pz are the momentum of the track, "SetMag(1.)" is meant to normalize so you get a direction vector.
1290  // if one of the two trajectory points somehow has zero momentum, avoid those cases.
1291  if (!(maxPointDir.Mag() > 0.)) {maxTrajIdx -= 1; maxPointDir=TVector3(mcParticle->Px(maxTrajIdx),mcParticle->Py(maxTrajIdx),mcParticle->Pz(maxTrajIdx));}
1292  if (!(minPointDir.Mag() > 0.)) {minTrajIdx -= 1; minPointDir=TVector3(mcParticle->Px(minTrajIdx),mcParticle->Py(minTrajIdx),mcParticle->Pz(minTrajIdx));}
1293  // if (!(maxPointDir.Mag() > 0.)) maxTrajIdx -= 1;
1294 
1295  //std:: cout << " minx: " << minPointDir.X() <<", miny: " << minPointDir.Y() <<", minz: " << minPointDir.Z() <<", minMag: " << minPointDir.Mag()
1296  // << ", maxx: " << maxPointDir.X() <<", maxy: " << maxPointDir.Y() <<", maxz: " << maxPointDir.Z() <<", maxMag: " << maxPointDir.Mag() << std::endl;
1297 
1298 
1299  minPointDir.SetMag(1.);
1300  maxPointDir.SetMag(1.);
1301 
1302  fmctstartx.push_back(mcParticle->Vx(minTrajIdx));
1303  fmctstarty.push_back(mcParticle->Vy(minTrajIdx));
1304  fmctstartz.push_back(mcParticle->Vz(minTrajIdx));
1305 
1306  fmctendx.push_back(mcParticle->Vx(maxTrajIdx));
1307  fmctendy.push_back(mcParticle->Vy(maxTrajIdx));
1308  fmctendz.push_back(mcParticle->Vz(maxTrajIdx));
1309 
1310  fmctstartdirx.push_back(minPointDir[0]);
1311  fmctstartdiry.push_back(minPointDir[1]);
1312  fmctstartdirz.push_back(minPointDir[2]);
1313 
1314  fmctenddirx.push_back(maxPointDir[0]);
1315  fmctenddiry.push_back(maxPointDir[1]);
1316  fmctenddirz.push_back(maxPointDir[2]);
1317 
1318  fmctThetaXZ.push_back( std::atan2(minPointDir.X(), minPointDir.Z()) );
1319  fmctThetaYZ.push_back( std::atan2(minPointDir.Y(), minPointDir.Z()) );
1320 
1321  fmctTheta. push_back( minPointDir.Theta() );
1322  fmctPhi. push_back( minPointDir.Phi() );
1323 
1324  //fmcstartx.push_back(mcParticle->Vx());
1325  //double earliest_time = 1e15;
1326  // TVector point;
1327  //std::cout << "Before: \t" <<partStartPos[0] <<"\t"<< partStartPos[1] <<"\t"<< partStartPos[2] << std::endl;
1328  std::vector<double> x, y, z, t;
1329  std::vector<double> px, py, pz, e;
1330  TLorentzVector dir_start;
1331  // std::vector<float> endx, endy, endz, endt;
1332  //int last = (int)trackIDToMCPartItr->second->NumberTrajectoryPoints() -1;
1333  //std::vector<double> vtime;
1334  for (unsigned int i=0; i<mcParticle->NumberTrajectoryPoints(); i++)
1335  {
1336  const TLorentzVector& pos = mcParticle->Position(i); // 4-position in World coordinates
1337  const TLorentzVector& mom = mcParticle->Momentum(i); // 4-position in World coordinates
1338  //point(pos.X(),pos.Y(),pos.Z()
1339  //const double point[3] = {pos.X(),pos.Y(),pos.Z()};
1340 
1341  //std::cout << "Mc Part--- Before: \t" << pos.X() <<"\t"<< pos.Y() <<"\t"<< pos.Z() << std::endl;
1342  // if (mcParticle->Vx() < -894.951 && mcParticle->Vx() > 894.951) std::cout <<"---------------->>>>>>>>>>>>>> " << mcParticle->Vx() << std::endl;
1343 
1344  geo::Point_t trackPoint(pos.X(),pos.Y(),pos.Z());
1345 
1346  const geo::TPCGeo* tpcGeo = fGeometry->PositionToTPCptr(trackPoint);
1347 
1348  if (!tpcGeo) continue;
1349  if (tpcGeo->ID() != C0id) continue; // point not in cryostat 0
1350  if (!tpcGeo->ActiveBoundingBox().ContainsPosition(trackPoint)) continue; // out of active volume
1351 
1352  // std::cout << "McPart--- After: " << tpcGeo->ID() << " - " << pos.X() <<"\t"<< pos.Y() <<"\t"<< pos.Z() << std::endl;
1353  //std::cout << "Mc Part--- After: \t" << pos.X() <<"\t"<< pos.Y() <<"\t"<< pos.Z() << std::endl;
1354 
1355  // x.push_back(mcParticle->Vx());
1356  x.push_back(pos.X());
1357  y.push_back(pos.Y());
1358  z.push_back(pos.Z());
1359  t.push_back(pos.T());
1360 
1361  px.push_back(mom.X());
1362  py.push_back(mom.Y());
1363  pz.push_back(mom.Z());
1364  e.push_back(mom.T());
1365 
1366  } // no. of trajectory loop
1367  // for (unsigned int i =0; i< x.size(); i++){
1368  //std::cout << x[i] << std::endl;
1369  //}
1370 
1371  // std::cout << " ==> filling overall MC arrays, size of x: " << x.size() << std::endl;
1372  //std::cout << "---- front -------- x: " << x.front() << ", y: " << y.front() << ",z: " << z.front() << std::endl;
1373  //std::cout << "---- back ---------- x: " << x.back() << ", y: " << y.back() << ",z: " << z.back() << std::endl;
1374 
1375  // Make sure something got added
1376  if (!x.empty())
1377  {
1378  fmcstartx.push_back(x.front());
1379  fmcstarty.push_back(y.front());
1380  fmcstartz.push_back(z.front());
1381 
1382  fmcendx.push_back(x.back());
1383  fmcendy.push_back(y.back());
1384  fmcendz.push_back(z.back());
1385 
1386  fmcstartdirx.push_back(px.front());
1387  fmcstartdiry.push_back(py.front());
1388  fmcstartdirz.push_back(pz.front());
1389 
1390  fmcenddirx.push_back(px.back());
1391  fmcenddiry.push_back(py.back());
1392  fmcenddirz.push_back(pz.back());
1393 
1394  dir_start.SetPxPyPzE(px.front(),py.front(),pz.front(),e.front());
1395 
1396  fmcThetaXZ.push_back( std::atan2(px.front(), pz.front()) );
1397  fmcThetaYZ.push_back( std::atan2(py.front(), pz.front()) );
1398 
1399  fmcTheta. push_back( dir_start.Theta() );
1400  fmcPhi. push_back( dir_start.Phi() );
1401  }
1402 
1403  x.clear(); y.clear(); z.clear(); t.clear();
1404  px.clear(); py.clear(); pz.clear(); e.clear();
1405  } // Done with loop over MCParticles associted to this track
1406  }
1407 
1408  // float diffex = std::min(std::abs((*recostartx)[0] - (*mcendx)[0]), std::abs((*recoendx)[0]-(*mcendx)[0]));
1409  //float diffey = std::min(std::abs(frecostarty[0] - fmcendy[0]), std::abs(frecoendy[0]-fmcendy[0]));
1410  // float diffez = std::min(std::abs((*recostartz)[0] - (*mcendz)[0]), std::abs((*recoendz)[0]-(*mcendz)[0]));
1411  //if (diffey < 34 && diffey > 32)
1412  //std::cout << diffey << std::endl;
1413  /*
1414  if (fNtracks_primary == 1){
1415  if (!fmcstartx.empty()){
1416  if (!frecostartx.empty()){
1417  // float diffx = std::min(std::abs(frecostartx[0] - fmcstartx[0]), std::abs(frecoendx[0]-fmcstartx[0]));
1418 
1419 
1420  float diffxer = std::abs(frecostartx[0] - fmcstartx[0]);
1421  float diffyer = std::abs(frecostarty[0] - fmcstarty[0]);
1422  float diffzer = std::abs(frecostartz[0] - fmcstartz[0]);
1423 
1424  float diffxsr = std::abs(frecoendx[0] - fmcendx[0]);
1425  float diffysr = std::abs(frecoendy[0] - fmcendy[0]);
1426  float diffzsr = std::abs(frecoendz[0] - fmcendz[0]);
1427 
1428  float diffdr = std::min( std::sqrt(std::pow(diffxsr,2)+std::pow(diffysr,2)+std::pow(diffzsr,2)) ,
1429  std::sqrt(std::pow(diffxer,2)+std::pow(diffyer,2)+std::pow(diffzer,2)) );
1430  if (diffdr <= 1.04){
1431  if ( frecostartx[0] < -268.49 || frecostartx[0] > -171.94){
1432  std::cout <<"------> outside fiducial volume ------->mcstartx: " << fmcstartx[0] <<" ; recoendx: "
1433  << frecoendx[0] << "; mcendx: "<< fmcendx[0] <<", recostartx: "<< frecostartx[0]
1434  << ", (reco - true) startx: "<< frecostartx[0] - fmcstartx[0] <<" ; mcstarty: " << fmcstarty[0] <<" ; recoendy: "
1435  << frecoendy[0] << "; mcendy: "<< fmcendy[0] <<", recostarty: "<< frecostarty[0]
1436  << ", (reco - true) starty: "<< frecostarty[0] - fmcstarty[0] << " ; mcstartz: " << fmcstartz[0] <<" ; recoendz: "
1437  << frecoendz[0] << "; mcendz: "<< fmcendz[0] <<", recostartz: "<< frecostartz[0]
1438  << ", (reco - true) startz: "<< frecostartz[0] - fmcstartz[0] << " ; checking minm: "<< diffdr << std::endl;
1439  // std::cout <<"-----------------------------> outside fiducial volume ........ mcstartx: " << (*mcstarty)[0] <<" ; recoendx: "
1440  // << (*recoendy)[0] << "; mcendx: "<< (*mcendy)[0] <<", recostartx: "<< (*recostarty)[0]
1441  // << ", (reco - true) startx: "<< (*recostarty)[0] - (*mcstarty)[0] << "; checking minm: "<< diffdr << std::endl;
1442  // std::cout <<"-----------------------------> outside fiducial volume ........ mcstartx: " << (*mcstartx)[0] <<" ; recoendx: "
1443  // << (*recoendx)[0] << "; mcendx: "<< (*mcendx)[0] <<", recostartx: "<< (*recostartx)[0]
1444  // << ", (reco - true) startx: "<< (*recostartx)[0] - (*mcstartx)[0] << "; checking minm: "<< diffdr << std::endl;
1445  }
1446  }
1447 
1448  // std::cout << "***************************************************************> everything work upto here before track loop" << std::endl;
1449  // if (!fmcstartx.empty()) //&& fmcendx[0] > 0)
1450  // std::cout <<"........................ mcstartx: " << fmcstartx[0] <<" , mcgenerationvx: " << fmcendx[0] <<", recostartx: "<< frecoendx[0] << std::endl;
1451 
1452  //std::cout <<"mcstartx: " << fmcstartx[0] <<" ; recoendx: " << frecoendx[0] << "; mcendx: "<< fmcendx[0] <<", recostartx: "<< frecostartx[0]
1453  // << ", (reco - true) startx: "<< frecostartx[0] - fmcstartx[0] << "; checking minm: "<< diffx << std::endl;
1454 
1455  //if ( frecostartx[0] < -268.49 || frecostartx[0] > -171.94){
1456  //std::cout <<"-----------------------------> outside fiducial volume ........ mcstartx: " << fmcstartx[0] <<" ; recoendx: "
1457  // << frecoendx[0] << "; mcendx: "<< fmcendx[0] <<", recostartx: "<< frecostartx[0]
1458  // << ", (reco - true) startx: "<< frecostartx[0] - fmcstartx[0] << "; checking minm: "<< diffx << std::endl;
1459  }
1460 
1461  }
1462  }
1463  }
1464  */
1465  //std::cout << x[0] << "\t" << y[0] << "\t" << z[0] << "\t" <<
1466  //x[last] << "\t" << y[last] << "\t" << z[last] << std::endl;
1467 
1468  std::cout << "-- final quick loop to fill summary hists" << std::endl;
1469 
1470  for(size_t idx = 0; idx < fGeometry->Nplanes();idx++)
1471  {
1472  if (nSimChannelHitVec[idx] > 10)
1473  {
1474  float hitEfficiency = float(nRecobHitVec[idx]) / float(nSimChannelHitVec[idx]);
1475 
1476  // fNSimChannelHitsVec.push_back(std::min(nSimChannelHitVec[idx],1999));
1477  //fNRecobHitVec.push_back(std::min(nRecobHitVec[idx],1999));
1478  fNSimChannelHitsVec.push_back(nSimChannelHitVec[idx]);
1479  fNRecobHitVec.push_back(nRecobHitVec[idx]);
1480  fNFakeHitVec.push_back(nFakeHitVec[idx]/(float)nSimulatedWiresVec[idx]);
1481  fHitEfficiencyVec.push_back(hitEfficiency);
1482  }
1483  }
1484 
1485 
1486  // std::cout << "# wire: \t "<< wireno << std::endl;
1487  //MF_LOG_DEBUG("ThroughgoingmuonAnalyzer")
1488  //<< "Event" << fEvent
1489  //<< "Run: " << fRun
1490  //<< "SubRun"<< fSubRun;
1491 
1492  std::cout << "-- fill the tree" << std::endl;
1493 
1494  fTree->Fill();
1495 
1496  std::cout << "-- done " << std::endl;
1497 
1498  return;
1499 }
1500 DEFINE_ART_MODULE(ThroughgoingmuonAnalyzer)
1501 } // end of namespace
geo::TPCGeo::ID
geo::TPCID const & ID() const
Returns the identifier of this TPC.
Definition: TPCGeo.h:333
recob::Hit::LocalIndex
short int LocalIndex() const
How well do we believe we know this hit?
Definition: Hit.h:227
LArG4Parameters.h
Store parameters for running LArG4.
z
process_name opflash particleana ie ie ie z
Definition: run_opflash_electron.fcl:92
thrugoingmuon::ThroughgoingmuonAnalyzer::tfs
art::ServiceHandle< art::TFileService > tfs
Definition: ThroughgoingmuonAnalyzer_module.cc:229
ServiceUtil.h
Utilities related to art service access.
Plane
BEGIN_PROLOG true icarus_rawdigitfilter FilterTools FilterPlane1 Plane
Definition: rawdigitfilter_icarus.fcl:8
CryostatGeo.h
Encapsulate the construction of a single cyostat.
thrugoingmuon::ThroughgoingmuonAnalyzer::ftotalelctroncut
float ftotalelctroncut
Threshold for electron (10k)
Definition: ThroughgoingmuonAnalyzer_module.cc:95
x
process_name opflash particleana ie x
Definition: run_opflash_electron.fcl:90
pz
BEGIN_PROLOG pz
Definition: filemuons.fcl:10
recob::Track::MomentumVectorAtPoint
Vector_t const & MomentumVectorAtPoint(size_t i) const
Definition: bj/lardataobj/RecoBase/Track.h:148
recob::Hit::WireID
geo::WireID WireID() const
Definition: Hit.h:233
recob::Hit::RMS
float RMS() const
RMS of the hit shape, in tick units.
Definition: Hit.h:220
Hit.h
Declaration of signal hit object.
recob::Track::LocationAtPoint
Point_t const & LocationAtPoint(size_t i) const
Definition: bj/lardataobj/RecoBase/Track.h:126
geo::TPCGeo
Geometry information for a single TPC.
Definition: TPCGeo.h:38
thrugoingmuon::ThroughgoingmuonAnalyzer::fGeometry
const geo::GeometryCore * fGeometry
pointer to Geometry service
Definition: ThroughgoingmuonAnalyzer_module.cc:239
recob::Track::NumberTrajectoryPoints
size_t NumberTrajectoryPoints() const
Various functions related to the presence and the number of (valid) points.
Definition: bj/lardataobj/RecoBase/Track.h:102
geo::TPCGeo::ActiveBoundingBox
geo::BoxBoundedGeo const & ActiveBoundingBox() const
Returns the box of the active volume of this TPC.
Definition: TPCGeo.h:320
geo::GeometryCore::PositionToTPCptr
geo::TPCGeo const * PositionToTPCptr(geo::Point_t const &point) const
Returns the TPC at specified location.
Definition: GeometryCore.cxx:410
recob::Hit::DegreesOfFreedom
int DegreesOfFreedom() const
Definition: Hit.h:229
geo::GeometryCore::ChannelToWire
std::vector< geo::WireID > ChannelToWire(raw::ChannelID_t const channel) const
Returns a list of wires connected to the specified TPC channel.
Definition: GeometryCore.cxx:1196
RawDigit.h
Definition of basic raw digits.
recob::Hit::Integral
float Integral() const
Integral under the calibrated signal waveform of the hit, in tick x ADC units.
Definition: Hit.h:224
track
process_name use argoneut_mc_hitfinder track
Definition: standard_reco.fcl:36
hit
process_name hit
Definition: cheaterreco.fcl:51
recob::Hit::GoodnessOfFit
float GoodnessOfFit() const
Degrees of freedom in the determination of the hit signal shape (-1 by default)
Definition: Hit.h:228
recob::Hit::Multiplicity
short int Multiplicity() const
How many hits could this one be shared with.
Definition: Hit.h:226
TPC
BEGIN_PROLOG TPC
Definition: sbnd_flashalgo.fcl:14
recob::Hit::PeakAmplitude
float PeakAmplitude() const
The estimated amplitude of the hit at its peak, in ADC units.
Definition: Hit.h:221
thrugoingmuon::ThroughgoingmuonAnalyzer::operator=
ThroughgoingmuonAnalyzer & operator=(ThroughgoingmuonAnalyzer const &)=delete
geo::GeometryCore::Nplanes
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
Definition: GeometryCore.h:2854
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: bj/lardataobj/RecoBase/Track.h:167
y
process_name opflash particleana ie ie y
Definition: run_opflash_electron.fcl:91
thrugoingmuon::ThroughgoingmuonAnalyzer::fSimTrackID
int fSimTrackID
GEANT ID of the particle being processed.
Definition: ThroughgoingmuonAnalyzer_module.cc:108
recob::Track::Theta
double Theta() const
Access to spherical or geographical angles at vertex or at any point.
Definition: bj/lardataobj/RecoBase/Track.h:176
thrugoingmuon::ThroughgoingmuonAnalyzer::ThroughgoingmuonAnalyzer
ThroughgoingmuonAnalyzer(fhicl::ParameterSet const &pset)
Definition: ThroughgoingmuonAnalyzer_module.cc:246
thrugoingmuon::ThroughgoingmuonAnalyzer::fSigmaVec
std::vector< float > fSigmaVec
Window size for matching to SimChannels.
Definition: ThroughgoingmuonAnalyzer_module.cc:93
util::quantities::tick
tick_as<> tick
Tick number, represented by std::ptrdiff_t.
Definition: electronics.h:75
lariov::ChannelStatusProvider::Status
virtual Status_t Status(raw::ChannelID_t channel) const
Returns a status integer with arbitrary meaning.
Definition: ChannelStatusProvider.h:89
lariov::ChannelStatusProvider
Class providing information about the quality of channels.
Definition: ChannelStatusProvider.h:49
geo::PlaneID::Plane
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:1476
thrugoingmuon::ThroughgoingmuonAnalyzer::fSimPDG
int fSimPDG
PDG ID of the particle being processed.
Definition: ThroughgoingmuonAnalyzer_module.cc:107
Track.h
Provides recob::Track data product.
geo::vect::norm
auto norm(Vector const &v)
Return norm of the specified vector.
Definition: geo_vectors_utils.h:1202
dropbox.dir
tuple dir
Definition: dropbox.py:28
px
BEGIN_PROLOG px
Definition: filemuons.fcl:10
recob::Hit::PeakTime
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
py
process_name physics producers generator physics producers generator physics producers generator py
Definition: prod_muon_2GeV_d45_Collection.fcl:21
icarus::trigger::to_string
std::string to_string(WindowPattern const &pattern)
Definition: WindowPattern.h:115
SimEnergyDeposit.h
contains information for a single step in the detector simulation
thrugoingmuon::ThroughgoingmuonAnalyzer::fMinAllowedChanStatus
int fMinAllowedChanStatus
Don&#39;t consider channels with lower status.
Definition: ThroughgoingmuonAnalyzer_module.cc:94
SimChannel.h
object containing MC truth information necessary for making RawDigits and doing back tracking ...
ChannelStatusProvider.h
Interface for experiment-specific channel quality info provider.
e
do i e
Definition: updateVersionList.sh:68
recob::Hit::SummedADC
float SummedADC() const
The sum of calibrated ADC counts of the hit (0. by default)
Definition: Hit.h:223
Wire.h
Declaration of basic channel signal object.
recob::Hit
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
tfs
art::ServiceHandle< art::TFileService > tfs
Definition: NuVertexChargeTree_module.cc:254
tca::evt
TCEvent evt
Definition: DataStructs.cxx:8
ChannelStatusService.h
Interface for experiment-specific service for channel quality info.
geo::Point_t
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double >, ROOT::Math::GlobalCoordinateSystemTag > Point_t
Type for representation of position in physical 3D space.
Definition: geo_vectors.h:184
thrugoingmuon::ThroughgoingmuonAnalyzer::fOffsetVec
std::vector< int > fOffsetVec
Allow offsets for each plane.
Definition: ThroughgoingmuonAnalyzer_module.cc:92
thrugoingmuon::ThroughgoingmuonAnalyzer::fEvent
int fEvent
number of the event being processed
Definition: ThroughgoingmuonAnalyzer_module.cc:100
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
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a &quot;fitted&quot; track:
Definition: bj/lardataobj/RecoBase/Track.h:49
geo::CryostatID
The data type to uniquely identify a cryostat.
Definition: geo_types.h:190
Generated by   doxygen 1.8.5