anatree_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 // Class: anatree
3 // Plugin Type: analyzer (art v2_11_02)
4 // File: anatree_module.cc
5 //
6 // Generated at Mon Aug 27 16:41:13 2018 by Thomas Junk using cetskelgen
7 // from cetlib version v3_03_01.
8 // Additions from Leo Bellantoni, 2019-21
9 ////////////////////////////////////////////////////////////////////////////////
10 
11 #include "art/Framework/Core/EDAnalyzer.h"
12 #include "art/Framework/Core/ModuleMacros.h"
13 #include "art/Framework/Principal/Event.h"
14 #include "art/Framework/Principal/Handle.h"
15 #include "art/Framework/Principal/Run.h"
16 #include "art/Framework/Principal/SubRun.h"
17 #include "art_root_io/TFileService.h"
18 #include "canvas/Utilities/InputTag.h"
19 #include "fhiclcpp/ParameterSet.h"
20 #include "messagefacility/MessageLogger/MessageLogger.h"
21 #include "art/Persistency/Common/PtrMaker.h"
22 #include "canvas/Persistency/Common/FindOne.h"
23 #include "canvas/Persistency/Common/FindOneP.h"
24 #include "canvas/Persistency/Common/FindMany.h"
25 #include "canvas/Persistency/Common/FindManyP.h"
26 
27 #include "SummaryDataProducts/RunData.h"
28 #include "SummaryDataProducts/POTSummary.h"
29 
30 #include "nusimdata/SimulationBase/GTruth.h"
31 #include "nusimdata/SimulationBase/MCTruth.h"
32 #include "nusimdata/SimulationBase/MCParticle.h"
33 #include "MCCheater/BackTracker.h"
34 #include "SimulationDataProducts/GenieParticle.h"
35 #include "SimulationDataProducts/EnergyDeposit.h"
36 #include "SimulationDataProducts/SimChannel.h"
37 #include "SimulationDataProducts/CaloDeposit.h"
38 
39 #include "ReconstructionDataProducts/TPCCluster.h"
40 #include "ReconstructionDataProducts/Hit.h"
41 #include "ReconstructionDataProducts/Track.h"
42 #include "ReconstructionDataProducts/TrackTrajectory.h"
43 #include "ReconstructionDataProducts/TrackIoniz.h"
44 #include "ReconstructionDataProducts/Vertex.h"
45 #include "ReconstructionDataProducts/Vee.h"
46 #include "ReconstructionDataProducts/CaloHit.h"
47 #include "ReconstructionDataProducts/Cluster.h"
48 
49 #include "RawDataProducts/CaloRawDigit.h"
50 
51 // nutools extensions
52 #include "nurandom/RandomUtils/NuRandomService.h"
53 
54 #include "CoreUtils/ServiceUtil.h"
55 #include "Geometry/GeometryGAr.h"
56 #include "Geometry/BitFieldCoder.h"
57 
58 #include "TTree.h"
59 #include "TDatabasePDG.h"
60 #include "TParticlePDG.h"
61 #include "TH2.h"
62 #include "TFile.h"
63 #include "TVector3.h"
64 
65 #include "CLHEP/Random/RandFlat.h"
66 
67 #include <string>
68 #include <vector>
69 #include <unordered_map>
70 
71 
72 
73 namespace gar {
74 
75  typedef std::pair<float, std::string> P;
76 
77  class anatree : public art::EDAnalyzer {
78  public:
79  explicit anatree(fhicl::ParameterSet const & p);
80  // The compiler-generated destructor is fine for non-base
81  // classes without bare pointers or other resource use.
82 
83  // Plugins should not be copied or assigned.
84  anatree(anatree const &) = delete;
85  anatree(anatree &&) = delete;
86  anatree & operator = (anatree const &) = delete;
87  anatree & operator = (anatree &&) = delete;
88 
89  virtual void endRun(art::Run const& run) override;
90  virtual void beginJob() override;
91 
92  // Required functions.
93  void analyze(art::Event const & e) override;
94 
95  private:
96  void ClearVectors();
97 
98  //Working Horse
99  cheat::BackTrackerCore* BackTrack;
100  void FillGeneratorMonteCarloInfo(art::Event const & e);
101  void FillRawInfo(art::Event const & e);
102  void FillRecoInfo(art::Event const & e);
103  void FillHighLevelRecoInfo(art::Event const & e);
104 
105  //Helpers
106  void processIonizationInfo(rec::TrackIoniz& ion, float ionizeTruncate,
107  float& forwardIonVal, float& backwardIonVal);
108  float processOneDirection(std::vector<std::pair<float,float>> SigData,
109  float ionizeTruncate);
110  // Helper method for processOneDirection
111  static bool lessThan_byE(std::pair<float,float> a, std::pair<float,float> b)
112  {return a.first < b.first;}
113 
114  // Compute T for coherent pion analysis
115  float computeT( simb::MCTruth theMCTruth );
116 
117  // Calculate the track PID based on reco momentum and Tom's parametrization
118  std::vector< std::pair<int, float> > processPIDInfo( float p );
119 
120  // Position of TPC from geometry service; 1 S Boston Ave.
121  float ItsInTulsa[3];
122  float xTPC;
123  float rTPC;
124 
125  // Input data labels
126  std::vector<std::string> fGeneratorLabels;
127  std::vector<std::string> fGENIEGeneratorLabels;
128  std::string fPOTtag;
129 
130  std::string fInstanceLabelCalo; ///< Instance name for ECAL
131  std::string fInstanceLabelMuID; ///< Instance name for MuID
132 
133  std::string fGeantLabel; ///< module label for geant4 simulated hits
134 
135  std::string fHitLabel; ///< module label for reco TPC hits rec::Hit
136  std::string fTPCClusterLabel; ///< module label for TPC Clusters rec::TPCCluster
137  std::string fTrackLabel; ///< module label for TPC Tracks rec:Track
138  std::string fTrackTrajectoryLabel; ///< module label for TPC Track Trajectories rec:TrackTrajectory
139  std::string fVertexLabel; ///< module label for vertexes rec:Vertex
140  std::string fVeeLabel; ///< module label for conversion/decay vertexes rec:Vee
141 
142  std::string fRawCaloHitLabel; ///< module label for digitized calo hits raw::CaloRawDigit
143  std::string fRawMuIDHitLabel;
144 
145  std::string fCaloHitLabel; ///< module label for reco calo hits rec::CaloHit
146  std::string fMuIDHitLabel;
147 
148  std::string fClusterLabel; ///< module label for calo clusters rec::Cluster
149  std::string fClusterMuIDLabel; ///< module label for calo clusters rec::Cluster in MuID
150  std::string fPFLabel; ///< module label for reco particles rec::PFParticle
151  std::string fECALAssnLabel; ///< module label for track-clusters associations
152 
153  // Optionally keep/drop parts of the analysis tree
154  bool fWriteMCinfo; ///< Info from MCTruth, GTruth Default=true
155  bool fWriteMCPTrajectory; ///< Write MCP Trajectory Default=true
156  bool fWriteMCPTrajMomenta; ///< Write Momenta associated with MCP Trajectory Default=false
157  bool fWriteMCCaloInfo; ///< Write MC info for calorimeter Default=true
158  float fMatchMCPtoVertDist; ///< MCParticle to MC vertex match Default=roundoff
159 
160  bool fWriteHits; ///< Write info about TPC Hits Default=false
161  bool fWriteTPCClusters; ///< Write TPCClusters info Default=true
162  bool fWriteTracks; ///< Start/end X, P for tracks Default=true
163  bool fWriteTrackTrajectories; ///< Point traj of reco tracks Default=false
164  bool fWriteVertices; ///< Reco vertexes & their tracks Default=true
165  bool fWriteVees; ///< Reco vees & their tracks Default=true
166 
167  bool fWriteMuID;
168  bool fWriteCaloDigits; ///< Raw digits for calorimetry. Default=false
169  bool fWriteCaloHits; ///< Write ECAL hits. Default=true
170  bool fWriteCaloClusters; ///< Write ECAL clusters. Default=true
171  bool fWriteMatchedTracks; ///< Write ECAL-track Assns Default=true
172 
173  // Truncation parameter for dE/dx (average this fraction of the lowest readings)
174  float fIonizTruncate; ///< Default=1.00;
175  bool fWriteCohInfo; ///< MC level t for coherent pi+. Default=false
176 
177  // the analysis tree
178  TTree *fTree;
179 
180  //Geometry
181  const geo::GeometryCore* fGeo; ///< pointer to the geometry
182  std::string fECALEncoding;
183  std::string fMuIDEncoding;
184  gar::geo::BitFieldCoder *fFieldDecoder_ECAL;
185  gar::geo::BitFieldCoder *fFieldDecoder_MuID;
186 
187  typedef int TrkId;
188  std::unordered_map<TrkId, Int_t> TrackIdToIndex;
189 
190 
191  // global event info
192  Int_t fEvent; ///< number of the event being processed
193  Int_t fRun; ///< number of the run being processed
194  Int_t fSubRun; ///< number of the sub-run being processed
195  Float_t fTPC_X; ///< center of TPC stored as per-event & compressed by root
196  Float_t fTPC_Y;
197  Float_t fTPC_Z;
198  Int_t fTotalPOT;
199  Int_t fNSpills;
200 
201  // MCTruth data.
202  // GENIE kinematics computed ignoring Fermi momentum and the off-shellness
203  // of the bound nucleon. Well, that's what the documentation says. Might
204  // not be true! But to get the right kinematics here, t has to be computed.
205  // Mode and InteractionType given in simb::MCNeutrino.h of the nusimdata
206  // product.
207  // Use Rtypes.h here, as these data get used by root
208 
209  std::vector<Int_t> fNeutrinoType;
210  std::vector<Int_t> fCCNC;
211  std::vector<Int_t> fMode;
212  std::vector<Int_t> fInteractionType;
213  std::vector<Float_t> fQ2;
214  std::vector<Float_t> fW;
215  std::vector<Float_t> fX;
216  std::vector<Float_t> fY;
217  std::vector<Float_t> fTheta;
218  std::vector<Float_t> fT;
219  std::vector<Float_t> fMCVertexX;
220  std::vector<Float_t> fMCVertexY;
221  std::vector<Float_t> fMCVertexZ;
222  std::vector<Float_t> fMCnuPx;
223  std::vector<Float_t> fMCnuPy;
224  std::vector<Float_t> fMCnuPz;
225 
226  // GTruth data
227  std::vector<Int_t> fGint;
228  std::vector<Int_t> fTgtPDG;
229  std::vector<Float_t> fWeight;
230  std::vector<Float_t> fgT;
231 
232  //GENIE particle list from event record
233  std::vector<Int_t> fnGPart;
234  std::vector<Int_t> fGPartIntIdx;
235  std::vector<Int_t> fGPartIdx;
236  std::vector<Int_t> fGPartPdg;
237  std::vector<Int_t> fGPartStatus;
238  std::vector<Int_t> fGPartFirstMom;
239  std::vector<Int_t> fGPartLastMom;
240  std::vector<Int_t> fGPartFirstDaugh;
241  std::vector<Int_t> fGPartLastDaugh;
242  std::vector<std::string> fGPartName;
243  std::vector<Float_t> fGPartPx;
244  std::vector<Float_t> fGPartPy;
245  std::vector<Float_t> fGPartPz;
246  std::vector<Float_t> fGPartE;
247  std::vector<Float_t> fGPartMass;
248 
249  // MCParticle data
250  std::vector<Int_t> fMCTrkID;
251  std::vector<Int_t> fMCPDG;
252  std::vector<Int_t> fMCMotherIndex;
253  std::vector<Int_t> fMCMotherTrkID;
254  std::vector<Int_t> fMCPDGMother;
255  std::vector<Float_t> fMCPStartX;
256  std::vector<Float_t> fMCPStartY;
257  std::vector<Float_t> fMCPStartZ;
258  std::vector<Float_t> fMCPTime;
259  std::vector<Float_t> fMCPStartPX;
260  std::vector<Float_t> fMCPStartPY;
261  std::vector<Float_t> fMCPStartPZ;
262  std::vector<Float_t> fMCPEndX;
263  std::vector<Float_t> fMCPEndY;
264  std::vector<Float_t> fMCPEndZ;
265  std::vector<Float_t> fMCPEndPX;
266  std::vector<Float_t> fMCPEndPY;
267  std::vector<Float_t> fMCPEndPZ;
268  std::vector<std::string> fMCPProc;
269  std::vector<std::string> fMCPEndProc;
270  std::vector<Int_t> fMCPVertIndex;
271 
272  // track trajectory of MCP
273  std::vector<Float_t> fTrajMCPX;
274  std::vector<Float_t> fTrajMCPY;
275  std::vector<Float_t> fTrajMCPZ;
276  std::vector<Float_t> fTrajMCPT;
277  std::vector<Float_t> fTrajMCPE;
278  std::vector<Int_t> fTrajMCPIndex;
279  std::vector<Int_t> fTrajMCPTrackID;
280  std::vector<Float_t> fTrajMCPPX;
281  std::vector<Float_t> fTrajMCPPY;
282  std::vector<Float_t> fTrajMCPPZ;
283 
284  // sim calo hit data
285  UInt_t fSimnHits;
286  std::vector<Float_t> fSimHitX;
287  std::vector<Float_t> fSimHitY;
288  std::vector<Float_t> fSimHitZ;
289  std::vector<Float_t> fSimHitTime;
290  std::vector<Float_t> fSimHitEnergy;
291  std::vector<Int_t> fSimHitTrackID;
292  std::vector<ULong64_t> fSimHitCellID; // ULong64_t is size_t on 64 bit machines
293  Float_t fSimEnergySum;
294 
295  //Muon system sim hits
296  UInt_t fSimnHits_MuID;
297  std::vector<Float_t> fSimHitX_MuID;
298  std::vector<Float_t> fSimHitY_MuID;
299  std::vector<Float_t> fSimHitZ_MuID;
300  std::vector<Float_t> fSimHitTime_MuID;
301  std::vector<Float_t> fSimHitEnergy_MuID;
302  std::vector<Int_t> fSimHitTrackID_MuID;
303  std::vector<ULong64_t> fSimHitCellID_MuID; // ULong64_t is size_t on 64 bit machines
304  Float_t fSimEnergySum_MuID;
305 
306  // Hit data
307  std::vector<Float_t> fHitX;
308  std::vector<Float_t> fHitY;
309  std::vector<Float_t> fHitZ;
310  std::vector<Float_t> fHitSig;
311  std::vector<Float_t> fHitRMS;
312  std::vector<UInt_t> fHitChan;
313 
314  // TPCCluster data
315  std::vector<Float_t> fTPCClusterX;
316  std::vector<Float_t> fTPCClusterY;
317  std::vector<Float_t> fTPCClusterZ;
318  std::vector<Float_t> fTPCClusterSig;
319  std::vector<Float_t> fTPCClusterRMS;
320  std::vector<ULong64_t> fTPCClusterTrkIDNumber;
321  std::vector<Float_t> fTPCClusterCovXX;
322  std::vector<Float_t> fTPCClusterCovXY;
323  std::vector<Float_t> fTPCClusterCovXZ;
324  std::vector<Float_t> fTPCClusterCovYY;
325  std::vector<Float_t> fTPCClusterCovYZ;
326  std::vector<Float_t> fTPCClusterCovZZ;
327  std::vector<Int_t> fTPCClusterMCindex; // Branch index (NOT the GEANT track ID) of MCParticle
328  std::vector<Float_t> fTPCClusterMCfrac; // that best matches & fraction of ionization therefrom
329 
330  // track data
331  std::vector<ULong64_t> fTrackIDNumber;
332  std::vector<Float_t> fTrackStartX;
333  std::vector<Float_t> fTrackStartY;
334  std::vector<Float_t> fTrackStartZ;
335  std::vector<Float_t> fTrackStartPX;
336  std::vector<Float_t> fTrackStartPY;
337  std::vector<Float_t> fTrackStartPZ;
338  std::vector<Int_t> fTrackStartQ;
339 
340  std::vector<Float_t> fTrackEndX;
341  std::vector<Float_t> fTrackEndY;
342  std::vector<Float_t> fTrackEndZ;
343  std::vector<Float_t> fTrackEndPX;
344  std::vector<Float_t> fTrackEndPY;
345  std::vector<Float_t> fTrackEndPZ;
346  std::vector<Int_t> fTrackEndQ;
347 
348  std::vector<Float_t> fTrackLenF; // from forward fit, from the Beg end to the End end
349  std::vector<Float_t> fTrackLenB;
350  std::vector<Float_t> fTrackChi2F; // from forward fit, from the Beg end to the End end
351  std::vector<Float_t> fTrackChi2B;
352  std::vector<Int_t> fNTPCClustersOnTrack;
353  std::vector<Float_t> fTrackAvgIonF; // from forward fit, from the Beg end to the End end
354  std::vector<Float_t> fTrackAvgIonB;
355 
356  std::vector<Int_t> fTrackPIDF;
357  std::vector<Float_t> fTrackPIDProbF;
358  std::vector<Int_t> fTrackPIDB;
359  std::vector<Float_t> fTrackPIDProbB;
360  std::vector<Int_t> fTrackMCindex; // Branch index (NOT the GEANT track ID) of MCParticle
361  std::vector<Float_t> fTrackMCfrac; // that best matches & fraction of ionization therefrom
362 
363  //TrackTrajectory
364  std::vector<Float_t> fTrackTrajectoryFWDX; //forward
365  std::vector<Float_t> fTrackTrajectoryFWDY; //forward
366  std::vector<Float_t> fTrackTrajectoryFWDZ; //forward
367  std::vector<Int_t> fTrackTrajectoryFWDID;
368 
369  std::vector<Float_t> fTrackTrajectoryBWDX; //backward
370  std::vector<Float_t> fTrackTrajectoryBWDY; //backward
371  std::vector<Float_t> fTrackTrajectoryBWDZ; //backward
372  std::vector<Int_t> fTrackTrajectoryBWDID;
373 
374  // vertex branches
375  std::vector<ULong64_t> fVertexIDNumber;
376  std::vector<Float_t> fVertexX;
377  std::vector<Float_t> fVertexY;
378  std::vector<Float_t> fVertexZ;
379  std::vector<ULong64_t> fVertexT;
380  std::vector<Int_t> fVertexN;
381  std::vector<Int_t> fVertexQ;
382 
383  std::vector<ULong64_t> fVTAssn_VertIDNumber; // Being the vertex which this Assn belongs to
384  std::vector<ULong64_t> fVTAssn_TrackIDNumber;
385  std::vector<gar::rec::TrackEnd> fVTAssn_TrackEnd;
386 
387  // Vee branches
388  std::vector<ULong64_t> fVeeIDNumber;
389  std::vector<Float_t> fVeeX;
390  std::vector<Float_t> fVeeY;
391  std::vector<Float_t> fVeeZ;
392  std::vector<ULong64_t> fVeeT;
393  std::vector<Float_t> fVeePXKpipi;
394  std::vector<Float_t> fVeePYKpipi;
395  std::vector<Float_t> fVeePZKpipi;
396  std::vector<Float_t> fVeeEKpipi;
397  std::vector<Float_t> fVeeMKpipi;
398  std::vector<Float_t> fVeePXLppi;
399  std::vector<Float_t> fVeePYLppi;
400  std::vector<Float_t> fVeePZLppi;
401  std::vector<Float_t> fVeeELppi;
402  std::vector<Float_t> fVeeMLppi;
403  std::vector<Float_t> fVeePXLpip;
404  std::vector<Float_t> fVeePYLpip;
405  std::vector<Float_t> fVeePZLpip;
406  std::vector<Float_t> fVeeELpip;
407  std::vector<Float_t> fVeeMLpip;
408  std::vector<ULong64_t> fVeeTAssn_VeeIDNumber; // Being the Vee which this Assn belongs to
409  std::vector<ULong64_t> fVeeTAssn_TrackIDNumber;
410  std::vector<gar::rec::TrackEnd> fVeeTAssn_TrackEnd;
411 
412  // raw calo digits data
413  UInt_t fDiginHits;
414  std::vector<Float_t> fDigiHitX;
415  std::vector<Float_t> fDigiHitY;
416  std::vector<Float_t> fDigiHitZ;
417  std::vector<Float_t> fDigiHitTime;
418  std::vector<UInt_t> fDigiHitADC; // UInt_t is unsigned 32 bit integer
419  std::vector<ULong64_t> fDigiHitCellID;
420 
421  //Muon system raw hits
422  UInt_t fDiginHits_MuID;
423  std::vector<Float_t> fDigiHitX_MuID;
424  std::vector<Float_t> fDigiHitY_MuID;
425  std::vector<Float_t> fDigiHitZ_MuID;
426  std::vector<Float_t> fDigiHitTime_MuID;
427  std::vector<UInt_t> fDigiHitADC_MuID; // UInt_t is unsigned 32 bit integer
428  std::vector<ULong64_t> fDigiHitCellID_MuID;
429 
430  // reco calo hit data
431  UInt_t fReconHits;
432  std::vector<ULong64_t> fReconHitIDNumber;
433  std::vector<Float_t> fRecoHitX;
434  std::vector<Float_t> fRecoHitY;
435  std::vector<Float_t> fRecoHitZ;
436  std::vector<Float_t> fRecoHitTime;
437  std::vector<Float_t> fRecoHitEnergy; // Hit energies have the sampling fraction
438  std::vector<ULong64_t> fRecoHitCellID; // correction factors in Reco/SiPMHitFinder.fcl
439  std::vector<Int_t> fRecoHitLayer; // These are propagated to clusters and the same
440  Float_t fRecoEnergySum; // is done for the NuID system.
441 
442  //Muon system reco hits
443  UInt_t fReconHits_MuID;
444  std::vector<ULong64_t> fReconHitIDNumber_MuID;
445  std::vector<Float_t> fRecoHitX_MuID;
446  std::vector<Float_t> fRecoHitY_MuID;
447  std::vector<Float_t> fRecoHitZ_MuID;
448  std::vector<Float_t> fRecoHitTime_MuID;
449  std::vector<Float_t> fRecoHitEnergy_MuID;
450  std::vector<ULong64_t> fRecoHitCellID_MuID;
451  std::vector<Int_t> fRecoHitLayer_MuID;
452  Float_t fRecoEnergySum_MuID;
453 
454  // calo cluster data
455  UInt_t fnCluster;
456  std::vector<ULong64_t> fClusterIDNumber;
457  std::vector<UInt_t> fClusterNhits;
458  std::vector<Float_t> fClusterEnergy;
459  std::vector<Float_t> fClusterTime;
460  std::vector<Float_t> fClusterTimeDiffFirstLast;
461  std::vector<Float_t> fClusterX;
462  std::vector<Float_t> fClusterY;
463  std::vector<Float_t> fClusterZ;
464  std::vector<Float_t> fClusterTheta;
465  std::vector<Float_t> fClusterPhi;
466  std::vector<Float_t> fClusterPID;
467  // std::vector<Float_t> fClusterShape;
468  std::vector<Float_t> fClusterMainAxisX;
469  std::vector<Float_t> fClusterMainAxisY;
470  std::vector<Float_t> fClusterMainAxisZ;
471  std::vector<Int_t> fClusterMCindex; // Branch index (NOT the GEANT track ID) of MCParticle
472  std::vector<Float_t> fClusterMCfrac; // that best matches & fraction of ionization therefrom
473 
474  // calo cluster data
475  UInt_t fnCluster_MuID;
476  std::vector<ULong64_t> fClusterIDNumber_MuID;
477  std::vector<UInt_t> fClusterNhits_MuID;
478  std::vector<Float_t> fClusterEnergy_MuID;
479  std::vector<Float_t> fClusterTime_MuID;
480  std::vector<Float_t> fClusterTimeDiffFirstLast_MuID;
481  std::vector<Float_t> fClusterX_MuID;
482  std::vector<Float_t> fClusterY_MuID;
483  std::vector<Float_t> fClusterZ_MuID;
484  std::vector<Float_t> fClusterTheta_MuID;
485  std::vector<Float_t> fClusterPhi_MuID;
486  std::vector<Float_t> fClusterPID_MuID;
487  // std::vector<Float_t> fClusterShape;
488  std::vector<Float_t> fClusterMainAxisX_MuID;
489  std::vector<Float_t> fClusterMainAxisY_MuID;
490  std::vector<Float_t> fClusterMainAxisZ_MuID;
491  std::vector<Int_t> fClusterMCindex_MuID; // Branch index (NOT the GEANT track ID) of MCParticle
492  std::vector<Float_t> fClusterMCfrac_MuID; // that best matches & fraction of ionization therefrom
493 
494  // ECAL cluster to ECAL hits association info
495  std::vector<std::vector<ULong64_t>> fClusterAssn_RecoHitIDNumber;
496 
497  // MuID cluster to MuID hits association info
498  std::vector<std::vector<ULong64_t>> fClusterMuIDAssn_MuIDHitIDNumber;
499 
500  // ECAL cluster to track association info
501  std::vector<ULong64_t> fCALAssn_ClusIDNumber; // Being the cluster which this Assn belongs to
502  std::vector<ULong64_t> fCALAssn_TrackIDNumber; // The rec::TrackEnd (see Track.h) that extrapolated to cluster
503  std::vector<gar::rec::TrackEnd> fCALAssn_TrackEnd;
504 
505  //map of PID TH2 per momentum value
506  CLHEP::HepRandomEngine &fEngine; ///< random engine
507  std::unordered_map<int, TH2F*> m_pidinterp;
508  };
509 }
510 
511 
512 
513 //==============================================================================
514 //==============================================================================
515 //==============================================================================
516 // constructor
517 gar::anatree::anatree(fhicl::ParameterSet const & p)
518  : EDAnalyzer(p),
519  fEngine(art::ServiceHandle<rndm::NuRandomService>()->createEngine(*this, p, "Seed")) {
520  fGeo = gar::providerFrom<geo::GeometryGAr>();
521 
522  fECALEncoding = fGeo->GetECALCellIDEncoding();
523  fFieldDecoder_ECAL = new gar::geo::BitFieldCoder( fECALEncoding );
524 
525  if(fGeo->HasMuonDetector()) {
526  fMuIDEncoding = fGeo->GetMuIDCellIDEncoding();
527  fFieldDecoder_MuID = new gar::geo::BitFieldCoder( fMuIDEncoding );
528  }
529 
530  bool usegenlabels =
531  p.get_if_present<std::vector<std::string> >("GeneratorLabels",fGeneratorLabels);
532  if (!usegenlabels) fGeneratorLabels.clear();
533 
534  bool usegeniegenlabels =
535  p.get_if_present<std::vector<std::string> >("GENIEGeneratorLabels",fGENIEGeneratorLabels);
536  if (!usegeniegenlabels) fGENIEGeneratorLabels.clear();
537 
538  fPOTtag = p.get<std::string>("SummaryDataLabel","generator");
539 
540  //Sim Hits
541  fGeantLabel = p.get<std::string>("GEANTLabel","geant");
542  fInstanceLabelCalo = p.get<std::string>("InstanceLabelCalo","ECAL");
543  fInstanceLabelMuID = p.get<std::string>("InstanceLabelMuID","MuID");
544 
545  fHitLabel = p.get<std::string>("HitLabel","hit");
546  fTPCClusterLabel = p.get<std::string>("TPCClusterLabel","tpccluster");
547  fTrackLabel = p.get<std::string>("TrackLabel","track");
548  fTrackTrajectoryLabel = p.get<std::string>("TrackTrajectoryLabel","track");
549  fVertexLabel = p.get<std::string>("VertexLabel","vertex");
550  fVeeLabel = p.get<std::string>("VeeLabel","veefinder1");
551 
552  //Calorimetric related ECAL/MuID
553  fRawCaloHitLabel = p.get<std::string>("RawCaloHitLabel","daqsipm");
554  fRawMuIDHitLabel = p.get<std::string>("RawMuIDHitLabel","daqsipmmuid");
555  fCaloHitLabel = p.get<std::string>("CaloHitLabel","sscalohit");
556  fMuIDHitLabel = p.get<std::string>("MuIDHitLabel","sscalohitmuid");
557 
558  fClusterLabel = p.get<std::string>("ClusterLabel","calocluster");
559  fClusterMuIDLabel = p.get<std::string>("MuIDClusterLabel","caloclustermuid");
560  fPFLabel = p.get<std::string>("PFLabel","pandora");
561  fECALAssnLabel = p.get<std::string>("ECALAssnLabel","trkecalassn");
562 
563  // What to write
564  fWriteMCinfo = p.get<bool>("WriteMCinfo", true);
565  fWriteMCPTrajectory = p.get<bool>("WriteMCPTrajectory", true);
566  fWriteMCPTrajMomenta = p.get<bool>("WriteMCPTrajMomenta",false);
567  fWriteMCCaloInfo = p.get<bool>("WriteMCCaloInfo", true);
568  float MCPtoVertDefault = 10.0*std::numeric_limits<Float_t>::epsilon();
569  fMatchMCPtoVertDist = p.get<float>("MatchMCPtoVertDist",MCPtoVertDefault);
570 
571  fWriteHits = p.get<bool>("WriteHits", false);
572  fWriteTPCClusters = p.get<bool>("WriteTPCClusters", true);
573  fWriteTracks = p.get<bool>("WriteTracks", true);
574  fWriteTrackTrajectories = p.get<bool>("WriteTrackTrajectories", false);
575  fWriteVertices = p.get<bool>("WriteVertices", true);
576  fWriteVees = p.get<bool>("WriteVees", true);
577 
578  fWriteMuID = p.get<bool>("WriteMuID", true);
579  fWriteCaloDigits = p.get<bool>("WriteCaloDigits", false);
580  fWriteCaloHits = p.get<bool>("WriteCaloHits", true);
581  fWriteCaloClusters = p.get<bool>("WriteCaloClusters", true);
582  fWriteMatchedTracks = p.get<bool>("WriteMatchedTracks",true);
583 
584  fIonizTruncate = p.get<float>("IonizTruncate", 0.70);
585  fWriteCohInfo = p.get<bool> ("WriteCohInfo", false);
586 
587  if (usegenlabels) {
588  for (size_t i=0; i<fGeneratorLabels.size(); ++i) {
589  consumes<std::vector<simb::MCTruth> >(fGeneratorLabels.at(i));
590  }
591  } else {
592  consumesMany<std::vector<simb::MCTruth> >();
593  }
594 
595  if (usegeniegenlabels) {
596  for (size_t i=0; i<fGENIEGeneratorLabels.size(); ++i) {
597  consumes<std::vector<simb::GTruth> >(fGENIEGeneratorLabels.at(i));
598  }
599  } else {
600  consumesMany<std::vector<simb::GTruth> >();
601  }
602 
603  consumesMany<std::vector<sdp::GenieParticle> >();
604  //consumes<art::Assns<simb::MCTruth, simb::MCParticle> >(fGeantLabel);
605  consumes<std::vector<simb::MCParticle> >(fGeantLabel);
606 
607  //TPC related
608  consumes<std::vector<sdp::EnergyDeposit> >(fGeantLabel);
609  consumes<std::vector<rec::TPCCluster> >(fTPCClusterLabel);
610  consumes<art::Assns<rec::Track, rec::TPCCluster> >(fTPCClusterLabel);
611  consumes<std::vector<rec::Hit> >(fHitLabel);
612  consumes<std::vector<rec::Track> >(fTrackLabel);
613  consumes<std::vector<rec::TrackTrajectory> >(fTrackTrajectoryLabel);
614  consumes<std::vector<rec::Vertex> >(fVertexLabel);
615  consumes<art::Assns<rec::Track, rec::Vertex> >(fVertexLabel);
616  consumes<std::vector<rec::Vee> >(fVeeLabel);
617  consumes<art::Assns<rec::Track, rec::Vee> >(fVeeLabel);
618  consumes<std::vector<rec::TrackIoniz>>(fTrackLabel);
619  consumes<art::Assns<rec::TrackIoniz, rec::Track>>(fTrackLabel);
620 
621  //Calorimetry related
622  art::InputTag ecalgeanttag(fGeantLabel, fInstanceLabelCalo);
623  consumes<std::vector<gar::sdp::CaloDeposit> >(ecalgeanttag);
624 
625  art::InputTag ecalrawtag(fRawCaloHitLabel, fInstanceLabelCalo);
626  consumes<std::vector<raw::CaloRawDigit> >(ecalrawtag);
627 
628  art::InputTag ecalhittag(fCaloHitLabel, fInstanceLabelCalo);
629  consumes<std::vector<rec::CaloHit> >(ecalhittag);
630 
631  art::InputTag ecalclustertag(fClusterLabel, fInstanceLabelCalo);
632  consumes<std::vector<rec::Cluster> >(ecalclustertag);
633  consumes<art::Assns<rec::Cluster, rec::CaloHit>>(ecalclustertag);
634 
635  //Muon system related
636  if (fGeo->HasMuonDetector() && fWriteMuID) {
637  art::InputTag muidgeanttag(fGeantLabel, fInstanceLabelMuID);
638  consumes<std::vector<gar::sdp::CaloDeposit> >(muidgeanttag);
639 
640  art::InputTag muidrawtag(fRawMuIDHitLabel, fInstanceLabelMuID);
641  consumes<std::vector<raw::CaloRawDigit> >(muidrawtag);
642 
643  art::InputTag muidhittag(fCaloHitLabel, fInstanceLabelMuID);
644  consumes<std::vector<rec::CaloHit> >(muidhittag);
645 
646  art::InputTag muidclustertag(fClusterMuIDLabel, fInstanceLabelMuID);
647  consumes<std::vector<rec::Cluster> >(muidclustertag);
648  consumes<art::Assns<rec::Cluster, rec::CaloHit>>(muidclustertag);
649  }
650 
651  consumes<art::Assns<rec::Cluster, rec::Track>>(fECALAssnLabel);
652 
653  return;
654 } // end constructor
655 
656 
657 
658 //==============================================================================
659 //==============================================================================
660 //==============================================================================
661 void gar::anatree::beginJob() {
662 
663  fTPC_X = ItsInTulsa[0] = fGeo->TPCXCent();
664  fTPC_Y = ItsInTulsa[1] = fGeo->TPCYCent();
665  fTPC_Z = ItsInTulsa[2] = fGeo->TPCZCent();
666 
667  xTPC = fGeo->TPCLength() / 2.;
668  rTPC = fGeo->TPCRadius();
669 
670  art::ServiceHandle<art::TFileService> tfs;
671  fTree = tfs->make<TTree>("GArAnaTree","GArAnaTree");
672 
673  fTree->Branch("Run", &fRun, "Run/I");
674  fTree->Branch("SubRun", &fSubRun, "SubRun/I");
675  fTree->Branch("Event", &fEvent, "Event/I");
676  fTree->Branch("TPC_X", &fTPC_X, "TPC_X/F");
677  fTree->Branch("TPC_Y", &fTPC_Y, "TPC_Y/F");
678  fTree->Branch("TPC_Z", &fTPC_Z, "TPC_Z/F");
679  fTree->Branch("POT", &fTotalPOT, "POT/I");
680  fTree->Branch("NSpills", &fNSpills, "NSpills/I");
681 
682  if (fWriteMCinfo) {
683  fTree->Branch("NType", &fNeutrinoType);
684  fTree->Branch("CCNC", &fCCNC);
685  fTree->Branch("Mode", &fMode);
686  fTree->Branch("InterT", &fInteractionType);
687  fTree->Branch("MC_Q2", &fQ2);
688  fTree->Branch("MC_W", &fW);
689  fTree->Branch("MC_X", &fX);
690  fTree->Branch("MC_Y", &fY);
691  fTree->Branch("MC_Theta", &fTheta);
692  if (fWriteCohInfo) fTree->Branch("MC_T", &fT);
693  fTree->Branch("MCVertX", &fMCVertexX);
694  fTree->Branch("MCVertY", &fMCVertexY);
695  fTree->Branch("MCVertZ", &fMCVertexZ);
696  fTree->Branch("MCNuPx", &fMCnuPx);
697  fTree->Branch("MCNuPy", &fMCnuPy);
698  fTree->Branch("MCNuPz", &fMCnuPz);
699 
700  fTree->Branch("Gint", &fGint);
701  fTree->Branch("TgtPDG", &fTgtPDG);
702  fTree->Branch("Weight", &fWeight);
703  fTree->Branch("GT_T", &fgT);
704 
705  //GENIE particle list from the event record
706  fTree->Branch("nGPart", &fnGPart);
707  fTree->Branch("GPartIntIdx", &fGPartIntIdx);
708  fTree->Branch("GPartIdx", &fGPartIdx);
709  fTree->Branch("GPartName", &fGPartName);
710  fTree->Branch("GPartPdg", &fGPartPdg);
711  fTree->Branch("GPartStatus", &fGPartStatus);
712  fTree->Branch("GPartFirstMom", &fGPartFirstMom);
713  fTree->Branch("GPartLastMom", &fGPartLastMom);
714  fTree->Branch("GPartFirstDaugh",&fGPartFirstDaugh);
715  fTree->Branch("GPartLastDaugh", &fGPartLastDaugh);
716  fTree->Branch("GPartPx", &fGPartPx);
717  fTree->Branch("GPartPy", &fGPartPy);
718  fTree->Branch("GPartPz", &fGPartPz);
719  fTree->Branch("GPartE", &fGPartE);
720  fTree->Branch("GPartMass", &fGPartMass);
721 
722  fTree->Branch("MCTrkID", &fMCTrkID);
723  fTree->Branch("PDG", &fMCPDG);
724  fTree->Branch("MotherIndex", &fMCMotherIndex); // Index into these vector branches
725  fTree->Branch("MotherTrkID", &fMCMotherTrkID); // trackid of the mother
726  fTree->Branch("PDGMother", &fMCPDGMother);
727  fTree->Branch("MCPStartX", &fMCPStartX);
728  fTree->Branch("MCPStartY", &fMCPStartY);
729  fTree->Branch("MCPStartZ", &fMCPStartZ);
730  fTree->Branch("MCPTime", &fMCPTime);
731  fTree->Branch("MCPStartPX", &fMCPStartPX);
732  fTree->Branch("MCPStartPY", &fMCPStartPY);
733  fTree->Branch("MCPStartPZ", &fMCPStartPZ);
734  fTree->Branch("MCPEndX", &fMCPEndX);
735  fTree->Branch("MCPEndY", &fMCPEndY);
736  fTree->Branch("MCPEndZ", &fMCPEndZ);
737  fTree->Branch("MCPEndPX", &fMCPEndPX);
738  fTree->Branch("MCPEndPY", &fMCPEndPY);
739  fTree->Branch("MCPEndPZ", &fMCPEndPZ);
740  fTree->Branch("MCPProc", &fMCPProc);
741  fTree->Branch("MCPEndProc", &fMCPEndProc);
742  fTree->Branch("MCPVertIndex",&fMCPVertIndex);
743  }
744 
745  if (fWriteMCPTrajectory) {
746  // Checking for parameter file validity only once to simplify
747  // later code in anatree::process etc.
748  if (!fWriteMCinfo) {
749  throw cet::exception("anatree")
750  << " fWriteMCPTrajectory, but !fWriteMCinfo."
751  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
752  } else {
753  //Write MCP Trajectory
754  fTree->Branch("TrajMCPX", &fTrajMCPX);
755  fTree->Branch("TrajMCPY", &fTrajMCPY);
756  fTree->Branch("TrajMCPZ", &fTrajMCPZ);
757  fTree->Branch("TrajMCPT", &fTrajMCPT);
758  if (fWriteMCPTrajMomenta) {
759  fTree->Branch("TrajMCPPX", &fTrajMCPPX);
760  fTree->Branch("TrajMCPPY", &fTrajMCPPY);
761  fTree->Branch("TrajMCPPZ", &fTrajMCPPZ);
762  }
763  fTree->Branch("TrajMCPE", &fTrajMCPE);
764  fTree->Branch("TrajMCPIndex", &fTrajMCPIndex);
765  fTree->Branch("TrajMCPTrackID", &fTrajMCPTrackID);
766  }
767  }
768 
769  if (fWriteMCCaloInfo) {
770  if (!fWriteMCinfo) {
771  throw cet::exception("anatree")
772  << " fWriteMCCaloInfo, but !fWriteMCinfo."
773  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
774  } else {
775  // Write calorimetry MC information
776  fTree->Branch("SimnHits", &fSimnHits);
777  fTree->Branch("SimHitX", &fSimHitX);
778  fTree->Branch("SimHitY", &fSimHitY);
779  fTree->Branch("SimHitZ", &fSimHitZ);
780  fTree->Branch("SimHitTime", &fSimHitTime);
781  fTree->Branch("SimHitEnergy", &fSimHitEnergy);
782  fTree->Branch("SimHitTrkID", &fSimHitTrackID);
783  fTree->Branch("SimHitCellID", &fSimHitCellID);
784  fTree->Branch("SimEnergySum", &fSimEnergySum);
785 
786  if(fGeo->HasMuonDetector() && fWriteMuID) {
787  fTree->Branch("SimnHits_MuID", &fSimnHits_MuID);
788  fTree->Branch("SimHitX_MuID", &fSimHitX_MuID);
789  fTree->Branch("SimHitY_MuID", &fSimHitY_MuID);
790  fTree->Branch("SimHitZ_MuID", &fSimHitZ_MuID);
791  fTree->Branch("SimHitTime_MuID", &fSimHitTime_MuID);
792  fTree->Branch("SimHitEnergy_MuID", &fSimHitEnergy_MuID);
793  fTree->Branch("SimHitTrkID_MuID", &fSimHitTrackID_MuID);
794  fTree->Branch("SimHitCellID_MuID", &fSimHitCellID_MuID);
795  fTree->Branch("SimEnergySum_MuID", &fSimEnergySum_MuID);
796  }
797  }
798  }
799 
800  if (fWriteHits) {
801  fTree->Branch("HitX", &fHitX);
802  fTree->Branch("HitY", &fHitY);
803  fTree->Branch("HitZ", &fHitZ);
804  fTree->Branch("HitSig", &fHitSig);
805  fTree->Branch("HitRMS", &fHitRMS);
806  fTree->Branch("HitChan", &fHitChan);
807  }
808 
809  if (fWriteTPCClusters) {
810  fTree->Branch("TPCClusterX", &fTPCClusterX);
811  fTree->Branch("TPCClusterY", &fTPCClusterY);
812  fTree->Branch("TPCClusterZ", &fTPCClusterZ);
813  fTree->Branch("TPCClusterSig", &fTPCClusterSig);
814  fTree->Branch("TPCClusterRMS", &fTPCClusterRMS);
815  fTree->Branch("TPCClusterTrkIDNumber", &fTPCClusterTrkIDNumber);
816  fTree->Branch("TPCClusterCovXX", &fTPCClusterCovXX);
817  fTree->Branch("TPCClusterCovXY", &fTPCClusterCovXY);
818  fTree->Branch("TPCClusterCovXZ", &fTPCClusterCovXZ);
819  fTree->Branch("TPCClusterCovYY", &fTPCClusterCovYY);
820  fTree->Branch("TPCClusterCovYZ", &fTPCClusterCovYZ);
821  fTree->Branch("TPCClusterCovZZ", &fTPCClusterCovZZ);
822  fTree->Branch("TPCClusterMCindex", &fTPCClusterMCindex);
823  fTree->Branch("TPCClusterMCfrac", &fTPCClusterMCfrac);
824  }
825 
826  // All position, momentum, etc
827  if (fWriteTracks) {
828  fTree->Branch("TrackIDNumber", &fTrackIDNumber);
829 
830  fTree->Branch("TrackStartX", &fTrackStartX);
831  fTree->Branch("TrackStartY", &fTrackStartY);
832  fTree->Branch("TrackStartZ", &fTrackStartZ);
833  fTree->Branch("TrackStartPX", &fTrackStartPX);
834  fTree->Branch("TrackStartPY", &fTrackStartPY);
835  fTree->Branch("TrackStartPZ", &fTrackStartPZ);
836  fTree->Branch("TrackStartQ", &fTrackStartQ);
837 
838  fTree->Branch("TrackEndX", &fTrackEndX);
839  fTree->Branch("TrackEndY", &fTrackEndY);
840  fTree->Branch("TrackEndZ", &fTrackEndZ);
841  fTree->Branch("TrackEndPX", &fTrackEndPX);
842  fTree->Branch("TrackEndPY", &fTrackEndPY);
843  fTree->Branch("TrackEndPZ", &fTrackEndPZ);
844  fTree->Branch("TrackEndQ", &fTrackEndQ);
845 
846  fTree->Branch("TrackLenF", &fTrackLenF);
847  fTree->Branch("TrackLenB", &fTrackLenB);
848  fTree->Branch("TrackChi2F", &fTrackChi2F);
849  fTree->Branch("TrackChi2B", &fTrackChi2B);
850  fTree->Branch("NTPCClustersOnTrack",&fNTPCClustersOnTrack);
851  fTree->Branch("TrackAvgIonF", &fTrackAvgIonF);
852  fTree->Branch("TrackAvgIonB", &fTrackAvgIonB);
853 
854  fTree->Branch("TrackPIDF", &fTrackPIDF);
855  fTree->Branch("TrackPIDProbF", &fTrackPIDProbF);
856  fTree->Branch("TrackPIDB", &fTrackPIDB);
857  fTree->Branch("TrackPIDProbB", &fTrackPIDProbB);
858  fTree->Branch("TrackMCindex", &fTrackMCindex);
859  fTree->Branch("TrackMCfrac", &fTrackMCfrac);
860 
861  //Track Trajectories
862  if(fWriteTrackTrajectories) {
863  fTree->Branch("TrackTrajectoryFWDX", &fTrackTrajectoryFWDX);
864  fTree->Branch("TrackTrajectoryFWDY", &fTrackTrajectoryFWDY);
865  fTree->Branch("TrackTrajectoryFWDZ", &fTrackTrajectoryFWDZ);
866  fTree->Branch("TrackTrajectoryFWDID", &fTrackTrajectoryFWDID);
867 
868  fTree->Branch("TrackTrajectoryBWDX", &fTrackTrajectoryBWDX);
869  fTree->Branch("TrackTrajectoryBWDY", &fTrackTrajectoryBWDY);
870  fTree->Branch("TrackTrajectoryBWDZ", &fTrackTrajectoryBWDZ);
871  fTree->Branch("TrackTrajectoryBWDID", &fTrackTrajectoryBWDID);
872  }
873  }
874 
875  // Reco'd verts & their track-ends
876  if (fWriteVertices) {
877  if (!fWriteTracks) {
878  throw cet::exception("anatree")
879  << " fWriteVertices, but !fWriteTracks."
880  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
881  } else {
882  fTree->Branch("VertIDNumber", &fVertexIDNumber);
883  fTree->Branch("VertX", &fVertexX);
884  fTree->Branch("VertY", &fVertexY);
885  fTree->Branch("VertZ", &fVertexZ);
886  fTree->Branch("VertT", &fVertexT);
887  fTree->Branch("VertN", &fVertexN);
888  fTree->Branch("VertQ", &fVertexQ);
889 
890  fTree->Branch("VT_VertIDNumber", &fVTAssn_VertIDNumber);
891  fTree->Branch("VT_TrackIDNumber",&fVTAssn_TrackIDNumber);
892  fTree->Branch("VT_TrackEnd", &fVTAssn_TrackEnd);
893  }
894  }
895 
896  // Reco'd vees & their track-ends
897  if (fWriteVees) {
898  if (!fWriteTracks) {
899  throw cet::exception("anatree")
900  << " fWriteVees, but !fWriteTracks."
901  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
902  } else {
903  fTree->Branch("VeeIDNumber", &fVeeIDNumber);
904  fTree->Branch("VeeX", &fVeeX);
905  fTree->Branch("VeeY", &fVeeY);
906  fTree->Branch("VeeZ", &fVeeZ);
907  fTree->Branch("VeeT", &fVeeT);
908  fTree->Branch("VeePXKpipi", &fVeePXKpipi);
909  fTree->Branch("VeePYKpipi", &fVeePYKpipi);
910  fTree->Branch("VeePZKpipi", &fVeePZKpipi);
911  fTree->Branch("VeeEKpipi", &fVeeEKpipi);
912  fTree->Branch("VeeMKpipi", &fVeeMKpipi);
913  fTree->Branch("VeePXLppi", &fVeePXLppi);
914  fTree->Branch("VeePYLppi", &fVeePYLppi);
915  fTree->Branch("VeePZLppi", &fVeePZLppi);
916  fTree->Branch("VeeELppi", &fVeeELppi);
917  fTree->Branch("VeeMLppi", &fVeeMLppi);
918  fTree->Branch("VeePXLpip", &fVeePXLpip);
919  fTree->Branch("VeePYLpip", &fVeePYLpip);
920  fTree->Branch("VeePZLpip", &fVeePZLpip);
921  fTree->Branch("VeeELpip", &fVeeELpip);
922  fTree->Branch("VeeMLpip", &fVeeMLpip);
923  fTree->Branch("VeeT_VertIDNumber", &fVeeTAssn_VeeIDNumber);
924  fTree->Branch("VeeT_TrackIDNumber",&fVeeTAssn_TrackIDNumber);
925  fTree->Branch("VeeT_TrackEnd", &fVeeTAssn_TrackEnd);
926  }
927  }
928 
929  // Write calorimetry digits
930  if (fWriteCaloDigits) {
931  fTree->Branch("DiginHits", &fDiginHits);
932  fTree->Branch("DigiHitX", &fDigiHitX);
933  fTree->Branch("DigiHitY", &fDigiHitY);
934  fTree->Branch("DigiHitZ", &fDigiHitZ);
935  fTree->Branch("DigiHitTime", &fDigiHitTime);
936  fTree->Branch("DigiHitADC", &fDigiHitADC);
937  fTree->Branch("DigiHitCellID", &fDigiHitCellID);
938 
939  if(fGeo->HasMuonDetector()) {
940  fTree->Branch("DiginHits_MuID", &fDiginHits_MuID);
941  fTree->Branch("DigiHitX_MuID", &fDigiHitX_MuID);
942  fTree->Branch("DigiHitY_MuID", &fDigiHitY_MuID);
943  fTree->Branch("DigiHitZ_MuID", &fDigiHitZ_MuID);
944  fTree->Branch("DigiHitTime_MuID", &fDigiHitTime_MuID);
945  fTree->Branch("DigiHitADC_MuID", &fDigiHitADC_MuID);
946  fTree->Branch("DigiHitCellID_MuID", &fDigiHitCellID_MuID);
947  }
948  }
949 
950  // Write calorimetry hits
951  if (fWriteCaloHits) {
952  fTree->Branch("ReconHits", &fReconHits);
953  fTree->Branch("ReconHitIDNumber", &fReconHitIDNumber);
954  fTree->Branch("RecoHitX", &fRecoHitX);
955  fTree->Branch("RecoHitY", &fRecoHitY);
956  fTree->Branch("RecoHitZ", &fRecoHitZ);
957  fTree->Branch("RecoHitTime", &fRecoHitTime);
958  fTree->Branch("RecoHitEnergy", &fRecoHitEnergy);
959  fTree->Branch("RecoHitCellID", &fRecoHitCellID);
960  fTree->Branch("RecoHitLayer", &fRecoHitLayer);
961  fTree->Branch("RecoEnergySum", &fRecoEnergySum);
962 
963  if(fGeo->HasMuonDetector() && fWriteMuID) {
964  fTree->Branch("ReconHits_MuID", &fReconHits_MuID);
965  fTree->Branch("ReconHitIDNumber_MuID", &fReconHitIDNumber_MuID);
966  fTree->Branch("RecoHitX_MuID", &fRecoHitX_MuID);
967  fTree->Branch("RecoHitY_MuID", &fRecoHitY_MuID);
968  fTree->Branch("RecoHitZ_MuID", &fRecoHitZ_MuID);
969  fTree->Branch("RecoHitTime_MuID", &fRecoHitTime_MuID);
970  fTree->Branch("RecoHitEnergy_MuID", &fRecoHitEnergy_MuID);
971  fTree->Branch("RecoHitCellID_MuID", &fRecoHitCellID_MuID);
972  fTree->Branch("RecoHitLayer_MuID", &fRecoHitLayer_MuID);
973  fTree->Branch("RecoEnergySum_MuID", &fRecoEnergySum_MuID);
974  }
975  }
976 
977  // Write calorimetry clusters
978  if (fWriteCaloClusters) {
979  fTree->Branch("nCluster", &fnCluster);
980  fTree->Branch("ClusterIDNumber", &fClusterIDNumber);
981  fTree->Branch("ClusterNhits", &fClusterNhits);
982  fTree->Branch("ClusterEnergy", &fClusterEnergy);
983  fTree->Branch("ClusterTime", &fClusterTime);
984  fTree->Branch("ClusterTimeDiffFirstLast", &fClusterTimeDiffFirstLast);
985  fTree->Branch("ClusterX", &fClusterX);
986  fTree->Branch("ClusterY", &fClusterY);
987  fTree->Branch("ClusterZ", &fClusterZ);
988  fTree->Branch("ClusterTheta", &fClusterTheta);
989  fTree->Branch("ClusterPhi", &fClusterPhi);
990  fTree->Branch("ClusterPID", &fClusterPID);
991  // fTree->Branch("ClusterShape", &fClusterShape);
992  fTree->Branch("ClusterMainAxisX", &fClusterMainAxisX);
993  fTree->Branch("ClusterMainAxisY", &fClusterMainAxisY);
994  fTree->Branch("ClusterMainAxisZ", &fClusterMainAxisZ);
995  fTree->Branch("ClusterMCindex", &fClusterMCindex);
996  fTree->Branch("ClusterMCfrac", &fClusterMCfrac);
997 
998  fTree->Branch("ClusterAssn_RecoHitIDNumber", &fClusterAssn_RecoHitIDNumber);
999 
1000  if(fGeo->HasMuonDetector() && fWriteMuID) {
1001  fTree->Branch("nCluster_MuID", &fnCluster_MuID);
1002  fTree->Branch("ClusterIDNumber_MuID", &fClusterIDNumber_MuID);
1003  fTree->Branch("ClusterNhits_MuID", &fClusterNhits_MuID);
1004  fTree->Branch("ClusterEnergy_MuID", &fClusterEnergy_MuID);
1005  fTree->Branch("ClusterTime_MuID", &fClusterTime_MuID);
1006  fTree->Branch("ClusterTimeDiffFirstLast_MuID", &fClusterTimeDiffFirstLast_MuID);
1007  fTree->Branch("ClusterX_MuID", &fClusterX_MuID);
1008  fTree->Branch("ClusterY_MuID", &fClusterY_MuID);
1009  fTree->Branch("ClusterZ_MuID", &fClusterZ_MuID);
1010  fTree->Branch("ClusterTheta_MuID", &fClusterTheta_MuID);
1011  fTree->Branch("ClusterPhi_MuID", &fClusterPhi_MuID);
1012  fTree->Branch("ClusterPID_MuID", &fClusterPID_MuID);
1013  // fTree->Branch("ClusterShape", &fClusterShape);
1014  fTree->Branch("ClusterMainAxisX_MuID", &fClusterMainAxisX_MuID);
1015  fTree->Branch("ClusterMainAxisY_MuID", &fClusterMainAxisY_MuID);
1016  fTree->Branch("ClusterMainAxisZ_MuID", &fClusterMainAxisZ_MuID);
1017  fTree->Branch("ClusterMCindex_MuID", &fClusterMCindex_MuID);
1018  fTree->Branch("ClusterMCfrac_MuID", &fClusterMCfrac_MuID);
1019 
1020  fTree->Branch("ClusterMuIDAssn_MuIDHitIDNumber", &fClusterMuIDAssn_MuIDHitIDNumber);
1021  }
1022  }
1023 
1024  if (fWriteMatchedTracks) {
1025  if (!fWriteTracks || !fWriteCaloClusters) {
1026  throw cet::exception("anatree")
1027  << " fWriteMatchedTracks, but (!fWriteTracks || !fWriteCaloClusters)."
1028  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1029  } else {
1030  fTree->Branch("ECALAssn_ClusIDNumber", &fCALAssn_ClusIDNumber);
1031  fTree->Branch("ECALAssn_TrackIDNumber", &fCALAssn_TrackIDNumber);
1032  fTree->Branch("ECALAssn_TrackEnd", &fCALAssn_TrackEnd);
1033  }
1034  }
1035 
1036  std::string filename = "${DUNE_PARDATA_DIR}/MPD/dedxPID/dedxpidmatrices8kevcm.root";
1037  TFile infile(filename.c_str(), "READ");
1038 
1039  m_pidinterp.clear();
1040  char str[11];
1041  for (int q = 0; q < 501; ++q) {
1042  sprintf(str, "%d", q);
1043  std::string s = "pidmatrix";
1044  s.append(str);
1045  // read the 500 histograms one by one; each histogram is a
1046  // 6 by 6 matrix of probabilities for a given momentum value
1047  m_pidinterp.insert( std::make_pair(q, (TH2F*) infile.Get(s.c_str())->Clone("pidinterp")) );
1048  }
1049 
1050  return;
1051 } // End of :anatree::beginJob
1052 
1053 //==============================================================================
1054 //==============================================================================
1055 //==============================================================================
1056 void gar::anatree::endRun(art::Run const& run) {
1057  auto const& ID = run.id();
1058 
1059  auto summaryHandle = run.getHandle<sumdata::POTSummary>(fPOTtag);
1060  if (!summaryHandle) {
1061  MF_LOG_DEBUG("anatree") << " No sumdata::POTSummary branch for run " << ID
1062  <<" Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1063  return;
1064  }
1065 
1066  sumdata::POTSummary const& RunPOT = *summaryHandle;
1067  fTotalPOT = RunPOT.TotalPOT();
1068  fNSpills = RunPOT.TotalSpills();
1069 
1070  MF_LOG_INFO("anatree") << "POT for this file is " << fTotalPOT
1071  << " The number of spills is " << fNSpills;
1072 }
1073 
1074 //==============================================================================
1075 //==============================================================================
1076 //==============================================================================
1077 void gar::anatree::analyze(art::Event const & e) {
1078 
1079  ClearVectors();
1080 
1081  fRun = e.run();
1082  fSubRun = e.subRun();
1083  fEvent = e.id().event();
1084 
1085  // Need a non-constant backtracker instance, for now, in analyze not beginJob
1086  cheat::BackTrackerCore const* const_bt = gar::providerFrom<cheat::BackTracker>();
1087  BackTrack = const_cast<cheat::BackTrackerCore*>(const_bt);
1088 
1089  //Fill generator and MC Information
1090  if (fWriteMCinfo) FillGeneratorMonteCarloInfo(e);
1091 
1092  //Fill raw calorimeter information
1093  if (fWriteCaloDigits) FillRawInfo(e);
1094 
1095  //Fill reco information
1096  FillRecoInfo(e);
1097 
1098  //Fill high-leve reco information
1099  FillHighLevelRecoInfo(e);
1100 
1101  fTree->Fill();
1102  return;
1103 }
1104 
1105 
1106 
1107 //==============================================================================
1108 //==============================================================================
1109 //==============================================================================
1110 void gar::anatree::ClearVectors() {
1111 
1112  // clear out all our vectors
1113  if (fWriteMCinfo) {
1114  fNeutrinoType.clear();
1115  fCCNC.clear();
1116  fMode.clear();
1117  fInteractionType.clear();
1118  fQ2.clear();
1119  fW.clear();
1120  fX.clear();
1121  fY.clear();
1122  fTheta.clear();
1123  if (fWriteCohInfo) fT.clear();
1124  fMCVertexX.clear();
1125  fMCVertexY.clear();
1126  fMCVertexZ.clear();
1127  fMCnuPx.clear();
1128  fMCnuPy.clear();
1129  fMCnuPz.clear();
1130 
1131  fGint.clear();
1132  fTgtPDG.clear();
1133  fWeight.clear();
1134  fgT.clear();
1135 
1136  fnGPart.clear();
1137  fGPartIntIdx.clear();
1138  fGPartIdx.clear();
1139  fGPartPdg.clear();
1140  fGPartStatus.clear();
1141  fGPartName.clear();
1142  fGPartFirstMom.clear();
1143  fGPartLastMom.clear();
1144  fGPartFirstDaugh.clear();
1145  fGPartLastDaugh.clear();
1146  fGPartPx.clear();
1147  fGPartPy.clear();
1148  fGPartPz.clear();
1149  fGPartE.clear();
1150  fGPartMass.clear();
1151 
1152  fMCTrkID.clear();
1153  fMCPDG.clear();
1154  fMCMotherIndex.clear();
1155  fMCMotherTrkID.clear();
1156  fMCPDGMother.clear();
1157  fMCPStartX.clear();
1158  fMCPStartY.clear();
1159  fMCPStartZ.clear();
1160  fMCPTime.clear();
1161  fMCPStartPX.clear();
1162  fMCPStartPY.clear();
1163  fMCPStartPZ.clear();
1164  fMCPEndX.clear();
1165  fMCPEndY.clear();
1166  fMCPEndZ.clear();
1167  fMCPEndPX.clear();
1168  fMCPEndPY.clear();
1169  fMCPEndPZ.clear();
1170  fMCPProc.clear();
1171  fMCPEndProc.clear();
1172  fMCPVertIndex.clear();
1173  }
1174 
1175  if (fWriteMCPTrajectory) {
1176  fTrajMCPX.clear();
1177  fTrajMCPY.clear();
1178  fTrajMCPZ.clear();
1179  fTrajMCPT.clear();
1180  if (fWriteMCPTrajMomenta) {
1181  fTrajMCPPX.clear();
1182  fTrajMCPPY.clear();
1183  fTrajMCPPZ.clear();
1184  }
1185  fTrajMCPE.clear();
1186  fTrajMCPIndex.clear();
1187  fTrajMCPTrackID.clear();
1188  }
1189 
1190  if (fWriteMCCaloInfo) {
1191  fSimnHits = 0;
1192  fSimHitX.clear();
1193  fSimHitY.clear();
1194  fSimHitZ.clear();
1195  fSimHitTime.clear();
1196  fSimHitEnergy.clear();
1197  fSimHitTrackID.clear();
1198  fSimHitCellID.clear();
1199  fSimEnergySum = 0.;
1200 
1201  if(fGeo->HasMuonDetector() && fWriteMuID) {
1202  fSimnHits_MuID = 0;
1203  fSimHitX_MuID.clear();
1204  fSimHitY_MuID.clear();
1205  fSimHitZ_MuID.clear();
1206  fSimHitTime_MuID.clear();
1207  fSimHitEnergy_MuID.clear();
1208  fSimHitTrackID_MuID.clear();
1209  fSimHitCellID_MuID.clear();
1210  fSimEnergySum_MuID = 0.;
1211  }
1212  }
1213 
1214  if (fWriteHits) {
1215  fHitX.clear();
1216  fHitY.clear();
1217  fHitZ.clear();
1218  fHitSig.clear();
1219  fHitRMS.clear();
1220  fHitChan.clear();
1221  }
1222 
1223  if (fWriteTPCClusters) {
1224  fTPCClusterX.clear();
1225  fTPCClusterY.clear();
1226  fTPCClusterZ.clear();
1227  fTPCClusterSig.clear();
1228  fTPCClusterRMS.clear();
1229  fTPCClusterTrkIDNumber.clear();
1230  fTPCClusterCovXX.clear();
1231  fTPCClusterCovXY.clear();
1232  fTPCClusterCovXZ.clear();
1233  fTPCClusterCovYY.clear();
1234  fTPCClusterCovYZ.clear();
1235  fTPCClusterCovZZ.clear();
1236  fTPCClusterMCindex.clear();
1237  fTPCClusterMCfrac.clear();
1238  }
1239 
1240  if (fWriteTracks) {
1241  fTrackIDNumber.clear();
1242  fTrackStartX.clear();
1243  fTrackStartY.clear();
1244  fTrackStartZ.clear();
1245  fTrackStartPX.clear();
1246  fTrackStartPY.clear();
1247  fTrackStartPZ.clear();
1248  fTrackStartQ.clear();
1249  fTrackEndX.clear();
1250  fTrackEndY.clear();
1251  fTrackEndZ.clear();
1252  fTrackEndPX.clear();
1253  fTrackEndPY.clear();
1254  fTrackEndPZ.clear();
1255  fTrackEndQ.clear();
1256  fTrackLenF.clear();
1257  fTrackLenB.clear();
1258  fTrackChi2F.clear();
1259  fTrackChi2B.clear();
1260  fNTPCClustersOnTrack.clear();
1261  fTrackAvgIonF.clear();
1262  fTrackAvgIonB.clear();
1263  fTrackPIDF.clear();
1264  fTrackPIDProbF.clear();
1265  fTrackPIDB.clear();
1266  fTrackPIDProbB.clear();
1267  fTrackMCindex.clear();
1268  fTrackMCfrac.clear();
1269 
1270  if(fWriteTrackTrajectories) {
1271  fTrackTrajectoryFWDX.clear();
1272  fTrackTrajectoryFWDY.clear();
1273  fTrackTrajectoryFWDZ.clear();
1274  fTrackTrajectoryFWDID.clear();
1275 
1276  fTrackTrajectoryBWDX.clear();
1277  fTrackTrajectoryBWDY.clear();
1278  fTrackTrajectoryBWDZ.clear();
1279  fTrackTrajectoryBWDID.clear();
1280  }
1281  }
1282 
1283  if (fWriteVertices) {
1284  fVertexIDNumber.clear();
1285  fVertexX.clear();
1286  fVertexY.clear();
1287  fVertexZ.clear();
1288  fVertexT.clear();
1289  fVertexN.clear();
1290  fVertexQ.clear();
1291 
1292  fVTAssn_VertIDNumber.clear();
1293  fVTAssn_TrackIDNumber.clear();
1294  fVTAssn_TrackEnd.clear();
1295  }
1296 
1297  if (fWriteVees) {
1298  fVeeIDNumber.clear();
1299  fVeeX.clear();
1300  fVeeY.clear();
1301  fVeeZ.clear();
1302  fVeeT.clear();
1303  fVeePXKpipi.clear();
1304  fVeePYKpipi.clear();
1305  fVeePZKpipi.clear();
1306  fVeeEKpipi.clear();
1307  fVeeMKpipi.clear();
1308  fVeePXLppi.clear();
1309  fVeePYLppi.clear();
1310  fVeePZLppi.clear();
1311  fVeeELppi.clear();
1312  fVeeMLppi.clear();
1313  fVeePXLpip.clear();
1314  fVeePYLpip.clear();
1315  fVeePZLpip.clear();
1316  fVeeELpip.clear();
1317  fVeeMLpip.clear();
1318  fVeeTAssn_VeeIDNumber.clear();
1319  fVeeTAssn_TrackIDNumber.clear();
1320  fVeeTAssn_TrackEnd.clear();
1321  }
1322 
1323  if (fWriteCaloDigits) {
1324  fDiginHits = 0;
1325  fDigiHitX.clear();
1326  fDigiHitY.clear();
1327  fDigiHitZ.clear();
1328  fDigiHitTime.clear();
1329  fDigiHitADC.clear();
1330  fDigiHitCellID.clear();
1331 
1332  if(fGeo->HasMuonDetector() && fWriteMuID) {
1333  fDiginHits_MuID = 0;
1334  fDigiHitX_MuID.clear();
1335  fDigiHitY_MuID.clear();
1336  fDigiHitZ_MuID.clear();
1337  fDigiHitTime_MuID.clear();
1338  fDigiHitADC_MuID.clear();
1339  fDigiHitCellID_MuID.clear();
1340  }
1341  }
1342 
1343  if (fWriteCaloHits) {
1344  fReconHits = 0;
1345  fReconHitIDNumber.clear();
1346  fRecoHitX.clear();
1347  fRecoHitY.clear();
1348  fRecoHitZ.clear();
1349  fRecoHitTime.clear();
1350  fRecoHitEnergy.clear();
1351  fRecoHitCellID.clear();
1352  fRecoHitLayer.clear();
1353  fRecoEnergySum = 0.;
1354 
1355  if(fGeo->HasMuonDetector() && fWriteMuID) {
1356  fReconHits_MuID = 0;
1357  fReconHitIDNumber_MuID.clear();
1358  fRecoHitX_MuID.clear();
1359  fRecoHitY_MuID.clear();
1360  fRecoHitZ_MuID.clear();
1361  fRecoHitTime_MuID.clear();
1362  fRecoHitEnergy_MuID.clear();
1363  fRecoHitCellID_MuID.clear();
1364  fRecoHitLayer_MuID.clear();
1365  fRecoEnergySum_MuID = 0.;
1366  }
1367  }
1368 
1369  if (fWriteCaloClusters) {
1370  fnCluster = 0;
1371  fClusterIDNumber.clear();
1372  fClusterNhits.clear();
1373  fClusterEnergy.clear();
1374  fClusterTime.clear();
1375  fClusterTimeDiffFirstLast.clear();
1376  fClusterX.clear();
1377  fClusterY.clear();
1378  fClusterZ.clear();
1379  fClusterTheta.clear();
1380  fClusterPhi.clear();
1381  fClusterPID.clear();
1382  // fClusterShape.clear();
1383  fClusterMainAxisX.clear();
1384  fClusterMainAxisY.clear();
1385  fClusterMainAxisZ.clear();
1386  fClusterMCindex.clear();
1387  fClusterMCfrac.clear();
1388 
1389  fClusterAssn_RecoHitIDNumber.clear();
1390 
1391  if(fGeo->HasMuonDetector() && fWriteMuID) {
1392  fnCluster_MuID = 0;
1393  fClusterIDNumber_MuID.clear();
1394  fClusterNhits_MuID.clear();
1395  fClusterEnergy_MuID.clear();
1396  fClusterTime_MuID.clear();
1397  fClusterTimeDiffFirstLast_MuID.clear();
1398  fClusterX_MuID.clear();
1399  fClusterY_MuID.clear();
1400  fClusterZ_MuID.clear();
1401  fClusterTheta_MuID.clear();
1402  fClusterPhi_MuID.clear();
1403  fClusterPID_MuID.clear();
1404  // fClusterShape.clear();
1405  fClusterMainAxisX_MuID.clear();
1406  fClusterMainAxisY_MuID.clear();
1407  fClusterMainAxisZ_MuID.clear();
1408  fClusterMCindex_MuID.clear();
1409  fClusterMCfrac_MuID.clear();
1410 
1411  fClusterMuIDAssn_MuIDHitIDNumber.clear();
1412  }
1413  }
1414 
1415  if (fWriteMatchedTracks) {
1416  fCALAssn_ClusIDNumber.clear();
1417  fCALAssn_TrackIDNumber.clear();
1418  fCALAssn_TrackEnd.clear();
1419  }
1420 
1421  return;
1422 } // end :anatree::ClearVectors
1423 
1424 
1425 
1426 //==============================================================================
1427 //==============================================================================
1428 //==============================================================================
1429 void gar::anatree::FillGeneratorMonteCarloInfo(art::Event const & e) {
1430 
1431  // ============= Get art handles ==========================================
1432  // Get handles for MCinfo, also good for MCPTrajectory
1433  std::vector< art::Handle< std::vector<simb::MCTruth> > > mcthandlelist;
1434  std::vector< art::Handle< std::vector<simb::GTruth> > > gthandlelist;
1435 
1436  if (fGeneratorLabels.size()<1) {
1437  mcthandlelist = e.getMany<std::vector<simb::MCTruth> >(); // get them all (even if there are none)
1438  } else {
1439  mcthandlelist.resize(fGeneratorLabels.size());
1440  for (size_t i=0; i< fGeneratorLabels.size(); ++i) {
1441  // complain if we wanted a specific one but didn't find it
1442  mcthandlelist.at(i) = e.getHandle<std::vector<simb::MCTruth> >(fGeneratorLabels.at(i));
1443  if (!mcthandlelist.at(i)) {
1444  throw cet::exception("anatree") << " No simb::MCTruth branch."
1445  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1446  }
1447  }
1448  }
1449 
1450  if (fGENIEGeneratorLabels.size()<1) {
1451  gthandlelist = e.getMany< std::vector<simb::GTruth> >(); // get them all (even if there are none)
1452  } else {
1453  gthandlelist.resize(fGENIEGeneratorLabels.size());
1454  for (size_t i=0; i< fGENIEGeneratorLabels.size(); ++i) {
1455  // complain if we wanted a specific one but didn't find it
1456  gthandlelist.at(i) = e.getHandle<std::vector<simb::GTruth> >(fGENIEGeneratorLabels.at(i));
1457  if (!gthandlelist.at(i)) {
1458  throw cet::exception("anatree") << " No simb::GTruth branch."
1459  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1460  }
1461  }
1462  }
1463 
1464  // save MCTruth info
1465  for (size_t imchl = 0; imchl < mcthandlelist.size(); ++imchl) {
1466  for ( auto const& mct : (*mcthandlelist.at(imchl)) ) {
1467  if (mct.NeutrinoSet()) {
1468  simb::MCNeutrino nuw = mct.GetNeutrino();
1469  fNeutrinoType.push_back(nuw.Nu().PdgCode());
1470  fCCNC.push_back(nuw.CCNC());
1471  fMode.push_back(nuw.Mode());
1472  fInteractionType.push_back(nuw.InteractionType());
1473  fQ2.push_back(nuw.QSqr());
1474  fW.push_back(nuw.W());
1475  fX.push_back(nuw.X());
1476  fY.push_back(nuw.Y());
1477  fTheta.push_back(nuw.Theta());
1478  if (fWriteCohInfo) {
1479  double getT = computeT(mct);
1480  fT.push_back( static_cast<Float_t>(getT) );
1481  }
1482  fMCVertexX.push_back(nuw.Nu().EndX());
1483  fMCVertexY.push_back(nuw.Nu().EndY());
1484  fMCVertexZ.push_back(nuw.Nu().EndZ());
1485  fMCnuPx.push_back(nuw.Nu().Px());
1486  fMCnuPy.push_back(nuw.Nu().Py());
1487  fMCnuPz.push_back(nuw.Nu().Pz());
1488  } // end MC info from MCTruth
1489  }
1490  }
1491 
1492  // save GTruth info
1493  for (size_t igthl = 0; igthl < gthandlelist.size(); ++igthl) {
1494  for ( auto const& gt : (*gthandlelist.at(igthl)) ) {
1495  fGint.push_back(gt.fGint);
1496  fTgtPDG.push_back(gt.ftgtPDG);
1497  fWeight.push_back(gt.fweight);
1498  fgT.push_back(gt.fgT);
1499  }
1500  }
1501 
1502  // Save the particle list from the GENIE event record
1503  auto gparthandlelist = e.getMany<std::vector<sdp::GenieParticle>>();
1504 
1505  for (size_t igphl = 0; igphl < gparthandlelist.size(); ++igphl) {
1506  unsigned int nGPart = 0;
1507  for ( auto const& gpart : (*gparthandlelist.at(igphl)) ) {
1508  fGPartIntIdx.push_back(gpart.InteractionIndex());
1509  fGPartIdx.push_back(gpart.Index());
1510  fGPartPdg.push_back(gpart.Pdg());
1511  fGPartStatus.push_back(gpart.Status());
1512  fGPartName.push_back(gpart.Name());
1513  fGPartFirstMom.push_back(gpart.FirstMother());
1514  fGPartLastMom.push_back(gpart.LastMother());
1515  fGPartFirstDaugh.push_back(gpart.FirstDaughter());
1516  fGPartLastDaugh.push_back(gpart.LastDaughter());
1517  fGPartPx.push_back(gpart.Px());
1518  fGPartPy.push_back(gpart.Py());
1519  fGPartPz.push_back(gpart.Pz());
1520  fGPartE.push_back(gpart.E());
1521  fGPartMass.push_back(gpart.Mass());
1522  nGPart++;
1523  }
1524  fnGPart.push_back(nGPart);
1525  }
1526 
1527  auto MCPHandle = e.getHandle<std::vector<simb::MCParticle> >(fGeantLabel);
1528  if (!MCPHandle) {
1529  throw cet::exception("anatree") << " No simb::MCParticle branch."
1530  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1531  }
1532 
1533  // Save MCParticle info (post-GEANT; for pre-GEANT, maybe try MCTruth?)
1534  // Helps to have a map from the TrackId from generator to index number
1535  // into *MCPHandle, which is a vector<MCParticle>. Not all TrackId value
1536  // appear in MCPHandle, although I think the missing ones are all from
1537  // the GEANT rather than the GENIE stage of GENIEGen. (*MCPHandle).size()
1538  // is the same as eg fMCPStartX.size() by the time this entry is written.
1539  // TracIdToIndex isa class variable for accessibility by per-Track and
1540  // per-Cluster code
1541  Int_t index = 0;
1542  for ( auto const& mcp : (*MCPHandle) ) {
1543  int TrackId = mcp.TrackId();
1544  TrackIdToIndex[TrackId] = index++;
1545  }
1546 
1547  for ( auto const& mcp : (*MCPHandle) ) {
1548  fMCTrkID.push_back(mcp.TrackId());
1549  fMCPDG.push_back(mcp.PdgCode());
1550 
1551  // If mcp.Mother() == 0, particle is from initial vertex;
1552  // Set MCMotherIndex to -1 in that case.
1553  TrkId momTrkId = mcp.Mother();
1554  Int_t momIndex = -1;
1555  int momPDG = 0;
1556  if (momTrkId>0) {
1557  //Check if it exists!
1558  if(TrackIdToIndex.find(momTrkId) != TrackIdToIndex.end()){
1559  momIndex = TrackIdToIndex[momTrkId];
1560  momPDG = (*MCPHandle).at(momIndex).PdgCode();
1561  } else {
1562  MF_LOG_DEBUG("Anatree_module")
1563  << " mcp trkid " << mcp.TrackId()
1564  << " pdg code " << mcp.PdgCode()
1565  << " could not find mother trk id " << momTrkId
1566  << " in the TrackIdToIndex map"
1567  << " creating process is [ " << mcp.Process() << " ]";
1568  }
1569  }
1570 
1571  fMCMotherIndex.push_back(momIndex);
1572  fMCMotherTrkID.push_back(mcp.Mother());//directly trackid not index
1573  fMCPDGMother.push_back(momPDG);
1574 
1575  const TLorentzVector& position = mcp.Position(0);
1576  const TLorentzVector& momentum = mcp.Momentum(0);
1577  fMCPStartX.push_back(position.X());
1578  fMCPStartY.push_back(position.Y());
1579  fMCPStartZ.push_back(position.Z());
1580  fMCPTime.push_back(mcp.T());
1581  fMCPStartPX.push_back(momentum.Px());
1582  fMCPStartPY.push_back(momentum.Py());
1583  fMCPStartPZ.push_back(momentum.Pz());
1584 
1585  const TLorentzVector& positionEnd = mcp.EndPosition();
1586  const TLorentzVector& momentumEnd = mcp.EndMomentum();
1587  fMCPEndX.push_back(positionEnd.X());
1588  fMCPEndY.push_back(positionEnd.Y());
1589  fMCPEndZ.push_back(positionEnd.Z());
1590  fMCPEndPX.push_back(momentumEnd.Px());
1591  fMCPEndPY.push_back(momentumEnd.Py());
1592  fMCPEndPZ.push_back(momentumEnd.Pz());
1593  fMCPProc.push_back(mcp.Process());
1594  fMCPEndProc.push_back(mcp.EndProcess());
1595  }
1596 
1597  // Get vertex for each MCParticle. In principle, the mct.GetParticle()
1598  // method should produce all of them; filter out the ones with a non-stable
1599  // status code and the GENIE pseudo-particles (Pid>= 2000000000) and you
1600  // should have the contents of (*MCPHandle). But you don't. On the 1st
1601  // event I (Leo Bellantoni) tested, with 40 vertices and 351 primary
1602  // particles, particle 143, a proton, evidently decayed or something and is
1603  // not in (*MCPHandle). So we use the following primitive earth technology.
1604  size_t nMCParticles = (*MCPHandle).size();
1605  size_t iMCParticle = 0;
1606  fMCPVertIndex.resize(nMCParticles);
1607  for (; iMCParticle<nMCParticles; ++iMCParticle) {
1608  // Assign noprimary to start with
1609  fMCPVertIndex[iMCParticle] = -1;
1610  // Do the primaries first
1611  if (fMCMotherIndex[iMCParticle]!=-1) break;
1612  Float_t trackX = fMCPStartX[iMCParticle];
1613  Float_t trackY = fMCPStartY[iMCParticle];
1614  Float_t trackZ = fMCPStartZ[iMCParticle];
1615  int vertexIndex = 0;
1616  for (size_t imchl = 0; imchl < mcthandlelist.size(); ++imchl) {
1617  for ( auto const& mct : (*mcthandlelist.at(imchl)) ) {
1618  if (mct.NeutrinoSet()) {
1619  simb::MCNeutrino nuw = mct.GetNeutrino();
1620  Float_t vertX = nuw.Nu().EndX();
1621  Float_t vertY = nuw.Nu().EndY();
1622  Float_t vertZ = nuw.Nu().EndZ();
1623  Float_t dist = std::hypot(trackX-vertX,trackY-vertY,trackZ-vertZ);
1624  if ( dist <= fMatchMCPtoVertDist ) {
1625  fMCPVertIndex[iMCParticle] = vertexIndex;
1626  goto foundMCvert; // break out of all inner loops
1627  }
1628  }
1629  ++vertexIndex;
1630  }
1631  }
1632  foundMCvert:
1633  ;
1634  }
1635 
1636  // Now the secondaries. As they are after the primaries, do not re-init iMCParticle
1637  for (; iMCParticle<nMCParticles; ++iMCParticle) {
1638  int momIndex = fMCMotherIndex[iMCParticle];
1639  int lastMCParticle = iMCParticle;
1640  while (momIndex != -1) {
1641  lastMCParticle = momIndex;
1642  momIndex = fMCMotherIndex[momIndex];
1643  }
1644  fMCPVertIndex[iMCParticle] = fMCPVertIndex[lastMCParticle];
1645  }
1646 
1647  if (fWriteMCPTrajectory) {
1648  // It's in the MCParticle table
1649  Int_t mcpIndex = 0;
1650  for ( auto const& mcp : (*MCPHandle) ) {
1651  const TDatabasePDG* databasePDG = TDatabasePDG::Instance();
1652  const TParticlePDG* definition = databasePDG->GetParticle( mcp.PdgCode() );
1653  //No charge don't store the trajectory
1654  // this test fails for alpha particles because they aren't in databasePDG
1655  // so skip it in this case.
1656  if (mcp.PdgCode() != 1000020040)
1657  {
1658  if (definition==nullptr || definition->Charge() == 0) continue;
1659  }
1660  //TrackID of the mcp to keep track to which mcp this trajectory is
1661  int trackId = mcp.TrackId();
1662  for(uint iTraj=0; iTraj < mcp.Trajectory().size(); iTraj++) {
1663  float xTraj = mcp.Trajectory().X(iTraj);
1664  float yTraj = mcp.Trajectory().Y(iTraj);
1665  float zTraj = mcp.Trajectory().Z(iTraj);
1666  float rTraj = std::hypot( yTraj - ItsInTulsa[1], zTraj - ItsInTulsa[2]);
1667 
1668  if (abs(xTraj - ItsInTulsa[0]) > xTPC) continue;
1669  if (rTraj > rTPC) continue;
1670 
1671  fTrajMCPX.push_back(xTraj);
1672  fTrajMCPY.push_back(yTraj);
1673  fTrajMCPZ.push_back(zTraj);
1674  fTrajMCPT.push_back(mcp.Trajectory().T(iTraj));
1675  if (fWriteMCPTrajMomenta) {
1676  fTrajMCPPX.push_back(mcp.Trajectory().Px(iTraj));
1677  fTrajMCPPY.push_back(mcp.Trajectory().Py(iTraj));
1678  fTrajMCPPZ.push_back(mcp.Trajectory().Pz(iTraj));
1679  }
1680  fTrajMCPE.push_back(mcp.Trajectory().E(iTraj));
1681  fTrajMCPIndex.push_back(mcpIndex);
1682  fTrajMCPTrackID.push_back(trackId);
1683  }
1684  mcpIndex++;
1685  }
1686  }
1687 
1688  // Get handles for MCCaloInfo
1689 
1690  if (fWriteMCCaloInfo) {
1691  art::InputTag ecalgeanttag(fGeantLabel, fInstanceLabelCalo);
1692  auto SimHitHandle = e.getHandle<std::vector<gar::sdp::CaloDeposit> >(ecalgeanttag);
1693  if (!SimHitHandle) {
1694  throw cet::exception("anatree") << " No gar::sdp::CaloDeposit branch."
1695  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1696  }
1697 
1698  // Save simulation ecal hit info
1699  for ( auto const& SimHit : (*SimHitHandle) ) {
1700  fSimnHits++;
1701  fSimHitX.push_back(SimHit.X());
1702  fSimHitY.push_back(SimHit.Y());
1703  fSimHitZ.push_back(SimHit.Z());
1704  fSimHitTime.push_back(SimHit.Time());
1705  fSimHitEnergy.push_back(SimHit.Energy());
1706  fSimHitTrackID.push_back(SimHit.TrackID());
1707  fSimHitCellID.push_back(SimHit.CellID());
1708  fSimEnergySum += SimHit.Energy();
1709  }
1710 
1711  if(fWriteMuID) {
1712  art::InputTag muidgeanttag(fGeantLabel, fInstanceLabelMuID);
1713  art::Handle< std::vector<gar::sdp::CaloDeposit> > MuIDSimHitHandle;
1714  if (fGeo->HasMuonDetector())
1715  {
1716  MuIDSimHitHandle = e.getHandle< std::vector<gar::sdp::CaloDeposit> >(muidgeanttag);
1717  if (!MuIDSimHitHandle) {
1718  throw cet::exception("anatree") << " No gar::sdp::CaloDeposit branch."
1719  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1720  }
1721  }
1722 
1723  // Save simulation muon system hit info
1724  for ( auto const& SimHit : (*MuIDSimHitHandle) ) {
1725  fSimnHits_MuID++;
1726  fSimHitX_MuID.push_back(SimHit.X());
1727  fSimHitY_MuID.push_back(SimHit.Y());
1728  fSimHitZ_MuID.push_back(SimHit.Z());
1729  fSimHitTime_MuID.push_back(SimHit.Time());
1730  fSimHitEnergy_MuID.push_back(SimHit.Energy());
1731  fSimHitTrackID_MuID.push_back(SimHit.TrackID());
1732  fSimHitCellID_MuID.push_back(SimHit.CellID());
1733  fSimEnergySum_MuID += SimHit.Energy();
1734  }
1735  }
1736  }
1737 }
1738 
1739 
1740 
1741 //==============================================================================
1742 //==============================================================================
1743 //==============================================================================
1744 void gar::anatree::FillRawInfo(art::Event const & e) {
1745 
1746  // save ecal raw digits info
1747  art::InputTag ecalrawtag(fRawCaloHitLabel, fInstanceLabelCalo);
1748  auto RawHitHandle = e.getHandle<std::vector<gar::raw::CaloRawDigit> >(ecalrawtag);
1749  if (!RawHitHandle) {
1750  throw cet::exception("anatree") << " No :raw::CaloRawDigit branch."
1751  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1752  }
1753  for ( auto const& DigiHit : (*RawHitHandle) ) {
1754  fDiginHits++;
1755  fDigiHitX.push_back(DigiHit.X());
1756  fDigiHitY.push_back(DigiHit.Y());
1757  fDigiHitZ.push_back(DigiHit.Z());
1758  fDigiHitTime.push_back( (DigiHit.Time().first + DigiHit.Time().second) / 2.0 );
1759  fDigiHitADC.push_back(DigiHit.ADC().first);
1760  fDigiHitCellID.push_back(DigiHit.CellID());
1761  }
1762 
1763  // save muon system raw digits info
1764  if(fWriteMuID) {
1765  art::InputTag muidrawtag(fRawMuIDHitLabel, fInstanceLabelMuID);
1766  art::Handle<std::vector<gar::raw::CaloRawDigit> > MuIDRawHitHandle;
1767  if (fGeo->HasMuonDetector()) {
1768  MuIDRawHitHandle = e.getHandle<std::vector<gar::raw::CaloRawDigit> >(muidrawtag);
1769  if (!MuIDRawHitHandle) {
1770  throw cet::exception("anatree") << " No :raw::CaloRawDigit branch."
1771  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1772  }
1773  }
1774 
1775  for ( auto const& DigiHit : (*MuIDRawHitHandle) ) {
1776  fDiginHits_MuID++;
1777  fDigiHitX_MuID.push_back(DigiHit.X());
1778  fDigiHitY_MuID.push_back(DigiHit.Y());
1779  fDigiHitZ_MuID.push_back(DigiHit.Z());
1780  fDigiHitTime_MuID.push_back( (DigiHit.Time().first + DigiHit.Time().second) / 2.0 );
1781  fDigiHitADC_MuID.push_back(DigiHit.ADC().first);
1782  fDigiHitCellID_MuID.push_back(DigiHit.CellID());
1783  }
1784  }
1785 }
1786 
1787 
1788 
1789 //==============================================================================
1790 //==============================================================================
1791 //==============================================================================
1792 void gar::anatree::FillRecoInfo(art::Event const & e) {
1793 
1794  // Get handle for TPC hit data
1795  if (fWriteHits) {
1796  auto HitHandle = e.getHandle<std::vector<rec::Hit> >(fHitLabel);
1797  if (!HitHandle) {
1798  throw cet::exception("anatree") << " No rec::Hit branch."
1799  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1800  }
1801 
1802  // save hits in the TPC
1803  for ( auto const& Hit : (*HitHandle) ) {
1804  fHitX.push_back(Hit.Position()[0]);
1805  fHitY.push_back(Hit.Position()[1]);
1806  fHitZ.push_back(Hit.Position()[2]);
1807  fHitSig.push_back(Hit.Signal());
1808  fHitRMS.push_back(Hit.RMS());
1809  fHitChan.push_back(Hit.Channel());
1810  }
1811  }
1812 
1813  if (fWriteCaloHits) {
1814  art::InputTag ecalrecotag(fCaloHitLabel, fInstanceLabelCalo);
1815  auto RecoHitHandle = e.getHandle<std::vector<rec::CaloHit> >(ecalrecotag);
1816  if (!RecoHitHandle) {
1817  throw cet::exception("anatree") << " No rec::CaloHit branch."
1818  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1819  }
1820  // save reco'd Calorimetry hits
1821  for ( auto const& Hit : (*RecoHitHandle) ) {
1822  fReconHits++;
1823  fReconHitIDNumber.push_back(Hit.getIDNumber());
1824  fRecoHitX.push_back(Hit.Position()[0]);
1825  fRecoHitY.push_back(Hit.Position()[1]);
1826  fRecoHitZ.push_back(Hit.Position()[2]);
1827  fRecoHitTime.push_back(Hit.Time().first);
1828  fRecoHitEnergy.push_back(Hit.Energy());
1829  fRecoHitCellID.push_back(Hit.CellID());
1830  fRecoHitLayer.push_back(fFieldDecoder_ECAL->get(Hit.CellID(), "layer"));
1831  fRecoEnergySum += Hit.Energy();
1832  }
1833 
1834  if(fWriteMuID) {
1835  art::InputTag muirecotag(fMuIDHitLabel, fInstanceLabelMuID);
1836  art::Handle<std::vector<rec::CaloHit> > MuIDRecoHitHandle;
1837  if (fGeo->HasMuonDetector()) {
1838  MuIDRecoHitHandle = e.getHandle<std::vector<rec::CaloHit> >(muirecotag);
1839  if (!MuIDRecoHitHandle) {
1840  throw cet::exception("anatree") << " No rec::CaloHit branch."
1841  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1842  }
1843  }
1844  for ( auto const& Hit : (*MuIDRecoHitHandle) ) {
1845  fReconHits_MuID++;
1846  fReconHitIDNumber_MuID.push_back(Hit.getIDNumber());
1847  fRecoHitX_MuID.push_back(Hit.Position()[0]);
1848  fRecoHitY_MuID.push_back(Hit.Position()[1]);
1849  fRecoHitZ_MuID.push_back(Hit.Position()[2]);
1850  fRecoHitTime_MuID.push_back(Hit.Time().first);
1851  fRecoHitEnergy_MuID.push_back(Hit.Energy());
1852  fRecoHitCellID_MuID.push_back(Hit.CellID());
1853  fRecoHitLayer_MuID.push_back(fFieldDecoder_MuID->get(Hit.CellID(), "layer"));
1854  fRecoEnergySum_MuID += Hit.Energy();
1855  }
1856  }
1857  }
1858 }
1859 
1860 //==============================================================================
1861 //==============================================================================
1862 //==============================================================================
1863 void gar::anatree::FillHighLevelRecoInfo(art::Event const & e) {
1864 
1865  // Get handle for TPCClusters
1866  art::Handle< std::vector<rec::TPCCluster> > TPCClusterHandle;
1867  art::FindManyP<rec::Hit>* findManyHits = NULL;
1868  if (fWriteTPCClusters) {
1869  TPCClusterHandle = e.getHandle< std::vector<rec::TPCCluster> >(fTPCClusterLabel);
1870  if (!TPCClusterHandle) {
1871  throw cet::exception("anatree") << " No rec::TPCCluster branch."
1872  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1873  }
1874  findManyHits = new art::FindManyP<rec::Hit>(TPCClusterHandle,e,fTPCClusterLabel);
1875  }
1876 
1877  // Get handles for Tracks and their ionizations; also Assn's to TPCClusters, TrackIoniz
1878  // null handles if we switch off reading in the data products
1879  art::Handle< std::vector<rec::Track> > TrackHandle;
1880  art::Handle< std::vector<rec::TrackIoniz> > TrackIonHandle;
1881  art::Handle< std::vector<rec::TrackTrajectory> > TrackTrajHandle;
1882  art::FindManyP<rec::TPCCluster>* findManyTPCClusters = NULL;
1883  art::FindOneP<rec::TrackIoniz>* findIonization = NULL;
1884  if(fWriteTracks) {
1885  TrackHandle = e.getHandle< std::vector<rec::Track> >(fTrackLabel);
1886  if (!TrackHandle) {
1887  throw cet::exception("anatree") << " No rec::Track branch."
1888  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1889  }
1890  TrackIonHandle = e.getHandle< std::vector<rec::TrackIoniz> >(fTrackLabel);
1891  if (!TrackIonHandle) {
1892  throw cet::exception("anatree") << " No rec::TrackIoniz branch."
1893  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1894  }
1895 
1896  findManyTPCClusters = new art::FindManyP<rec::TPCCluster>(TrackHandle,e,fTrackLabel);
1897  findIonization = new art::FindOneP<rec::TrackIoniz>(TrackHandle,e,fTrackLabel);
1898 
1899  if(fWriteTrackTrajectories) {
1900  TrackTrajHandle = e.getHandle< std::vector<rec::TrackTrajectory> >(fTrackTrajectoryLabel);
1901  if (!TrackTrajHandle) {
1902  throw cet::exception("anatree") << " No rec::TrackTrajectory branch."
1903  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1904  }
1905  }
1906  }
1907 
1908  // Get handle for Vertices; also Assn's to Tracks
1909  art::Handle< std::vector<rec::Vertex> > VertexHandle;
1910  art::FindManyP<rec::Track, rec::TrackEnd>* findManyTrackEnd = NULL;
1911  if (fWriteVertices) {
1912  VertexHandle = e.getHandle< std::vector<rec::Vertex> >(fVertexLabel);
1913  if (!VertexHandle) {
1914  throw cet::exception("anatree") << " No rec::Vertex branch."
1915  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1916  }
1917  findManyTrackEnd = new art::FindManyP<rec::Track, rec::TrackEnd>(VertexHandle,e,fVertexLabel);
1918  }
1919 
1920  // Get handle for Vees; also Assn's to Tracks
1921  art::Handle< std::vector<rec::Vee> > VeeHandle;
1922  art::FindManyP<rec::Track, rec::TrackEnd>* findManyVeeTrackEnd = NULL;
1923  if (fWriteVees) {
1924  VeeHandle = e.getHandle< std::vector<rec::Vee> >(fVeeLabel);
1925  if (!VeeHandle) {
1926  throw cet::exception("anatree") << " No rec::Vee branch."
1927  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1928  }
1929  findManyVeeTrackEnd = new art::FindManyP<rec::Track, rec::TrackEnd>(VeeHandle,e,fVeeLabel);
1930  }
1931 
1932  // Get handle for CaloClusters; also Assn for matching tracks
1933  art::InputTag ecalclustertag(fClusterLabel, fInstanceLabelCalo);
1934  art::InputTag muidclustertag(fClusterMuIDLabel, fInstanceLabelMuID);
1935  art::Handle< std::vector<rec::Cluster> > RecoClusterHandle;
1936  art::Handle< std::vector<rec::Cluster> > RecoClusterMuIDHandle;
1937  art::FindManyP<rec::Track, rec::TrackEnd>* findManyCALTrackEnd = NULL;
1938  art::FindManyP<gar::rec::CaloHit>* findManyClusterRecoHit = NULL;
1939  art::FindManyP<gar::rec::CaloHit>* findManyClusterMuIDHit = NULL;
1940 
1941  if (fWriteCaloClusters) {
1942  RecoClusterHandle = e.getHandle< std::vector<rec::Cluster> >(ecalclustertag);
1943  if (!RecoClusterHandle) {
1944  throw cet::exception("anatree") << " No rec::Cluster branch."
1945  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1946  }
1947 
1948  if(fWriteMuID) {
1949  if (fGeo->HasMuonDetector()){
1950  RecoClusterMuIDHandle = e.getHandle< std::vector<rec::Cluster> >(muidclustertag);
1951  if (!RecoClusterMuIDHandle) {
1952  throw cet::exception("anatree") << " No rec::Cluster (MuID) branch."
1953  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
1954  }
1955  }
1956  }
1957 
1958  findManyClusterRecoHit = new art::FindManyP<gar::rec::CaloHit>(RecoClusterHandle,e,ecalclustertag);
1959 
1960  if (fGeo->HasMuonDetector() && fWriteMuID)
1961  findManyClusterMuIDHit = new art::FindManyP<gar::rec::CaloHit>(RecoClusterMuIDHandle,e,muidclustertag);
1962 
1963 
1964  if (fWriteTracks)
1965  findManyCALTrackEnd = new art::FindManyP<rec::Track, rec::TrackEnd>(RecoClusterHandle,e,fECALAssnLabel);
1966 
1967  }
1968 
1969  // save clusters in the TPC. For some reason, can't get FindOneP<rec::Track> or
1970  // FindManyP<rec::Track> to work; seems the underlying Assn isn't found. Have
1971  // to FindManyP<TPCCluster> instead and iterate if (fWriteTracks). :(
1972  if (fWriteTPCClusters) {
1973  size_t iTPCCluster = 0;
1974  for ( auto const& TPCCluster : (*TPCClusterHandle) ) {
1975  fTPCClusterX.push_back(TPCCluster.Position()[0]);
1976  fTPCClusterY.push_back(TPCCluster.Position()[1]);
1977  fTPCClusterZ.push_back(TPCCluster.Position()[2]);
1978  fTPCClusterSig.push_back(TPCCluster.Signal());
1979  fTPCClusterRMS.push_back(TPCCluster.RMS());
1980  const float* cov;
1981  cov = TPCCluster.CovMatPacked();
1982  fTPCClusterCovXX.push_back(cov[0]);
1983  fTPCClusterCovXY.push_back(cov[1]);
1984  fTPCClusterCovXZ.push_back(cov[2]);
1985  fTPCClusterCovYY.push_back(cov[3]);
1986  fTPCClusterCovYZ.push_back(cov[4]);
1987  fTPCClusterCovZZ.push_back(cov[5]);
1988 
1989  // To get MC matching info from TPCClusters, 1st get the associated Hits
1990  int indexToPush = -1; float valueToPush = 0;
1991  if ( findManyHits->isValid() ) {
1992  std::map<int,float> sumEforTrkID;
1993  float eTotCluster = 0;
1994  auto const& hitsInTPCCluster = findManyHits->at(iTPCCluster);
1995  for (size_t iHits = 0; iHits<hitsInTPCCluster.size(); ++iHits) {
1996  std::vector<cheat::HitIDE> IDEs = BackTrack->HitToHitIDEs( hitsInTPCCluster[iHits] );
1997  for (size_t iIDE = 0; iIDE<IDEs.size(); ++iIDE) {
1998  int trackID = IDEs[iIDE].trackID;
1999  float thisEdepE = IDEs[iIDE].energyTot;
2000  if ( sumEforTrkID.find(trackID) == sumEforTrkID.end() ) {
2001  sumEforTrkID[trackID] = 0;
2002  }
2003  sumEforTrkID[trackID] += thisEdepE;
2004  eTotCluster += thisEdepE;
2005  }
2006  }
2007  if (sumEforTrkID.size()!=0) {
2008  // Sort the map by value. Start by declaring the type of the sorting predicate
2009  typedef std::function<bool(std::pair<int,float>, std::pair<int,float>)> Comparator;
2010  // Declare a set that will store the pairs using above comparison logic
2011  std::set<std::pair<int,float>, Comparator> setOfTrkIDs(
2012  sumEforTrkID.begin(), sumEforTrkID.end(),
2013  [](std::pair<int,float> a ,std::pair<int,float> b) {
2014  return a.second > b.second;
2015  }
2016  );
2017  auto iReturnSet = setOfTrkIDs.begin();
2018  indexToPush = TrackIdToIndex[iReturnSet->first];
2019  valueToPush = iReturnSet->second/eTotCluster;
2020  }
2021  }
2022  fTPCClusterMCindex.push_back(indexToPush);
2023  fTPCClusterMCfrac.push_back(valueToPush);
2024 
2025  Int_t trackForThisTPCluster = -1;
2026  if (fWriteTracks) {
2027  size_t iTrack = 0;
2028  for ( auto const& track : (*TrackHandle) ) {
2029  for (size_t iCluster=0; iCluster<track.NHits(); iCluster++) {
2030  auto const& trackedCluster =
2031  *(findManyTPCClusters->at(iTrack).at(iCluster));
2032  if (TPCCluster==trackedCluster) {
2033  trackForThisTPCluster = track.getIDNumber();
2034  // No cluster is in 2 tracks (don't mess up dE/dx!)
2035  goto pushit; // break 2 loops
2036  }
2037  }
2038  iTrack++;
2039  }
2040  }
2041  pushit:
2042  fTPCClusterTrkIDNumber.push_back(trackForThisTPCluster);
2043  iTPCCluster++;
2044  }
2045  }
2046 
2047  // save per-track info
2048  if (fWriteTracks) {
2049  size_t iTrack = 0;
2050  for ( auto const& track : (*TrackHandle) ) {
2051  // track is a rec::Track, not a rec::TrackPar
2052  fTrackIDNumber.push_back(track.getIDNumber());
2053 
2054  fTrackStartX.push_back(track.Vertex()[0]);
2055  fTrackStartY.push_back(track.Vertex()[1]);
2056  fTrackStartZ.push_back(track.Vertex()[2]);
2057  fTrackStartPX.push_back(track.Momentum_beg()*track.VtxDir()[0]);
2058  fTrackStartPY.push_back(track.Momentum_beg()*track.VtxDir()[1]);
2059  fTrackStartPZ.push_back(track.Momentum_beg()*track.VtxDir()[2]);
2060  fTrackStartQ.push_back(track.ChargeBeg());
2061 
2062  fTrackEndX.push_back(track.End()[0]);
2063  fTrackEndY.push_back(track.End()[1]);
2064  fTrackEndZ.push_back(track.End()[2]);
2065  fTrackEndPX.push_back(track.Momentum_end()*track.EndDir()[0]);
2066  fTrackEndPY.push_back(track.Momentum_end()*track.EndDir()[1]);
2067  fTrackEndPZ.push_back(track.Momentum_end()*track.EndDir()[2]);
2068  fTrackEndQ.push_back(track.ChargeEnd());
2069 
2070  fTrackLenF.push_back(track.LengthForward());
2071  fTrackLenB.push_back(track.LengthBackward());
2072  fTrackChi2F.push_back(track.ChisqForward());
2073  fTrackChi2B.push_back(track.ChisqBackward());
2074  fNTPCClustersOnTrack.push_back(track.NHits());
2075 
2076  //Add the PID information based on Tom's parametrization
2077  TVector3 momF(track.Momentum_beg()*track.VtxDir()[0], track.Momentum_beg()*track.VtxDir()[1], track.Momentum_beg()*track.VtxDir()[2]);
2078  TVector3 momB(track.Momentum_end()*track.EndDir()[0], track.Momentum_end()*track.EndDir()[1], track.Momentum_end()*track.EndDir()[2]);
2079 
2080  //Reconstructed momentum forward and backward
2081  float pF = momF.Mag();
2082  float pB = momB.Mag();
2083  std::vector< std::pair<int, float> > pidF = processPIDInfo( pF );
2084  std::vector< std::pair<int, float> > pidB = processPIDInfo( pB );
2085 
2086  //Fill the pid and its probability
2087  for(size_t ipid = 0; ipid < pidF.size(); ipid++) {
2088  fTrackPIDF.push_back( pidF.at(ipid).first );
2089  fTrackPIDProbF.push_back( pidF.at(ipid).second );
2090  }
2091  for(size_t ipid = 0; ipid < pidB.size(); ipid++) {
2092  fTrackPIDB.push_back( pidB.at(ipid).first );
2093  fTrackPIDProbB.push_back( pidB.at(ipid).second );
2094  }
2095 
2096  // Matching MCParticle info
2097  std::vector<std::pair<simb::MCParticle*,float>> trakt;
2098  trakt = BackTrack->TrackToMCParticles( const_cast<rec::Track*>(&track) );
2099  int eileen = -1;
2100  if (trakt.size()>0 && TrackIdToIndex.size()!=0) {
2101  eileen = TrackIdToIndex[trakt[0].first->TrackId()];
2102  }
2103  fTrackMCindex.push_back(eileen);
2104  if (eileen > -1) {
2105  fTrackMCfrac.push_back(trakt[0].second);
2106  } else {
2107  fTrackMCfrac.push_back(0.0);
2108  }
2109 
2110  if (findIonization->isValid()) {
2111  // No calibration for now. Someday this should all be in reco
2112  rec::TrackIoniz ionization = *(findIonization->at(iTrack));
2113  float avgIonF, avgIonB;
2114  processIonizationInfo(ionization, fIonizTruncate, avgIonF, avgIonB);
2115  fTrackAvgIonF.push_back( avgIonF );
2116  fTrackAvgIonB.push_back( avgIonB );
2117  } else {
2118  // must push_back something so that fTrackAvgIonF,B are of correct size.
2119  fTrackAvgIonF.push_back( 0.0 );
2120  fTrackAvgIonB.push_back( 0.0 );
2121  }
2122  iTrack++;
2123  } // end loop over TrackHandle
2124  }
2125 
2126  //TrackTrajectories
2127  if(fWriteTrackTrajectories) {
2128  size_t iTrackTraj = 0;
2129  for ( auto const& tracktraj : (*TrackTrajHandle) ) {
2130 
2131  std::vector<TVector3> temp = tracktraj.getFWDTrajectory();
2132  for(size_t i = 0; i < temp.size(); i++) {
2133  fTrackTrajectoryFWDX.push_back(temp.at(i).X());
2134  fTrackTrajectoryFWDY.push_back(temp.at(i).Y());
2135  fTrackTrajectoryFWDZ.push_back(temp.at(i).Z());
2136  fTrackTrajectoryFWDID.push_back(iTrackTraj);
2137  }
2138 
2139  temp = tracktraj.getBAKTrajectory();
2140  for(size_t i = 0; i < temp.size(); i++) {
2141  fTrackTrajectoryBWDX.push_back(temp.at(i).X());
2142  fTrackTrajectoryBWDY.push_back(temp.at(i).Y());
2143  fTrackTrajectoryBWDZ.push_back(temp.at(i).Z());
2144  fTrackTrajectoryBWDID.push_back(iTrackTraj);
2145  }
2146  iTrackTraj++;
2147  }
2148  }
2149 
2150  // save Vertex and Track-Vertex association info
2151  if (fWriteVertices) {
2152  size_t iVertex = 0;
2153  for ( auto const& vertex : (*VertexHandle) ) {
2154  fVertexIDNumber.push_back(vertex.getIDNumber());
2155  fVertexX.push_back(vertex.Position()[0]);
2156  fVertexY.push_back(vertex.Position()[1]);
2157  fVertexZ.push_back(vertex.Position()[2]);
2158  fVertexT.push_back(vertex.Time());
2159 
2160  int nVertexedTracks = 0;
2161  if ( findManyTrackEnd->isValid() ) {
2162  nVertexedTracks = findManyTrackEnd->at(iVertex).size();
2163  }
2164  fVertexN.push_back(nVertexedTracks);
2165 
2166  int vertexCharge = 0;
2167  for (int iVertexedTrack=0; iVertexedTrack<nVertexedTracks; ++iVertexedTrack) {
2168  fVTAssn_VertIDNumber.push_back(vertex.getIDNumber());
2169 
2170  // Get this vertexed track.
2171  rec::Track track = *(findManyTrackEnd->at(iVertex).at(iVertexedTrack));
2172  fVTAssn_TrackIDNumber.push_back(track.getIDNumber());
2173 
2174  // Get the end of the track in the vertex. It isn't that odd for the end
2175  // of the track not used in the vertex to be closer to the vertex than the
2176  // one actually used; you might have a very short stub track in a 3 track
2177  // vertex with small opening angles and the other 2 tracks might pull the
2178  // vertex some distance towards the far end of the stub track
2179  rec::TrackEnd fee = *(findManyTrackEnd->data(iVertex).at(iVertexedTrack));
2180  // TrackEnd is defined in Track.h; 1 means use Beg values, 0 means use End
2181  fVTAssn_TrackEnd.push_back(fee);
2182 
2183  if (fee==rec::TrackEndBeg) {
2184  vertexCharge += track.ChargeBeg();
2185  } else {
2186  vertexCharge += track.ChargeEnd();
2187  }
2188  }
2189  fVertexQ.push_back(vertexCharge);
2190  ++iVertex;
2191  } // end loop over VertexHandle
2192  }
2193 
2194  // save Vee and Track-Vee association info
2195  if (fWriteVees) {
2196  size_t iVee = 0;
2197  for ( auto const& vee : (*VeeHandle) ) {
2198  fVeeIDNumber.push_back(vee.getIDNumber());
2199  fVeeX.push_back(vee.Position()[0]);
2200  fVeeY.push_back(vee.Position()[1]);
2201  fVeeZ.push_back(vee.Position()[2]);
2202  fVeeT.push_back(vee.Time());
2203  fVeePXKpipi.push_back(vee.FourMomentum(0).X());
2204  fVeePYKpipi.push_back(vee.FourMomentum(0).Y());
2205  fVeePZKpipi.push_back(vee.FourMomentum(0).Z());
2206  fVeeEKpipi.push_back(vee.FourMomentum(0).E());
2207  fVeeMKpipi.push_back(vee.FourMomentum(0).M());
2208  fVeePXLppi.push_back(vee.FourMomentum(1).X());
2209  fVeePYLppi.push_back(vee.FourMomentum(1).Y());
2210  fVeePZLppi.push_back(vee.FourMomentum(1).Z());
2211  fVeeELppi.push_back(vee.FourMomentum(1).E());
2212  fVeeMLppi.push_back(vee.FourMomentum(1).M());
2213  fVeePXLpip.push_back(vee.FourMomentum(2).X());
2214  fVeePYLpip.push_back(vee.FourMomentum(2).Y());
2215  fVeePZLpip.push_back(vee.FourMomentum(2).Z());
2216  fVeeELpip.push_back(vee.FourMomentum(2).E());
2217  fVeeMLpip.push_back(vee.FourMomentum(2).M());
2218 
2219  int nVeeTracks = 0;
2220  if ( findManyVeeTrackEnd->isValid() ) {
2221  nVeeTracks = findManyVeeTrackEnd->at(iVee).size();
2222  }
2223 
2224  for (int iVeeTrack=0; iVeeTrack<nVeeTracks; ++iVeeTrack) {
2225  fVeeTAssn_VeeIDNumber.push_back(vee.getIDNumber());
2226 
2227  // Get this vertexed track.
2228  rec::Track track = *(findManyVeeTrackEnd->at(iVee).at(iVeeTrack));
2229  fVeeTAssn_TrackIDNumber.push_back(track.getIDNumber());
2230 
2231  rec::TrackEnd fee = *(findManyVeeTrackEnd->data(iVee).at(iVeeTrack));
2232  // TrackEnd is defined in Track.h; 1 means use Beg values, 0 means use End
2233  fVeeTAssn_TrackEnd.push_back(fee);
2234  }
2235  ++iVee;
2236  } // end loop over VeeHandle
2237  }
2238 
2239  // save Cluster info
2240  if (fWriteCaloClusters) {
2241  size_t iCluster = 0;
2242  for ( auto const& cluster : (*RecoClusterHandle) ) {
2243  fnCluster++;
2244  fClusterIDNumber.push_back(cluster.getIDNumber());
2245  fClusterNhits.push_back(cluster.CalorimeterHits().size());
2246  fClusterEnergy.push_back(cluster.Energy());
2247  fClusterTime.push_back(cluster.Time());
2248  fClusterTimeDiffFirstLast.push_back(cluster.TimeDiffFirstLast());
2249  fClusterX.push_back(cluster.Position()[0]);
2250  fClusterY.push_back(cluster.Position()[1]);
2251  fClusterZ.push_back(cluster.Position()[2]);
2252  fClusterTheta.push_back(cluster.ITheta());
2253  fClusterPhi.push_back(cluster.IPhi());
2254  fClusterPID.push_back(cluster.ParticleID());
2255  // fClusterShape.push_back(cluster.Shape());
2256  fClusterMainAxisX.push_back(cluster.EigenVectors()[0]);
2257  fClusterMainAxisY.push_back(cluster.EigenVectors()[1]);
2258  fClusterMainAxisZ.push_back(cluster.EigenVectors()[2]);
2259 
2260  //Get the associated reco hit
2261  std::vector<ULong64_t> fVecHitIDs = {};
2262  if (findManyClusterRecoHit->isValid()) {
2263  int nClusterHit = findManyClusterRecoHit->at(iCluster).size();
2264  for (int iClusterHit=0; iClusterHit<nClusterHit; ++iClusterHit) {
2265  rec::CaloHit hit = *(findManyClusterRecoHit->at(iCluster).at(iClusterHit));
2266  fVecHitIDs.push_back(hit.getIDNumber());
2267  }
2268  }
2269 
2270  fClusterAssn_RecoHitIDNumber.push_back(fVecHitIDs);
2271 
2272  // Matching MCParticle info
2273  std::vector<std::pair<simb::MCParticle*,float>> trakt;
2274  trakt = BackTrack->ClusterToMCParticles( const_cast<rec::Cluster*>(&cluster) );
2275  int eileen = -1;
2276  if (trakt.size()>0 && TrackIdToIndex.size()!=0) {
2277  eileen = TrackIdToIndex[trakt[0].first->TrackId()];
2278  }
2279  fClusterMCindex.push_back(eileen);
2280  if (eileen > -1) {
2281  fClusterMCfrac.push_back(trakt[0].second);
2282  } else {
2283  fClusterMCfrac.push_back(0.0);
2284  }
2285  iCluster++;
2286  }
2287 
2288  if(fGeo->HasMuonDetector() && fWriteMuID) {
2289  size_t iCluster_local = 0;
2290  for ( auto const& cluster : (*RecoClusterMuIDHandle) ) {
2291  fnCluster_MuID++;
2292  fClusterIDNumber_MuID.push_back(cluster.getIDNumber());
2293  fClusterNhits_MuID.push_back(cluster.CalorimeterHits().size());
2294  fClusterEnergy_MuID.push_back(cluster.Energy());
2295  fClusterTime_MuID.push_back(cluster.Time());
2296  fClusterTimeDiffFirstLast_MuID.push_back(cluster.TimeDiffFirstLast());
2297  fClusterX_MuID.push_back(cluster.Position()[0]);
2298  fClusterY_MuID.push_back(cluster.Position()[1]);
2299  fClusterZ_MuID.push_back(cluster.Position()[2]);
2300  fClusterTheta_MuID.push_back(cluster.ITheta());
2301  fClusterPhi_MuID.push_back(cluster.IPhi());
2302  fClusterPID_MuID.push_back(cluster.ParticleID());
2303  // fClusterShape.push_back(cluster.Shape());
2304  fClusterMainAxisX_MuID.push_back(cluster.EigenVectors()[0]);
2305  fClusterMainAxisY_MuID.push_back(cluster.EigenVectors()[1]);
2306  fClusterMainAxisZ_MuID.push_back(cluster.EigenVectors()[2]);
2307 
2308  //Get the associated reco hit
2309  std::vector<ULong64_t> fVecHitIDs = {};
2310  if (findManyClusterMuIDHit->isValid()) {
2311  int nClusterHit = findManyClusterMuIDHit->at(iCluster_local).size();
2312  for (int iClusterHit=0; iClusterHit<nClusterHit; ++iClusterHit) {
2313  rec::CaloHit hit = *(findManyClusterMuIDHit->at(iCluster_local).at(iClusterHit));
2314  fVecHitIDs.push_back(hit.getIDNumber());
2315  }
2316  }
2317 
2318  fClusterMuIDAssn_MuIDHitIDNumber.push_back(fVecHitIDs);
2319 
2320  /*** TO DO ***/
2321  // Matching MCParticle info
2322  // std::vector<std::pair<simb::MCParticle*,float>> trakt;
2323  // trakt = BackTrack->ClusterToMCParticles( const_cast<rec::Cluster*>(&cluster) );
2324  // int eileen = -1;
2325  // if (trakt.size()>0 && TrackIdToIndex.size()!=0) {
2326  // eileen = TrackIdToIndex[trakt[0].first->TrackId()];
2327  // }
2328  // fClusterMCindex_MuID.push_back(eileen);
2329  // if (eileen > -1) {
2330  // fClusterMCfrac_MuID.push_back(trakt[0].second);
2331  // } else {
2332  // fClusterMCfrac_MuID.push_back(0.0);
2333  // }
2334  fClusterMCindex_MuID.push_back(0);
2335  fClusterMCfrac_MuID.push_back(0.0);
2336  iCluster_local++;
2337  }
2338  }
2339  }
2340 
2341  // Write info for ECAL-matched tracks
2342  if (fWriteMatchedTracks) {
2343  size_t iCluster = 0;
2344  for ( auto const& cluster : (*RecoClusterHandle) ) {
2345  int nCALedTracks(0);
2346  if ( findManyCALTrackEnd->isValid() ) {
2347  nCALedTracks = findManyCALTrackEnd->at(iCluster).size();
2348  }
2349  for (int iCALedTrack=0; iCALedTrack<nCALedTracks; ++iCALedTrack) {
2350  fCALAssn_ClusIDNumber.push_back(cluster.getIDNumber());
2351  rec::Track track = *(findManyCALTrackEnd->at(iCluster).at(iCALedTrack));
2352  fCALAssn_TrackIDNumber.push_back( track.getIDNumber() );
2353 
2354  rec::TrackEnd fee = *(findManyCALTrackEnd->data(iCluster).at(iCALedTrack));
2355  fCALAssn_TrackEnd.push_back(fee); // The rec::TrackEnd (see Track.h) that extrapolated to cluster
2356  }
2357  iCluster++;
2358  }
2359  } // end branch on fWriteCaloInfo
2360 
2361  return;
2362 } // end :anatree::FillVectors
2363 
2364 
2365 
2366  //==============================================================================
2367  //==============================================================================
2368  //==============================================================================
2369  // Process ionization. Eventually this moves into the reco code.
2370 void gar::anatree::processIonizationInfo(rec::TrackIoniz& ion, float ionizeTruncate,
2371  float& forwardIonVal, float& backwardIonVal) {
2372 
2373  // NO CALIBRATION SERVICE FOR NOW
2374 
2375  std::vector<std::pair<float,float>> SigData = ion.getFWD_dSigdXs();
2376  forwardIonVal = processOneDirection(SigData, ionizeTruncate);
2377 
2378  SigData = ion.getBAK_dSigdXs();
2379  backwardIonVal = processOneDirection(SigData, ionizeTruncate);
2380 
2381  return;
2382 }
2383 
2384 
2385 
2386 float gar::anatree::processOneDirection(std::vector<std::pair<float,float>> SigData, float ionizeTruncate) {
2387 
2388  std::vector<std::pair<float,float>> dEvsX; // Ionization vs distance along track
2389 
2390  // The first hit on the track never had its ionization info stored. Not a problem
2391  // really. Each pair is a hit and the step along the track that ends at the hit
2392  // For the last hit, just take the step from the n-1 hit; don't guess some distance
2393  // to (nonexistant!) n+1 hit. Using pointer arithmetic because you are a real K&R
2394  // C nerd! Except that C++ doesn't know you are such a nerd and if
2395  // SigData.size()==0, then SigData.end()-1 is 0xFFFFFFFFFFFFFFF8.
2396  if (SigData.size()==0) return 0.0;
2397  float distAlongTrack = 0;
2398  std::vector<std::pair<float,float>>::iterator littlebit = SigData.begin();
2399  for (; littlebit<(SigData.end()-1); ++littlebit) {
2400  float dE = std::get<0>(*littlebit);
2401  // tpctrackfit2_module.cc fills the TrackIoniz data product so that
2402  // this quantity is really dL > 0 not dX, a coordinate on the drift axis
2403  float dX = std::get<1>(*littlebit);
2404  distAlongTrack += dX; // But count full step to get hit position on track
2405  // Take dX to be 1/2 the previous + last segment
2406  dX += std::get<1>(*(littlebit+1));
2407  float dEdX = dE/(0.5*dX);
2408 
2409  std::pair pushme = std::make_pair(dEdX,distAlongTrack);
2410  dEvsX.push_back( pushme );
2411  }
2412 
2413  // Get the truncated mean; first sort then take mean
2414  std::sort(dEvsX.begin(),dEvsX.end(), lessThan_byE);
2415 
2416  // Get the dEdX vs length data, truncated.
2417  int goUpTo = ionizeTruncate * dEvsX.size() +0.5;
2418  if (goUpTo > (int)dEvsX.size()) goUpTo = dEvsX.size();
2419  int i = 1; float returnvalue = 0;
2420  littlebit = dEvsX.begin();
2421  for (; littlebit<dEvsX.end(); ++littlebit) {
2422  returnvalue += std::get<0>(*littlebit);
2423  ++i;
2424  if (i>goUpTo) break;
2425  }
2426  returnvalue /= goUpTo;
2427  return returnvalue;
2428 }
2429 
2430 
2431 
2432 //==============================================================================
2433 //==============================================================================
2434 //==============================================================================
2435 // Coherent pion analysis specific code
2436 float gar::anatree::computeT( simb::MCTruth theMCTruth ) {
2437  // Warning. You probably want the absolute value of t, not t.
2438  int nPart = theMCTruth.NParticles();
2439  enum { nu, mu, pi};
2440  float E[3], Px[3], Py[3], Pz[3];
2441  E[nu] = E[mu] = E[pi] = -1e42;
2442 
2443  for (int i=0; i<3;++i) {
2444  Px[i] = 0;
2445  Py[i] = 0;
2446  Pz[i] = 0;
2447  E[i] = 0;
2448  }
2449  // Find t from the MCParticles via the
2450  for (int iPart=0; iPart<nPart; iPart++) {
2451  simb::MCParticle Part = theMCTruth.GetParticle(iPart);
2452  int code = Part.PdgCode();
2453  int mom = Part.Mother();
2454 
2455  // get the neutrino
2456  if ( abs(code) == 12 || abs(code) == 14 || abs(code) == 16 ) {
2457  if (mom == -1) {
2458  E[nu] = Part.E(); Px[nu] = Part.Px(); Py[nu] = Part.Py(); Pz[nu] = Part.Pz();
2459  }
2460  }
2461 
2462  // get the lepton
2463  if ( abs(code) == 11 || abs(code) == 13 || abs(code) == 15 ) {
2464  if (mom == 0) {
2465  E[mu] = Part.E(); Px[mu] = Part.Px(); Py[mu] = Part.Py(); Pz[mu] = Part.Pz();
2466  }
2467  }
2468 
2469  // get the pion
2470  if ( code==111 || abs(code)==211 ) {
2471  if (mom == 1) {
2472  E[pi] = Part.E(); Px[pi] = Part.Px(); Py[pi] = Part.Py(); Pz[pi] = Part.Pz();
2473  }
2474  }
2475 
2476  // get outa here
2477  if ( E[nu]!=0 && E[mu]!=0 && E[pi]!=0) break;
2478 
2479  }
2480 
2481  // Compute t; reuse nu 4-vector to get first q, then t.
2482  E[nu] -= E[mu]; Px[nu] -= Px[mu]; Py[nu] -= Py[mu]; Pz[nu] -= Pz[mu];
2483  E[nu] -= E[pi]; Px[nu] -= Px[pi]; Py[nu] -= Py[pi]; Pz[nu] -= Pz[pi];
2484  float t = E[nu]*E[nu] -Px[nu]*Px[nu] -Py[nu]*Py[nu] -Pz[nu]*Pz[nu];
2485  return t;
2486 }
2487 
2488 
2489 
2490 //==============================================================================
2491 //==============================================================================
2492 //==============================================================================
2493 std::vector< std::pair<int, float> > gar::anatree::processPIDInfo( float p ) {
2494 
2495  std::vector<std::string> recopnamelist = {"#pi", "#mu", "p", "K", "d", "e"};
2496  std::vector<int> pdg_charged = {211, 13, 2212, 321, 1000010020, 11};
2497  std::vector< std::pair<int, float> > pid;
2498  pid.resize(6);
2499 
2500  int qclosest = 0;
2501  float dist = 100000000.;
2502  CLHEP::RandFlat FlatRand(fEngine);
2503 
2504  for (int q = 0; q < 501; ++q) {
2505  //Check the title and the reco momentum take only the one that fits
2506  std::string fulltitle = m_pidinterp[q]->GetTitle();
2507  unsigned first = fulltitle.find("=");
2508  unsigned last = fulltitle.find("GeV");
2509  std::string substr = fulltitle.substr(first+1, last - first-1);
2510  float pidinterp_mom = std::atof(substr.c_str());
2511  //calculate the distance between the bin and mom, store the q the closest
2512  float disttemp = std::abs(pidinterp_mom - p);
2513 
2514  if( disttemp < dist ) {
2515  dist = disttemp;
2516  qclosest = q;
2517  }
2518  } // closes the "pidmatrix" loop
2519 
2520  //Compute all the probabities for each type of true to reco
2521  //loop over the columns (true pid)
2522  for (int pidm = 0; pidm < 6; ++pidm) {
2523 
2524  //loop over the columns (true pid)
2525  std::vector< std::pair<float, std::string> > v_prob;
2526 
2527  //loop over the rows (reco pid)
2528  for (int pidr = 0; pidr < 6; ++pidr) {
2529  std::string recoparticlename = m_pidinterp[qclosest]->GetYaxis()->GetBinLabel(pidr+1);
2530  float prob = m_pidinterp[qclosest]->GetBinContent(pidm+1, pidr+1);
2531  //Need to check random number value and prob value then associate the recopdg to the reco prob
2532  v_prob.push_back( std::make_pair(prob, recoparticlename) );
2533  }
2534 
2535  //Compute the pid from it
2536  if (v_prob.size() > 1) {
2537  //Order the vector of prob
2538  std::sort(v_prob.begin(), v_prob.end());
2539  //Throw a random number between 0 and 1
2540  float random_number = FlatRand.fire();
2541  //Make cumulative sum to get the range
2542  std::partial_sum(v_prob.begin(), v_prob.end(), v_prob.begin(), [](const P& _x, const P& _y){return P(_x.first + _y.first, _y.second);});
2543 
2544  for(size_t ivec = 0; ivec < v_prob.size()-1; ivec++) {
2545  if( random_number < v_prob.at(ivec+1).first && random_number >= v_prob.at(ivec).first ) {
2546  pid.push_back( std::make_pair(pdg_charged.at( std::distance( recopnamelist.begin(), std::find(recopnamelist.begin(), recopnamelist.end(), v_prob.at(ivec+1).second) ) ), v_prob.at(ivec+1).first) );
2547  }
2548  }
2549  } else {
2550  pid.push_back( std::make_pair(pdg_charged.at( std::distance( recopnamelist.begin(), std::find(recopnamelist.begin(), recopnamelist.end(), v_prob.at(0).second) ) ), v_prob.at(0).first) );
2551  }
2552  }
2553 
2554  //return a vector of pid and prob
2555  return pid;
2556 }
2557 
2558 DEFINE_ART_MODULE(gar::anatree)
gar::anatree::fRecoHitCellID
std::vector< ULong64_t > fRecoHitCellID
Definition: anatree_module.cc:438
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
art::detail::EngineCreator::createEngine
base_engine_t & createEngine(seed_t seed)
Definition: EngineCreator.cc:20
gar::anatree::fWriteHits
bool fWriteHits
Write info about TPC Hits Default=false.
Definition: anatree_module.cc:160
gar::anatree::fVeeY
std::vector< Float_t > fVeeY
Definition: anatree_module.cc:390
gar::anatree::fIonizTruncate
float fIonizTruncate
Default=1.00;.
Definition: anatree_module.cc:174
gar::anatree::fClusterPhi
std::vector< Float_t > fClusterPhi
Definition: anatree_module.cc:465
gar::anatree::fRecoHitY_MuID
std::vector< Float_t > fRecoHitY_MuID
Definition: anatree_module.cc:446
gar::anatree::fSimHitTime_MuID
std::vector< Float_t > fSimHitTime_MuID
Definition: anatree_module.cc:300
gar::anatree::fGENIEGeneratorLabels
std::vector< std::string > fGENIEGeneratorLabels
Definition: anatree_module.cc:127
gar::anatree::fVertexY
std::vector< Float_t > fVertexY
Definition: anatree_module.cc:377
gar::anatree::fMCPDGMother
std::vector< Int_t > fMCPDGMother
Definition: anatree_module.cc:254
gar::anatree::fMCPStartX
std::vector< Float_t > fMCPStartX
Definition: anatree_module.cc:255
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
simb::MCNeutrino::CCNC
int CCNC() const
Definition: MCNeutrino.h:148
gar::anatree::fTrackChi2F
std::vector< Float_t > fTrackChi2F
Definition: anatree_module.cc:350
gar::anatree::fGPartFirstDaugh
std::vector< Int_t > fGPartFirstDaugh
Definition: anatree_module.cc:240
gar::anatree::fVTAssn_TrackIDNumber
std::vector< ULong64_t > fVTAssn_TrackIDNumber
Definition: anatree_module.cc:384
simb::MCNeutrino::QSqr
double QSqr() const
Definition: MCNeutrino.h:157
gar::anatree::fGPartFirstMom
std::vector< Int_t > fGPartFirstMom
Definition: anatree_module.cc:238
simb::MCNeutrino::Theta
double Theta() const
angle between incoming and outgoing leptons, in radians
Definition: MCNeutrino.cxx:63
gar::rec::TrackEnd
int TrackEnd
Definition: Track.h:32
gar::anatree::fY
std::vector< Float_t > fY
Definition: anatree_module.cc:216
simb::MCParticle::Py
double Py(const int i=0) const
Definition: MCParticle.h:231
gar::anatree::fRun
Int_t fRun
number of the run being processed
Definition: anatree_module.cc:193
gar::anatree::processPIDInfo
std::vector< std::pair< int, float > > processPIDInfo(float p)
Definition: anatree_module.cc:2493
gar::anatree::fRecoHitZ
std::vector< Float_t > fRecoHitZ
Definition: anatree_module.cc:435
gar::anatree::fDigiHitADC_MuID
std::vector< UInt_t > fDigiHitADC_MuID
Definition: anatree_module.cc:427
gar::anatree::fMCPStartPZ
std::vector< Float_t > fMCPStartPZ
Definition: anatree_module.cc:261
art::Run::id
RunID id() const
Definition: Run.cc:17
gar::anatree::fWriteTPCClusters
bool fWriteTPCClusters
Write TPCClusters info Default=true.
Definition: anatree_module.cc:161
gar::anatree::fDigiHitADC
std::vector< UInt_t > fDigiHitADC
Definition: anatree_module.cc:418
simb::MCParticle::EndZ
double EndZ() const
Definition: MCParticle.h:228
gar::anatree::fClusterZ
std::vector< Float_t > fClusterZ
Definition: anatree_module.cc:463
gar::anatree::fTPC_X
Float_t fTPC_X
center of TPC stored as per-event & compressed by root
Definition: anatree_module.cc:195
gar::anatree::fTrackLenB
std::vector< Float_t > fTrackLenB
Definition: anatree_module.cc:349
gar::anatree::fGint
std::vector< Int_t > fGint
Definition: anatree_module.cc:227
gar::anatree::fWriteMCCaloInfo
bool fWriteMCCaloInfo
Write MC info for calorimeter Default=true.
Definition: anatree_module.cc:157
gar::anatree::fWriteVees
bool fWriteVees
Reco vees & their tracks Default=true.
Definition: anatree_module.cc:165
gar::geo::GeometryCore::GetMuIDCellIDEncoding
const std::string GetMuIDCellIDEncoding() const
Definition: GeometryCore.cxx:856
gar::anatree::fClusterX
std::vector< Float_t > fClusterX
Definition: anatree_module.cc:461
gar::anatree::fMCPVertIndex
std::vector< Int_t > fMCPVertIndex
Definition: anatree_module.cc:270
gar::anatree::fSimHitX
std::vector< Float_t > fSimHitX
Definition: anatree_module.cc:286
gar::anatree::fX
std::vector< Float_t > fX
Definition: anatree_module.cc:215
gar::anatree::fDigiHitZ_MuID
std::vector< Float_t > fDigiHitZ_MuID
Definition: anatree_module.cc:425
gar::anatree::fVeePXLpip
std::vector< Float_t > fVeePXLpip
Definition: anatree_module.cc:403
gar::anatree::fVeeIDNumber
std::vector< ULong64_t > fVeeIDNumber
Definition: anatree_module.cc:388
art::DataViewImpl::getHandle
Handle< PROD > getHandle(SelectorBase const &) const
Definition: DataViewImpl.h:382
string
std::string string
Definition: nybbler.cc:12
Hardware::ID
unsigned int ID
Definition: HardwareElements.h:23
gar::anatree::lessThan_byE
static bool lessThan_byE(std::pair< float, float > a, std::pair< float, float > b)
Definition: anatree_module.cc:111
gar::anatree::fClusterY_MuID
std::vector< Float_t > fClusterY_MuID
Definition: anatree_module.cc:482
gar::anatree::fWriteMCinfo
bool fWriteMCinfo
Info from MCTruth, GTruth Default=true.
Definition: anatree_module.cc:154
gar::anatree::fTrackEndX
std::vector< Float_t > fTrackEndX
Definition: anatree_module.cc:340
gar::anatree::fTPCClusterCovYZ
std::vector< Float_t > fTPCClusterCovYZ
Definition: anatree_module.cc:325
gar::anatree::fTPCClusterMCindex
std::vector< Int_t > fTPCClusterMCindex
Definition: anatree_module.cc:327
gar::anatree::fNeutrinoType
std::vector< Int_t > fNeutrinoType
Definition: anatree_module.cc:209
simb::MCParticle::Mother
int Mother() const
Definition: MCParticle.h:213
gar::anatree::fTPCClusterLabel
std::string fTPCClusterLabel
module label for TPC Clusters rec::TPCCluster
Definition: anatree_module.cc:136
gar::anatree::fWriteTrackTrajectories
bool fWriteTrackTrajectories
Point traj of reco tracks Default=false.
Definition: anatree_module.cc:163
gar::anatree::fWriteVertices
bool fWriteVertices
Reco vertexes & their tracks Default=true.
Definition: anatree_module.cc:164
gar::anatree::fTrackEndPZ
std::vector< Float_t > fTrackEndPZ
Definition: anatree_module.cc:345
gar::anatree::fClusterMuIDAssn_MuIDHitIDNumber
std::vector< std::vector< ULong64_t > > fClusterMuIDAssn_MuIDHitIDNumber
Definition: anatree_module.cc:498
gar::anatree::fVeeTAssn_VeeIDNumber
std::vector< ULong64_t > fVeeTAssn_VeeIDNumber
Definition: anatree_module.cc:408
gar::anatree::fHitZ
std::vector< Float_t > fHitZ
Definition: anatree_module.cc:309
gar::anatree::fClusterTheta
std::vector< Float_t > fClusterTheta
Definition: anatree_module.cc:464
gar::anatree::fSimHitTrackID
std::vector< Int_t > fSimHitTrackID
Definition: anatree_module.cc:291
simb::MCNeutrino::Nu
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
gar::anatree::fVeeTAssn_TrackEnd
std::vector< gar::rec::TrackEnd > fVeeTAssn_TrackEnd
Definition: anatree_module.cc:410
gar::anatree::fDigiHitZ
std::vector< Float_t > fDigiHitZ
Definition: anatree_module.cc:416
gar::anatree::fClusterNhits_MuID
std::vector< UInt_t > fClusterNhits_MuID
Definition: anatree_module.cc:477
gar::anatree::fWeight
std::vector< Float_t > fWeight
Definition: anatree_module.cc:229
gar::anatree::beginJob
virtual void beginJob() override
Definition: anatree_module.cc:661
simb::MCParticle::Px
double Px(const int i=0) const
Definition: MCParticle.h:230
gar::anatree::fClusterZ_MuID
std::vector< Float_t > fClusterZ_MuID
Definition: anatree_module.cc:483
gar::anatree::computeT
float computeT(simb::MCTruth theMCTruth)
Definition: anatree_module.cc:2436
gar::anatree::fGeo
const geo::GeometryCore * fGeo
pointer to the geometry
Definition: anatree_module.cc:181
gar::anatree::fMCPProc
std::vector< std::string > fMCPProc
Definition: anatree_module.cc:268
gar::anatree::fReconHitIDNumber_MuID
std::vector< ULong64_t > fReconHitIDNumber_MuID
Definition: anatree_module.cc:444
vector
struct vector vector
gar::anatree::fClusterTimeDiffFirstLast
std::vector< Float_t > fClusterTimeDiffFirstLast
Definition: anatree_module.cc:460
gar::anatree::fT
std::vector< Float_t > fT
Definition: anatree_module.cc:218
gar::anatree::FillRawInfo
void FillRawInfo(art::Event const &e)
Definition: anatree_module.cc:1744
gar::anatree::fSimHitCellID
std::vector< ULong64_t > fSimHitCellID
Definition: anatree_module.cc:292
gar::anatree::fTrackTrajectoryFWDZ
std::vector< Float_t > fTrackTrajectoryFWDZ
Definition: anatree_module.cc:366
gar::anatree::fVertexX
std::vector< Float_t > fVertexX
Definition: anatree_module.cc:376
gar::anatree::fDigiHitX_MuID
std::vector< Float_t > fDigiHitX_MuID
Definition: anatree_module.cc:423
gar::anatree::fVeePZLppi
std::vector< Float_t > fVeePZLppi
Definition: anatree_module.cc:400
gar::anatree::fTrajMCPTrackID
std::vector< Int_t > fTrajMCPTrackID
Definition: anatree_module.cc:279
gar::geo::GeometryCore::TPCXCent
float TPCXCent() const
Returns the X location of the center of the TPC in cm.
Definition: GeometryCore.h:778
gar::anatree::fTrackStartPZ
std::vector< Float_t > fTrackStartPZ
Definition: anatree_module.cc:337
gar::P
std::pair< float, std::string > P
Definition: anatree_module.cc:75
gar::anatree::fTPCClusterY
std::vector< Float_t > fTPCClusterY
Definition: anatree_module.cc:316
gar::geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:436
gar::anatree::fNTPCClustersOnTrack
std::vector< Int_t > fNTPCClustersOnTrack
Definition: anatree_module.cc:352
gar::anatree::fTrackTrajectoryLabel
std::string fTrackTrajectoryLabel
module label for TPC Track Trajectories rec:TrackTrajectory
Definition: anatree_module.cc:138
gar::anatree::fTrajMCPT
std::vector< Float_t > fTrajMCPT
Definition: anatree_module.cc:276
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:75
gar::anatree::fVeePZKpipi
std::vector< Float_t > fVeePZKpipi
Definition: anatree_module.cc:395
gar::anatree::ItsInTulsa
float ItsInTulsa[3]
Definition: anatree_module.cc:121
gar::anatree::fMCnuPx
std::vector< Float_t > fMCnuPx
Definition: anatree_module.cc:222
cluster
Cluster finding and building.
Definition: SmallClusterFilter_module.cc:34
gar::anatree::fVertexN
std::vector< Int_t > fVertexN
Definition: anatree_module.cc:380
gar::anatree::fMCPStartPX
std::vector< Float_t > fMCPStartPX
Definition: anatree_module.cc:259
gar::anatree::fGPartPdg
std::vector< Int_t > fGPartPdg
Definition: anatree_module.cc:236
gar::anatree::FillHighLevelRecoInfo
void FillHighLevelRecoInfo(art::Event const &e)
Definition: anatree_module.cc:1863
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
gar::anatree::fnCluster_MuID
UInt_t fnCluster_MuID
Definition: anatree_module.cc:475
gar::anatree::fTrajMCPPZ
std::vector< Float_t > fTrajMCPPZ
Definition: anatree_module.cc:282
gar::anatree::fClusterPhi_MuID
std::vector< Float_t > fClusterPhi_MuID
Definition: anatree_module.cc:485
MCParticle.h
Particle class.
simb::MCParticle::EndY
double EndY() const
Definition: MCParticle.h:227
gar::anatree::anatree
anatree(fhicl::ParameterSet const &p)
Definition: anatree_module.cc:517
train.filename
string filename
Definition: train.py:213
gar::anatree::fTrackTrajectoryBWDZ
std::vector< Float_t > fTrackTrajectoryBWDZ
Definition: anatree_module.cc:371
gar::anatree::fMatchMCPtoVertDist
float fMatchMCPtoVertDist
MCParticle to MC vertex match Default=roundoff.
Definition: anatree_module.cc:158
gar::anatree::fTPCClusterTrkIDNumber
std::vector< ULong64_t > fTPCClusterTrkIDNumber
Definition: anatree_module.cc:320
gar::anatree::fTrackEndPY
std::vector< Float_t > fTrackEndPY
Definition: anatree_module.cc:344
gar::anatree::fWriteCaloHits
bool fWriteCaloHits
Write ECAL hits. Default=true.
Definition: anatree_module.cc:169
gar::anatree::fMCnuPz
std::vector< Float_t > fMCnuPz
Definition: anatree_module.cc:224
gar::anatree::fVertexZ
std::vector< Float_t > fVertexZ
Definition: anatree_module.cc:378
gar::anatree::fClusterMCindex_MuID
std::vector< Int_t > fClusterMCindex_MuID
Definition: anatree_module.cc:491
gar::anatree::fClusterMainAxisZ_MuID
std::vector< Float_t > fClusterMainAxisZ_MuID
Definition: anatree_module.cc:490
pionana::dEdX
double dEdX(double KE, const simb::MCParticle *part)
Definition: TruthAnalyzer_module.cc:67
gar::anatree::fMCPStartY
std::vector< Float_t > fMCPStartY
Definition: anatree_module.cc:256
gar::anatree::fSimHitY_MuID
std::vector< Float_t > fSimHitY_MuID
Definition: anatree_module.cc:298
cet::hypot
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
Definition: hypot.h:60
art::Run
Definition: Run.h:17
gar::anatree::fQ2
std::vector< Float_t > fQ2
Definition: anatree_module.cc:213
gar::anatree::fClusterLabel
std::string fClusterLabel
module label for calo clusters rec::Cluster
Definition: anatree_module.cc:148
gar::anatree::fSimHitTrackID_MuID
std::vector< Int_t > fSimHitTrackID_MuID
Definition: anatree_module.cc:302
gar::sumdata::POTSummary::TotalPOT
double const & TotalPOT() const
Definition: POTSummary.h:45
gar::anatree::fInteractionType
std::vector< Int_t > fInteractionType
Definition: anatree_module.cc:212
gar::anatree::fVeeMLppi
std::vector< Float_t > fVeeMLppi
Definition: anatree_module.cc:402
gar::anatree::fTrackEndY
std::vector< Float_t > fTrackEndY
Definition: anatree_module.cc:341
gar::anatree::fGPartLastMom
std::vector< Int_t > fGPartLastMom
Definition: anatree_module.cc:239
gar::anatree::fMCPEndPY
std::vector< Float_t > fMCPEndPY
Definition: anatree_module.cc:266
gar::anatree::fTPCClusterSig
std::vector< Float_t > fTPCClusterSig
Definition: anatree_module.cc:318
gar::anatree::fVeePXLppi
std::vector< Float_t > fVeePXLppi
Definition: anatree_module.cc:398
gar::anatree::fTrackPIDProbF
std::vector< Float_t > fTrackPIDProbF
Definition: anatree_module.cc:357
gar::anatree::fTrackIDNumber
std::vector< ULong64_t > fTrackIDNumber
Definition: anatree_module.cc:331
gar::anatree::fTrackStartX
std::vector< Float_t > fTrackStartX
Definition: anatree_module.cc:332
gar::anatree::FillGeneratorMonteCarloInfo
void FillGeneratorMonteCarloInfo(art::Event const &e)
Definition: anatree_module.cc:1429
gar::anatree::fTrackMCindex
std::vector< Int_t > fTrackMCindex
Definition: anatree_module.cc:360
gar::anatree::fMCMotherIndex
std::vector< Int_t > fMCMotherIndex
Definition: anatree_module.cc:252
gar::rec::CaloHit::getIDNumber
gar::rec::IDNumber getIDNumber() const
Definition: CaloHit.cxx:40
simb::MCNeutrino::InteractionType
int InteractionType() const
Definition: MCNeutrino.h:150
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
gar::anatree::fMCPEndPZ
std::vector< Float_t > fMCPEndPZ
Definition: anatree_module.cc:267
gar::anatree::fTrackTrajectoryFWDID
std::vector< Int_t > fTrackTrajectoryFWDID
Definition: anatree_module.cc:367
gar::cheat::BackTrackerCore::ClusterToMCParticles
std::vector< std::pair< simb::MCParticle *, float > > ClusterToMCParticles(rec::Cluster *const c)
Definition: BackTrackerCore.cxx:933
gar::anatree::fGPartLastDaugh
std::vector< Int_t > fGPartLastDaugh
Definition: anatree_module.cc:241
gar::anatree::fMCPEndZ
std::vector< Float_t > fMCPEndZ
Definition: anatree_module.cc:264
gar::anatree::fSimHitX_MuID
std::vector< Float_t > fSimHitX_MuID
Definition: anatree_module.cc:297
gar::anatree::fVeeELpip
std::vector< Float_t > fVeeELpip
Definition: anatree_module.cc:406
gar::anatree::fTrackEndZ
std::vector< Float_t > fTrackEndZ
Definition: anatree_module.cc:342
gar::anatree::fSimHitEnergy
std::vector< Float_t > fSimHitEnergy
Definition: anatree_module.cc:290
gar::anatree::fGPartPy
std::vector< Float_t > fGPartPy
Definition: anatree_module.cc:244
gar::anatree::fVeeTAssn_TrackIDNumber
std::vector< ULong64_t > fVeeTAssn_TrackIDNumber
Definition: anatree_module.cc:409
e
const double e
Definition: gUpMuFluxGen.cxx:165
gar::anatree::fGeneratorLabels
std::vector< std::string > fGeneratorLabels
Definition: anatree_module.cc:126
gar::anatree::fRecoHitX
std::vector< Float_t > fRecoHitX
Definition: anatree_module.cc:433
simb::MCNeutrino::W
double W() const
Definition: MCNeutrino.h:154
gar::anatree::fTrackEndPX
std::vector< Float_t > fTrackEndPX
Definition: anatree_module.cc:343
gar::anatree::fRecoHitLayer_MuID
std::vector< Int_t > fRecoHitLayer_MuID
Definition: anatree_module.cc:451
gar::rec::Track::ChargeEnd
int ChargeEnd() const
Definition: Track.cxx:236
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
infile
string infile
Definition: INukeNucleonCorr.cxx:62
gar::anatree::fTrackAvgIonF
std::vector< Float_t > fTrackAvgIonF
Definition: anatree_module.cc:353
gar::anatree::fRecoEnergySum_MuID
Float_t fRecoEnergySum_MuID
Definition: anatree_module.cc:452
gar::anatree::fRecoHitLayer
std::vector< Int_t > fRecoHitLayer
Definition: anatree_module.cc:439
gar::anatree::fVeePYLppi
std::vector< Float_t > fVeePYLppi
Definition: anatree_module.cc:399
gar::anatree::fVeePYLpip
std::vector< Float_t > fVeePYLpip
Definition: anatree_module.cc:404
gar::geo::GeometryCore::TPCRadius
float TPCRadius() const
Returns the radius of the TPC (y or z direction)
Definition: GeometryCore.h:755
gar::anatree::fTrackTrajectoryBWDX
std::vector< Float_t > fTrackTrajectoryBWDX
Definition: anatree_module.cc:369
simb::MCNeutrino::Y
double Y() const
Definition: MCNeutrino.h:156
gar::anatree::fGPartStatus
std::vector< Int_t > fGPartStatus
Definition: anatree_module.cc:237
gar::anatree::fDigiHitCellID
std::vector< ULong64_t > fDigiHitCellID
Definition: anatree_module.cc:419
gar::anatree::fFieldDecoder_ECAL
gar::geo::BitFieldCoder * fFieldDecoder_ECAL
Definition: anatree_module.cc:184
gar::anatree::fnGPart
std::vector< Int_t > fnGPart
Definition: anatree_module.cc:233
gar::anatree::fSimHitTime
std::vector< Float_t > fSimHitTime
Definition: anatree_module.cc:289
gar::anatree::fTPCClusterCovXZ
std::vector< Float_t > fTPCClusterCovXZ
Definition: anatree_module.cc:323
gar::anatree::fHitX
std::vector< Float_t > fHitX
Definition: anatree_module.cc:307
submit_ppfxjobs.temp
temp
Definition: submit_ppfxjobs.py:110
gar::geo::BitFieldCoder
Helper class for decoding and encoding a bit field of 64bits for convenient declaration.
Definition: BitFieldCoder.h:130
a
const double a
Definition: gUpMuFluxGen.cxx:164
gar::anatree::fClusterPID_MuID
std::vector< Float_t > fClusterPID_MuID
Definition: anatree_module.cc:486
reco_momentum_tuples.t
t
Definition: reco_momentum_tuples.py:25
art::DataViewImpl::getMany
std::vector< Handle< PROD > > getMany(SelectorBase const &selector=MatchAllSelector{}) const
Definition: DataViewImpl.h:479
gar::rec::Track::ChargeBeg
int ChargeBeg() const
Definition: Track.cxx:229
gar::anatree::fGeantLabel
std::string fGeantLabel
module label for geant4 simulated hits
Definition: anatree_module.cc:133
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:271
gar::anatree::fClusterEnergy_MuID
std::vector< Float_t > fClusterEnergy_MuID
Definition: anatree_module.cc:478
gar::anatree::fVeeLabel
std::string fVeeLabel
module label for conversion/decay vertexes rec:Vee
Definition: anatree_module.cc:140
gar::anatree::fDigiHitCellID_MuID
std::vector< ULong64_t > fDigiHitCellID_MuID
Definition: anatree_module.cc:428
gar::anatree::fClusterMCindex
std::vector< Int_t > fClusterMCindex
Definition: anatree_module.cc:471
gar::anatree::fTrackAvgIonB
std::vector< Float_t > fTrackAvgIonB
Definition: anatree_module.cc:354
gar::anatree::fCALAssn_TrackIDNumber
std::vector< ULong64_t > fCALAssn_TrackIDNumber
Definition: anatree_module.cc:502
gar::anatree::fReconHits_MuID
UInt_t fReconHits_MuID
Definition: anatree_module.cc:443
gar::anatree::fRecoEnergySum
Float_t fRecoEnergySum
Definition: anatree_module.cc:440
gar::anatree::operator=
anatree & operator=(anatree const &)=delete
test.p
p
Definition: test.py:223
gar::anatree::fClusterIDNumber_MuID
std::vector< ULong64_t > fClusterIDNumber_MuID
Definition: anatree_module.cc:476
gar::geo::GeometryCore::TPCZCent
float TPCZCent() const
Returns the Z location of the center of the TPC in cm.
Definition: GeometryCore.h:792
gar::anatree::fClusterIDNumber
std::vector< ULong64_t > fClusterIDNumber
Definition: anatree_module.cc:456
gar::anatree::fClusterTime_MuID
std::vector< Float_t > fClusterTime_MuID
Definition: anatree_module.cc:479
art::DataViewImpl::subRun
SubRunNumber_t subRun() const
Definition: DataViewImpl.cc:78
gar::anatree::fClusterAssn_RecoHitIDNumber
std::vector< std::vector< ULong64_t > > fClusterAssn_RecoHitIDNumber
Definition: anatree_module.cc:495
simb::MCNeutrino::X
double X() const
Definition: MCNeutrino.h:155
gar::anatree::fMuIDHitLabel
std::string fMuIDHitLabel
Definition: anatree_module.cc:146
gar::anatree::fMCPStartPY
std::vector< Float_t > fMCPStartPY
Definition: anatree_module.cc:260
code
CodeOutputInterface * code
Definition: tclscanner.cpp:1114
gar::anatree::fRawCaloHitLabel
std::string fRawCaloHitLabel
module label for digitized calo hits raw::CaloRawDigit
Definition: anatree_module.cc:142
gar::anatree::fPFLabel
std::string fPFLabel
module label for reco particles rec::PFParticle
Definition: anatree_module.cc:150
gar::cheat::BackTrackerCore::HitToHitIDEs
std::vector< HitIDE > HitToHitIDEs(art::Ptr< rec::Hit > const &hit) const
Definition: BackTrackerCore.cxx:225
CaloRawDigit.h
Definition of basic calo raw digits.
protoana::distance
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
Definition: truepionXsection_module.cc:217
gar::anatree::fGPartName
std::vector< std::string > fGPartName
Definition: anatree_module.cc:242
gar::anatree::fTrackStartY
std::vector< Float_t > fTrackStartY
Definition: anatree_module.cc:333
art::DataViewImpl::run
RunNumber_t run() const
Definition: DataViewImpl.cc:71
MF_LOG_INFO
#define MF_LOG_INFO(category)
Definition: MessageLogger.h:225
gar::anatree::fSimEnergySum
Float_t fSimEnergySum
Definition: anatree_module.cc:293
gar::anatree::fSimHitZ
std::vector< Float_t > fSimHitZ
Definition: anatree_module.cc:288
gar::anatree::fClusterMCfrac_MuID
std::vector< Float_t > fClusterMCfrac_MuID
Definition: anatree_module.cc:492
gar::anatree::fTPCClusterX
std::vector< Float_t > fTPCClusterX
Definition: anatree_module.cc:315
gar::anatree::fMode
std::vector< Int_t > fMode
Definition: anatree_module.cc:211
gar::anatree::fMCPDG
std::vector< Int_t > fMCPDG
Definition: anatree_module.cc:251
gar::anatree::fClusterTimeDiffFirstLast_MuID
std::vector< Float_t > fClusterTimeDiffFirstLast_MuID
Definition: anatree_module.cc:480
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
gar::anatree::fClusterMuIDLabel
std::string fClusterMuIDLabel
module label for calo clusters rec::Cluster in MuID
Definition: anatree_module.cc:149
gar::anatree::fRecoHitEnergy_MuID
std::vector< Float_t > fRecoHitEnergy_MuID
Definition: anatree_module.cc:449
gar::anatree::fMuIDEncoding
std::string fMuIDEncoding
Definition: anatree_module.cc:183
gar::anatree::fTrackTrajectoryFWDY
std::vector< Float_t > fTrackTrajectoryFWDY
Definition: anatree_module.cc:365
simb::MCTruth::GetParticle
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
gar::anatree::fTrajMCPPX
std::vector< Float_t > fTrajMCPPX
Definition: anatree_module.cc:280
gar::anatree::fVeeT
std::vector< ULong64_t > fVeeT
Definition: anatree_module.cc:392
gar::anatree::fVeePYKpipi
std::vector< Float_t > fVeePYKpipi
Definition: anatree_module.cc:394
gar::anatree::fVertexIDNumber
std::vector< ULong64_t > fVertexIDNumber
Definition: anatree_module.cc:375
gar::anatree::fTrackEndQ
std::vector< Int_t > fTrackEndQ
Definition: anatree_module.cc:346
gar::anatree::fDiginHits_MuID
UInt_t fDiginHits_MuID
Definition: anatree_module.cc:422
gar::anatree::fWriteCaloDigits
bool fWriteCaloDigits
Raw digits for calorimetry. Default=false.
Definition: anatree_module.cc:168
gar::anatree::fSimHitY
std::vector< Float_t > fSimHitY
Definition: anatree_module.cc:287
gar::anatree::fVeeMKpipi
std::vector< Float_t > fVeeMKpipi
Definition: anatree_module.cc:397
gar::anatree::fReconHitIDNumber
std::vector< ULong64_t > fReconHitIDNumber
Definition: anatree_module.cc:432
gar::anatree::fEngine
CLHEP::HepRandomEngine & fEngine
random engine
Definition: anatree_module.cc:506
gar
General GArSoft Utilities.
Definition: anatest_module.cc:57
gar::anatree::fVTAssn_TrackEnd
std::vector< gar::rec::TrackEnd > fVTAssn_TrackEnd
Definition: anatree_module.cc:385
gar::anatree::fMCTrkID
std::vector< Int_t > fMCTrkID
Definition: anatree_module.cc:250
gar::geo::BitFieldCoder::get
long64 get(long64 bitfield, size_t index) const
Definition: BitFieldCoder.h:181
gar::anatree::fVTAssn_VertIDNumber
std::vector< ULong64_t > fVTAssn_VertIDNumber
Definition: anatree_module.cc:383
gar::anatree::fPOTtag
std::string fPOTtag
Definition: anatree_module.cc:128
gar::anatree::fRecoHitY
std::vector< Float_t > fRecoHitY
Definition: anatree_module.cc:434
gar::rec::TrackEndBeg
TrackEnd const TrackEndBeg
Definition: Track.h:33
gar::anatree::fVeeEKpipi
std::vector< Float_t > fVeeEKpipi
Definition: anatree_module.cc:396
gar::anatree::fClusterMainAxisZ
std::vector< Float_t > fClusterMainAxisZ
Definition: anatree_module.cc:470
gar::anatree::fClusterMCfrac
std::vector< Float_t > fClusterMCfrac
Definition: anatree_module.cc:472
gar::anatree::fECALAssnLabel
std::string fECALAssnLabel
module label for track-clusters associations
Definition: anatree_module.cc:151
gar::anatree::fgT
std::vector< Float_t > fgT
Definition: anatree_module.cc:230
gar::anatree::fTrajMCPE
std::vector< Float_t > fTrajMCPE
Definition: anatree_module.cc:277
gar::anatree::fTrajMCPX
std::vector< Float_t > fTrajMCPX
Definition: anatree_module.cc:273
gar::anatree::fDigiHitY
std::vector< Float_t > fDigiHitY
Definition: anatree_module.cc:415
gar::anatree::BackTrack
cheat::BackTrackerCore * BackTrack
Definition: anatree_module.cc:99
gar::anatree::fTPCClusterCovYY
std::vector< Float_t > fTPCClusterCovYY
Definition: anatree_module.cc:324
gar::anatree::fVertexT
std::vector< ULong64_t > fVertexT
Definition: anatree_module.cc:379
gar::anatree::fClusterEnergy
std::vector< Float_t > fClusterEnergy
Definition: anatree_module.cc:458
gar::anatree::fRecoHitEnergy
std::vector< Float_t > fRecoHitEnergy
Definition: anatree_module.cc:437
gar::anatree::fCaloHitLabel
std::string fCaloHitLabel
module label for reco calo hits rec::CaloHit
Definition: anatree_module.cc:145
gar::anatree::fWriteMCPTrajectory
bool fWriteMCPTrajectory
Write MCP Trajectory Default=true.
Definition: anatree_module.cc:155
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:563
gar::anatree::fMCPStartZ
std::vector< Float_t > fMCPStartZ
Definition: anatree_module.cc:257
gar::geo::GeometryCore::TPCYCent
float TPCYCent() const
Returns the Y location of the center of the TPC in cm.
Definition: GeometryCore.h:785
gar::anatree::fW
std::vector< Float_t > fW
Definition: anatree_module.cc:214
gar::anatree::fVeeELppi
std::vector< Float_t > fVeeELppi
Definition: anatree_module.cc:401
gar::anatree::fSimHitZ_MuID
std::vector< Float_t > fSimHitZ_MuID
Definition: anatree_module.cc:299
gar::anatree::fSubRun
Int_t fSubRun
number of the sub-run being processed
Definition: anatree_module.cc:194
gar::anatree::fTrackPIDB
std::vector< Int_t > fTrackPIDB
Definition: anatree_module.cc:358
gar::anatree::fRecoHitZ_MuID
std::vector< Float_t > fRecoHitZ_MuID
Definition: anatree_module.cc:447
gar::anatree::fTrackStartPY
std::vector< Float_t > fTrackStartPY
Definition: anatree_module.cc:336
gar::anatree::fClusterTime
std::vector< Float_t > fClusterTime
Definition: anatree_module.cc:459
gar::cheat::BackTrackerCore::TrackToMCParticles
std::vector< std::pair< simb::MCParticle *, float > > TrackToMCParticles(rec::Track *const t)
Definition: BackTrackerCore.cxx:823
gar::anatree::fTPCClusterCovZZ
std::vector< Float_t > fTPCClusterCovZZ
Definition: anatree_module.cc:326
fhicl::ParameterSet::get_if_present
std::optional< T > get_if_present(std::string const &key) const
Definition: ParameterSet.h:224
gar::anatree::fHitChan
std::vector< UInt_t > fHitChan
Definition: anatree_module.cc:312
gar::anatree::fTrackPIDProbB
std::vector< Float_t > fTrackPIDProbB
Definition: anatree_module.cc:359
gar::anatree::processOneDirection
float processOneDirection(std::vector< std::pair< float, float >> SigData, float ionizeTruncate)
Definition: anatree_module.cc:2386
gar::geo::GeometryCore::HasMuonDetector
bool HasMuonDetector() const
Definition: GeometryCore.h:999
gar::anatree::fClusterMainAxisY
std::vector< Float_t > fClusterMainAxisY
Definition: anatree_module.cc:469
gar::geo::GeometryCore::GetECALCellIDEncoding
const std::string GetECALCellIDEncoding() const
Definition: GeometryCore.cxx:838
MF_LOG_DEBUG
#define MF_LOG_DEBUG(id)
Definition: MessageLogger.h:268
gar::anatree::endRun
virtual void endRun(art::Run const &run) override
Definition: anatree_module.cc:1056
simb::MCParticle::Pz
double Pz(const int i=0) const
Definition: MCParticle.h:232
gar::anatree::fMCVertexZ
std::vector< Float_t > fMCVertexZ
Definition: anatree_module.cc:221
gar::anatree::fGPartMass
std::vector< Float_t > fGPartMass
Definition: anatree_module.cc:247
gar::anatree::fClusterPID
std::vector< Float_t > fClusterPID
Definition: anatree_module.cc:466
gar::anatree::fTPCClusterMCfrac
std::vector< Float_t > fTPCClusterMCfrac
Definition: anatree_module.cc:328
gar::anatree::fMCPEndPX
std::vector< Float_t > fMCPEndPX
Definition: anatree_module.cc:265
gar::anatree::fDigiHitTime_MuID
std::vector< Float_t > fDigiHitTime_MuID
Definition: anatree_module.cc:426
gar::anatree::fTPCClusterRMS
std::vector< Float_t > fTPCClusterRMS
Definition: anatree_module.cc:319
gar::anatree::fTgtPDG
std::vector< Int_t > fTgtPDG
Definition: anatree_module.cc:228
gar::anatree::fVeeX
std::vector< Float_t > fVeeX
Definition: anatree_module.cc:389
gar::anatree::fTrajMCPY
std::vector< Float_t > fTrajMCPY
Definition: anatree_module.cc:274
b
static bool * b
Definition: config.cpp:1043
gar::anatree::fVeePXKpipi
std::vector< Float_t > fVeePXKpipi
Definition: anatree_module.cc:393
gar::anatree::fWriteTracks
bool fWriteTracks
Start/end X, P for tracks Default=true.
Definition: anatree_module.cc:162
gar::anatree::fVeeZ
std::vector< Float_t > fVeeZ
Definition: anatree_module.cc:391
gar::anatree::TrackIdToIndex
std::unordered_map< TrkId, Int_t > TrackIdToIndex
Definition: anatree_module.cc:188
art::EventID::event
EventNumber_t event() const
Definition: EventID.h:116
gar::anatree::fDigiHitTime
std::vector< Float_t > fDigiHitTime
Definition: anatree_module.cc:417
gar::anatree::fMCPTime
std::vector< Float_t > fMCPTime
Definition: anatree_module.cc:258
gar::anatree::fSimHitEnergy_MuID
std::vector< Float_t > fSimHitEnergy_MuID
Definition: anatree_module.cc:301
gar::anatree::fGPartPx
std::vector< Float_t > fGPartPx
Definition: anatree_module.cc:243
gar::anatree::fTrajMCPZ
std::vector< Float_t > fTrajMCPZ
Definition: anatree_module.cc:275
gar::anatree::fInstanceLabelCalo
std::string fInstanceLabelCalo
Instance name for ECAL.
Definition: anatree_module.cc:130
gar::anatree::fTrackChi2B
std::vector< Float_t > fTrackChi2B
Definition: anatree_module.cc:351
gar::anatree::fDigiHitX
std::vector< Float_t > fDigiHitX
Definition: anatree_module.cc:414
gar::anatree::fHitSig
std::vector< Float_t > fHitSig
Definition: anatree_module.cc:310
gar::geo::GeometryCore::TPCLength
float TPCLength() const
Returns the length of the TPC (x direction)
Definition: GeometryCore.h:763
gar::anatree::fWriteMatchedTracks
bool fWriteMatchedTracks
Write ECAL-track Assns Default=true.
Definition: anatree_module.cc:171
gar::anatree::fFieldDecoder_MuID
gar::geo::BitFieldCoder * fFieldDecoder_MuID
Definition: anatree_module.cc:185
gar::rec::Track::getIDNumber
gar::rec::IDNumber getIDNumber() const
Definition: Track.cxx:42
gar::anatree::fWriteMCPTrajMomenta
bool fWriteMCPTrajMomenta
Write Momenta associated with MCP Trajectory Default=false.
Definition: anatree_module.cc:156
gar::anatree::fMCnuPy
std::vector< Float_t > fMCnuPy
Definition: anatree_module.cc:223
wirecell.units.pi
float pi
Definition: units.py:11
gar::anatree::fClusterMainAxisY_MuID
std::vector< Float_t > fClusterMainAxisY_MuID
Definition: anatree_module.cc:489
gar::anatree::fMCPEndY
std::vector< Float_t > fMCPEndY
Definition: anatree_module.cc:263
gar::anatree::fSimEnergySum_MuID
Float_t fSimEnergySum_MuID
Definition: anatree_module.cc:304
gar::anatree::fCALAssn_TrackEnd
std::vector< gar::rec::TrackEnd > fCALAssn_TrackEnd
Definition: anatree_module.cc:503
gar::anatree::fHitLabel
std::string fHitLabel
module label for reco TPC hits rec::Hit
Definition: anatree_module.cc:135
simb::MCTruth
Event generator information.
Definition: MCTruth.h:32
gar::anatree::fTPCClusterCovXY
std::vector< Float_t > fTPCClusterCovXY
Definition: anatree_module.cc:322
gar::anatree::fTrackStartZ
std::vector< Float_t > fTrackStartZ
Definition: anatree_module.cc:334
gar::anatree::fRecoHitCellID_MuID
std::vector< ULong64_t > fRecoHitCellID_MuID
Definition: anatree_module.cc:450
gar::rec::TrackIoniz::getBAK_dSigdXs
const std::vector< std::pair< float, float > > getBAK_dSigdXs() const
Definition: TrackIoniz.h:23
reco_momentum_tuples.momentum
def momentum(x1, x2, x3, scale=1.)
Definition: reco_momentum_tuples.py:198
gar::anatree::fWriteCohInfo
bool fWriteCohInfo
MC level t for coherent pi+. Default=false.
Definition: anatree_module.cc:175
simb::MCParticle::EndX
double EndX() const
Definition: MCParticle.h:226
gar::anatree::fSimnHits_MuID
UInt_t fSimnHits_MuID
Definition: anatree_module.cc:296
gar::anatree::fMCPEndProc
std::vector< std::string > fMCPEndProc
Definition: anatree_module.cc:269
gar::anatree::fTrackStartQ
std::vector< Int_t > fTrackStartQ
Definition: anatree_module.cc:338
gar::anatree::fEvent
Int_t fEvent
number of the event being processed
Definition: anatree_module.cc:192
gar::anatree::fTrackMCfrac
std::vector< Float_t > fTrackMCfrac
Definition: anatree_module.cc:361
gar::anatree::fTrajMCPIndex
std::vector< Int_t > fTrajMCPIndex
Definition: anatree_module.cc:278
gar::anatree::fTrackTrajectoryBWDY
std::vector< Float_t > fTrackTrajectoryBWDY
Definition: anatree_module.cc:370
util::quantities::second
second_as<> second
Type of time stored in seconds, in double precision.
Definition: spacetime.h:85
gar::anatree::FillRecoInfo
void FillRecoInfo(art::Event const &e)
Definition: anatree_module.cc:1792
gar::anatree::fVeePZLpip
std::vector< Float_t > fVeePZLpip
Definition: anatree_module.cc:405
gar::rec::TrackIoniz::getFWD_dSigdXs
const std::vector< std::pair< float, float > > getFWD_dSigdXs() const
Definition: TrackIoniz.h:22
simb::MCNeutrino
Event generator information.
Definition: MCNeutrino.h:18
gar::anatree::fTrackTrajectoryFWDX
std::vector< Float_t > fTrackTrajectoryFWDX
Definition: anatree_module.cc:364
GeometryGAr.h
art framework interface to geometry description
gar::sumdata::POTSummary::TotalSpills
unsigned int const & TotalSpills() const
Definition: POTSummary.h:46
uint
unsigned uint
Definition: qglobal.h:351
gar::anatree::fCCNC
std::vector< Int_t > fCCNC
Definition: anatree_module.cc:210
gar::anatree::fGPartPz
std::vector< Float_t > fGPartPz
Definition: anatree_module.cc:245
simb::MCNeutrino::Mode
int Mode() const
Definition: MCNeutrino.h:149
gar::anatree::fRecoHitX_MuID
std::vector< Float_t > fRecoHitX_MuID
Definition: anatree_module.cc:445
gar::anatree::fTrackPIDF
std::vector< Int_t > fTrackPIDF
Definition: anatree_module.cc:356
gar::anatree::fTrackStartPX
std::vector< Float_t > fTrackStartPX
Definition: anatree_module.cc:335
gar::anatree::fGPartE
std::vector< Float_t > fGPartE
Definition: anatree_module.cc:246
gar::anatree::fGPartIdx
std::vector< Int_t > fGPartIdx
Definition: anatree_module.cc:235
gar::anatree::fTrackLabel
std::string fTrackLabel
module label for TPC Tracks rec:Track
Definition: anatree_module.cc:137
gar::anatree::fClusterX_MuID
std::vector< Float_t > fClusterX_MuID
Definition: anatree_module.cc:481
gar::anatree::fVertexLabel
std::string fVertexLabel
module label for vertexes rec:Vertex
Definition: anatree_module.cc:139
gar::anatree::m_pidinterp
std::unordered_map< int, TH2F * > m_pidinterp
Definition: anatree_module.cc:507
gar::anatree::fInstanceLabelMuID
std::string fInstanceLabelMuID
Instance name for MuID.
Definition: anatree_module.cc:131
gar::anatree::fCALAssn_ClusIDNumber
std::vector< ULong64_t > fCALAssn_ClusIDNumber
Definition: anatree_module.cc:501
gar::anatree::processIonizationInfo
void processIonizationInfo(rec::TrackIoniz &ion, float ionizeTruncate, float &forwardIonVal, float &backwardIonVal)
Definition: anatree_module.cc:2370
gar::anatree::fTPCClusterZ
std::vector< Float_t > fTPCClusterZ
Definition: anatree_module.cc:317
gar::anatree::fHitRMS
std::vector< Float_t > fHitRMS
Definition: anatree_module.cc:311
gar::anatree::fMCPEndX
std::vector< Float_t > fMCPEndX
Definition: anatree_module.cc:262
s
static QCString * s
Definition: config.cpp:1042
gar::anatree::fWriteCaloClusters
bool fWriteCaloClusters
Write ECAL clusters. Default=true.
Definition: anatree_module.cc:170
gar::anatree::fMCVertexY
std::vector< Float_t > fMCVertexY
Definition: anatree_module.cc:220
gar::anatree::fVertexQ
std::vector< Int_t > fVertexQ
Definition: anatree_module.cc:381
art::Event::id
EventID id() const
Definition: Event.cc:34
gar::anatree::fClusterMainAxisX
std::vector< Float_t > fClusterMainAxisX
Definition: anatree_module.cc:468
str
static QCString str
Definition: fortrancode.cpp:27098
gar::anatree::fTrackTrajectoryBWDID
std::vector< Int_t > fTrackTrajectoryBWDID
Definition: anatree_module.cc:372
gar::anatree::fTPCClusterCovXX
std::vector< Float_t > fTPCClusterCovXX
Definition: anatree_module.cc:321
gar::anatree::fDigiHitY_MuID
std::vector< Float_t > fDigiHitY_MuID
Definition: anatree_module.cc:424
gar::anatree::fRawMuIDHitLabel
std::string fRawMuIDHitLabel
Definition: anatree_module.cc:143
gar::anatree::fClusterMainAxisX_MuID
std::vector< Float_t > fClusterMainAxisX_MuID
Definition: anatree_module.cc:488
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
gar::anatree::fRecoHitTime_MuID
std::vector< Float_t > fRecoHitTime_MuID
Definition: anatree_module.cc:448
gar::anatree::fRecoHitTime
std::vector< Float_t > fRecoHitTime
Definition: anatree_module.cc:436
gar::anatree::fHitY
std::vector< Float_t > fHitY
Definition: anatree_module.cc:308
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
QCString::append
QCString & append(const char *s)
Definition: qcstring.cpp:383
gar::anatree::fMCVertexX
std::vector< Float_t > fMCVertexX
Definition: anatree_module.cc:219
gar::anatree::fSimHitCellID_MuID
std::vector< ULong64_t > fSimHitCellID_MuID
Definition: anatree_module.cc:303
gar::anatree::fClusterNhits
std::vector< UInt_t > fClusterNhits
Definition: anatree_module.cc:457
gar::anatree::fMCMotherTrkID
std::vector< Int_t > fMCMotherTrkID
Definition: anatree_module.cc:253
gar::anatree::fTheta
std::vector< Float_t > fTheta
Definition: anatree_module.cc:217
gar::anatree::fTrackLenF
std::vector< Float_t > fTrackLenF
Definition: anatree_module.cc:348
gar::anatree::fVeeMLpip
std::vector< Float_t > fVeeMLpip
Definition: anatree_module.cc:407
gar::anatree::fECALEncoding
std::string fECALEncoding
Definition: anatree_module.cc:182
gar::anatree::fClusterTheta_MuID
std::vector< Float_t > fClusterTheta_MuID
Definition: anatree_module.cc:484
gar::anatree::fTrajMCPPY
std::vector< Float_t > fTrajMCPPY
Definition: anatree_module.cc:281
gar::anatree::fGPartIntIdx
std::vector< Int_t > fGPartIntIdx
Definition: anatree_module.cc:234
gar::anatree::analyze
void analyze(art::Event const &e) override
Definition: anatree_module.cc:1077
gar::anatree::fClusterY
std::vector< Float_t > fClusterY
Definition: anatree_module.cc:462
vertex
vertex reconstruction
Definition: AggregateVertexAna_module.cc:32
Generated by   doxygen 1.8.11