LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
MVAAlg.cxx
Go to the documentation of this file.
1 // \fileMVAAlg.cxx
3 // m.haigh@warwick.ac.uk
5 
16 
21 #include "cetlib/search_path.h"
22 #include "cetlib_except/exception.h"
23 #include "fhiclcpp/ParameterSet.h"
24 
25 #include "Fit/Fitter.h"
26 #include "Math/Functor.h"
27 #include "TPrincipal.h"
28 
29 #include <cmath>
30 #include <string>
31 #include <vector>
32 
34  : fCaloAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg")), fReader("")
35 {
36  fHitLabel = pset.get<std::string>("HitLabel");
37  fTrackLabel = pset.get<std::string>("TrackLabel");
38  fShowerLabel = pset.get<std::string>("ShowerLabel");
39  fSpacePointLabel = pset.get<std::string>("SpacePointLabel");
40  fTrackingLabel = pset.get<std::string>("TrackingLabel", "");
41 
42  fCheatVertex = pset.get<bool>("CheatVertex", false);
43 
44  fReader.AddVariable("evalRatio", &fResHolder.evalRatio);
45  fReader.AddVariable("coreHaloRatio", &fResHolder.coreHaloRatio);
46  fReader.AddVariable("concentration", &fResHolder.concentration);
47  fReader.AddVariable("conicalness", &fResHolder.conicalness);
48  fReader.AddVariable("dEdxStart", &fResHolder.dEdxStart);
49  fReader.AddVariable("dEdxEnd", &fResHolder.dEdxEnd);
50  fReader.AddVariable("dEdxEndRatio", &fResHolder.dEdxEndRatio);
51 
52  fMVAMethods = pset.get<std::vector<std::string>>("MVAMethods");
53  std::vector<std::string> weightFileBnames = pset.get<std::vector<std::string>>("WeightFiles");
54 
55  cet::search_path searchPath("FW_SEARCH_PATH");
56  for (auto fileIter = weightFileBnames.begin(); fileIter != weightFileBnames.end(); ++fileIter) {
57  std::string fileWithPath;
58  if (!searchPath.find_file(*fileIter, fileWithPath)) {
59  fWeightFiles.clear();
60  fMVAMethods.clear();
61  throw cet::exception("MVAPID") << "Unable to find weight file " << *fileIter
62  << " in search path " << searchPath.to_string() << std::endl;
63  }
64  fWeightFiles.push_back(fileWithPath);
65  }
66 
67  if (fMVAMethods.size() != fWeightFiles.size()) {
68  std::cerr << "Mismatch in number of MVA methods and weight files!" << std::endl;
69  exit(1);
70  }
71 
72  for (unsigned int iMethod = 0; iMethod != fMVAMethods.size(); ++iMethod) {
73  fReader.BookMVA(fMVAMethods[iMethod], fWeightFiles[iMethod]);
74  }
75 }
76 
77 int mvapid::MVAAlg::IsInActiveVol(const TVector3& pos)
78 {
79  const double fiducialDist = 5.0;
80 
81  if (pos.X() > (fDetMinX + fiducialDist) && pos.X() < (fDetMaxX - fiducialDist) &&
82  pos.Y() > (fDetMinY + fiducialDist) && pos.Y() < (fDetMaxY - fiducialDist) &&
83  pos.Z() > (fDetMinZ + fiducialDist) && pos.Z() < (fDetMaxZ - fiducialDist))
84  return 1;
85  else
86  return 0;
87 }
88 
90 {
91 
93 
94  fDetMinX = 999999.9;
95  fDetMaxX = -999999.9;
96  fDetMinY = 999999.9;
97  fDetMaxY = -999999.9;
98  fDetMinZ = 999999.9;
99  fDetMaxZ = -999999.9;
100 
101  for (auto const& tpc : geom->Iterate<geo::TPCGeo>()) {
102  if (tpc.MinX() < fDetMinX) fDetMinX = tpc.MinX();
103  if (tpc.MaxX() > fDetMaxX) fDetMaxX = tpc.MaxX();
104  if (tpc.MinY() < fDetMinY) fDetMinY = tpc.MinY();
105  if (tpc.MaxY() > fDetMaxY) fDetMaxY = tpc.MaxY();
106  if (tpc.MinZ() < fDetMinZ) fDetMinZ = tpc.MinZ();
107  if (tpc.MaxZ() > fDetMaxZ) fDetMaxZ = tpc.MaxZ();
108  }
109 }
110 
112 {
113 
115 
116  fNormToWiresY.clear();
117  fNormToWiresZ.clear();
118 
119  int planeKey;
120 
121  //Get normals to wires for each plane in the detector
122  //This assumes equal numbers of TPCs in each cryostat and equal numbers of planes in each TPC
123  for (auto const& plane : geom->Iterate<geo::PlaneGeo>()) {
124  std::string id = std::string(plane.ID());
125  int pcryo = id.find("C");
126  int ptpc = id.find("T");
127  int pplane = id.find("P");
128  std::string scryo = id.substr(pcryo + 2, 2);
129  std::string stpc = id.substr(ptpc + 2, 2);
130  std::string splane = id.substr(pplane + 2, 2);
131  int icryo = std::stoi(scryo);
132  int itpc = std::stoi(stpc);
133  int iplane = std::stoi(splane);
134  planeKey = icryo * geom->NTPC() * geom->Nplanes() + itpc * geom->Nplanes() +
135  iplane; //single index for all planes in detector
136  fNormToWiresY.insert(
137  std::make_pair(planeKey, -plane.Wire(0).Direction().Z())); //y component of normal
138  fNormToWiresZ.insert(
139  std::make_pair(planeKey, plane.Wire(0).Direction().Y())); //z component of normal
140  }
141 }
142 
144  std::vector<anab::MVAPIDResult>& result,
147 {
148  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
149  auto const detProp =
151  this->PrepareEvent(evt, clockData);
152 
153  for (auto trackIter = fTracks.begin(); trackIter != fTracks.end(); ++trackIter) {
154  mvapid::MVAAlg::SortedObj sortedObj;
155 
156  std::vector<double> eVals, eVecs;
157  int isStoppingReco;
158  this->RunPCA(fTracksToHits[*trackIter], eVals, eVecs);
159  double evalRatio;
160  if (eVals[0] < 0.0001)
161  evalRatio = 0.0;
162  else
163  evalRatio = std::sqrt(eVals[1] * eVals[1] + eVals[2] * eVals[2]) / eVals[0];
164  this->FitAndSortTrack(*trackIter, isStoppingReco, sortedObj);
165  double coreHaloRatio, concentration, conicalness;
166  this->_Var_Shape(sortedObj, coreHaloRatio, concentration, conicalness);
167  double dEdxStart = CalcSegmentdEdxFrac(clockData, detProp, sortedObj, 0., 0.05);
168  double dEdxEnd = CalcSegmentdEdxFrac(clockData, detProp, sortedObj, 0.9, 1.0);
169  double dEdxPenultimate = CalcSegmentdEdxFrac(clockData, detProp, sortedObj, 0.8, 0.9);
170 
171  fResHolder.isTrack = 1;
172  fResHolder.isStoppingReco = isStoppingReco;
173  fResHolder.nSpacePoints = sortedObj.hitMap.size();
174  fResHolder.trackID = (*trackIter)->ID();
175  fResHolder.evalRatio = evalRatio;
176  fResHolder.concentration = concentration;
177  fResHolder.coreHaloRatio = coreHaloRatio;
178  fResHolder.conicalness = conicalness;
179  fResHolder.dEdxStart = dEdxStart;
180  fResHolder.dEdxEnd = dEdxEnd;
181  if (dEdxPenultimate < 0.1)
182  fResHolder.dEdxEndRatio = 1.0;
183  else
184  fResHolder.dEdxEndRatio = dEdxEnd / dEdxPenultimate;
185  fResHolder.length = sortedObj.length;
186 
187  for (auto methodIter = fMVAMethods.begin(); methodIter != fMVAMethods.end(); ++methodIter) {
188  fResHolder.mvaOutput[*methodIter] = fReader.EvaluateMVA(*methodIter);
189  }
190  result.push_back(fResHolder);
191  util::CreateAssn(evt, result, *trackIter, trackAssns);
192  }
193 
194  for (auto showerIter = fShowers.begin(); showerIter != fShowers.end(); ++showerIter) {
195  mvapid::MVAAlg::SortedObj sortedObj;
196 
197  std::vector<double> eVals, eVecs;
198  int isStoppingReco;
199 
200  this->RunPCA(fShowersToHits[*showerIter], eVals, eVecs);
201 
202  double evalRatio;
203  if (eVals[0] < 0.0001)
204  evalRatio = 0.0;
205  else
206  evalRatio = std::sqrt(eVals[1] * eVals[1] + eVals[2] * eVals[2]) / eVals[0];
207 
208  this->SortShower(*showerIter, isStoppingReco, sortedObj);
209 
210  double coreHaloRatio, concentration, conicalness;
211  this->_Var_Shape(sortedObj, coreHaloRatio, concentration, conicalness);
212  double dEdxStart = CalcSegmentdEdxFrac(clockData, detProp, sortedObj, 0., 0.05);
213  double dEdxEnd = CalcSegmentdEdxFrac(clockData, detProp, sortedObj, 0.9, 1.0);
214  double dEdxPenultimate = CalcSegmentdEdxFrac(clockData, detProp, sortedObj, 0.8, 0.9);
215 
216  fResHolder.isTrack = 0;
217  fResHolder.isStoppingReco = isStoppingReco;
218  fResHolder.nSpacePoints = sortedObj.hitMap.size();
220  (*showerIter)->ID() + 1000; //For the moment label showers by adding 1000 to ID
221 
222  fResHolder.evalRatio = evalRatio;
223  fResHolder.concentration = concentration;
224  fResHolder.coreHaloRatio = coreHaloRatio;
225  fResHolder.conicalness = conicalness;
226  fResHolder.dEdxStart = dEdxStart;
227  fResHolder.dEdxEnd = dEdxEnd;
228  if (dEdxPenultimate < 0.1)
229  fResHolder.dEdxEndRatio = 1.0;
230  else
231  fResHolder.dEdxEndRatio = dEdxEnd / dEdxPenultimate;
232  fResHolder.length = sortedObj.length;
233 
234  for (auto methodIter = fMVAMethods.begin(); methodIter != fMVAMethods.end(); ++methodIter) {
235  fResHolder.mvaOutput[*methodIter] = fReader.EvaluateMVA(*methodIter);
236  }
237  result.push_back(fResHolder);
238  util::CreateAssn(evt, result, *showerIter, showerAssns);
239  }
240 }
241 
243  const detinfo::DetectorClocksData& clockData)
244 {
245 
246  fHits.clear();
247  fSpacePoints.clear();
248  fTracks.clear();
249  fShowers.clear();
250  fSpacePointsToHits.clear();
251  fHitsToSpacePoints.clear();
252  fTracksToHits.clear();
253  fTracksToSpacePoints.clear();
254  fShowersToHits.clear();
255  fShowersToSpacePoints.clear();
256 
257  fEventT0 = trigger_offset(clockData);
258 
260  evt.getByLabel(fHitLabel, hitsHandle);
261 
262  for (unsigned int iHit = 0; iHit < hitsHandle->size(); ++iHit) {
263  const art::Ptr<recob::Hit> hit(hitsHandle, iHit);
264  fHits.push_back(hit);
265  }
266 
268  evt.getByLabel(fTrackLabel, tracksHandle);
269 
270  for (unsigned int iTrack = 0; iTrack < tracksHandle->size(); ++iTrack) {
271  const art::Ptr<recob::Track> track(tracksHandle, iTrack);
272  fTracks.push_back(track);
273  }
274 
276  evt.getByLabel(fShowerLabel, showersHandle);
277 
278  for (unsigned int iShower = 0; iShower < showersHandle->size(); ++iShower) {
279  const art::Ptr<recob::Shower> shower(showersHandle, iShower);
280  fShowers.push_back(shower);
281  }
282 
284  evt.getByLabel(fSpacePointLabel, spHandle);
285 
286  for (unsigned int iSP = 0; iSP < spHandle->size(); ++iSP) {
287  const art::Ptr<recob::SpacePoint> spacePoint(spHandle, iSP);
288  fSpacePoints.push_back(spacePoint);
289  }
290 
291  art::FindManyP<recob::Hit> findTracksToHits(fTracks, evt, fTrackLabel);
292  art::FindManyP<recob::Hit> findShowersToHits(fShowers, evt, fShowerLabel);
294 
295  for (unsigned int iSP = 0; iSP < fSpacePoints.size(); ++iSP) {
296  const art::Ptr<recob::SpacePoint> spacePoint = fSpacePoints.at(iSP);
297 
298  const art::Ptr<recob::Hit> hit = findSPToHits.at(iSP);
299  fSpacePointsToHits[spacePoint] = hit;
300  fHitsToSpacePoints[hit] = spacePoint;
301  }
302 
303  for (unsigned int iTrack = 0; iTrack < fTracks.size(); ++iTrack) {
304  const art::Ptr<recob::Track> track = fTracks.at(iTrack);
305 
306  const std::vector<art::Ptr<recob::Hit>> trackHits = findTracksToHits.at(iTrack);
307 
308  for (unsigned int iHit = 0; iHit < trackHits.size(); ++iHit) {
309  const art::Ptr<recob::Hit> hit = trackHits.at(iHit);
310  fTracksToHits[track].push_back(hit);
311  if (fHitsToSpacePoints.count(hit)) {
312  fTracksToSpacePoints[track].push_back(fHitsToSpacePoints.at(hit));
313  }
314  }
315  }
316 
317  for (unsigned int iShower = 0; iShower < fShowers.size(); ++iShower) {
318  const art::Ptr<recob::Shower> shower = fShowers.at(iShower);
319  const std::vector<art::Ptr<recob::Hit>> showerHits = findShowersToHits.at(iShower);
320 
321  for (unsigned int iHit = 0; iHit < showerHits.size(); ++iHit) {
322  const art::Ptr<recob::Hit> hit = showerHits.at(iHit);
323  fShowersToHits[shower].push_back(hit);
324  if (fHitsToSpacePoints.count(hit)) {
325  fShowersToSpacePoints[shower].push_back(fHitsToSpacePoints.at(hit));
326  }
327  }
328  }
329 
330  if (fCheatVertex) {
332  evt.getByLabel(fTrackingLabel, partHandle);
333 
334  if (partHandle->size() == 0 || partHandle->at(0).TrackId() != 1) {
335  std::cout << "Error, ID of first track in largeant list is not 0" << std::endl;
336  exit(1);
337  }
338  fVertex4Vect = partHandle->at(0).Position();
339  }
340 }
341 
343  int& isStoppingReco,
344  mvapid::MVAAlg::SortedObj& sortedTrack)
345 {
346 
347  sortedTrack.hitMap.clear();
348  TVector3 trackPoint, trackDir;
349  this->LinFit(track, trackPoint, trackDir);
350 
351  TVector3 nearestPointStart, nearestPointEnd;
352 
353  //For single-particle events can opt to cheat vertex from start of primary trajectory.
354  //Ok since in real events it should be possible to identify the true vertex.
355  if (fCheatVertex) {
356  if ((track->End<TVector3>() - fVertex4Vect.Vect()).Mag() >
357  (track->Vertex<TVector3>() - fVertex4Vect.Vect()).Mag()) {
358  nearestPointStart =
359  trackPoint +
360  trackDir * (trackDir.Dot(track->Vertex<TVector3>() - trackPoint) / trackDir.Mag2());
361  nearestPointEnd = trackPoint + trackDir * (trackDir.Dot(track->End<TVector3>() - trackPoint) /
362  trackDir.Mag2());
363  isStoppingReco = this->IsInActiveVol(track->End<TVector3>());
364  }
365  else {
366  nearestPointStart =
367  trackPoint +
368  trackDir * (trackDir.Dot(track->End<TVector3>() - trackPoint) / trackDir.Mag2());
369  nearestPointEnd =
370  trackPoint +
371  trackDir * (trackDir.Dot(track->Vertex<TVector3>() - trackPoint) / trackDir.Mag2());
372  isStoppingReco = this->IsInActiveVol(track->Vertex<TVector3>());
373  trackDir *= -1.;
374  }
375  }
376  else {
377  if (track->End<TVector3>().Z() >=
378  track->Vertex<TVector3>().Z()) { //Otherwise assume particle is forward-going for now...
379  nearestPointStart =
380  trackPoint +
381  trackDir * (trackDir.Dot(track->Vertex<TVector3>() - trackPoint) / trackDir.Mag2());
382  nearestPointEnd = trackPoint + trackDir * (trackDir.Dot(track->End<TVector3>() - trackPoint) /
383  trackDir.Mag2());
384  isStoppingReco = this->IsInActiveVol(track->End<TVector3>());
385  }
386  else {
387  nearestPointStart =
388  trackPoint +
389  trackDir * (trackDir.Dot(track->End<TVector3>() - trackPoint) / trackDir.Mag2());
390  nearestPointEnd =
391  trackPoint +
392  trackDir * (trackDir.Dot(track->Vertex<TVector3>() - trackPoint) / trackDir.Mag2());
393  isStoppingReco = this->IsInActiveVol(track->Vertex<TVector3>());
394  }
395 
396  if (trackDir.Z() <= 0) {
397  trackDir.SetX(-trackDir.X());
398  trackDir.SetY(-trackDir.Y());
399  trackDir.SetZ(-trackDir.Z());
400  }
401  }
402 
403  sortedTrack.start = nearestPointStart;
404  sortedTrack.end = nearestPointEnd;
405  sortedTrack.dir = trackDir;
406  sortedTrack.length = (nearestPointEnd - nearestPointStart).Mag();
407 
408  std::vector<art::Ptr<recob::Hit>> hits = fTracksToHits[track];
409 
410  for (auto hitIter = hits.begin(); hitIter != hits.end(); ++hitIter) {
411 
412  if (!fHitsToSpacePoints.count(*hitIter)) continue;
414 
415  TVector3 nearestPoint =
416  trackPoint + trackDir * (trackDir.Dot(TVector3(sp->XYZ()) - trackPoint) / trackDir.Mag2());
417  double lengthAlongTrack = (nearestPointStart - nearestPoint).Mag();
418  sortedTrack.hitMap.insert(std::pair<double, art::Ptr<recob::Hit>>(lengthAlongTrack, *hitIter));
419  }
420 }
421 
422 //void mvapid::MVAAlg::SortShower(art::Ptr<recob::Shower> shower,TVector3 dir,int& isStoppingReco,
423 // mvapid::MVAAlg::SortedObj& sortedShower){
425  int& isStoppingReco,
426  mvapid::MVAAlg::SortedObj& sortedShower)
427 {
428  sortedShower.hitMap.clear();
429 
430  std::vector<art::Ptr<recob::Hit>> hits = fShowersToHits[shower];
431 
432  TVector3 showerEnd(0, 0, 0);
433  double furthestHitFromStart = -999.9;
434  for (auto hitIter = hits.begin(); hitIter != hits.end(); ++hitIter) {
435 
436  if (!fHitsToSpacePoints.count(*hitIter)) continue;
438  if ((TVector3(sp->XYZ()) - shower->ShowerStart()).Mag() > furthestHitFromStart) {
439  showerEnd = TVector3(sp->XYZ());
440  furthestHitFromStart = (TVector3(sp->XYZ()) - shower->ShowerStart()).Mag();
441  }
442  }
443 
444  TVector3 showerPoint, showerDir;
445  this->LinFitShower(shower, showerPoint, showerDir);
446 
447  TVector3 nearestPointStart, nearestPointEnd;
448 
449  //Ensure that shower is fitted in correct direction (assuming for now that particle moves in +z direction)
450 
451  if (fCheatVertex) {
452  if ((showerEnd - fVertex4Vect.Vect()).Mag() >
453  (shower->ShowerStart() - fVertex4Vect.Vect()).Mag()) {
454  nearestPointStart =
455  showerPoint +
456  showerDir * (showerDir.Dot(shower->ShowerStart() - showerPoint) / showerDir.Mag2());
457  nearestPointEnd =
458  showerPoint + showerDir * (showerDir.Dot(showerEnd - showerPoint) / showerDir.Mag2());
459  isStoppingReco = this->IsInActiveVol(showerEnd);
460  }
461  else {
462  nearestPointStart =
463  showerPoint + showerDir * (showerDir.Dot(showerEnd - showerPoint) / showerDir.Mag2());
464  nearestPointEnd =
465  showerPoint +
466  showerDir * (showerDir.Dot(shower->ShowerStart() - showerPoint) / showerDir.Mag2());
467  isStoppingReco = this->IsInActiveVol(shower->ShowerStart());
468  showerDir *= -1.;
469  }
470  }
471  else {
472  if (showerEnd.Z() >= shower->ShowerStart().Z()) {
473  nearestPointStart =
474  showerPoint +
475  showerDir * (showerDir.Dot(shower->ShowerStart() - showerPoint) / showerDir.Mag2());
476  nearestPointEnd =
477  showerPoint + showerDir * (showerDir.Dot(showerEnd - showerPoint) / showerDir.Mag2());
478  isStoppingReco = this->IsInActiveVol(showerEnd);
479  }
480  else {
481  nearestPointStart =
482  showerPoint + showerDir * (showerDir.Dot(showerEnd - showerPoint) / showerDir.Mag2());
483  nearestPointEnd =
484  showerPoint +
485  showerDir * (showerDir.Dot(shower->ShowerStart() - showerPoint) / showerDir.Mag2());
486  isStoppingReco = this->IsInActiveVol(shower->ShowerStart());
487  }
488 
489  if (showerDir.Z() <= 0) {
490  showerDir.SetX(-showerDir.X());
491  showerDir.SetY(-showerDir.Y());
492  showerDir.SetZ(-showerDir.Z());
493  }
494  }
495 
496  sortedShower.start = nearestPointStart;
497  sortedShower.end = nearestPointEnd;
498  sortedShower.dir = showerDir;
499  sortedShower.length = (nearestPointEnd - nearestPointStart).Mag();
500 
501  for (auto hitIter = hits.begin(); hitIter != hits.end(); ++hitIter) {
502 
503  if (!fHitsToSpacePoints.count(*hitIter)) continue;
505 
506  TVector3 nearestPoint =
507  showerPoint +
508  showerDir * (showerDir.Dot(TVector3(sp->XYZ()) - showerPoint) / showerDir.Mag2());
509  double lengthAlongShower = (nearestPointStart - nearestPoint).Mag();
510  sortedShower.hitMap.insert(
511  std::pair<double, art::Ptr<recob::Hit>>(lengthAlongShower, *hitIter));
512  }
513 }
515  std::vector<double>& eVals,
516  std::vector<double>& eVecs)
517 {
518  TPrincipal* principal = new TPrincipal(3, "D");
519 
520  for (auto hitIter = hits.begin(); hitIter != hits.end(); ++hitIter) {
521 
522  if (fHitsToSpacePoints.count(*hitIter)) {
523  principal->AddRow(fHitsToSpacePoints.at(*hitIter)->XYZ());
524  }
525  }
526 
527  // PERFORM PCA
528  principal->MakePrincipals();
529  // GET EIGENVALUES AND EIGENVECTORS
530  for (unsigned int i = 0; i < 3; ++i) {
531  eVals.push_back(principal->GetEigenValues()->GetMatrixArray()[i]);
532  }
533 
534  for (unsigned int i = 0; i < 9; ++i) {
535  eVecs.push_back(principal->GetEigenVectors()->GetMatrixArray()[i]);
536  }
537 }
539  double& coreHaloRatio,
540  double& concentration,
541  double& conicalness)
542 {
543 
544  static const unsigned int conMinHits = 3;
545  static const double minCharge = 0.1;
546  static const double conFracRange = 0.2;
547  static const double MoliereRadius = 10.1;
548  static const double MoliereRadiusFraction = 0.2;
549 
550  double totalCharge = 0;
551  double totalChargeStart = 0;
552  double totalChargeEnd = 0;
553 
554  double chargeCore = 0;
555  double chargeHalo = 0;
556  double chargeCon = 0;
557  unsigned int nHits = 0;
558 
559  //stuff for conicalness
560  double chargeConStart = 0;
561  double chargeConEnd = 0;
562  unsigned int nHitsConStart = 0;
563  unsigned int nHitsConEnd = 0;
564 
565  for (auto hitIter = track.hitMap.begin(); hitIter != track.hitMap.end(); ++hitIter) {
566  if (fHitsToSpacePoints.count(hitIter->second)) {
567  art::Ptr<recob::SpacePoint> sp = fHitsToSpacePoints.at(hitIter->second);
568 
569  double distFromTrackFit = ((TVector3(sp->XYZ()) - track.start).Cross(track.dir)).Mag();
570 
571  ++nHits;
572 
573  if (distFromTrackFit < MoliereRadiusFraction * MoliereRadius)
574  chargeCore += hitIter->second->Integral();
575  else
576  chargeHalo += hitIter->second->Integral();
577 
578  totalCharge += hitIter->second->Integral();
579 
580  chargeCon += hitIter->second->Integral() / std::max(1.E-2, distFromTrackFit);
581  if (hitIter->first / track.length < conFracRange) {
582  chargeConStart += distFromTrackFit * distFromTrackFit * hitIter->second->Integral();
583  ++nHitsConStart;
584  totalChargeStart += hitIter->second->Integral();
585  }
586  else if (1. - hitIter->first / track.length < conFracRange) {
587  chargeConEnd += distFromTrackFit * distFromTrackFit * hitIter->second->Integral();
588  ++nHitsConEnd;
589  totalChargeEnd += hitIter->second->Integral();
590  }
591  }
592  }
593 
594  coreHaloRatio = chargeHalo / TMath::Max(1.0E-3, chargeCore);
595  coreHaloRatio = TMath::Min(100.0, coreHaloRatio);
596  concentration = chargeCon / totalCharge;
597  if (nHitsConStart >= conMinHits && nHitsConEnd >= conMinHits && totalChargeEnd > minCharge &&
598  sqrt(chargeConStart) > minCharge && totalChargeStart > minCharge) {
599  conicalness = (sqrt(chargeConEnd) / totalChargeEnd) / (sqrt(chargeConStart) / totalChargeStart);
600  }
601  else {
602  conicalness = 1.;
603  }
604 }
605 
607  const detinfo::DetectorPropertiesData& det_prop,
609  double start,
610  double end)
611 {
612 
613  double trackLength = (track.end - track.start).Mag();
614  return CalcSegmentdEdxDist(clock_data, det_prop, track, start * trackLength, end * trackLength);
615 }
616 
618  const detinfo::DetectorPropertiesData& det_prop,
620  double distAtEnd)
621 {
622 
623  double trackLength = (track.end - track.start).Mag();
624  return CalcSegmentdEdxDist(clock_data, det_prop, track, trackLength - distAtEnd, trackLength);
625 }
626 
628  const detinfo::DetectorPropertiesData& det_prop,
630  double start,
631  double end)
632 {
634 
635  double totaldEdx = 0;
636  unsigned int nHits = 0;
637 
638  //Loop over hits again to calculate average dE/dx and shape variables
639  for (auto hitIter = track.hitMap.begin(); hitIter != track.hitMap.end(); ++hitIter) {
640 
641  if (hitIter->first < start) continue;
642  if (hitIter->first >= end) break;
643 
644  art::Ptr<recob::Hit> hit = hitIter->second;
645 
646  //Pitch to use in dEdx calculation
647  double yzPitch = geom->WirePitch(
648  hit->WireID().asPlaneID()); //pitch not taking into account angle of track or shower
649  double xComponent, pitch3D;
650 
651  TVector3 dir = track.dir;
652 
653  //This assumes equal numbers of TPCs in each cryostat and equal numbers of planes in each TPC
654  int planeKey = hit->WireID().Cryostat * geom->NTPC() * geom->Nplanes() +
655  hit->WireID().TPC * geom->Nplanes() + hit->WireID().Plane;
656 
657  if (fNormToWiresY.count(planeKey) && fNormToWiresZ.count(planeKey)) {
658  TVector3 normToWires(0.0, fNormToWiresY.at(planeKey), fNormToWiresZ.at(planeKey));
659  yzPitch = geom->WirePitch(hit->WireID().asPlaneID()) / fabs(dir.Dot(normToWires));
660  }
661 
662  xComponent = yzPitch * dir[0] / sqrt(dir[1] * dir[1] + dir[2] * dir[2]);
663  pitch3D = sqrt(xComponent * xComponent + yzPitch * yzPitch);
664 
665  double dEdx = fCaloAlg.dEdx_AREA(clock_data, det_prop, *hit, pitch3D, fEventT0);
666  if (dEdx < 50.) {
667  ++nHits;
668  totaldEdx += dEdx;
669  }
670  }
671 
672  return nHits ? totaldEdx / nHits : 0;
673 }
674 
676  TVector3& trackPoint,
677  TVector3& trackDir)
678 {
679 
680  const std::vector<art::Ptr<recob::SpacePoint>>& sp = fTracksToSpacePoints.at(track);
681 
682  TGraph2D grFit(1);
683  unsigned int iPt = 0;
684  for (auto spIter = sp.begin(); spIter != sp.end(); ++spIter) {
685  TVector3 point = (*spIter)->XYZ();
686  grFit.SetPoint(iPt++, point.X(), point.Y(), point.Z());
687  }
688 
689  //Lift from the ROOT line3Dfit.C tutorial
690  ROOT::Fit::Fitter fitter;
691  // make the functor object
692  mvapid::MVAAlg::SumDistance2 sdist(&grFit);
693 
694  ROOT::Math::Functor fcn(sdist, 6);
695 
696  //Initial fit parameters from track start and end...
697  TVector3 trackStart = track->Vertex<TVector3>();
698  TVector3 trackEnd = track->End<TVector3>();
699  trackDir = (trackEnd - trackStart).Unit();
700 
701  TVector3 x0 = trackStart - trackDir;
702  TVector3 u = trackDir;
703 
704  double pStart[6] = {x0.X(), u.X(), x0.Y(), u.Y(), x0.Z(), u.Z()};
705 
706  fitter.SetFCN(fcn, pStart);
707 
708  bool ok = fitter.FitFCN();
709  if (!ok) {
710  trackPoint.SetXYZ(x0.X(), x0.Y(), x0.Z());
711  trackDir.SetXYZ(u.X(), u.Y(), u.Z());
712  trackDir = trackDir.Unit();
713  return 1;
714  }
715  else {
716  const ROOT::Fit::FitResult& result = fitter.Result();
717  const double* parFit = result.GetParams();
718  trackPoint.SetXYZ(parFit[0], parFit[2], parFit[4]);
719  trackDir.SetXYZ(parFit[1], parFit[3], parFit[5]);
720  trackDir = trackDir.Unit();
721  return 0;
722  }
723 }
724 
726  TVector3& showerPoint,
727  TVector3& showerDir)
728 {
729 
730  const std::vector<art::Ptr<recob::SpacePoint>>& sp = fShowersToSpacePoints.at(shower);
731 
732  TGraph2D grFit(1);
733  unsigned int iPt = 0;
734  for (auto spIter = sp.begin(); spIter != sp.end(); ++spIter) {
735  TVector3 point = (*spIter)->XYZ();
736  grFit.SetPoint(iPt++, point.X(), point.Y(), point.Z());
737  }
738 
739  //Lift from the ROOT line3Dfit.C tutorial
740  ROOT::Fit::Fitter fitter;
741  // make the functor object
742  mvapid::MVAAlg::SumDistance2 sdist(&grFit);
743 
744  ROOT::Math::Functor fcn(sdist, 6);
745 
746  //Initial fit parameters from shower start and end...
747  TVector3 showerStart = shower->ShowerStart();
748  showerDir = shower->Direction().Unit();
749 
750  TVector3 x0 = showerStart - showerDir;
751  TVector3 u = showerDir;
752 
753  double pStart[6] = {x0.X(), u.X(), x0.Y(), u.Y(), x0.Z(), u.Z()};
754 
755  fitter.SetFCN(fcn, pStart);
756 
757  bool ok = fitter.FitFCN();
758  if (!ok) {
759  showerPoint.SetXYZ(x0.X(), x0.Y(), x0.Z());
760  showerDir.SetXYZ(u.X(), u.Y(), u.Z());
761  showerDir = showerDir.Unit();
762  return 1;
763  }
764  else {
765  const ROOT::Fit::FitResult& result = fitter.Result();
766  const double* parFit = result.GetParams();
767  showerPoint.SetXYZ(parFit[0], parFit[2], parFit[4]);
768  showerDir.SetXYZ(parFit[1], parFit[3], parFit[5]);
769  showerDir = showerDir.Unit();
770  return 0;
771  }
772 }
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
std::vector< art::Ptr< recob::Shower > > fShowers
Definition: MVAAlg.h:151
const TVector3 & ShowerStart() const
Definition: Shower.h:197
unsigned int NTPC(CryostatID const &cryoid=cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
Definition: GeometryCore.h:686
anab::MVAPIDResult fResHolder
Definition: MVAAlg.h:163
double CalcSegmentdEdxDistAtEnd(const detinfo::DetectorClocksData &clock_data, const detinfo::DetectorPropertiesData &det_prop, const mvapid::MVAAlg::SortedObj &track, double distAtEnd)
Definition: MVAAlg.cxx:617
const calo::CalorimetryAlg fCaloAlg
Definition: MVAAlg.h:135
int LinFit(const art::Ptr< recob::Track > track, TVector3 &trackPoint, TVector3 &trackDir)
Definition: MVAAlg.cxx:675
void RunPCA(std::vector< art::Ptr< recob::Hit >> &hits, std::vector< double > &eVals, std::vector< double > &eVecs)
Definition: MVAAlg.cxx:514
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::SpacePoint > > fHitsToSpacePoints
Definition: MVAAlg.h:160
double fDetMaxX
Definition: MVAAlg.h:139
std::map< std::string, double > mvaOutput
Definition: MVAPIDResult.h:27
std::vector< std::string > fMVAMethods
Definition: MVAAlg.h:167
Declaration of signal hit object.
Geometry information for a single TPC.
Definition: TPCGeo.h:36
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:211
int IsInActiveVol(const TVector3 &pos)
Definition: MVAAlg.cxx:77
std::map< art::Ptr< recob::Track >, std::vector< art::Ptr< recob::Hit > > > fTracksToHits
Definition: MVAAlg.h:155
geo::WireID const & WireID() const
Initial tdc tick for hit.
Definition: Hit.h:280
std::string fSpacePointLabel
Definition: MVAAlg.h:147
void PrepareEvent(const art::Event &event, const detinfo::DetectorClocksData &clockData)
Definition: MVAAlg.cxx:242
MVAAlg(fhicl::ParameterSet const &pset)
Definition: MVAAlg.cxx:33
TMVA::Reader fReader
Definition: MVAAlg.h:165
Particle class.
bool fCheatVertex
Definition: MVAAlg.h:170
void RunPID(art::Event &evt, std::vector< anab::MVAPIDResult > &result, art::Assns< recob::Track, anab::MVAPIDResult, void > &trackAssns, art::Assns< recob::Shower, anab::MVAPIDResult, void > &showerAssns)
Definition: MVAAlg.cxx:143
std::vector< art::Ptr< recob::Hit > > fHits
Definition: MVAAlg.h:153
int LinFitShower(const art::Ptr< recob::Shower > shower, TVector3 &showerPoint, TVector3 &showerDir)
Definition: MVAAlg.cxx:725
Float_t Z
Definition: plot.C:37
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
void SortShower(art::Ptr< recob::Shower > shower, int &isStoppingReco, mvapid::MVAAlg::SortedObj &sortedShower)
Definition: MVAAlg.cxx:424
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
void GetDetectorEdges()
Definition: MVAAlg.cxx:89
void hits()
Definition: readHits.C:15
Float_t E
Definition: plot.C:20
double fDetMinY
Definition: MVAAlg.h:139
parameter set interface
std::map< double, const art::Ptr< recob::Hit > > hitMap
Definition: MVAAlg.h:49
std::map< int, double > fNormToWiresY
Definition: MVAAlg.h:141
double fDetMinZ
Definition: MVAAlg.h:139
T get(std::string const &key) const
Definition: ParameterSet.h:314
void GetWireNormals()
Definition: MVAAlg.cxx:111
std::vector< art::Ptr< recob::SpacePoint > > fSpacePoints
Definition: MVAAlg.h:152
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:78
Point_t const & Vertex() const
Access to track position at different points.
Definition: Track.h:158
double fDetMinX
Definition: MVAAlg.h:139
Provides recob::Track data product.
const Double32_t * XYZ() const
Definition: SpacePoint.h:78
TLorentzVector fVertex4Vect
Definition: MVAAlg.h:172
const TVector3 & Direction() const
Definition: Shower.h:188
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
Definition: geo_types.h:520
std::map< art::Ptr< recob::Shower >, std::vector< art::Ptr< recob::Hit > > > fShowersToHits
Definition: MVAAlg.h:157
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:481
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Definition: PFPUtils.cxx:2675
Detector simulation of raw signals on wires.
std::string fTrackingLabel
Definition: MVAAlg.h:148
double CalcSegmentdEdxFrac(const detinfo::DetectorClocksData &clock_data, const detinfo::DetectorPropertiesData &det_prop, const SortedObj &track, double start, double end)
Definition: MVAAlg.cxx:606
std::string fTrackLabel
Definition: MVAAlg.h:144
bool CreateAssn(art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t index=UINT_MAX)
Creates a single one-to-one association.
std::string fHitLabel
Definition: MVAAlg.h:146
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
Definition: StdUtils.h:120
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Utility object to perform functions of association.
Contains all timing reference information for the detector.
TDirectory * dir
Definition: macro.C:5
std::map< art::Ptr< recob::Track >, std::vector< art::Ptr< recob::SpacePoint > > > fTracksToSpacePoints
Definition: MVAAlg.h:156
double dEdx_AREA(detinfo::DetectorClocksData const &clock_data, detinfo::DetectorPropertiesData const &det_prop, recob::Hit const &hit, double pitch, double T0=0) const
void _Var_Shape(const SortedObj &track, double &coreHaloRatio, double &concentration, double &conicalness)
Definition: MVAAlg.cxx:538
std::vector< art::Ptr< recob::Track > > fTracks
Definition: MVAAlg.h:150
int trigger_offset(DetectorClocksData const &data)
Point_t const & End() const
Access to track position at different points.
Definition: Track.h:159
std::string fShowerLabel
Definition: MVAAlg.h:145
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
Definition: GeometryCore.h:977
double fDetMaxZ
Definition: MVAAlg.h:139
std::map< art::Ptr< recob::Shower >, std::vector< art::Ptr< recob::SpacePoint > > > fShowersToSpacePoints
Definition: MVAAlg.h:159
TCEvent evt
Definition: DataStructs.cxx:8
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:399
std::map< art::Ptr< recob::SpacePoint >, art::Ptr< recob::Hit > > fSpacePointsToHits
Definition: MVAAlg.h:161
double fDetMaxY
Definition: MVAAlg.h:139
void FitAndSortTrack(art::Ptr< recob::Track > track, int &isStoppingReco, SortedObj &sortedObj)
Definition: MVAAlg.cxx:342
unsigned int trackID
Definition: MVAPIDResult.h:21
Float_t track
Definition: plot.C:35
std::map< int, double > fNormToWiresZ
Definition: MVAAlg.h:142
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
double fEventT0
Definition: MVAAlg.h:137
double CalcSegmentdEdxDist(const detinfo::DetectorClocksData &clock_data, const detinfo::DetectorPropertiesData &det_prop, const SortedObj &track, double start, double end)
Definition: MVAAlg.cxx:627
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< std::string > fWeightFiles
Definition: MVAAlg.h:168