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