HNLAnalysis_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 // Module for heavy neutral lepton analysis
3 //
4 // File: HNLAnalysis_module.cc
5 //
6 // Original version by Georgios Christodoulou, November 2020
7 //
8 // Some code is borrowed from anatree_module.cc
9 //
10 ////////////////////////////////////////////////////////////////////////////////
11 
12 #include "art/Framework/Core/EDAnalyzer.h"
13 #include "art/Framework/Core/ModuleMacros.h"
14 #include "art/Framework/Principal/Event.h"
15 #include "art/Framework/Principal/Handle.h"
16 #include "art/Framework/Principal/Run.h"
17 #include "art/Framework/Principal/SubRun.h"
18 #include "art_root_io/TFileService.h"
19 #include "canvas/Utilities/InputTag.h"
20 #include "fhiclcpp/ParameterSet.h"
21 #include "messagefacility/MessageLogger/MessageLogger.h"
22 #include "art/Persistency/Common/PtrMaker.h"
23 #include "canvas/Persistency/Common/FindOne.h"
24 #include "canvas/Persistency/Common/FindOneP.h"
25 #include "canvas/Persistency/Common/FindMany.h"
26 #include "canvas/Persistency/Common/FindManyP.h"
27 
28 #include "nusimdata/SimulationBase/MCTruth.h"
29 #include "nusimdata/SimulationBase/MCParticle.h"
30 #include "MCCheater/BackTracker.h"
31 
32 #include "ReconstructionDataProducts/Track.h"
33 #include "ReconstructionDataProducts/Cluster.h"
34 #include "ReconstructionDataProducts/TrackIoniz.h"
35 
36 #include "CoreUtils/ServiceUtil.h"
37 #include "Geometry/GeometryGAr.h"
38 
39 // nutools extensions
40 #include "nurandom/RandomUtils/NuRandomService.h"
41 
42 #include "TTree.h"
43 #include "TH2.h"
44 #include "TFile.h"
45 #include "TDatabasePDG.h"
46 #include "TParticlePDG.h"
47 #include "TRandom3.h"
48 
49 #include "CLHEP/Random/RandFlat.h"
50 
51 #include <string>
52 #include <vector>
53 #include <unordered_map>
54 
55 namespace gar {
56 
57  typedef std::pair<float, std::string> P;
58 
59  class HNLAnalysis : public art::EDAnalyzer {
60  public:
61  explicit HNLAnalysis(fhicl::ParameterSet const & p);
62  // The compiler-generated destructor is fine for non-base
63  // classes without bare pointers or other resource use.
64 
65  // Plugins should not be copied or assigned.
66  HNLAnalysis(HNLAnalysis const &) = delete;
67  HNLAnalysis(HNLAnalysis &&) = delete;
68  HNLAnalysis & operator = (HNLAnalysis const &) = delete;
69  HNLAnalysis & operator = (HNLAnalysis &&) = delete;
70 
71  virtual void beginJob() override;
72 
73  // Required functions.
74  void analyze(art::Event const & e) override;
75 
76  private:
77  void ClearVectors();
78  void FillVectors(art::Event const & e);
79 
80  // Helpers copied from anatree_module
81  void processIonizationInfo(rec::TrackIoniz& ion, float ionizeTruncate, float& forwardIonVal, float& backwardIonVal);
82  float processOneDirection(std::vector<std::pair<float,float>> SigData, float ionizeTruncate);
83  // Helper method for processOneDirection
84  static bool lessThan_byE(std::pair<float,float> a, std::pair<float,float> b) {return a.first < b.first;}
85 
86  // Get parameterized TPC PID
87  std::vector< std::pair<int, float> > processPIDInfo(float p);
88 
89  // Get the parameterized ECal PID - this from CAF.cpp
90  int GetECalPIDParameterized(int pdg, float ptrue);
91 
92  // Position of TPC from geometry service
93  double ftpccentre[3];
94  // TPC length
95  double ftpclength;
96  // TPC radius
97  double ftpcradius;
98  // ECal geometry
99  double fECALBarrelInnerRadius;
100  double fECALBarrelOuterRadius;
101  double fECALEndcapInnerRadius;
102  double fECALEndcapOuterRadius;
103  double fECALStartX;
104  double fECALEndX;
105 
106  // Input data labels
107  std::string fGeantLabel;
108  //std::string fGeantInstanceCalo; ///< Instance name ECAL for sdp::CaloDeposit
109  //std::string fGeantInstanceMuID; ///< Instance name MuID for sdp::CaloDeposit
110 
111  std::string fTrackLabel;
112  std::string fClusterLabel;
113  std::string fECALAssnLabel;
114 
115  // Truncation parameter for dE/dx (average this fraction lowest readings)
116  float fIonizTruncate; ///< Default=1.00;
117  // MCParticle to MC vertex match
118  float fMatchMCPtoVertDist; ///< Default=roundoff
119 
120  // The analysis tree
121  TTree *fTree;
122  // The neutrino truth tree
123  TTree *fTruthTree;
124  // Configuration tree
125  TTree *fconfig;
126 
127  // Geometry service
128  const gar::geo::GeometryCore* fGeo;
129  // Backtracker service
130  cheat::BackTrackerCore* fBack;
131  // PDG database from ROOT
132  TDatabasePDG* pdgInstance;
133 
134  // Map of PID TH2 per momentum value
135  CLHEP::HepRandomEngine &fEngine; ///< random engine
136  std::unordered_map<int, TH2F*> m_pidinterp;
137 
138  // The analysis tree. Use Rtypes.h here, as these data get used by root global event info
139  Int_t fRun;
140  Int_t fSubRun;
141  Int_t fEvent;
142  Float_t fWeight;
143 
144  // Number of TPC tracks
145  Int_t fNTPCTracks;
146  Int_t fNTPCECalTracks;
147 
148  // This analysis starts from reco'd tracks
149  std::vector<ULong64_t> fTrackIDNumber;
150  std::vector<Float_t> fTrackX;
151  std::vector<Float_t> fTrackY;
152  std::vector<Float_t> fTrackZ;
153  std::vector<Float_t> fTrackPX;
154  std::vector<Float_t> fTrackPY;
155  std::vector<Float_t> fTrackPZ;
156  std::vector<Float_t> fTrackP;
157  std::vector<Int_t> fTrackQ;
158  std::vector<Float_t> fTrackLen;
159  std::vector<Float_t> fTrackChi2;
160  std::vector<Int_t> fNTPCClustersOnTrack;
161  std::vector<Double_t> fTrackTime;
162  std::vector<Float_t> fTrackAvgIon;
163  std::vector<Float_t> fTrackEndX;
164  std::vector<Float_t> fTrackEndY;
165  std::vector<Float_t> fTrackEndZ;
166  std::vector<Float_t> fTrackEndPX;
167  std::vector<Float_t> fTrackEndPY;
168  std::vector<Float_t> fTrackEndPZ;
169  std::vector<Float_t> fTrackEndP;
170  std::vector<Int_t> fTrackEndQ;
171  std::vector<Int_t> fTrackPID;
172  // Which match to qualified MCParticles
173  std::vector<Int_t> fMCPTrackID;
174  std::vector<Int_t> fMCPDG;
175  std::vector<Int_t> fMCPDGMother;
176  std::vector<Float_t> fMCPStartX;
177  std::vector<Float_t> fMCPStartY;
178  std::vector<Float_t> fMCPStartZ;
179  std::vector<Float_t> fMCPStartPX;
180  std::vector<Float_t> fMCPStartPY;
181  std::vector<Float_t> fMCPStartPZ;
182  std::vector<Float_t> fMCPStartP;
183  std::vector<Float_t> fMCPStartEndP;
184  std::vector<std::string> fMCPProc;
185  std::vector<std::string> fMCPEndProc;
186  std::vector<Float_t> fMCPTime;
187  // MC particle enters the ECal
188  std::vector<Int_t> fMCPInECal;
189  // Matched vertex id
190  //std::vector<Int_t> fMCPFromVertexID;
191 
192  // ECal variables
193  std::vector<Int_t> fNECalClustersForTrack;
194  std::vector<Int_t> fClusterTrackIDMatched;
195  std::vector<Int_t> fClusterNhits;
196  std::vector<Float_t> fClusterEnergy;
197  std::vector<Float_t> fClusterTime;
198  std::vector<Float_t> fClusterTimeDiffFirstLast;
199  std::vector<Float_t> fClusterX;
200  std::vector<Float_t> fClusterY;
201  std::vector<Float_t> fClusterZ;
202  std::vector<Float_t> fClusterTheta;
203  std::vector<Float_t> fClusterPhi;
204  std::vector<Int_t> fClusterPID;
205  std::vector<Float_t> fClusterMainAxisX;
206  std::vector<Float_t> fClusterMainAxisY;
207  std::vector<Float_t> fClusterMainAxisZ;
208 
209  // Neutrino truth
210  std::vector<Int_t> f_NeutrinoType;
211  std::vector<Int_t> f_CCNC;
212  std::vector<Int_t> f_Mode;
213  std::vector<Int_t> f_InteractionType;
214  std::vector<Float_t> f_Q2;
215  std::vector<Float_t> f_W;
216  std::vector<Float_t> f_X;
217  std::vector<Float_t> f_Y;
218  std::vector<Float_t> f_Theta;
219  std::vector<Float_t> f_T;
220  std::vector<Float_t> f_MCVertexX;
221  std::vector<Float_t> f_MCVertexY;
222  std::vector<Float_t> f_MCVertexZ;
223  std::vector<Float_t> f_MCnuE;
224  std::vector<Int_t> f_TgtPDG;
225  std::vector<Int_t> f_MCVertexInGasTPC;
226 
227  // MC truth
228  std::vector<Int_t> fNeutrinoType;
229  std::vector<Int_t> fCCNC;
230  std::vector<Int_t> fMode;
231  std::vector<Int_t> fInteractionType;
232  std::vector<Float_t> fQ2;
233  std::vector<Float_t> fW;
234  std::vector<Float_t> fX;
235  std::vector<Float_t> fY;
236  std::vector<Float_t> fTheta;
237  std::vector<Float_t> fT;
238  std::vector<Float_t> fMCVertexX;
239  std::vector<Float_t> fMCVertexY;
240  std::vector<Float_t> fMCVertexZ;
241  std::vector<Float_t> fMCnuPx;
242  std::vector<Float_t> fMCnuPy;
243  std::vector<Float_t> fMCnuPz;
244  std::vector<Float_t> fMCnuE;
245  std::vector<Int_t> fTgtPDG;
246  std::vector<Int_t> fMCVertexInGasTPC;
247  //std::vector<Int_t> fMCVertexID;
248 
249  // MC HNL truth decay products
250  std::vector<Int_t> fMCP_TrackID;
251  std::vector<Int_t> fMCP_PDG;
252  std::vector<Float_t> fMCP_X;
253  std::vector<Float_t> fMCP_Y;
254  std::vector<Float_t> fMCP_Z;
255  std::vector<Float_t> fMCP_T;
256  std::vector<Float_t> fMCP_PX;
257  std::vector<Float_t> fMCP_PY;
258  std::vector<Float_t> fMCP_PZ;
259  std::vector<Float_t> fMCP_E;
260  std::vector<Int_t> fMCP_InGasTPC;
261  //std::vector<Int_t> fMCP_VertexID;
262 
263  };
264 
265 }
266 
267 //==============================================================================
268 // constructor
269 gar::HNLAnalysis::HNLAnalysis(fhicl::ParameterSet const & p) : EDAnalyzer(p), fEngine(art::ServiceHandle<rndm::NuRandomService>()->createEngine(*this, p, "Seed")) {
270 
271  fGeantLabel = p.get<std::string>("GEANTLabel", "geant");
272  //fGeantInstanceCalo = p.get<std::string>("GEANTInstanceCalo","ECAL");
273  //fGeantInstanceMuID = p.get<std::string>("GEANTInstanceMuID","MuID");
274 
275  fTrackLabel = p.get<std::string>("TrackLabel", "track");
276  fClusterLabel = p.get<std::string>("ClusterLabel", "calocluster");
277  fECALAssnLabel = p.get<std::string>("ECALAssnLabel","trkecalassn");
278 
279  fIonizTruncate = p.get<float>("IonizTruncate", 1.00);
280 
281  float MCPtoVertDefault = 10.0*std::numeric_limits<Float_t>::epsilon();
282  fMatchMCPtoVertDist = p.get<float>("MatchMCPtoVertDist",MCPtoVertDefault);
283 
284  pdgInstance = TDatabasePDG::Instance();
285 
286  consumesMany<std::vector<simb::MCTruth> >();
287  consumes<std::vector<simb::MCParticle> >(fGeantLabel);
288  consumes<std::vector<rec::Track> >(fTrackLabel);
289  consumes<std::vector<rec::Cluster> >(fClusterLabel);
290  consumes<art::Assns<rec::Cluster, rec::Track>>(fECALAssnLabel);
291 
292  // Calorimetry related
293  //art::InputTag ecalgeanttag(fGeantLabel, fGeantInstanceCalo);
294  //consumes<std::vector<gar::sdp::CaloDeposit> >(ecalgeanttag);
295 
296  // Muon system related
297  //if (fGeo->HasMuonDetector()) {
298  //art::InputTag muidgeanttag(fGeantLabel, fGeantInstanceMuID);
299  //consumes<std::vector<gar::sdp::CaloDeposit> >(muidgeanttag);
300  //}
301 
302 
303 } // end constructor
304 
305 //==============================================================================
306 void gar::HNLAnalysis::beginJob() {
307 
308  // Backtracker
309  cheat::BackTrackerCore const* const_bt = gar::providerFrom<cheat::BackTracker>();
310  fBack = const_cast<cheat::BackTrackerCore*>(const_bt);
311 
312  // Geometry
313  fGeo = gar::providerFrom<geo::GeometryGAr>();
314  ftpccentre[0] = fGeo->TPCXCent();
315  ftpccentre[1] = fGeo->TPCYCent();
316  ftpccentre[2] = fGeo->TPCZCent();
317  ftpcradius = fGeo->TPCRadius();
318  ftpclength = fGeo->TPCLength();
319  fECALBarrelInnerRadius = fGeo->GetECALInnerBarrelRadius();
320  fECALBarrelOuterRadius = fGeo->GetECALOuterBarrelRadius();
321  fECALEndcapInnerRadius = fGeo->GetECALInnerEndcapRadius();
322  fECALEndcapOuterRadius = fGeo->GetECALOuterEndcapRadius();
323  fECALStartX = fGeo->GetECALEndcapStartX();
324  fECALEndX = fGeo->GetECALEndcapOuterX();
325 
326  art::ServiceHandle<art::TFileService> tfs;
327 
328  // Save TPC/ECal geometry
329  fconfig = tfs->make<TTree>("config","config");
330  fconfig->Branch("tpccentre", &ftpccentre, "tpccentre[3]/D");
331  fconfig->Branch("tpcradius", &ftpcradius, "tpcradius/D");
332  fconfig->Branch("tpclength", &ftpclength, "tpclength/D");
333  fconfig->Branch("ECALBarrelInnerRadius", &fECALBarrelInnerRadius, "ECALBarrelInnerRadius/D");
334  fconfig->Branch("ECALBarrelOuterRadius", &fECALBarrelOuterRadius, "ECALBarrelOuterRadius/D");
335  fconfig->Branch("ECALEndcapInnerRadius", &fECALEndcapInnerRadius, "ECALEndcapInnerRadius/D");
336  fconfig->Branch("ECALEndcapOuterRadius", &fECALEndcapOuterRadius, "ECALEndcapOuterRadius/D");
337  fconfig->Branch("ECALStartX", &fECALStartX, "ECALStartX/D");
338  fconfig->Branch("ECALEndX", &fECALEndX, "ECALEndX/D");
339  fconfig->Fill();
340 
341  // Analysis tree with reco
342  fTree = tfs->make<TTree>("GArAnaTree","GArAnaTree");
343 
344  fTree->Branch("Run", &fRun, "Run/I");
345  fTree->Branch("SubRun", &fSubRun, "SubRun/I");
346  fTree->Branch("Event", &fEvent, "Event/I");
347  fTree->Branch("Weight", &fWeight, "Weight/F");
348  fTree->Branch("NTPCTracks", &fNTPCTracks, "NTPCTracks/I");
349  fTree->Branch("NTPCECalTracks", &fNTPCECalTracks, "NTPCECalTracks/I");
350 
351  fTree->Branch("TrackIDNumber", &fTrackIDNumber);
352  fTree->Branch("TrackX", &fTrackX);
353  fTree->Branch("TrackY", &fTrackY);
354  fTree->Branch("TrackZ", &fTrackZ);
355  fTree->Branch("TrackPX", &fTrackPX);
356  fTree->Branch("TrackPY", &fTrackPY);
357  fTree->Branch("TrackPZ", &fTrackPZ);
358  fTree->Branch("TrackP", &fTrackP);
359  fTree->Branch("TrackQ", &fTrackQ);
360  fTree->Branch("TrackLen", &fTrackLen);
361  fTree->Branch("TrackChi2", &fTrackChi2);
362  fTree->Branch("NTPCClustersOnTrack", &fNTPCClustersOnTrack);
363  fTree->Branch("TrackTime", &fTrackTime);
364  fTree->Branch("TrackAvgIon", &fTrackAvgIon);
365  fTree->Branch("TrackEndX", &fTrackEndX);
366  fTree->Branch("TrackEndY", &fTrackEndY);
367  fTree->Branch("TrackEndZ", &fTrackEndZ);
368  fTree->Branch("TrackEndPX", &fTrackEndPX);
369  fTree->Branch("TrackEndPY", &fTrackEndPY);
370  fTree->Branch("TrackEndPZ", &fTrackEndPZ);
371  fTree->Branch("TrackEndP", &fTrackEndP);
372  fTree->Branch("TrackEndQ", &fTrackEndQ);
373  fTree->Branch("TrackPID", &fTrackPID);
374 
375  fTree->Branch("NECalClustersForTrack", &fNECalClustersForTrack);
376  fTree->Branch("ClusterTrackIDMatched", &fClusterTrackIDMatched);
377  fTree->Branch("ClusterNhits", &fClusterNhits);
378  fTree->Branch("ClusterEnergy", &fClusterEnergy);
379  fTree->Branch("ClusterTime", &fClusterTime);
380  fTree->Branch("ClusterTimeDiffFirstLast", &fClusterTimeDiffFirstLast);
381  fTree->Branch("ClusterX", &fClusterX);
382  fTree->Branch("ClusterY", &fClusterY);
383  fTree->Branch("ClusterZ", &fClusterZ);
384  fTree->Branch("ClusterTheta", &fClusterTheta);
385  fTree->Branch("ClusterPhi", &fClusterPhi);
386  fTree->Branch("ClusterPID", &fClusterPID);
387  fTree->Branch("ClusterMainAxisX", &fClusterMainAxisX);
388  fTree->Branch("ClusterMainAxisY", &fClusterMainAxisY);
389  fTree->Branch("ClusterMainAxisZ", &fClusterMainAxisZ);
390 
391  fTree->Branch("MCPTrackID", &fMCPTrackID);
392  fTree->Branch("PDG", &fMCPDG);
393  fTree->Branch("PDGMother", &fMCPDGMother);
394  fTree->Branch("MCPStartX", &fMCPStartX);
395  fTree->Branch("MCPStartY", &fMCPStartY);
396  fTree->Branch("MCPStartZ", &fMCPStartZ);
397  fTree->Branch("MCPStartPX", &fMCPStartPX);
398  fTree->Branch("MCPStartPY", &fMCPStartPY);
399  fTree->Branch("MCPStartPZ", &fMCPStartPZ);
400  fTree->Branch("MCPStartP", &fMCPStartP);
401  fTree->Branch("MCPStartEndP", &fMCPStartEndP);
402  fTree->Branch("MCPProc", &fMCPProc);
403  fTree->Branch("MCPEndProc", &fMCPEndProc);
404  fTree->Branch("MCPTime", &fMCPTime);
405  fTree->Branch("MCPInECal", &fMCPInECal);
406  //fTree->Branch("MCPFromVertexID", &fMCPFromVertexID);
407 
408  // Neutrino truth
409  fTree->Branch("Nu_Type", &f_NeutrinoType);
410  fTree->Branch("Nu_CCNC", &f_CCNC);
411  fTree->Branch("Nu_Mode", &f_Mode);
412  fTree->Branch("Nu_InterT", &f_InteractionType);
413  fTree->Branch("Nu_Q2", &f_Q2);
414  fTree->Branch("Nu_W", &f_W);
415  fTree->Branch("Nu_X", &f_X);
416  fTree->Branch("Nu_Y", &f_Y);
417  fTree->Branch("Nu_nuTheta", &f_Theta);
418  fTree->Branch("Nu_nuVertX", &f_MCVertexX);
419  fTree->Branch("Nu_nuVertY", &f_MCVertexY);
420  fTree->Branch("Nu_nuVertZ", &f_MCVertexZ);
421  fTree->Branch("Nu_E", &f_MCnuE);
422  fTree->Branch("Nu_TgtPDG", &f_TgtPDG);
423  fTree->Branch("Nu_VertexInGasTPC", &f_MCVertexInGasTPC);
424 
425  // Truth tree
426  fTruthTree = tfs->make<TTree>("GArTruthTree","GArTruthTree");
427  fTruthTree->Branch("Run", &fRun, "Run/I");
428  fTruthTree->Branch("SubRun", &fSubRun, "SubRun/I");
429  fTruthTree->Branch("Event", &fEvent, "Event/I");
430  fTruthTree->Branch("Weight", &fWeight, "Weight/F");
431 
432  fTruthTree->Branch("Nu_Type", &fNeutrinoType);
433  fTruthTree->Branch("Nu_CCNC", &fCCNC);
434  fTruthTree->Branch("Nu_Mode", &fMode);
435  fTruthTree->Branch("Nu_InterT", &fInteractionType);
436  fTruthTree->Branch("Nu_Q2", &fQ2);
437  fTruthTree->Branch("Nu_W", &fW);
438  fTruthTree->Branch("Nu_X", &fX);
439  fTruthTree->Branch("Nu_Y", &fY);
440  fTruthTree->Branch("Nu_Theta", &fTheta);
441  fTruthTree->Branch("Nu_VertX", &fMCVertexX);
442  fTruthTree->Branch("Nu_VertY", &fMCVertexY);
443  fTruthTree->Branch("Nu_VertZ", &fMCVertexZ);
444  fTruthTree->Branch("Nu_Px", &fMCnuPx);
445  fTruthTree->Branch("Nu_Py", &fMCnuPy);
446  fTruthTree->Branch("Nu_Pz", &fMCnuPz);
447  fTruthTree->Branch("Nu_E", &fMCnuE);
448  fTruthTree->Branch("Nu_TgtPDG", &fTgtPDG);
449  fTruthTree->Branch("Nu_VertexInGasTPC", &fMCVertexInGasTPC);
450  //fTruthTree->Branch("MCVertexID", &fMCVertexID);
451 
452  fTruthTree->Branch("MCP_TrackID", &fMCP_TrackID);
453  fTruthTree->Branch("MCP_PDG", &fMCP_PDG);
454  fTruthTree->Branch("MCP_X", &fMCP_X);
455  fTruthTree->Branch("MCP_Y", &fMCP_Y);
456  fTruthTree->Branch("MCP_Z", &fMCP_Z);
457  fTruthTree->Branch("MCP_T", &fMCP_T);
458  fTruthTree->Branch("MCP_PX", &fMCP_PX);
459  fTruthTree->Branch("MCP_PY", &fMCP_PY);
460  fTruthTree->Branch("MCP_PZ", &fMCP_PZ);
461  fTruthTree->Branch("MCP_E", &fMCP_E);
462  fTruthTree->Branch("MCP_InGasTPC", &fMCP_InGasTPC);
463  //fTruthTree->Branch("MCP_VertexID", &fMCP_VertexID);
464 
465  // Get histograms for TPC PID parameterization
466  std::string filename = "${DUNE_PARDATA_DIR}/MPD/dedxPID/dedxpidmatrices8kevcm.root";
467  TFile infile(filename.c_str(), "READ");
468 
469  m_pidinterp.clear();
470  char str[11];
471  for(int q = 0; q < 501; ++q){
472  sprintf(str, "%d", q);
473  std::string s = "pidmatrix";
474  s.append(str);
475  // read the 500 histograms one by one; each histogram is a
476  // 6 by 6 matrix of probabilities for a given momentum value
477  m_pidinterp.insert( std::make_pair(q, (TH2F*) infile.Get(s.c_str())->Clone("pidinterp")) );
478  }
479 
480  return;
481 
482 } // End of :HNLAnalysis::beginJob
483 
484 //==============================================================================
485 void gar::HNLAnalysis::analyze(art::Event const & event) {
486 
487  /// Is the backtracker good to go?
488  //cheat::BackTrackerCore const* const_bt = gar::providerFrom<cheat::BackTracker>();
489  //fBack = const_cast<cheat::BackTrackerCore*>(const_bt);
490  //if ( !fBack->HasTracks() || !fBack->HasClusters() ) {
491  //throw cet::exception("HNLAnalysis") << " Not enough backtracker info"
492  //<< " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
493  //}
494 
495  ClearVectors();
496  FillVectors(event);
497 
498  return;
499 
500 }
501 
502 //==============================================================================
503 void gar::HNLAnalysis::ClearVectors() {
504 
505  fNTPCECalTracks = 0;
506 
507  // clear out all our vectors
508  fTrackIDNumber.clear();
509  fTrackX.clear();
510  fTrackY.clear();
511  fTrackZ.clear();
512  fTrackPX.clear();
513  fTrackPY.clear();
514  fTrackPZ.clear();
515  fTrackP.clear();
516  fTrackQ.clear();
517  fTrackLen.clear();
518  fTrackChi2.clear();
519  fNTPCClustersOnTrack.clear();
520  fTrackTime.clear();
521  fTrackAvgIon.clear();
522  fTrackEndX.clear();
523  fTrackEndY.clear();
524  fTrackEndZ.clear();
525  fTrackEndPX.clear();
526  fTrackEndPY.clear();
527  fTrackEndPZ.clear();
528  fTrackEndP.clear();
529  fTrackEndQ.clear();
530  fTrackPID.clear();
531 
532  fNECalClustersForTrack.clear();
533  fClusterTrackIDMatched.clear();
534  fClusterNhits.clear();
535  fClusterEnergy.clear();
536  fClusterTime.clear();
537  fClusterTimeDiffFirstLast.clear();
538  fClusterX.clear();
539  fClusterY.clear();
540  fClusterZ.clear();
541  fClusterTheta.clear();
542  fClusterPhi.clear();
543  fClusterPID.clear();
544  fClusterMainAxisX.clear();
545  fClusterMainAxisY.clear();
546  fClusterMainAxisZ.clear();
547 
548  fMCPTrackID.clear();
549  fMCPDG.clear();
550  fMCPDGMother.clear();
551  fMCPStartX.clear();
552  fMCPStartY.clear();
553  fMCPStartZ.clear();
554  fMCPStartPX.clear();
555  fMCPStartPY.clear();
556  fMCPStartPZ.clear();
557  fMCPStartP.clear();
558  fMCPStartEndP.clear();
559  fMCPProc.clear();
560  fMCPEndProc.clear();
561  fMCPTime.clear();
562  fMCPInECal.clear();
563  //fMCPFromVertexID.clear();
564 
565  f_NeutrinoType.clear();
566  f_CCNC.clear();
567  f_Mode.clear();
568  f_InteractionType.clear();
569  f_Q2.clear();
570  f_W.clear();
571  f_X.clear();
572  f_Y.clear();
573  f_Theta.clear();
574  f_MCVertexX.clear();
575  f_MCVertexY.clear();
576  f_MCVertexZ.clear();
577  f_MCnuE.clear();
578  f_TgtPDG.clear();
579  f_MCVertexInGasTPC.clear();
580 
581  fNeutrinoType.clear();
582  fCCNC.clear();
583  fMode.clear();
584  fInteractionType.clear();
585  fQ2.clear();
586  fW.clear();
587  fX.clear();
588  fY.clear();
589  fTheta.clear();
590  fMCVertexX.clear();
591  fMCVertexY.clear();
592  fMCVertexZ.clear();
593  fMCnuPx.clear();
594  fMCnuPy.clear();
595  fMCnuPz.clear();
596  fMCnuE.clear();
597  fTgtPDG.clear();
598  fMCVertexInGasTPC.clear();
599  //fMCVertexID.clear();
600 
601  fMCP_TrackID.clear();
602  fMCP_PDG.clear();
603  fMCP_X.clear();
604  fMCP_Y.clear();
605  fMCP_Z.clear();
606  fMCP_T.clear();
607  fMCP_PX.clear();
608  fMCP_PY.clear();
609  fMCP_PZ.clear();
610  fMCP_E.clear();
611  fMCP_InGasTPC.clear();
612  //fMCP_VertexID.clear();
613 
614  return;
615 
616 } // end :HNLAnalysis::ClearVectors
617 
618 //==============================================================================
619 void gar::HNLAnalysis::FillVectors(art::Event const& event) {
620 
621  // ============= Get art handles ==========================================
622  // Get handles for Tracks, clusters, associations between them
623  auto TrackHandle = event.getHandle< std::vector<rec::Track> >(fTrackLabel);
624  if (!TrackHandle) {
625  throw cet::exception("HNLAnalysis") << " No rec::Track"
626  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
627  }
628  auto RecoClusterHandle = event.getHandle< std::vector<rec::Cluster> >(fClusterLabel);
629  if (!RecoClusterHandle) {
630  throw cet::exception("HNLAnalysis") << " No rec::Cluster branch."
631  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
632  }
633 
634  auto TrackIonHandle = event.getHandle< std::vector<rec::TrackIoniz> >(fTrackLabel);
635  if (!TrackIonHandle) {
636  throw cet::exception("HNLAnalysis") << " No rec::TrackIoniz branch."
637  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
638  }
639 
640  art::FindOneP<rec::TrackIoniz>* findIonization = NULL;
641  findIonization = new art::FindOneP<rec::TrackIoniz>(TrackHandle,event,fTrackLabel);
642 
643  art::FindManyP<rec::Track, rec::TrackEnd>* findManyCALTrackEnd = NULL;
644  findManyCALTrackEnd = new art::FindManyP<rec::Track, rec::TrackEnd>(RecoClusterHandle,event,fECALAssnLabel);
645 
646  // Get handles for MCinfo, also good for MCPTrajectory. Want to get all MCTruths, regardless of generator label.
647 
648  auto mctruthHandles = event.getMany<std::vector<simb::MCTruth> >();
649  if (mctruthHandles.size()!=1) {
650  throw cet::exception("HNLAnalysis") << " Need just 1 simb::MCTruth"
651  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
652  }
653 
654  auto MCPHandle = event.getHandle< std::vector<simb::MCParticle> >(fGeantLabel);
655  if (!MCPHandle) {
656  throw cet::exception("HNLAnalysis") << " No simb::MCParticle"
657  << " Line " << __LINE__ << " in file " << __FILE__ << std::endl;
658  }
659 
660  // ============= Start analysis =========================================
661  fRun = event.run();
662  fSubRun = event.subRun();
663  fEvent = event.id().event();
664  fWeight = 1.0;
665 
666  // Save all neutrino vertices
667  //int nuvertexid = 0;
668  for (size_t imchl = 0; imchl < mctruthHandles.size(); ++imchl) {
669  for ( auto const& mct : (*mctruthHandles.at(imchl)) ) {
670  if(mct.NeutrinoSet()){ // Standard neutrino
671  simb::MCNeutrino nuw = mct.GetNeutrino();
672  fNeutrinoType.push_back(nuw.Nu().PdgCode());
673  fCCNC.push_back(nuw.CCNC());
674  fMode.push_back(nuw.Mode());
675  fInteractionType.push_back(nuw.InteractionType());
676  fQ2.push_back(nuw.QSqr());
677  fW.push_back(nuw.W());
678  fX.push_back(nuw.X());
679  fY.push_back(nuw.Y());
680  fTheta.push_back(nuw.Theta());
681  fMCVertexX.push_back(nuw.Nu().EndX());
682  fMCVertexY.push_back(nuw.Nu().EndY());
683  fMCVertexZ.push_back(nuw.Nu().EndZ());
684  fMCnuPx.push_back(nuw.Nu().Px());
685  fMCnuPy.push_back(nuw.Nu().Py());
686  fMCnuPz.push_back(nuw.Nu().Pz());
687  fMCnuE.push_back( sqrt(nuw.Nu().Px()*nuw.Nu().Px() + nuw.Nu().Py()*nuw.Nu().Py() + nuw.Nu().Pz()*nuw.Nu().Pz()) );
688  fTgtPDG.push_back(nuw.Target());
689  TVector3 nupos(nuw.Nu().EndX(), nuw.Nu().EndY(), nuw.Nu().EndZ());
690  fMCVertexInGasTPC.push_back(fGeo->PointInGArTPC(nupos));
691  //fMCVertexID.push_back(nuvertexid++);
692  }
693  else{
694  fNeutrinoType.push_back(2000020000);
695  fCCNC.push_back(0);
696  fMode.push_back(-1);
697  fInteractionType.push_back(-1);
698  fQ2.push_back(0);
699  fW.push_back(0);
700  fX.push_back(0);
701  fY.push_back(0);
702  fTheta.push_back(0);
703  fMCVertexX.push_back(0);
704  fMCVertexY.push_back(0);
705  fMCVertexZ.push_back(0);
706  fMCnuPx.push_back(0);
707  fMCnuPy.push_back(0);
708  fMCnuPz.push_back(0);
709  fMCnuE.push_back(0);
710  fTgtPDG.push_back(0);
711  fMCVertexInGasTPC.push_back(1);
712  //fMCVertexID.push_back(nuvertexid++);
713  }
714 
715  int nparticles = mct.NParticles();
716  for(int i = 0; i < nparticles; i++){
717  simb::MCParticle mcpart = mct.GetParticle(i);
718 
719  if(mcpart.Weight() != fWeight) fWeight = mcpart.Weight();
720 
721  const TLorentzVector& position = mcpart.Position(0);
722  const TLorentzVector& momentum = mcpart.Momentum(0);
723 
724  fMCP_PDG.push_back(mcpart.PdgCode());
725  fMCP_TrackID.push_back(mcpart.TrackId());
726  fMCP_X.push_back(position.X());
727  fMCP_Y.push_back(position.Y());
728  fMCP_Z.push_back(position.Z());
729  fMCP_T.push_back(mcpart.T());
730  fMCP_PX.push_back(momentum.Px());
731  fMCP_PY.push_back(momentum.Py());
732  fMCP_PZ.push_back(momentum.Pz());
733  fMCP_E.push_back(mcpart.E());
734  TVector3 nupos(position.X(), position.Y(), position.Z());
735  fMCP_InGasTPC.push_back(fGeo->PointInGArTPC(nupos));
736  //fMCP_VertexID.push_back(nuvertexid);
737  }
738  }
739  }
740 
741  // Fill the truth tree
742  fTruthTree->Fill();
743 
744  // Make the usual map for the MC info
745  typedef int TrkId;
746  std::unordered_map<TrkId, Int_t> TrackIdToIndex;
747  Int_t index = 0;
748  for ( auto const& mcp : (*MCPHandle) ) {
749  int TrackId = mcp.TrackId();
750  TrackIdToIndex[TrackId] = index++;
751  }
752 
753  // Loop over TPC tracks
754  size_t iTrack = 0;
755  for ( auto const& track : (*TrackHandle) ) {
756 
757  std::vector<std::pair<simb::MCParticle*,float>> MCsfromTrack;
758  gar::rec::Track* nonconst_track = const_cast<gar::rec::Track*>(&track);
759  MCsfromTrack = fBack->TrackToMCParticles(nonconst_track);
760  int nMCsfromTrack = MCsfromTrack.size();
761  if(nMCsfromTrack==0){
762  continue;
763  }
764 
765  simb::MCParticle theMCPart;
766  for (int iMCfromTrack =0; iMCfromTrack<nMCsfromTrack; ++iMCfromTrack) {
767  // Plausible MCParticle to make this track?
768  theMCPart = *(MCsfromTrack[iMCfromTrack].first);
769  TParticlePDG* particle = pdgInstance->GetParticle(theMCPart.PdgCode());
770  if ( particle->Charge() == 0.0 ) continue;
771  if ( particle->Stable() == 0 ) continue;
772  //if ( particle->PdgCode() == 22 || particle->PdgCode() == 2112) continue;
773  break; // Take highest fraction candidate
774  }
775 
776  // Does theMCParticle go out to the ECAL?
777  int nTraj = theMCPart.NumberTrajectoryPoints();
778  TVector3 doink;
779  bool whackThatECAL = false;
780  for (int iTraj=0; iTraj<nTraj; ++iTraj) {
781  doink.SetXYZ(theMCPart.Vx(iTraj),theMCPart.Vy(iTraj),theMCPart.Vz(iTraj));
782  if ( fGeo->PointInECALBarrel(doink) || fGeo->PointInECALEndcap(doink) ) {
783  whackThatECAL = true;
784  break;
785  }
786  }
787 
788  fMCPInECal.push_back(whackThatECAL);
789 
790  // Which end of this track corresponds to the start of the reco track?
791  const TLorentzVector& positionMCP = theMCPart.Position(0);
792  float distStart = std::hypot(track.Vertex()[0] -positionMCP[0],
793  track.Vertex()[1] -positionMCP[1],
794  track.Vertex()[2] -positionMCP[2]);
795  float distEnd = std::hypot(track.End()[0] -positionMCP[0],
796  track.End()[1] -positionMCP[1],
797  track.End()[2] -positionMCP[2]);
798  rec::TrackEnd kate = (distStart<distEnd) ? rec::TrackEndBeg : rec::TrackEndEnd;
799 
800  // To store parameterized TPC PID
801  std::vector< std::pair<int, float> > pid;
802 
803  // Ready to record some info about the track & matching MC info
804  fTrackIDNumber.push_back(track.getIDNumber());
805  if(kate==rec::TrackEndBeg){
806  fTrackX.push_back (track.Vertex()[0]);
807  fTrackY.push_back (track.Vertex()[1]);
808  fTrackZ.push_back (track.Vertex()[2]);
809  fTrackPX.push_back (track.Momentum_beg()*track.VtxDir()[0]);
810  fTrackPY.push_back (track.Momentum_beg()*track.VtxDir()[1]);
811  fTrackPZ.push_back (track.Momentum_beg()*track.VtxDir()[2]);
812  fTrackP.push_back (track.Momentum_beg());
813  fTrackQ.push_back (track.ChargeBeg());
814  fTrackLen.push_back (track.LengthForward());
815  fTrackChi2.push_back (track.ChisqForward());
816 
817  fTrackEndX.push_back (track.End()[0]);
818  fTrackEndY.push_back (track.End()[1]);
819  fTrackEndZ.push_back (track.End()[2]);
820  fTrackEndPX.push_back (track.Momentum_end()*track.EndDir()[0]);
821  fTrackEndPY.push_back (track.Momentum_end()*track.EndDir()[1]);
822  fTrackEndPZ.push_back (track.Momentum_end()*track.EndDir()[2]);
823  fTrackEndP.push_back (track.Momentum_end());
824  fTrackEndQ.push_back (track.ChargeEnd());
825 
826  // Add the TPC PID information based on Tom's parametrization
827  pid = processPIDInfo( track.Momentum_beg() );
828  }
829  else{
830  fTrackX.push_back (track.End()[0]);
831  fTrackY.push_back (track.End()[1]);
832  fTrackZ.push_back (track.End()[2]);
833  fTrackPX.push_back (track.Momentum_end()*track.EndDir()[0]);
834  fTrackPY.push_back (track.Momentum_end()*track.EndDir()[1]);
835  fTrackPZ.push_back (track.Momentum_end()*track.EndDir()[2]);
836  fTrackP.push_back (track.Momentum_end());
837  fTrackQ.push_back (track.ChargeEnd());
838  fTrackLen.push_back (track.LengthBackward());
839  fTrackChi2.push_back (track.ChisqBackward());
840 
841  fTrackEndX.push_back (track.Vertex()[0]);
842  fTrackEndY.push_back (track.Vertex()[1]);
843  fTrackEndZ.push_back (track.Vertex()[2]);
844  fTrackEndPX.push_back (track.Momentum_beg()*track.VtxDir()[0]);
845  fTrackEndPY.push_back (track.Momentum_beg()*track.VtxDir()[1]);
846  fTrackEndPZ.push_back (track.Momentum_beg()*track.VtxDir()[2]);
847  fTrackEndP.push_back (track.Momentum_beg());
848  fTrackEndQ.push_back (track.ChargeBeg());
849 
850  // Add the TPC PID information based on Tom's parametrization
851  pid = processPIDInfo( track.Momentum_end() );
852  }
853  fNTPCClustersOnTrack.push_back(track.NHits());
854  fTrackTime.push_back(track.Time());
855 
856  // Check TPC PID
857  float pidhypo[5] = {5.0, 5.0, 5.0, 5.0, 5.0};
858  for(size_t ipid = 0; ipid < pid.size(); ipid++){
859  if(pid.at(ipid).first == 1000010020) continue;
860  if (pid.at(ipid).first == 13 && pid.at(ipid).second < pidhypo[0]) pidhypo[0] = pid.at(ipid).second;
861  else if(pid.at(ipid).first == 211 && pid.at(ipid).second < pidhypo[1]) pidhypo[1] = pid.at(ipid).second;
862  else if(pid.at(ipid).first == 321 && pid.at(ipid).second < pidhypo[2]) pidhypo[2] = pid.at(ipid).second;
863  else if(pid.at(ipid).first == 2212 && pid.at(ipid).second < pidhypo[3]) pidhypo[3] = pid.at(ipid).second;
864  else if(pid.at(ipid).first == 11 && pid.at(ipid).second < pidhypo[4]) pidhypo[4] = pid.at(ipid).second;
865  }
866 
867  for(int ihypo = 0; ihypo < 5; ihypo++){
868  if(pidhypo[ihypo] == 5.0) pidhypo[ihypo] = 0.0;
869  }
870 
871  // Backtracker is not perfect
872  int trpdg = abs(theMCPart.PdgCode());
873  if(trpdg == 22) trpdg = 11;
874  else if(trpdg == 2112) trpdg = 2212;
875 
876  CLHEP::RandFlat FlatRand(fEngine);
877  if( (trpdg == 11 and pidhypo[4] == 1) ||
878  (trpdg == 13 and pidhypo[0] == 1) ||
879  (trpdg == 211 and pidhypo[1] == 1) ||
880  (trpdg == 321 and pidhypo[2] == 1) ||
881  (trpdg == 2212 and pidhypo[3] == 1) ) {
882  fTrackPID.push_back(trpdg);
883  }
884  else{
885  int pidpdg = -1;
886  for(int ihypo = 0; ihypo < 5; ihypo++){
887  if(pidhypo[ihypo] == 1 || pidhypo[ihypo] == 0) continue;
888  // Throw a random number between 0 and 1
889  float random_number = FlatRand.fire();
890  if(random_number < pidhypo[ihypo]){
891  pidpdg = ihypo;
892  break;
893  }
894  }
895 
896  if (pidpdg == 0) fTrackPID.push_back(13);
897  else if(pidpdg == 1) fTrackPID.push_back(211);
898  else if(pidpdg == 2) fTrackPID.push_back(321);
899  else if(pidpdg == 3) fTrackPID.push_back(2212);
900  else if(pidpdg == 4) fTrackPID.push_back(11);
901  else fTrackPID.push_back(0);
902  }
903 
904  // Check if the track goes in the ECal and save the ECal info
905  std::vector<ULong64_t> clusterID;
906  size_t iCluster = 0;
907  for ( auto const& cluster : (*RecoClusterHandle) ) {
908  int nCALedTracks(0);
909  if ( findManyCALTrackEnd->isValid() ) {
910  nCALedTracks = findManyCALTrackEnd->at(iCluster).size();
911  }
912 
913  for (int iCALedTrack=0; iCALedTrack<nCALedTracks; ++iCALedTrack) {
914  rec::Track temptrack = *(findManyCALTrackEnd->at(iCluster).at(iCALedTrack));
915  if(temptrack.getIDNumber() == track.getIDNumber()){
916  clusterID.push_back(cluster.getIDNumber());
917  }
918  }
919  iCluster++;
920  }
921 
922  //std::sort(clusterID.begin(), clusterID.end());
923  //auto iter = std::unique(clusterID.begin(), clusterID.end());
924  //clusterID.erase(iter, clusterID.end());
925 
926  fNECalClustersForTrack.push_back(clusterID.size());
927  if(clusterID.size() > 0){
928  fNTPCECalTracks++;
929  for ( auto const& cluster : (*RecoClusterHandle) ) {
930  for(long unsigned int i = 0; i < clusterID.size(); i++){
931  if(clusterID.at(i) == cluster.getIDNumber()){
932  fClusterTrackIDMatched.push_back(track.getIDNumber());
933  fClusterNhits.push_back(cluster.CalorimeterHits().size());
934  fClusterEnergy.push_back(cluster.Energy());
935  fClusterTime.push_back(cluster.Time());
936  fClusterTimeDiffFirstLast.push_back(cluster.TimeDiffFirstLast());
937  fClusterX.push_back(cluster.Position()[0]);
938  fClusterY.push_back(cluster.Position()[1]);
939  fClusterZ.push_back(cluster.Position()[2]);
940  fClusterTheta.push_back(cluster.ITheta());
941  fClusterPhi.push_back(cluster.IPhi());
942  //fClusterPID.push_back(cluster.ParticleID());
943  fClusterPID.push_back(GetECalPIDParameterized( abs(theMCPart.PdgCode()), sqrt(theMCPart.EndPx()*theMCPart.EndPx() + theMCPart.EndPy()*theMCPart.EndPy() + theMCPart.EndPz()*theMCPart.EndPz())) );
944  fClusterMainAxisX.push_back(cluster.EigenVectors()[0]);
945  fClusterMainAxisY.push_back(cluster.EigenVectors()[1]);
946  fClusterMainAxisZ.push_back(cluster.EigenVectors()[2]);
947  }
948  }
949  }
950  }
951 
952  // Save dE/dx
953  if (findIonization->isValid()) {
954  // No calibration for now. Someday this should all be in reco
955  rec::TrackIoniz ionization = *(findIonization->at(iTrack));
956  float avgIonF, avgIonB;
957  processIonizationInfo(ionization, fIonizTruncate, avgIonF, avgIonB);
958  if (kate==rec::TrackEndBeg) {
959  fTrackAvgIon.push_back( avgIonF );
960  }
961  else{
962  fTrackAvgIon.push_back( avgIonB );
963  }
964  }
965  else {
966  // must push_back something so that fTrackAvgIon is of correct size.
967  fTrackAvgIon.push_back( 0.0 );
968  }
969  iTrack++;
970 
971  fMCPDG.push_back(theMCPart.PdgCode());
972  TrkId momTrkId = theMCPart.Mother();
973  int momPDG = 0;
974  if (momTrkId>0) {
975  int momIndex = TrackIdToIndex[momTrkId];
976  momPDG = (*MCPHandle).at(momIndex).PdgCode();
977  }
978  fMCPDGMother.push_back(momPDG);
979 
980  const TLorentzVector& position = theMCPart.Position(0);
981  fMCPStartX.push_back(position.X());
982  fMCPStartY.push_back(position.Y());
983  fMCPStartZ.push_back(position.Z());
984  const TLorentzVector& momentum = theMCPart.Momentum(0);
985  fMCPStartPX.push_back(momentum.Px());
986  fMCPStartPY.push_back(momentum.Py());
987  fMCPStartPZ.push_back(momentum.Pz());
988  fMCPStartP.push_back(theMCPart.P());
989  fMCPStartEndP.push_back(sqrt(theMCPart.EndPx()*theMCPart.EndPx() + theMCPart.EndPy()*theMCPart.EndPy() + theMCPart.EndPz()*theMCPart.EndPz()));
990 
991  fMCPProc.push_back(theMCPart.Process());
992  fMCPEndProc.push_back(theMCPart.EndProcess());
993  fMCPTime.push_back(theMCPart.T());
994  fMCPTrackID.push_back(theMCPart.TrackId());
995 
996  // Find the vertex that this MCParticle is coming from
997  bool vfound = false;
998  for (size_t imchl = 0; imchl < mctruthHandles.size(); ++imchl) {
999  for ( auto const& mct : (*mctruthHandles.at(imchl)) ) {
1000  if (mct.NeutrinoSet()) {
1001  simb::MCNeutrino nuw = mct.GetNeutrino();
1002  Float_t vertX = nuw.Nu().EndX();
1003  Float_t vertY = nuw.Nu().EndY();
1004  Float_t vertZ = nuw.Nu().EndZ();
1005  Float_t dist = std::hypot(position.X()-vertX,position.Y()-vertY,position.Z()-vertZ);
1006 
1007  if ( dist <= fMatchMCPtoVertDist ) { // First neutrino found
1008  f_NeutrinoType.push_back(nuw.Nu().PdgCode());
1009  f_CCNC.push_back(nuw.CCNC());
1010  f_Mode.push_back(nuw.Mode());
1011  f_InteractionType.push_back(nuw.InteractionType());
1012  f_Q2.push_back(nuw.QSqr());
1013  f_W.push_back(nuw.W());
1014  f_X.push_back(nuw.X());
1015  f_Y.push_back(nuw.Y());
1016  f_Theta.push_back(nuw.Theta());
1017  f_MCVertexX.push_back(nuw.Nu().EndX());
1018  f_MCVertexY.push_back(nuw.Nu().EndY());
1019  f_MCVertexZ.push_back(nuw.Nu().EndZ());
1020  f_MCnuE.push_back( sqrt(nuw.Nu().Px()*nuw.Nu().Px() + nuw.Nu().Py()*nuw.Nu().Py() + nuw.Nu().Pz()*nuw.Nu().Pz()) );
1021  f_TgtPDG.push_back(nuw.Target());
1022  TVector3 nupos(nuw.Nu().EndX(), nuw.Nu().EndY(), nuw.Nu().EndZ());
1023  f_MCVertexInGasTPC.push_back(fGeo->PointInGArTPC(nupos));
1024  vfound = true;
1025 
1026  break;
1027  }
1028  }
1029  }
1030  }
1031 
1032  // Fill something if no neutrino is found
1033  if(!vfound){
1034  f_NeutrinoType.push_back(0);
1035  f_CCNC.push_back(0);
1036  f_Mode.push_back(-1);
1037  f_InteractionType.push_back(-1);
1038  f_Q2.push_back(0);
1039  f_W.push_back(0);
1040  f_X.push_back(0);
1041  f_Y.push_back(0);
1042  f_Theta.push_back(0);
1043  f_MCVertexX.push_back(0);
1044  f_MCVertexY.push_back(0);
1045  f_MCVertexZ.push_back(0);
1046  f_MCnuE.push_back(0);
1047  f_TgtPDG.push_back(0);
1048  f_MCVertexInGasTPC.push_back(0);
1049  }
1050 
1051  /*
1052  // Find the vertex that this MCParticle is coming from
1053  nuvertexid = 0; int hnlvertexid = 0; int count = -1;
1054  for (size_t imchl = 0; imchl < mctruthHandles.size(); ++imchl) {
1055  for ( auto const& mct : (*mctruthHandles.at(imchl)) ) {
1056  if (mct.NeutrinoSet()) {
1057  simb::MCNeutrino nuw = mct.GetNeutrino();
1058  Float_t vertX = nuw.Nu().EndX();
1059  Float_t vertY = nuw.Nu().EndY();
1060  Float_t vertZ = nuw.Nu().EndZ();
1061  Float_t dist = std::hypot(position.X()-vertX,position.Y()-vertY,position.Z()-vertZ);
1062  nuvertexid++;
1063  if ( dist <= fMatchMCPtoVertDist ) {
1064  count = nuvertexid;
1065  break;
1066  }
1067  }
1068  else{
1069  int nparticles = mct.NParticles();
1070  hnlvertexid++;
1071  for(int i = 0; i < nparticles; i++){
1072  simb::MCParticle mcpart = mct.GetParticle(i);
1073 
1074  if(mcpart.TrackId() == theMCPart.TrackId()){
1075  count = hnlvertexid;
1076  break;
1077  }
1078 
1079  }
1080  } // else
1081 
1082  }
1083  }
1084 
1085  fMCPFromVertexID.push_back(count);
1086  */
1087  } // end loop over TrackHandle
1088 
1089  fNTPCTracks = (Int_t)iTrack;
1090  if(iTrack > 0) fTree->Fill();
1091 
1092  return;
1093 
1094 } // end HNLAnalysis::FillVectors
1095 
1096 //==============================================================================
1097 // Process ionization. Eventually this moves into the reco code.
1098 void gar::HNLAnalysis::processIonizationInfo(rec::TrackIoniz& ion, float ionizeTruncate, float& forwardIonVal, float& backwardIonVal) {
1099 
1100  // NO CALIBRATION SERVICE FOR NOW
1101 
1102  std::vector<std::pair<float,float>> SigData = ion.getFWD_dSigdXs();
1103  forwardIonVal = processOneDirection(SigData, ionizeTruncate);
1104 
1105  SigData = ion.getBAK_dSigdXs();
1106  backwardIonVal = processOneDirection(SigData, ionizeTruncate);
1107 
1108  return;
1109 
1110 }
1111 
1112 //==============================================================================
1113 float gar::HNLAnalysis::processOneDirection(std::vector<std::pair<float,float>> SigData, float ionizeTruncate) {
1114  //==============================================================================
1115 
1116  std::vector<std::pair<float,float>> dEvsX; // Ionization vs distance along track
1117 
1118  // The first hit on the track never had its ionization info stored. Not a problem
1119  // really. Each pair is a hit and the step along the track that ends at the hit
1120  // For the last hit, just take the step from the n-1 hit; don't guess some distance
1121  // to (nonexistant!) n+1 hit. Using pointer arithmetic because you are a real K&R
1122  // C nerd! Except that C++ doesn't know you are such a nerd and if
1123  // SigData.size()==0, then SigData.end()-1 is 0xFFFFFFFFFFFFFFF8.
1124 
1125  if (SigData.size()==0) return 0.0;
1126  float distAlongTrack = 0;
1127  std::vector<std::pair<float,float>>::iterator littlebit = SigData.begin();
1128  for (; littlebit<(SigData.end()-1); ++littlebit) {
1129  float dE = std::get<0>(*littlebit);
1130  // tpctrackfit2_module.cc fills the TrackIoniz data product so that
1131  // this quantity is really dL > 0 not dX, a coordinate on the drift axis
1132  float dX = std::get<1>(*littlebit);
1133  distAlongTrack += dX; // But count full step to get hit position on track
1134  // Take dX to be 1/2 the previous + last segment
1135  dX += std::get<1>(*(littlebit+1));
1136  float dEdX = dE/(0.5*dX);
1137 
1138  std::pair pushme = std::make_pair(dEdX,distAlongTrack);
1139  dEvsX.push_back( pushme );
1140  }
1141 
1142  // Get the truncated mean; first sort then take mean
1143  std::sort(dEvsX.begin(),dEvsX.end(), lessThan_byE);
1144 
1145  // Get the dEdX vs length data, truncated.
1146  int goUpTo = ionizeTruncate * dEvsX.size() +0.5;
1147  if (goUpTo > (int)dEvsX.size()) goUpTo = dEvsX.size();
1148  int i = 1; float returnvalue = 0;
1149  littlebit = dEvsX.begin();
1150  for (; littlebit<dEvsX.end(); ++littlebit) {
1151  returnvalue += std::get<0>(*littlebit);
1152  ++i;
1153  if (i>goUpTo) break;
1154  }
1155  returnvalue /= goUpTo;
1156 
1157  return returnvalue;
1158 
1159 }
1160 
1161 //==============================================================================
1162 std::vector< std::pair<int, float> > gar::HNLAnalysis::processPIDInfo(float p){
1163  //==============================================================================
1164 
1165  std::vector<std::string> recopnamelist = {"#pi", "#mu", "p", "K", "d", "e"};
1166  std::vector<int> pdg_charged = {211, 13, 2212, 321, 1000010020, 11};
1167  std::vector< std::pair<int, float> > pid;
1168  pid.resize(6);
1169 
1170  int qclosest = 0;
1171  float dist = 100000000.;
1172  CLHEP::RandFlat FlatRand(fEngine);
1173 
1174  for(int q = 0; q < 501; ++q){
1175  //Check the title and the reco momentum take only the one that fits
1176  std::string fulltitle = m_pidinterp[q]->GetTitle();
1177  unsigned first = fulltitle.find("=");
1178  unsigned last = fulltitle.find("GeV");
1179  std::string substr = fulltitle.substr(first+1, last - first-1);
1180  float pidinterp_mom = std::atof(substr.c_str());
1181  //calculate the distance between the bin and mom, store the q the closest
1182  float disttemp = std::abs(pidinterp_mom - p);
1183 
1184  if( disttemp < dist ){
1185  dist = disttemp;
1186  qclosest = q;
1187  }
1188  } // closes the "pidmatrix" loop
1189 
1190  //Compute all the probabities for each type of true to reco
1191  //loop over the columns (true pid)
1192  for(int pidm = 0; pidm < 6; ++pidm){
1193  //loop over the columns (true pid)
1194  std::vector< std::pair<float, std::string> > v_prob;
1195 
1196  //loop over the rows (reco pid)
1197  for(int pidr = 0; pidr < 6; ++pidr){
1198  std::string recoparticlename = m_pidinterp[qclosest]->GetYaxis()->GetBinLabel(pidr+1);
1199  float prob = m_pidinterp[qclosest]->GetBinContent(pidm+1, pidr+1);
1200  //Need to check random number value and prob value then associate the recopdg to the reco prob
1201  v_prob.push_back( std::make_pair(prob, recoparticlename) );
1202  }
1203 
1204  //Compute the pid from it
1205  if(v_prob.size() > 1){
1206  //Order the vector of prob
1207  std::sort(v_prob.begin(), v_prob.end());
1208  //Throw a random number between 0 and 1
1209  float random_number = FlatRand.fire();
1210  //Make cumulative sum to get the range
1211  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);});
1212 
1213  for(size_t ivec = 0; ivec < v_prob.size()-1; ivec++){
1214  if( random_number < v_prob.at(ivec+1).first && random_number >= v_prob.at(ivec).first ){
1215  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) );
1216  }
1217  }
1218  }
1219  else{
1220  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) );
1221  }
1222  }
1223 
1224  //return a vector of pid and prob
1225  return pid;
1226 
1227 }
1228 
1229 //==============================================================================
1230 int gar::HNLAnalysis::GetECalPIDParameterized(int pdg, float ptrue){
1231  //==============================================================================
1232 
1233  CLHEP::RandFlat FlatRand(fEngine);
1234  float ranNum = FlatRand.fire();
1235 
1236  if(abs(pdg) == 11 || pdg == 22) return pdg;
1237  else if(abs(pdg) == 2212) return pdg;
1238  else if(abs(pdg) == 13 || abs(pdg) == 211){
1239  Int_t pdg2 = 13;
1240  if(abs(pdg) == 13) pdg2 = 211;
1241 
1242  if(ptrue < 0.48) return pdg;
1243  if(ptrue >= 0.48 && ptrue < 0.75 && ranNum < 0.8) return pdg;
1244  if(ptrue >= 0.75 && ptrue < 0.90 && ranNum < 0.9) return pdg;
1245  if(ptrue >= 0.90 && ranNum < 0.95) return pdg;
1246 
1247  if(ptrue >= 0.48 && ptrue < 0.75 && ranNum >= 0.8) return pdg2;
1248  if(ptrue >= 0.75 && ptrue < 0.90 && ranNum >= 0.9) return pdg2;
1249  if(ptrue >= 0.90 && ranNum >= 0.95) return pdg2;
1250  }
1251  else{// for any other case return a random PID
1252  if(ranNum < 0.25) return 11;
1253  else if(ranNum >= 0.25 && ranNum < 0.5) return 13;
1254  else if(ranNum >= 0.5 && ranNum < 0.75) return 211;
1255  else return 2212;
1256  }
1257 
1258  return pdg;
1259 
1260 }
1261 
1262 DEFINE_ART_MODULE(gar::HNLAnalysis)
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
gar::HNLAnalysis::fMCPStartPX
std::vector< Float_t > fMCPStartPX
Definition: HNLAnalysis_module.cc:179
gar::HNLAnalysis::fY
std::vector< Float_t > fY
Definition: HNLAnalysis_module.cc:235
art::detail::EngineCreator::createEngine
base_engine_t & createEngine(seed_t seed)
Definition: EngineCreator.cc:20
gar::HNLAnalysis::m_pidinterp
std::unordered_map< int, TH2F * > m_pidinterp
Definition: HNLAnalysis_module.cc:136
gar::HNLAnalysis::fTrackEndPX
std::vector< Float_t > fTrackEndPX
Definition: HNLAnalysis_module.cc:166
simb::MCParticle::NumberTrajectoryPoints
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:218
gar::HNLAnalysis::fClusterNhits
std::vector< Int_t > fClusterNhits
Definition: HNLAnalysis_module.cc:195
gar::HNLAnalysis::fTrackPZ
std::vector< Float_t > fTrackPZ
Definition: HNLAnalysis_module.cc:155
simb::MCParticle::Position
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:219
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
simb::MCNeutrino::CCNC
int CCNC() const
Definition: MCNeutrino.h:148
gar::HNLAnalysis::fMCVertexZ
std::vector< Float_t > fMCVertexZ
Definition: HNLAnalysis_module.cc:240
simb::MCNeutrino::QSqr
double QSqr() const
Definition: MCNeutrino.h:157
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
simb::MCParticle::Py
double Py(const int i=0) const
Definition: MCParticle.h:231
gar::HNLAnalysis::fTheta
std::vector< Float_t > fTheta
Definition: HNLAnalysis_module.cc:236
gar::HNLAnalysis::processPIDInfo
std::vector< std::pair< int, float > > processPIDInfo(float p)
Definition: HNLAnalysis_module.cc:1162
gar::HNLAnalysis::f_Mode
std::vector< Int_t > f_Mode
Definition: HNLAnalysis_module.cc:212
gar::HNLAnalysis::fQ2
std::vector< Float_t > fQ2
Definition: HNLAnalysis_module.cc:232
gar::HNLAnalysis::fMCnuPz
std::vector< Float_t > fMCnuPz
Definition: HNLAnalysis_module.cc:243
gar::HNLAnalysis::analyze
void analyze(art::Event const &e) override
Definition: HNLAnalysis_module.cc:485
simb::MCParticle::EndZ
double EndZ() const
Definition: MCParticle.h:228
simb::MCParticle::Weight
double Weight() const
Definition: MCParticle.h:254
gar::HNLAnalysis::fMCPStartPY
std::vector< Float_t > fMCPStartPY
Definition: HNLAnalysis_module.cc:180
gar::HNLAnalysis::fMCP_PZ
std::vector< Float_t > fMCP_PZ
Definition: HNLAnalysis_module.cc:258
gar::HNLAnalysis::fMCPStartZ
std::vector< Float_t > fMCPStartZ
Definition: HNLAnalysis_module.cc:178
gar::HNLAnalysis::operator=
HNLAnalysis & operator=(HNLAnalysis const &)=delete
gar::HNLAnalysis::fTrackPX
std::vector< Float_t > fTrackPX
Definition: HNLAnalysis_module.cc:153
string
std::string string
Definition: nybbler.cc:12
gar::geo::GeometryCore::PointInECALEndcap
bool PointInECALEndcap(TVector3 const &point) const
Definition: GeometryCore.cxx:624
simb::MCParticle::Mother
int Mother() const
Definition: MCParticle.h:213
gar::HNLAnalysis::fW
std::vector< Float_t > fW
Definition: HNLAnalysis_module.cc:233
gar::HNLAnalysis::fBack
cheat::BackTrackerCore * fBack
Definition: HNLAnalysis_module.cc:130
gar::HNLAnalysis::f_NeutrinoType
std::vector< Int_t > f_NeutrinoType
Definition: HNLAnalysis_module.cc:210
simb::MCNeutrino::Nu
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
gar::HNLAnalysis::fMCP_PDG
std::vector< Int_t > fMCP_PDG
Definition: HNLAnalysis_module.cc:251
gar::HNLAnalysis::HNLAnalysis
HNLAnalysis(fhicl::ParameterSet const &p)
Definition: HNLAnalysis_module.cc:269
gar::geo::GeometryCore::PointInGArTPC
bool PointInGArTPC(TVector3 const &point) const
Definition: GeometryCore.cxx:553
gar::HNLAnalysis::fTrackTime
std::vector< Double_t > fTrackTime
Definition: HNLAnalysis_module.cc:161
simb::MCParticle::Px
double Px(const int i=0) const
Definition: MCParticle.h:230
gar::HNLAnalysis::fMCP_PY
std::vector< Float_t > fMCP_PY
Definition: HNLAnalysis_module.cc:257
vector
struct vector vector
gar::HNLAnalysis::fX
std::vector< Float_t > fX
Definition: HNLAnalysis_module.cc:234
gar::HNLAnalysis::fT
std::vector< Float_t > fT
Definition: HNLAnalysis_module.cc:237
gar::HNLAnalysis::fClusterTheta
std::vector< Float_t > fClusterTheta
Definition: HNLAnalysis_module.cc:202
gar::HNLAnalysis::fMCVertexInGasTPC
std::vector< Int_t > fMCVertexInGasTPC
Definition: HNLAnalysis_module.cc:246
gar::geo::GeometryCore::TPCXCent
float TPCXCent() const
Returns the X location of the center of the TPC in cm.
Definition: GeometryCore.h:778
gar::P
std::pair< float, std::string > P
Definition: anatree_module.cc:75
gar::HNLAnalysis::fMCnuPx
std::vector< Float_t > fMCnuPx
Definition: HNLAnalysis_module.cc:241
gar::geo::GeometryCore
Description of geometry of one entire detector.
Definition: GeometryCore.h:436
gar::HNLAnalysis::pdgInstance
TDatabasePDG * pdgInstance
Definition: HNLAnalysis_module.cc:132
gar::HNLAnalysis::processOneDirection
float processOneDirection(std::vector< std::pair< float, float >> SigData, float ionizeTruncate)
Definition: HNLAnalysis_module.cc:1113
gar::geo::GeometryCore::GetECALOuterEndcapRadius
float GetECALOuterEndcapRadius() const
Definition: GeometryCore.h:951
cluster
Cluster finding and building.
Definition: SmallClusterFilter_module.cc:34
gar::HNLAnalysis::f_MCnuE
std::vector< Float_t > f_MCnuE
Definition: HNLAnalysis_module.cc:223
gar::rec::TrackEndEnd
TrackEnd const TrackEndEnd
Definition: Track.h:33
gar::HNLAnalysis::fMCnuE
std::vector< Float_t > fMCnuE
Definition: HNLAnalysis_module.cc:244
art::EDAnalyzer::EDAnalyzer
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:25
simb::MCParticle::Process
std::string Process() const
Definition: MCParticle.h:215
MCParticle.h
Particle class.
simb::MCParticle::EndY
double EndY() const
Definition: MCParticle.h:227
gar::HNLAnalysis::fMCP_X
std::vector< Float_t > fMCP_X
Definition: HNLAnalysis_module.cc:252
gar::HNLAnalysis::fMCP_InGasTPC
std::vector< Int_t > fMCP_InGasTPC
Definition: HNLAnalysis_module.cc:260
train.filename
string filename
Definition: train.py:213
gar::HNLAnalysis::fCCNC
std::vector< Int_t > fCCNC
Definition: HNLAnalysis_module.cc:229
gar::HNLAnalysis::fMCPDG
std::vector< Int_t > fMCPDG
Definition: HNLAnalysis_module.cc:174
gar::HNLAnalysis::fClusterZ
std::vector< Float_t > fClusterZ
Definition: HNLAnalysis_module.cc:201
gar::HNLAnalysis::fTrackEndP
std::vector< Float_t > fTrackEndP
Definition: HNLAnalysis_module.cc:169
gar::HNLAnalysis::f_Q2
std::vector< Float_t > f_Q2
Definition: HNLAnalysis_module.cc:214
gar::HNLAnalysis::fMCVertexX
std::vector< Float_t > fMCVertexX
Definition: HNLAnalysis_module.cc:238
gar::HNLAnalysis::fMCP_TrackID
std::vector< Int_t > fMCP_TrackID
Definition: HNLAnalysis_module.cc:250
pionana::dEdX
double dEdX(double KE, const simb::MCParticle *part)
Definition: TruthAnalyzer_module.cc:67
gar::HNLAnalysis::fTrackEndQ
std::vector< Int_t > fTrackEndQ
Definition: HNLAnalysis_module.cc:170
cet::hypot
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
Definition: hypot.h:60
gar::HNLAnalysis::f_W
std::vector< Float_t > f_W
Definition: HNLAnalysis_module.cc:215
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:210
gar::HNLAnalysis::fTrackEndX
std::vector< Float_t > fTrackEndX
Definition: HNLAnalysis_module.cc:163
gar::HNLAnalysis::f_Theta
std::vector< Float_t > f_Theta
Definition: HNLAnalysis_module.cc:218
gar::HNLAnalysis::fTrackEndZ
std::vector< Float_t > fTrackEndZ
Definition: HNLAnalysis_module.cc:165
gar::HNLAnalysis::fMCPEndProc
std::vector< std::string > fMCPEndProc
Definition: HNLAnalysis_module.cc:185
gar::HNLAnalysis::fMode
std::vector< Int_t > fMode
Definition: HNLAnalysis_module.cc:230
gar::HNLAnalysis::fClusterX
std::vector< Float_t > fClusterX
Definition: HNLAnalysis_module.cc:199
gar::HNLAnalysis::fTrackZ
std::vector< Float_t > fTrackZ
Definition: HNLAnalysis_module.cc:152
simb::MCNeutrino::InteractionType
int InteractionType() const
Definition: MCNeutrino.h:150
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
gar::HNLAnalysis::fMCPStartEndP
std::vector< Float_t > fMCPStartEndP
Definition: HNLAnalysis_module.cc:183
gar::HNLAnalysis::fMCPStartPZ
std::vector< Float_t > fMCPStartPZ
Definition: HNLAnalysis_module.cc:181
gar::HNLAnalysis::fEngine
CLHEP::HepRandomEngine & fEngine
random engine
Definition: HNLAnalysis_module.cc:135
e
const double e
Definition: gUpMuFluxGen.cxx:165
simb::MCNeutrino::W
double W() const
Definition: MCNeutrino.h:154
DEFINE_ART_MODULE
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
gar::HNLAnalysis::fClusterTimeDiffFirstLast
std::vector< Float_t > fClusterTimeDiffFirstLast
Definition: HNLAnalysis_module.cc:198
infile
string infile
Definition: INukeNucleonCorr.cxx:62
gar::HNLAnalysis::processIonizationInfo
void processIonizationInfo(rec::TrackIoniz &ion, float ionizeTruncate, float &forwardIonVal, float &backwardIonVal)
Definition: HNLAnalysis_module.cc:1098
simb::MCParticle::EndProcess
std::string EndProcess() const
Definition: MCParticle.h:216
gar::HNLAnalysis::beginJob
virtual void beginJob() override
Definition: HNLAnalysis_module.cc:306
gar::geo::GeometryCore::TPCRadius
float TPCRadius() const
Returns the radius of the TPC (y or z direction)
Definition: GeometryCore.h:755
simb::MCNeutrino::Y
double Y() const
Definition: MCNeutrino.h:156
gar::HNLAnalysis::fTrackEndY
std::vector< Float_t > fTrackEndY
Definition: HNLAnalysis_module.cc:164
simb::MCParticle::EndPz
double EndPz() const
Definition: MCParticle.h:243
gar::geo::GeometryCore::GetECALInnerEndcapRadius
float GetECALInnerEndcapRadius() const
Definition: GeometryCore.h:948
a
const double a
Definition: gUpMuFluxGen.cxx:164
gar::HNLAnalysis::fMCPDGMother
std::vector< Int_t > fMCPDGMother
Definition: HNLAnalysis_module.cc:175
simb::MCParticle::P
double P(const int i=0) const
Definition: MCParticle.h:234
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:271
simb::MCParticle::T
double T(const int i=0) const
Definition: MCParticle.h:224
gar::HNLAnalysis::fMCPStartY
std::vector< Float_t > fMCPStartY
Definition: HNLAnalysis_module.cc:177
gar::HNLAnalysis::fClusterPhi
std::vector< Float_t > fClusterPhi
Definition: HNLAnalysis_module.cc:203
gar::HNLAnalysis::fMCPTrackID
std::vector< Int_t > fMCPTrackID
Definition: HNLAnalysis_module.cc:173
test.p
p
Definition: test.py:223
gar::HNLAnalysis::fTrackLen
std::vector< Float_t > fTrackLen
Definition: HNLAnalysis_module.cc:158
gar::geo::GeometryCore::TPCZCent
float TPCZCent() const
Returns the Z location of the center of the TPC in cm.
Definition: GeometryCore.h:792
gar::HNLAnalysis::fClusterPID
std::vector< Int_t > fClusterPID
Definition: HNLAnalysis_module.cc:204
simb::MCNeutrino::X
double X() const
Definition: MCNeutrino.h:155
gar::HNLAnalysis::fIonizTruncate
float fIonizTruncate
Default=1.00;.
Definition: HNLAnalysis_module.cc:116
gar::geo::GeometryCore::PointInECALBarrel
bool PointInECALBarrel(TVector3 const &point) const
Definition: GeometryCore.cxx:606
protoana::distance
double distance(double x1, double y1, double z1, double x2, double y2, double z2)
Definition: truepionXsection_module.cc:217
gar::HNLAnalysis::fTrackPY
std::vector< Float_t > fTrackPY
Definition: HNLAnalysis_module.cc:154
gar::HNLAnalysis::fMatchMCPtoVertDist
float fMatchMCPtoVertDist
Default=roundoff.
Definition: HNLAnalysis_module.cc:118
gar::HNLAnalysis::f_MCVertexZ
std::vector< Float_t > f_MCVertexZ
Definition: HNLAnalysis_module.cc:222
gar::HNLAnalysis::fClusterMainAxisX
std::vector< Float_t > fClusterMainAxisX
Definition: HNLAnalysis_module.cc:205
gar::HNLAnalysis::fClusterTime
std::vector< Float_t > fClusterTime
Definition: HNLAnalysis_module.cc:197
gar::HNLAnalysis::f_TgtPDG
std::vector< Int_t > f_TgtPDG
Definition: HNLAnalysis_module.cc:224
gar::HNLAnalysis::fMCP_PX
std::vector< Float_t > fMCP_PX
Definition: HNLAnalysis_module.cc:256
gar::HNLAnalysis::fTrackY
std::vector< Float_t > fTrackY
Definition: HNLAnalysis_module.cc:151
gar::geo::GeometryCore::GetECALEndcapStartX
float GetECALEndcapStartX() const
Definition: GeometryCore.h:975
gar::HNLAnalysis::fTrackQ
std::vector< Int_t > fTrackQ
Definition: HNLAnalysis_module.cc:157
gar::HNLAnalysis::fTrackEndPY
std::vector< Float_t > fTrackEndPY
Definition: HNLAnalysis_module.cc:167
gar::HNLAnalysis::fMCPInECal
std::vector< Int_t > fMCPInECal
Definition: HNLAnalysis_module.cc:188
gar::HNLAnalysis::f_X
std::vector< Float_t > f_X
Definition: HNLAnalysis_module.cc:216
gar::geo::GeometryCore::GetECALOuterBarrelRadius
float GetECALOuterBarrelRadius() const
Definition: GeometryCore.h:945
gar::HNLAnalysis::fECALAssnLabel
std::string fECALAssnLabel
Definition: HNLAnalysis_module.cc:113
gar
General GArSoft Utilities.
Definition: anatest_module.cc:57
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:221
gar::HNLAnalysis::fMCPStartP
std::vector< Float_t > fMCPStartP
Definition: HNLAnalysis_module.cc:182
gar::HNLAnalysis::fTrackEndPZ
std::vector< Float_t > fTrackEndPZ
Definition: HNLAnalysis_module.cc:168
simb::MCNeutrino::Target
int Target() const
Definition: MCNeutrino.h:151
gar::rec::TrackEndBeg
TrackEnd const TrackEndBeg
Definition: Track.h:33
gar::HNLAnalysis::fTrackChi2
std::vector< Float_t > fTrackChi2
Definition: HNLAnalysis_module.cc:159
gar::HNLAnalysis::fNECalClustersForTrack
std::vector< Int_t > fNECalClustersForTrack
Definition: HNLAnalysis_module.cc:193
gar::HNLAnalysis::f_MCVertexInGasTPC
std::vector< Int_t > f_MCVertexInGasTPC
Definition: HNLAnalysis_module.cc:225
gar::HNLAnalysis::f_T
std::vector< Float_t > f_T
Definition: HNLAnalysis_module.cc:219
gar::HNLAnalysis::fTrackIDNumber
std::vector< ULong64_t > fTrackIDNumber
Definition: HNLAnalysis_module.cc:149
simb::MCParticle::EndPy
double EndPy() const
Definition: MCParticle.h:242
gar::HNLAnalysis::f_CCNC
std::vector< Int_t > f_CCNC
Definition: HNLAnalysis_module.cc:211
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:563
gar::geo::GeometryCore::TPCYCent
float TPCYCent() const
Returns the Y location of the center of the TPC in cm.
Definition: GeometryCore.h:785
gar::cheat::BackTrackerCore::TrackToMCParticles
std::vector< std::pair< simb::MCParticle *, float > > TrackToMCParticles(rec::Track *const t)
Definition: BackTrackerCore.cxx:823
gar::HNLAnalysis::fTrackAvgIon
std::vector< Float_t > fTrackAvgIon
Definition: HNLAnalysis_module.cc:162
gar::HNLAnalysis::f_MCVertexX
std::vector< Float_t > f_MCVertexX
Definition: HNLAnalysis_module.cc:220
gar::HNLAnalysis::f_MCVertexY
std::vector< Float_t > f_MCVertexY
Definition: HNLAnalysis_module.cc:221
gar::HNLAnalysis::fNeutrinoType
std::vector< Int_t > fNeutrinoType
Definition: HNLAnalysis_module.cc:228
simb::MCParticle::Momentum
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:220
gar::HNLAnalysis::f_InteractionType
std::vector< Int_t > f_InteractionType
Definition: HNLAnalysis_module.cc:213
simb::MCParticle::Pz
double Pz(const int i=0) const
Definition: MCParticle.h:232
gar::HNLAnalysis::fMCnuPy
std::vector< Float_t > fMCnuPy
Definition: HNLAnalysis_module.cc:242
gar::HNLAnalysis::fClusterMainAxisZ
std::vector< Float_t > fClusterMainAxisZ
Definition: HNLAnalysis_module.cc:207
gar::HNLAnalysis::fNTPCClustersOnTrack
std::vector< Int_t > fNTPCClustersOnTrack
Definition: HNLAnalysis_module.cc:160
gar::HNLAnalysis::fMCVertexY
std::vector< Float_t > fMCVertexY
Definition: HNLAnalysis_module.cc:239
gar::HNLAnalysis::fInteractionType
std::vector< Int_t > fInteractionType
Definition: HNLAnalysis_module.cc:231
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:223
gar::HNLAnalysis::lessThan_byE
static bool lessThan_byE(std::pair< float, float > a, std::pair< float, float > b)
Definition: HNLAnalysis_module.cc:84
b
static bool * b
Definition: config.cpp:1043
gar::geo::GeometryCore::GetECALInnerBarrelRadius
float GetECALInnerBarrelRadius() const
Definition: GeometryCore.h:942
gar::HNLAnalysis::fTrackX
std::vector< Float_t > fTrackX
Definition: HNLAnalysis_module.cc:150
gar::HNLAnalysis::f_Y
std::vector< Float_t > f_Y
Definition: HNLAnalysis_module.cc:217
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
gar::HNLAnalysis::FillVectors
void FillVectors(art::Event const &e)
Definition: HNLAnalysis_module.cc:619
gar::geo::GeometryCore::TPCLength
float TPCLength() const
Returns the length of the TPC (x direction)
Definition: GeometryCore.h:763
simb::MCParticle::EndPx
double EndPx() const
Definition: MCParticle.h:241
gar::rec::Track::getIDNumber
gar::rec::IDNumber getIDNumber() const
Definition: Track.cxx:42
gar::HNLAnalysis::fMCP_E
std::vector< Float_t > fMCP_E
Definition: HNLAnalysis_module.cc:259
gar::HNLAnalysis::GetECalPIDParameterized
int GetECalPIDParameterized(int pdg, float ptrue)
Definition: HNLAnalysis_module.cc:1230
gar::HNLAnalysis::fMCPStartX
std::vector< Float_t > fMCPStartX
Definition: HNLAnalysis_module.cc:176
gar::HNLAnalysis::fTgtPDG
std::vector< Int_t > fTgtPDG
Definition: HNLAnalysis_module.cc:245
gar::HNLAnalysis::fMCPProc
std::vector< std::string > fMCPProc
Definition: HNLAnalysis_module.cc:184
gar::HNLAnalysis::fMCPTime
std::vector< Float_t > fMCPTime
Definition: HNLAnalysis_module.cc:186
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
simb::MCParticle::EndX
double EndX() const
Definition: MCParticle.h:226
gar::HNLAnalysis::fMCP_Z
std::vector< Float_t > fMCP_Z
Definition: HNLAnalysis_module.cc:254
gar::HNLAnalysis::fGeo
const gar::geo::GeometryCore * fGeo
Definition: HNLAnalysis_module.cc:128
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
GeometryGAr.h
art framework interface to geometry description
gar::HNLAnalysis::fMCP_T
std::vector< Float_t > fMCP_T
Definition: HNLAnalysis_module.cc:255
simb::MCNeutrino::Mode
int Mode() const
Definition: MCNeutrino.h:149
gar::HNLAnalysis::fMCP_Y
std::vector< Float_t > fMCP_Y
Definition: HNLAnalysis_module.cc:253
gar::HNLAnalysis::fTrackPID
std::vector< Int_t > fTrackPID
Definition: HNLAnalysis_module.cc:171
gar::HNLAnalysis::fClusterTrackIDMatched
std::vector< Int_t > fClusterTrackIDMatched
Definition: HNLAnalysis_module.cc:194
s
static QCString * s
Definition: config.cpp:1042
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:222
str
static QCString str
Definition: fortrancode.cpp:27098
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
endl
QTextStream & endl(QTextStream &s)
Definition: qtextstream.cpp:2030
QCString::append
QCString & append(const char *s)
Definition: qcstring.cpp:383
event
Event finding and building.
Definition: EventCheater_module.cc:32
gar::geo::GeometryCore::GetECALEndcapOuterX
float GetECALEndcapOuterX() const
Definition: GeometryCore.h:978
gar::HNLAnalysis::fClusterEnergy
std::vector< Float_t > fClusterEnergy
Definition: HNLAnalysis_module.cc:196
gar::HNLAnalysis::fClusterY
std::vector< Float_t > fClusterY
Definition: HNLAnalysis_module.cc:200
gar::HNLAnalysis::fClusterMainAxisY
std::vector< Float_t > fClusterMainAxisY
Definition: HNLAnalysis_module.cc:206
gar::HNLAnalysis::fTrackP
std::vector< Float_t > fTrackP
Definition: HNLAnalysis_module.cc:156
Generated by   doxygen 1.8.11