RecoMuon_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file RecoMuon_module.cc
3 // \brief A module for creating reconstructed muon id objects
4 // \version $Id: RecoMuon_module.cxx,v 1.2 2012-10-06 01:56:18 gsdavies Exp $
5 // \author Nicholas Raddatz - raddatz@physics.umn.edu
6 ////////////////////////////////////////////////////////////////////////
7 #include <vector>
8 #include <string>
9 #include <limits>
10 
11 // ROOT includes
12 #include "TH2D.h"
13 #include "TVector3.h"
14 #include "TFile.h"
15 #include "TMath.h"
16 #include "TMVA/Factory.h"
17 #include "TMVA/Tools.h"
18 #include "TMVA/Reader.h"
19 
20 // Framework includes
24 #include "fhiclcpp/ParameterSet.h"
29 #include "Utilities/AssociationUtil.h"
31 #include "cetlib/search_path.h"
32 #include "cetlib/filesystem.h"
33 
34 // NOvA includes
35 #include "ReMId/ReMId.h"
36 #include "Calibrator/Calibrator.h"
37 #include "Geometry/Geometry.h"
38 #include "GeometryObjects/Geo.h"
41 #include "RecoBase/Track.h"
42 #include "RecoBase/Cluster.h"
43 #include "RecoBase/CellHit.h"
44 #include "RecoBase/RecoHit.h"
46 #include "Utilities/func/MathUtil.h"
48 #include "NumuEnergy/NumuEAlg.h"
49 #include "SummaryData/SpillData.h"
51 
52 
53 
54 /// A PID for muons
55 namespace remid {
56  class RecoMuon : public art::EDProducer {
57  public:
58 
59  explicit RecoMuon(fhicl::ParameterSet const &pset);
60 
61  virtual ~RecoMuon();
62 
63  void produce(art::Event &evt);
64 
65  private:
66  void Init(art::Event &evt);
67 
68  std::string fHistLibPath; ///< path to the REMID library histogram files listed below (usually set to $REMID_LIB_PATH)
69  std::string fPIDLibPath; ///< path to the REMID library PID files listed below (usually set to $REMID_LIB_PATH)
70 
71  std::string fDedxSignalFile; ///<Path to library samples
72  std::string fDedxBackgroundFile; ///<Path to library samples
73  std::string fMuonScatterFile; ///<Path to library samples
74  std::string fPionScatterFile; ///<Path to library samples
75  std::string fPIDMethodFD; ///< name of the TMVA method for the FD PID, fileName = "tmva_"+method+".xml"
76  std::string fPIDMethodND; ///< name of the TMVA method for the ND PID, fileName = "tmva_"+method+".xml"
77  std::string fSliceLabel; ///<Where to find reconstructed slices
78  std::string fTrackLabel; ///<Where to find reconstructed tracks
79  std::string fGeneratorLabel; ///<Used for the FHC/RHC switch with simulated data
80  std::string fNuMIBeamLabel; ///<Used for the FHC/RHC switch with real data
81 
82  fhicl::ParameterSet fTrackCleanUpAlgPSet; ///< Parameter set to configure the TrackCleanUpAlg object
83  fhicl::ParameterSet fNumuEAlgPSet; ///< Parameter set to configure the NumuEAlg object
84 
87 
90 
91  std::vector<double> fmuondedxnorm;
92  std::vector<double> fpiondedxnorm;
93  std::vector<double> fmuonscatnorm;
94  std::vector<double> fpionscatnorm;
95 
97  TMVA::Reader *fReader;
98  Float_t LLE, LLScat, length, measfrac;
99 
102 
105 
107 
108  void Get2DHistoBins(TH2D *histogram, double xval, double yval, int &xbin, int &ybin);
109 
110  };
111 
112  //....................................................................................
113 
115  trainfileloaded(false),
116  fReader(nullptr),
117  fTrackCleanUpAlg(0),
118  fNumuEAlg(0) {
119  fDedxSignalFile = pset.get<std::string>("DedxSignalFile");
120  fDedxBackgroundFile = pset.get<std::string>("DedxBackgroundFile");
121  fMuonScatterFile = pset.get<std::string>("MuonScatterFile");
122  fPionScatterFile = pset.get<std::string>("PionScatterFile");
123 
124  fHistLibPath = pset.get<std::string>("HistLibPath");
125  fPIDLibPath = pset.get<std::string>("PIDLibPath");
126 
127  fGeneratorLabel = pset.get<std::string>("GeneratorLabel");
128  fNuMIBeamLabel = pset.get<std::string>("NuMIBeamLabel");
129 
130  fPIDMethodFD = pset.get<std::string>("PIDMethodFD");
131  fPIDMethodND = pset.get<std::string>("PIDMethodND");
132 
133  fSliceLabel = pset.get<std::string>("SliceLabel");
134  fTrackLabel = pset.get<std::string>("TrackLabel");
135  fTrackCleanUpAlgPSet = pset.get<fhicl::ParameterSet>("TrackCleanUpAlgPSet");
136  fNumuEAlgPSet = pset.get<fhicl::ParameterSet>("NumuEAlgPSet");
137 
138  produces < std::vector < remid::ReMId > > ();
139  produces < art::Assns < remid::ReMId, rb::Track > > ();
140 
141  }
142 
143  //....................................................................................
144 
146  if (fTrackCleanUpAlg) { delete fTrackCleanUpAlg; }
147  if (fNumuEAlg) { delete fNumuEAlg; }
148  if (fReader) { delete fReader; }
149  }
150 
151  //....................................................................................
152 
154  // One of many signs that would indicate this function has already done
155  // its thing.
156  if (fTrackCleanUpAlg != 0) return;
157 
159 
160  trainfileloaded = false;
161 
164 
165 
167 
168  if (!evt.isRealData()) {
169  evt.getByLabel(fGeneratorLabel, spillPot);
170  } else {
171  evt.getByLabel(fNuMIBeamLabel, spillPot);
172  }
173  bool isRHC = false;
174  if (spillPot.failedToGet()) {
175  throw cet::exception("RecoMuon") << "Spill Data not found, aborting without horn current information" << std::endl;
176  } else {
177  isRHC = spillPot->isRHC;
178  }
179 
180  std::string hornSuffix;
181 
182 
183  // Modifying histogram names based on FHC/RHC
184  if (isRHC) {
185  hornSuffix = "_RHC";
186  } else {
187  hornSuffix = "_FHC";
188  }
189 
190  // Grab the Dedx histograms - signal
191  fileName = libpath + fDedxSignalFile + hornSuffix + ".root";
192  if (!cet::file_exists(fileName)) {
193  throw cet::exception("RecoMuon")
194  << "dEdX signal file: " << fileName << " can't be located" << std::endl;
195  }
196  TFile muondedxfile(fileName.c_str(), "READ");
197 
198  dEdxVxMuonSample = (TH2D *) muondedxfile.Get("DedxVxS");
199 
200  // Grab the Dedx histograms - background
201  fileName = libpath + fDedxBackgroundFile + hornSuffix + ".root";
202  if (!cet::file_exists(fileName)) {
203  throw cet::exception("RecoMuon")
204  << "dEdX background file: " << fileName << " can't be located" << std::endl;
205  }
206  TFile piondedxfile(fileName.c_str(), "READ");
207 
208 // std::string dEdxVxPionName = "DedxVxB" + hornSuffix;
209  dEdxVxPionSample = (TH2D *) piondedxfile.Get("DedxVxB");
210 
211 
212  // normalize these histograms to give probabilities in each bin
213  for (int i = 1; i <= dEdxVxMuonSample->GetNbinsX(); ++i) {
214  double xbintot = 0;
215  // JOSHNOTE: Check if these two loops can be merged into one.
216  for (int j = 0; j <= dEdxVxMuonSample->GetNbinsY() + 1; ++j) {
217  xbintot += dEdxVxMuonSample->GetBinContent(i, j);
218  }
219  //scale the y bins in this x bin column to total number in this x bin, i.e. normalize to total probability of 1
220  for (int j = 0; j <= dEdxVxMuonSample->GetNbinsY() + 1; ++j) {
221  if (xbintot != 0) {
222  dEdxVxMuonSample->SetBinContent(i, j, dEdxVxMuonSample->GetBinContent(i, j) / xbintot);
223  } else { dEdxVxMuonSample->SetBinContent(i, j, -1); }
224  }
225  fmuondedxnorm.push_back(xbintot);
226  }
227 
228  // normalize these histograms to give probabilities in each bin
229  for (int i = 1; i <= dEdxVxPionSample->GetNbinsX(); ++i) {
230  double xbintot = 0;
231  for (int j = 0; j <= dEdxVxPionSample->GetNbinsY() + 1; ++j) {
232  xbintot += dEdxVxPionSample->GetBinContent(i, j);
233  }
234  //scale the y bins in this x bin column to total number in this x bin, i.e. normalize to total probability of 1
235  for (int j = 0; j <= dEdxVxPionSample->GetNbinsY() + 1; ++j) {
236  if (xbintot != 0) {
237  dEdxVxPionSample->SetBinContent(i, j, dEdxVxPionSample->GetBinContent(i, j) / xbintot);
238  } else { dEdxVxPionSample->SetBinContent(i, j, -1); }
239  }
240  fpiondedxnorm.push_back(xbintot);
241  }
242  maxpionlength = 550.0; // Maximum length to sample the pion ll dedx histogram
243 
244  // Grab scattering histograms - signal
245  fileName = libpath + fMuonScatterFile + hornSuffix + ".root";
246  if (!cet::file_exists(fileName)) {
247  throw cet::exception("RecoMuon")
248  << "scattering signal file: " << fileName << " can't be located" << std::endl;
249  }
250  TFile muonscatfile(fileName.c_str(), "READ");
251 // std::string ScatVxMuonName = "ScatVxS" + hornSuffix;
252  ScatVxMuonSample = (TH2D *) muonscatfile.Get("ScatVxS");
253 
254  // Grab scattering histograms - background
255  fileName = libpath + fPionScatterFile + hornSuffix + ".root";
256  if (!cet::file_exists(fileName)) {
257  throw cet::exception("RecoMuon")
258  << "scattering background file: " << fileName << " can't be located" << std::endl;
259  }
260  TFile pionscatfile(fileName.c_str(), "READ");
261 // std::string ScatVxPionName = "ScatVxB" + hornSuffix;
262  ScatVxPionSample = (TH2D *) pionscatfile.Get("ScatVxB");
263 
264  // normalize these histograms to give probabilities in each bin - signal
265  for (int i = 1; i <= ScatVxMuonSample->GetNbinsX(); ++i) {
266  double xbintot = 0;
267  for (int j = 0; j <= ScatVxMuonSample->GetNbinsY() + 1; ++j) {
268  xbintot += ScatVxMuonSample->GetBinContent(i, j);
269  }
270  //scale the y bins in this x bin column to total number in this x bin, i.e. normalize to total probability of 1
271  for (int j = 0; j <= ScatVxMuonSample->GetNbinsY() + 1; ++j) {
272  if (xbintot != 0) {
273  ScatVxMuonSample->SetBinContent(i, j, ScatVxMuonSample->GetBinContent(i, j) / xbintot);
274  } else { ScatVxMuonSample->SetBinContent(i, j, -1); }
275  }
276  fmuonscatnorm.push_back(xbintot);
277  }
278 
279  // normalize these histograms to give probabilities in each bin - background
280  for (int i = 1; i <= ScatVxPionSample->GetNbinsX(); ++i) {
281  double xbintot = 0;
282  for (int j = 0; j <= ScatVxPionSample->GetNbinsY() + 1; ++j) {
283  xbintot += ScatVxPionSample->GetBinContent(i, j);
284  }
285  //scale the y bins in this x bin column to total number in this x bin, i.e. normalize to total probability of 1
286  for (int j = 0; j <= ScatVxPionSample->GetNbinsY() + 1; ++j) {
287  if (xbintot != 0) {
288  ScatVxPionSample->SetBinContent(i, j, ScatVxPionSample->GetBinContent(i, j) / xbintot);
289  } else { ScatVxPionSample->SetBinContent(i, j, -1); }
290  }
291  fpionscatnorm.push_back(xbintot);
292  }
293 
294 
296 
298 
299  std::string xmlFileName;
300 
301  if (geom->DetId() == novadaq::cnv::kFARDET) {
302  xmlFileName = libpath + "tmva_" + fPIDMethodFD + hornSuffix + ".fd" + ".xml";
303  } else if (geom->DetId() == novadaq::cnv::kNEARDET) {
304  xmlFileName = libpath + "tmva_" + fPIDMethodND + hornSuffix + ".nd" + ".xml";
305  } else {
306  std::cerr << "Unknown detector ID. Expecting FarDet or NearDet, but got geom->DetID(): \"" << geom->DetId()
307  << "\"" << std::endl;
308  }
309 
310  if (!cet::file_exists(xmlFileName)) {
311  throw cet::exception("RecoMuon")
312  << "TMVA weight file \"" << (xmlFileName) << "\" can't be found. Got the correct ReMId UPS version?"
313  << std::endl;
314  }
315 
316  // Book Method
317  if (!trainfileloaded) {
318  // Set up the TMVA::Reader to read in MVA weight file
319  fReader = new TMVA::Reader();
320  // Add the pid variables
321  fReader->AddVariable("trackLength", &length);
322  fReader->AddVariable("dedxSep", &LLE);
323  fReader->AddVariable("scatSep", &LLScat);
324  fReader->AddVariable("measFrac", &measfrac);
325  if (geom->DetId() == novadaq::cnv::kFARDET) fReader->BookMVA(fPIDMethodFD + " Method", xmlFileName);
326  else if (geom->DetId() == novadaq::cnv::kNEARDET) fReader->BookMVA(fPIDMethodND + " Method", xmlFileName);
327  else
328  throw cet::exception("RecoMuon") << "Unknown detector ID at TMVA Book time. Got \"" << geom->DetId()
329  << "\"" << std::endl;
330  trainfileloaded = true;
331  }
332  }
333  //....................................................................................
334 
335 
337 // art::Event &evnt;
338 // art::Event evnt = evt;
339  Init(evt);
340 
341  std::unique_ptr <std::vector<remid::ReMId>> pidcol(new std::vector <remid::ReMId>);
342  std::unique_ptr <art::Assns<remid::ReMId, rb::Track>> assncol(new art::Assns <remid::ReMId, rb::Track>);
343 
344  // get a geometry handle
347 
348  fDetWidth = geom->DetHalfWidth();
349  fDetHeight = geom->DetHalfHeight();
350  fDetLength = geom->DetLength();
351 
352  double tranZPos = -1.0;
353  if (geom->DetId() == novadaq::cnv::kNEARDET) { tranZPos = fNumuEAlg->NDTransitionZPos(); }
354 
356  evt.getByLabel(fSliceLabel, slicecol);
358  for (unsigned int i = 0; i < slicecol->size(); ++i) {
359  art::Ptr <rb::Cluster> slice(slicecol, i);
360  slicelist.push_back(slice);
361  }
362 
363 
365 
366  // loop over all of the slices
367  for (unsigned int iSlice = 0; iSlice < slicelist.size(); ++iSlice) {
368  // get all of the tracks in this slice
369  const std::vector <art::Ptr<rb::Track>> slicetrks = fmt.at(iSlice);
370  for (unsigned int iTrack = 0; iTrack < slicetrks.size(); ++iTrack) {
371  // make a pid object for this track
372  ReMId pid;
373  pid.SetVal(-1); // Set default pid value
374  pid.SetPdg(13); // Set that this is a pid for muons
375 
376  art::Ptr <rb::Track> track = slicetrks[iTrack];
377 
378  // store some useful track info
379  double trklen = track->TotalLength();
380 
381  unsigned int ntraj = track->NTrajectoryPoints();
382  // vector holds effective active track length from the end of the track
383  std::vector<double> lengths(ntraj);
384  // vector holds catcher track length from the end of the track
385  std::vector<double> catlengths(ntraj);
386 
387  // loop over trajectory points starting from the end
388  for (int itrj = ntraj - 1; itrj >= 0; --itrj) {
389  TVector3 trajpoint = track->TrajectoryPoint(itrj);
390  catlengths[itrj] = 0.0;
391  if (itrj == (int) ntraj - 1) {
392  double lend = 0.0;
393  if (geom->DetId() == novadaq::cnv::kNEARDET && trajpoint.Z() > tranZPos) {
394  double catlength = 0.0;
395  // convert the catcher track length to an energy
396  double catenergy = fNumuEAlg->NDMuonCatcherEForActiveAndCatcher(catlength);
397  // get the effective active length of this catcher energy
398  lend = fNumuEAlg->ActiveTrackLengthFromMuonE(catenergy);
399  catlengths[itrj] = catlength;
400  }
401  lengths[itrj] = lend;
402  } else {
403  TVector3 lasttrajpoint = track->TrajectoryPoint(itrj + 1);
404  lengths[itrj] = (trajpoint - lasttrajpoint).Mag() + lengths[itrj + 1];
405  catlengths[itrj] = catlengths[itrj + 1];
406  if (geom->DetId() == novadaq::cnv::kNEARDET) {
407  double catlength = 0.0;
408  double actlength = 0.0;
409  // check if this is all or partially in the muon catcher
410  if (lasttrajpoint.Z() > tranZPos) {
411  // all in the muon catcher
412  if (trajpoint.Z() > tranZPos) { catlength = (trajpoint - lasttrajpoint).Mag(); }
413  else {
414  // only part of this length is in the muon catcher,
415  // split it up appropriately between active and catcher length
416  TVector3 diff = lasttrajpoint - trajpoint;
417  catlength = fabs(tranZPos - lasttrajpoint.Z()) * diff.Mag() / fabs(diff.Z());
418  actlength = fabs(tranZPos - trajpoint.Z()) * diff.Mag() / fabs(diff.Z());
419  }
420  } else { actlength = (trajpoint - lasttrajpoint).Mag(); }
421  // convert the catlength to an active det length
422  double effactlen = 0.0;
423  if (catlength > 0.0) {
424  double totcatlength = catlength + catlengths[itrj + 1];
425  catlengths[itrj] = totcatlength;
426  // convert the total catcher track length to an energy
427  double catenergy = fNumuEAlg->NDMuonCatcherEForActiveAndCatcher(totcatlength);
428  // get the effective active length of this catcher energy
429  effactlen = fNumuEAlg->ActiveTrackLengthFromMuonE(catenergy);
430  } else { effactlen = lengths[itrj + 1]; }
431  lengths[itrj] = (effactlen + actlength);
432  } // end of if ND
433  }
434  } // end over loop over trajectory points
435 
436  trklen = lengths[0];
437 
438  unsigned int minPlane = track->MinPlane();
439  unsigned int maxPlane = track->MaxPlane();
440 
441  int maxXPlane = -1;
442  int maxYPlane = -1;
443 
444  double hitCosx[track->NXCell()];
445  double hitCosy[track->NYCell()];
446  for (unsigned int ihit = 0; ihit < track->NXCell(); ++ihit) {
447  if (track->XCell(ihit)->Plane() > maxXPlane) { maxXPlane = track->XCell(ihit)->Plane(); }
448  // get direction with this x hit
449  double xyzplane[3];
450  geom->Plane(track->XCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
451  double dz = geom->Plane(track->XCell(ihit)->Plane())->Cell(0)->HalfD();
452  double minz = xyzplane[2] - dz;
453  double maxz = xyzplane[2] + dz;
454  unsigned int currtrj = 0;
455  // find the matching trajectory point for this plane;
456  for (unsigned int itraj = 0; itraj < track->NTrajectoryPoints(); ++itraj) {
457  if (track->TrajectoryPoint(itraj).Z() > minz &&
458  track->TrajectoryPoint(itraj).Z() < maxz) { currtrj = itraj; }
459  }
460  TVector3 direction = track->Dir();
461  if (currtrj > 0 && currtrj < track->NTrajectoryPoints() - 2) {
462  direction = track->TrajectoryPoint(currtrj) - track->TrajectoryPoint(currtrj - 1);
463  }
464  direction = direction.Unit();
465  hitCosx[ihit] =
466  sqrt(1.0 / (direction.X() * direction.X() + direction.Z() * direction.Z())) * direction.Z();
467  }
468 
469  for (unsigned int ihit = 0; ihit < track->NYCell(); ++ihit) {
470  if (track->YCell(ihit)->Plane() > maxYPlane) { maxYPlane = track->YCell(ihit)->Plane(); }
471  // get direction with this y hit
472  double xyzplane[3];
473  geom->Plane(track->YCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
474  double dz = geom->Plane(track->YCell(ihit)->Plane())->Cell(0)->HalfD();
475  double minz = xyzplane[2] - dz;
476  double maxz = xyzplane[2] + dz;
477  unsigned int currtrj = 0;
478  // find the matching trajectory point for this plane;
479  for (unsigned int itraj = 0; itraj < track->NTrajectoryPoints(); ++itraj) {
480  if (track->TrajectoryPoint(itraj).Z() > minz &&
481  track->TrajectoryPoint(itraj).Z() < maxz) { currtrj = itraj; }
482  }
483  TVector3 direction = track->Dir();
484  if (currtrj > 0 && currtrj < track->NTrajectoryPoints() - 2) {
485  direction = track->TrajectoryPoint(currtrj) - track->TrajectoryPoint(currtrj - 1);
486  }
487  direction = direction.Unit();
488  hitCosy[ihit] =
489  sqrt(1.0 / (direction.Y() * direction.Y() + direction.Z() * direction.Z())) * direction.Z();
490  }
491 
492  // calculate dedx info for this track
493  std::vector<double> dedxs;
494  std::vector<double> xdedxs;
495  std::vector<unsigned int> dedxpls;
496 
497  for (unsigned int iPlane = minPlane; iPlane <= maxPlane; ++iPlane) {
498  double planeGeV = 0;
499  double totcosinview = 0;
500  double avgCosOtherView = 1.0;
501  bool useOtherView = true;
502  double nhits = 0;
503  double xyzplane[3];
504  geom->Plane(iPlane)->Cell(0)->GetCenter(xyzplane);
505  double dz = geom->Plane(iPlane)->Cell(0)->HalfD();
506  double minz = xyzplane[2] - dz;
507  double maxz = xyzplane[2] + dz;
508  int matchpt = 0;
509  int minCell = 10000;
510  int maxCell = 0;
511  double nhitsother = 0;
512  int view = 0;
513  // find the matching trajectory point for this plane;
514  for (unsigned int itraj = 0; itraj < track->NTrajectoryPoints(); ++itraj) {
515  if (track->TrajectoryPoint(itraj).Z() > minz &&
516  track->TrajectoryPoint(itraj).Z() < maxz) { matchpt = itraj; }
517  }
518 
519  // find the track length up to this trajectory point
520  double xlength = trklen - lengths[matchpt];
521 
522  if (minPlane == maxPlane) { xlength = trklen / 2; }
523  if (track->Stop().Z() < track->Start().Z()) { xlength = (trklen - xlength); }
524  if (geom->Plane(iPlane)->View() == geo::kX) {
525  for (unsigned int iHit = 0; iHit < track->NXCell(); ++iHit) {
526  if (track->XCell(iHit)->Plane() == iPlane && track->XCell(iHit)->View() != geo::kY) {
527  rb::RecoHit rhit = track->RecoHit(track->XCell(iHit));
528  if (minPlane == maxPlane) {
529  rhit = cal->MakeRecoHit(*(track->XCell(iHit)), 0.25);
530  }
531  if (rhit.IsCalibrated()) {
532  planeGeV += rhit.GeV();
533  totcosinview += hitCosx[iHit];
534  nhits += 1.0;
535  if (track->XCell(iHit)->Cell() < minCell) { minCell = track->XCell(iHit)->Cell(); }
536  if (track->XCell(iHit)->Cell() > maxCell) { maxCell = track->XCell(iHit)->Cell(); }
537  }
538  }
539  }
540  // try to get the avgcos from the other view
541  double nhitsbefore = 0;
542  double nhitsafter = 0;
543  double totCosBefore = 0;
544  double totCosAfter = 0;
545  for (unsigned int iHit = 0; iHit < track->NYCell(); ++iHit) {
546  if (track->YCell(iHit)->Plane() == iPlane - 1) {
547  ++nhitsbefore;
548  totCosBefore += hitCosy[iHit];
549  } else if (track->YCell(iHit)->Plane() == iPlane + 1) {
550  ++nhitsafter;
551  totCosAfter += hitCosy[iHit];
552  }
553  }
554  if (nhitsbefore > 0 && nhitsafter > 0) {
555  avgCosOtherView = ((totCosAfter + totCosBefore) / (nhitsbefore + nhitsafter));
556  nhitsother = (nhitsbefore + nhitsafter) / 2;
557  } else if (nhitsbefore > 0) {
558  avgCosOtherView = totCosBefore / nhitsbefore;
559  nhitsother = nhitsbefore;
560  } else if (nhitsafter > 0) {
561  avgCosOtherView = totCosAfter / nhitsafter;
562  nhitsother = nhitsafter;
563  } else { useOtherView = false; }
564  } else if (geom->Plane(iPlane)->View() == geo::kY) {
565  view = 1;
566  for (unsigned int iHit = 0; iHit < track->NYCell(); ++iHit) {
567  if (track->YCell(iHit)->Plane() == iPlane && track->YCell(iHit)->View() != geo::kX
568  && track->YCell(iHit)->View() != geo::kXorY) {
569  rb::RecoHit rhit = track->RecoHit(track->YCell(iHit));
570  if (minPlane == maxPlane) {
571  rhit = cal->MakeRecoHit(*(track->YCell(iHit)), 0.25);
572  }
573  if (rhit.IsCalibrated()) {
574  planeGeV += rhit.GeV();
575  totcosinview += hitCosy[iHit];
576  nhits += 1.0;
577  if (track->YCell(iHit)->Cell() < minCell) { minCell = track->YCell(iHit)->Cell(); }
578  if (track->YCell(iHit)->Cell() > maxCell) { maxCell = track->YCell(iHit)->Cell(); }
579  }
580  }
581  }
582  // try to get the avgcos from the other view
583  double nhitsbefore = 0;
584  double nhitsafter = 0;
585  double totCosBefore = 0;
586  double totCosAfter = 0;
587  for (unsigned int iHit = 0; iHit < track->NXCell(); ++iHit) {
588  if (track->XCell(iHit)->Plane() == iPlane - 1) {
589  ++nhitsbefore;
590  totCosBefore += hitCosx[iHit];
591  } else if (track->XCell(iHit)->Plane() == iPlane + 1) {
592  ++nhitsafter;
593  totCosAfter += hitCosx[iHit];
594  }
595  }
596  if (nhitsbefore > 0 && nhitsafter > 0) {
597  avgCosOtherView = TMath::Abs(((totCosAfter + totCosBefore) / (nhitsbefore + nhitsafter)));
598  nhitsother = (nhitsbefore + nhitsafter) / 2;
599  } else if (nhitsbefore > 0) {
600  avgCosOtherView = TMath::Abs(totCosBefore / nhitsbefore);
601  nhitsother = nhitsbefore;
602  } else if (nhitsafter > 0) {
603  avgCosOtherView = TMath::Abs(totCosAfter / nhitsafter);
604  nhitsother = nhitsbefore;
605  } else { useOtherView = false; }
606  }
607 
608  if (nhits > 0) {
609  double avgCos = TMath::Abs(totcosinview / nhits);
610  double xyzp[3];
611  geom->Plane(iPlane)->Cell(minCell)->GetCenter(xyzp);
612  double xyzp1[3];
613  geom->Plane(iPlane)->Cell(maxCell)->GetCenter(xyzp1);
614  double planewidth = 2 * geom->Plane(iPlane)->Cell(1)->HalfD();
615 
616  // calculate the amount of dead material to subtract off
617  // first add up live material
618  double liveLength = 0;
619  for (int i = minCell; i < maxCell + 1; ++i) {
620  // get the cell half width
621  const geo::PlaneGeo *p = geom->Plane(iPlane);
622  const geo::CellGeo *c = p->Cell(i);
623  if (i == minCell || i == maxCell) { liveLength += c->HalfW(); }
624  else { liveLength += (2.0 * c->HalfW()); }
625  }
626  // dead material is total length - live length only if min and max cell are not the same
627  double deadLength = 0;
628  if (minCell != maxCell) { deadLength = (TMath::Abs(xyzp1[view] - xyzp[view]) - liveLength); }
629 
630  double planeheight = planewidth * TMath::Tan(TMath::ACos(avgCos)) - deadLength;
631  double planePathLength = TMath::Sqrt(planewidth * planewidth + planeheight * planeheight);
632  // in the case of completely vertical track be a bit careful
633  if (avgCos == 0 || avgCos != avgCos) { planePathLength = liveLength; }
634  if (useOtherView) {
635  double deltaz = planewidth;
636  double deltav = planewidth * TMath::Tan(TMath::ACos(avgCos)) - deadLength;
637  if (avgCos == 0 || avgCos != avgCos) {
638  deltav = liveLength;
639  // if the track looks completely vertical in both views,
640  // than no part of the track length should be in z direction
641  if (avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView) { deltaz = 0; }
642  }
643  // in opposite view dead length
644  double deltaov = (planewidth - deadLength) * TMath::Tan(TMath::ACos(avgCosOtherView));
645  // in the case that the other view looks completely vertical scale up
646  // the length between cells in this view by the number of hits in the other view
647  // this case should really never happen
648  if (avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView) {
649  deltaov = nhitsother * (TMath::Abs(xyzp1[view] - xyzp[view])) / nhits;
650  }
651  planePathLength = TMath::Sqrt(deltaz * deltaz + deltav * deltav + deltaov * deltaov);
652  }
653  // also the case of completely vertical track
654  if (minPlane == maxPlane) { planePathLength = liveLength; }
655  assert(planeGeV / (planePathLength) >= 0);
656 
657  dedxs.push_back(planeGeV / (planePathLength));
658  xdedxs.push_back(xlength);
659  dedxpls.push_back(iPlane);
660 
661  } // nhits > 0
662 
663  } // end of loop over planes of the track
664 
665  // check if this track is contained
666  bool trackcont = IsTrackContained(track);
667  pid.SetContained(trackcont);
668 
669  // dedx log-likelihoods for muon and pion
670  double llmuon = 0;
671  double nllmuon = 0;
672  double llpion = 0;
673  double nllpion = 0;
674 
675  // get the vertex region planes of this track
676  int vertexPlanes[2];
677  fTrackCleanUpAlg->LastVertexPlanes(*(track), vertexPlanes);
678 
679  // is there any vertex free planes?
680  int minVert[2] = {0, 0};
681  if (vertexPlanes[0] < maxXPlane || vertexPlanes[1] < maxYPlane) {
682  minVert[0] = vertexPlanes[0];
683  minVert[1] = vertexPlanes[1];
684  }
685 
686  // calculate log-likelihood for every dedx measurement made not in vertex planes
687  for (unsigned int idedx = 0; idedx < dedxs.size(); ++idedx) {
688 
689  // check if in vertex region
690  if (geom->Plane(dedxpls[idedx])->View() == geo::kX) {
691  if ((int) dedxpls[idedx] < minVert[0]) { continue; }
692  } else {
693  if ((int) dedxpls[idedx] < minVert[1]) { continue; }
694  }
695  if (dedxpls[idedx] == maxPlane && maxPlane != minPlane) { continue; }
696 
697  // get the distance from track end of this dedx measurement
698  double distFromEnd = trklen - xdedxs[idedx];
699 
700  int xbin = 0;
701  int dedxbin = 0;
702  Get2DHistoBins(dEdxVxMuonSample, distFromEnd, dedxs[idedx], xbin, dedxbin);
703 
704  // compute the probability from the muon sample histogram
705  double prob = dEdxVxMuonSample->GetBinContent(xbin, dedxbin);
706  if (prob <= 0) {
707  // find the normalization of this distance bin
708  double norm = fmuondedxnorm[xbin - 1];
709  prob = 1 / norm;
710  }
711  llmuon += TMath::Log(prob);
712  ++nllmuon;
713 
714  } // end of loop over dedxs
715 
716  // get pion dedx log-likelihood
717  int maxpionbin = dEdxVxPionSample->GetXaxis()->FindFixBin(maxpionlength);
718  for (unsigned int idedx = 0; idedx < dedxs.size(); ++idedx) {
719 
720  pid.AddDedx(dedxs[idedx], xdedxs[idedx] / trklen);
721 
722  // check if in vertex region
723  if (geom->Plane(dedxpls[idedx])->View() == geo::kX) {
724  if ((int) dedxpls[idedx] < minVert[0]) {
725  pid.SetDedxVertex(idedx, true);
726  pid.SetDedxUsed(idedx, false);
727  continue;
728  }
729  } else {
730  if ((int) dedxpls[idedx] < minVert[1]) {
731  pid.SetDedxVertex(idedx, true);
732  pid.SetDedxUsed(idedx, false);
733  continue;
734  }
735  }
736 
737  if (dedxpls[idedx] == maxPlane && maxPlane != minPlane) {
738  pid.SetDedxUsed(idedx, false);
739  continue;
740  }
741 
742  // get the distance from track end of this dedx measurement
743  double distFromEnd = trklen - xdedxs[idedx];
744  int xbin = 0;
745  int dedxbin = 0;
746  Get2DHistoBins(dEdxVxPionSample, distFromEnd, dedxs[idedx], xbin, dedxbin);
747 
748  if (xbin > maxpionbin) { xbin = maxpionbin; }
749  // compute the probability from the pion sample histogram
750  double prob = dEdxVxPionSample->GetBinContent(xbin, dedxbin);
751  if (prob <= 0) {
752  // find the normalization of this distance bin
753  double norm = fmuondedxnorm[xbin - 1];
754  prob = 1 / norm;
755  }
756  llpion += TMath::Log(prob);
757  ++nllpion;
758  }
759 
760 
761  double llscatmuon = 0;
762  double nllscatmuon = 0;
763  double llscatpion = 0;
764  double nllscatpion = 0;
765 
766  for (unsigned int idir = 1; idir < track->NTrajectoryPoints() - 1; ++idir) {
767  // get the scattering angle
768  TVector3 initDir = track->TrajectoryPoint(idir) - track->TrajectoryPoint(idir - 1);
769  TVector3 currDir = track->TrajectoryPoint(idir + 1) - track->TrajectoryPoint(idir);
770  double cosScat = initDir.Dot(currDir) / (initDir.Mag() * currDir.Mag());
771 
772  // get position where this scatter occurs
773  double lastdist = fabs(lengths[idir] - lengths[idir - 1]);
774 
775  pid.AddScatter(track->TrajectoryPoint(idir), cosScat);
776 
777  double scatVar = TMath::ACos(cosScat) * TMath::ACos(cosScat) / lastdist;
778 
779  double distFromEnd = lengths[idir];
780 
781  // get the muon likelihood of this scat at this position
782  bool useMu = true;
783  int xbin = 0;
784  int scatbin = 0;
785  Get2DHistoBins(ScatVxMuonSample, distFromEnd, scatVar, xbin, scatbin);
786  double prob = ScatVxMuonSample->GetBinContent(xbin, scatbin);
787  if (useMu) {
788  if (prob <= 0) {
789  // find the normalization of this distance bin
790  double norm = fmuonscatnorm[xbin - 1];
791  prob = 1 / norm;
792  }
793  llscatmuon += TMath::Log(prob);
794  ++nllscatmuon;
795  }
796 
797  // get the pion likelihood of this scat at this position
798  bool usePi = true;
799  xbin = 0;
800  scatbin = 0;
801  Get2DHistoBins(ScatVxPionSample, distFromEnd, scatVar, xbin, scatbin);
802  prob = ScatVxPionSample->GetBinContent(xbin, scatbin);
803  if (usePi) {
804  if (prob <= 0) {
805  // find the normalization of this distance bin
806  double norm = fmuonscatnorm[xbin - 1];
807  prob = 1 / norm;
808  }
809  llscatpion += TMath::Log(prob);
810  ++nllscatpion;
811  }
812  }
813 
814  // Set the Muon and the Pion Log-Likelihoods
815  // First default values
816  pid.SetMuonDedxLL(-100);
817  pid.SetPionDedxLL(-100);
818  pid.SetMuonScatLL(0);
819  pid.SetPionScatLL(0);
820 
821  if (nllmuon > 0 && nllpion > 0 && nllscatmuon > 0 && nllscatpion > 0) {
822  // If non-zero measurement planes, then get real values of Log-Likelihoods
823  pid.SetMuonDedxLL(llmuon / nllmuon);
824  pid.SetPionDedxLL(llpion / nllpion);
825  pid.SetMuonScatLL(llscatmuon / nllscatmuon);
826  pid.SetPionScatLL(llscatpion / nllscatpion);
827 
828  LLScat = pid.ScatSeparation();
829  LLE = pid.DedxSeparation();
830  measfrac = ((float) pid.NMeasurementPlanes()) /
831  ((float) pid.NMeasurementPlanes() + (float) pid.NVertexPlanes());
832 
833  // get the PID output value for this track
834  std::vector<float> mvaInput(4);
835  mvaInput[0] = trklen;
836  mvaInput[1] = LLE;
837  mvaInput[2] = LLScat;
838  mvaInput[3] = measfrac;
839 
840  // get pid value
841  double score = std::numeric_limits<double>::signaling_NaN();
842  if (geom->DetId() == novadaq::cnv::kFARDET)
843  score = fReader->EvaluateMVA(mvaInput, fPIDMethodFD + " Method");
844  else if (geom->DetId() == novadaq::cnv::kNEARDET)
845  score = fReader->EvaluateMVA(mvaInput, fPIDMethodND + " Method");
846 
847  // Rescale the BDTG output to be [0,1] (i.e. - the same as the original kNN)
848  // NOTE: The BDT{A,D} have a different range and have NOT been rescaled here.
849  // That is a low-priority to-do in the future...
850  if (geom->DetId() == novadaq::cnv::kNEARDET && fPIDMethodND == "BDTG") {
851  score = (score + 1.0) / 2.0;
852  }
853  if (geom->DetId() == novadaq::cnv::kFARDET && fPIDMethodFD == "BDTG") {
854  score = (score + 1.0) / 2.0;
855  }
856 
857  pid.SetVal(score);
858 
859  }
860 
861  pidcol->push_back(pid);
862  util::CreateAssn(*this, evt, *(pidcol.get()), track, *(assncol.get()));
863 
864  } // end of loop over tracks on slice
865  } // end of loop over slice
866 
867  evt.put(std::move(pidcol));
868  evt.put(std::move(assncol));
869  }
870 
871  //..........................................................................
873 
874  // only check end of track containment, asssume uncontained until proven otherwise
875  bool stopcontained = false;
876 
878 
879  // FD track end containment file
880  if (geom->DetId() == novadaq::cnv::kFARDET) {
881  double xmin, xmax, ymin, ymax, zmin, zmax;
882  double voffset = 25.0;
883  xmin = -fDetWidth + voffset;
884  xmax = fDetWidth - voffset;
885  ymin = -fDetHeight + voffset;
886  ymax = fDetHeight - voffset;
887  zmin = 0 + voffset;
888  zmax = fDetLength - voffset;
889  TVector3 mins(xmin, ymin, zmin);
890  TVector3 maxs(xmax, ymax, zmax);
891 
892  // track stop containment
893  TVector3 stop = track->Stop();
894  // if outside of detector bounds, not contained
895  if (!(stop.X() > maxs.X() || stop.X() < mins.X() ||
896  stop.Y() > maxs.Y() || stop.Y() < mins.Y() ||
897  stop.Z() > maxs.Z() || stop.Z() < mins.Z())) {
898  if (track->NTrajectoryPoints() > 2) {
899  double end[3] = {stop.X(), stop.Y(), stop.Z()};
900  TVector3 traj = track->TrajectoryPoint(track->NTrajectoryPoints() - 2);
901  TVector3 enddir = (stop - traj).Unit();
902  double dir[3] = {enddir.X(), enddir.Y(), enddir.Z()};
903  double edge[3];
905  edge);
906  // check distance from track end to detector edge
907  TVector3 projedge(edge[0], edge[1], edge[2]);
908  double distToEdge = (projedge - stop).Mag();
909  if (distToEdge > 50) { stopcontained = true; }
910  }
911  }
912  } else if (geom->DetId() == novadaq::cnv::kNEARDET) {
913  // ND track end containment
914  TVector3 stop = track->Stop();
915  double xmin, xmax, ymin, ymax, zmin, zmax;
916  double voffset = 25.0;
917  // check if z position places us in the muon catcher
918  double muoncxyz[3];
919  geom->Plane(geom->FirstPlaneInMuonCatcher() - 1)->Cell(0)->GetCenter(muoncxyz);
920  if (stop.Z() > muoncxyz[2]) {
921  // stop is in the muon catcher
922  voffset = 50.0;
923  xmin = -fDetWidth + voffset;
924  xmax = fDetWidth - voffset;
925  ymin = -fDetHeight + voffset;
926  ymax = (fDetHeight / 3.0) - voffset;
927  zmin = 0 + voffset;
928  zmax = fDetLength - voffset;
929  } else {
930  // end is not in the muon catcher
931  xmin = -fDetWidth + voffset;
932  xmax = fDetWidth - voffset;
933  ymin = -fDetHeight + voffset;
934  ymax = fDetHeight - voffset;
935  zmin = 0 + voffset;
936  zmax = fDetLength - voffset;
937  }
938  TVector3 mins(xmin, ymin, zmin);
939  TVector3 maxs(xmax, ymax, zmax);
940  // if outside of detector bounds, not contained
941  if (!(stop.X() > maxs.X() || stop.X() < mins.X() ||
942  stop.Y() > maxs.Y() || stop.Y() < mins.Y() ||
943  stop.Z() > maxs.Z() || stop.Z() < mins.Z())) {
944  if (track->NTrajectoryPoints() > 2) {
945  double end[3] = {stop.X(), stop.Y(), stop.Z()};
946  TVector3 traj = track->TrajectoryPoint(track->NTrajectoryPoints() - 2);
947  TVector3 enddir = (stop - traj).Unit();
948  double dir[3] = {enddir.X(), enddir.Y(), enddir.Z()};
949  double edge[3];
950  if (stop.Z() < muoncxyz[2]) {
952  muoncxyz[2], edge);
953  } else {
955  edge);
956  }
957  // check distance from track end to detector edge
958  TVector3 projedge(edge[0], edge[1], edge[2]);
959  double distToEdge = (projedge - stop).Mag();
960  if (distToEdge > 50) { stopcontained = true; }
961  }
962  }
963  }
964 
965  return stopcontained;
966  } // End of IsTrackContained function
967 
968  //..........................................................................
969 
970  void RecoMuon::Get2DHistoBins(TH2D *histogram, double xval, double yval, int &xbin, int &ybin) {
971 
972  // get the x bin for the histogram
973  TAxis *xaxis = histogram->GetXaxis();
974  if (xval < xaxis->GetXmin()) { xbin = 1; }
975  else if (xval > xaxis->GetXmax()) { xbin = histogram->GetNbinsX(); }
976  else { xbin = xaxis->FindFixBin(xval); }
977 
978  // get the y bin for the histogram
979  TAxis *yaxis = histogram->GetYaxis();
980  if (yval < yaxis->GetXmin()) { ybin = 1; }
981  else if (yval > yaxis->GetXmax()) { ybin = histogram->GetNbinsY(); }
982  else { ybin = yaxis->FindFixBin(yval); }
983 
984  } // End of Get2DHistoBins
985 
986  //..........................................................................
987 
989 
990 }
std::string fPIDMethodND
name of the TMVA method for the ND PID, fileName = "tmva_"+method+".xml"
std::vector< double > fpiondedxnorm
size_t NTrajectoryPoints() const
Definition: Track.h:83
bool isRHC
is the beam in antineutrino mode, aka RHC
Definition: SpillData.h:28
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
double NDMuonCatcherEForActiveAndCatcher(double catTrkLen) const
Function that, given muon track length in cm, returns muon energy in GeV. This is only for ND muons t...
Definition: NumuEAlg.cxx:1098
fileName
Definition: plotROC.py:78
void SetPionScatLL(double ll)
Definition: ReMId.cxx:162
double HalfD() const
Definition: CellGeo.cxx:205
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
std::map< std::string, double > xmax
void Init(art::Event &evt)
std::string fGeneratorLabel
Used for the FHC/RHC switch with simulated data.
numue::NumuEAlg * fNumuEAlg
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::vector< double > fmuondedxnorm
X or Y views.
Definition: PlaneGeo.h:30
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
double ScatSeparation() const
Return the scattering separation variable used as an input to the kNN that determines the pid value...
Definition: ReMId.cxx:206
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
std::vector< double > fpionscatnorm
std::string EnvExpansion(const std::string &inString)
Function to expand environment variables.
Definition: EnvExpand.cxx:8
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::string yaxis
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
void SetDedxVertex(unsigned int idx, bool vert)
Definition: ReMId.cxx:174
Definition: event.h:19
double HalfW() const
Definition: CellGeo.cxx:191
double DetLength() const
TMVA::Reader * fReader
std::string fMuonScatterFile
Path to library samples.
OStream cerr
Definition: OStream.cxx:7
std::string fSliceLabel
Where to find reconstructed slices.
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void SetDedxUsed(unsigned int idx, bool used)
Definition: ReMId.cxx:168
const PlaneGeo * Plane(unsigned int i) const
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
void SetContained(bool cont)
Definition: ReMId.cxx:180
std::string fDedxBackgroundFile
Path to library samples.
std::string xaxis
bool IsTrackContained(art::Ptr< rb::Track > track)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
std::string fDedxSignalFile
Path to library samples.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
Double_t ymax
Definition: plot.C:25
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
void Get2DHistoBins(TH2D *histogram, double xval, double yval, int &xbin, int &ybin)
void SetMuonDedxLL(double ll)
Definition: ReMId.cxx:56
unsigned short Cell() const
Definition: CellHit.h:40
std::string fPIDMethodFD
name of the TMVA method for the FD PID, fileName = "tmva_"+method+".xml"
fhicl::ParameterSet fTrackCleanUpAlgPSet
Parameter set to configure the TrackCleanUpAlg object.
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
murem::TrackCleanUpAlg * fTrackCleanUpAlg
Far Detector at Ash River, MN.
fhicl::ParameterSet fNumuEAlgPSet
Parameter set to configure the NumuEAlg object.
double dz[NP][NC]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
std::string fHistLibPath
path to the REMID library histogram files listed below (usually set to $REMID_LIB_PATH) ...
void AddDedx(double dedx, double distFromStart, bool used=true, bool vert=false)
Definition: ReMId.cxx:47
int NMeasurementPlanes() const
Definition: ReMId.cxx:223
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
RecoMuon(fhicl::ParameterSet const &pset)
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
int evt
An algorithm to determine the energy of a slice.
Definition: NumuEAlg.h:28
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
Near Detector in the NuMI cavern.
Collect Geo headers and supply basic geometry functions.
void SetPionDedxLL(double ll)
Definition: ReMId.cxx:62
int NVertexPlanes() const
Definition: ReMId.cxx:233
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
const double j
Definition: BetheBloch.cxx:29
double DetHalfHeight() const
double NDTransitionZPos() const
Function that returns the z position of the middle of the transition plane - the last active plane be...
Definition: NumuEAlg.cxx:1193
void AddScatter(TVector3 pos, double cosScat)
Definition: ReMId.cxx:40
TVector3 Unit() const
std::string fTrackLabel
Where to find reconstructed tracks.
void SetMuonScatLL(double ll)
Definition: ReMId.cxx:156
size_type size() const
Definition: PtrVector.h:308
std::string fPIDLibPath
path to the REMID library PID files listed below (usually set to $REMID_LIB_PATH) ...
double ActiveTrackLengthFromMuonE(double muonE) const
Function that, given muon energy in GeV, returns an effective track length in active material in cm...
Definition: NumuEAlg.cxx:951
void SetPdg(int pdg)
Definition: PID.h:23
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
float GeV() const
Definition: RecoHit.cxx:69
double DetHalfWidth() const
std::string fPionScatterFile
Path to library samples.
Float_t norm
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
void geom(int which=0)
Definition: geom.C:163
float Mag() const
void produce(art::Event &evt)
assert(nhit_max >=nhit_nbins)
Double_t ymin
Definition: plot.C:24
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
bool file_exists(std::string const &qualified_filename)
void ProjectToBoxEdge(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi, double xyzout[])
Project along a direction from a particular starting point to the edge of a box.
Definition: Geo.cxx:38
A PID for muons.
Definition: FillPIDs.h:10
std::vector< double > fmuonscatnorm
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
void SetVal(double val)
Definition: PID.h:24
Encapsulate the cell geometry.
Definition: CellGeo.h:25
void LastVertexPlanes(const rb::Track &muonTrack, int *vertexPlanes, double *vertexRegionDeDxCutOff=NULL)
A function to decide how large a vertex region should be considered for clean up. vertexRegion is an ...
double DedxSeparation() const
Return the dE/dx separation variable used as an input to the kNN that determines the pid value...
Definition: ReMId.cxx:199
Encapsulate the geometry of one entire detector (near, far, ndos)
std::string fNuMIBeamLabel
Used for the FHC/RHC switch with real data.
bool failedToGet() const
Definition: Handle.h:196
const unsigned int FirstPlaneInMuonCatcher() const
Returns the index of the first plane contained in the muon catcher.
enum BeamMode string