21 #include "cetlib/search_path.h" 22 #include "cetlib_except/exception.h" 25 #include "Fit/Fitter.h" 26 #include "Math/Functor.h" 27 #include "TPrincipal.h" 34 : fCaloAlg(pset.
get<
fhicl::ParameterSet>(
"CalorimetryAlg")), fReader(
"")
53 std::vector<std::string> weightFileBnames = pset.
get<std::vector<std::string>>(
"WeightFiles");
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)) {
61 throw cet::exception(
"MVAPID") <<
"Unable to find weight file " << *fileIter
62 <<
" in search path " << searchPath.to_string() << std::endl;
68 std::cerr <<
"Mismatch in number of MVA methods and weight files!" << std::endl;
72 for (
unsigned int iMethod = 0; iMethod !=
fMVAMethods.size(); ++iMethod) {
79 const double fiducialDist = 5.0;
81 if (pos.X() > (
fDetMinX + fiducialDist) && pos.X() < (
fDetMaxX - fiducialDist) &&
82 pos.Y() > (
fDetMinY + fiducialDist) && pos.Y() < (
fDetMaxY - fiducialDist) &&
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);
137 std::make_pair(planeKey, -plane.Wire(0).Direction().Z()));
139 std::make_pair(planeKey, plane.Wire(0).Direction().Y()));
144 std::vector<anab::MVAPIDResult>& result,
153 for (
auto trackIter =
fTracks.begin(); trackIter !=
fTracks.end(); ++trackIter) {
156 std::vector<double> eVals, eVecs;
160 if (eVals[0] < 0.0001)
163 evalRatio = std::sqrt(eVals[1] * eVals[1] + eVals[2] * eVals[2]) / eVals[0];
165 double coreHaloRatio, concentration, conicalness;
166 this->
_Var_Shape(sortedObj, coreHaloRatio, concentration, conicalness);
181 if (dEdxPenultimate < 0.1)
194 for (
auto showerIter =
fShowers.begin(); showerIter !=
fShowers.end(); ++showerIter) {
197 std::vector<double> eVals, eVecs;
203 if (eVals[0] < 0.0001)
206 evalRatio = std::sqrt(eVals[1] * eVals[1] + eVals[2] * eVals[2]) / eVals[0];
208 this->
SortShower(*showerIter, isStoppingReco, sortedObj);
210 double coreHaloRatio, concentration, conicalness;
211 this->
_Var_Shape(sortedObj, coreHaloRatio, concentration, conicalness);
220 (*showerIter)->ID() + 1000;
228 if (dEdxPenultimate < 0.1)
262 for (
unsigned int iHit = 0; iHit < hitsHandle->size(); ++iHit) {
264 fHits.push_back(hit);
270 for (
unsigned int iTrack = 0; iTrack < tracksHandle->size(); ++iTrack) {
278 for (
unsigned int iShower = 0; iShower < showersHandle->size(); ++iShower) {
286 for (
unsigned int iSP = 0; iSP < spHandle->size(); ++iSP) {
295 for (
unsigned int iSP = 0; iSP <
fSpacePoints.size(); ++iSP) {
303 for (
unsigned int iTrack = 0; iTrack <
fTracks.size(); ++iTrack) {
306 const std::vector<art::Ptr<recob::Hit>> trackHits = findTracksToHits.at(iTrack);
308 for (
unsigned int iHit = 0; iHit < trackHits.size(); ++iHit) {
317 for (
unsigned int iShower = 0; iShower <
fShowers.size(); ++iShower) {
319 const std::vector<art::Ptr<recob::Hit>> showerHits = findShowersToHits.at(iShower);
321 for (
unsigned int iHit = 0; iHit < showerHits.size(); ++iHit) {
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;
347 sortedTrack.
hitMap.clear();
348 TVector3 trackPoint, trackDir;
349 this->
LinFit(track, trackPoint, trackDir);
351 TVector3 nearestPointStart, nearestPointEnd;
360 trackDir * (trackDir.Dot(track->
Vertex<TVector3>() - trackPoint) / trackDir.Mag2());
361 nearestPointEnd = trackPoint + trackDir * (trackDir.Dot(track->
End<TVector3>() - trackPoint) /
368 trackDir * (trackDir.Dot(track->
End<TVector3>() - trackPoint) / trackDir.Mag2());
371 trackDir * (trackDir.Dot(track->
Vertex<TVector3>() - trackPoint) / trackDir.Mag2());
377 if (track->
End<TVector3>().Z() >=
378 track->
Vertex<TVector3>().
Z()) {
381 trackDir * (trackDir.Dot(track->
Vertex<TVector3>() - trackPoint) / trackDir.Mag2());
382 nearestPointEnd = trackPoint + trackDir * (trackDir.Dot(track->
End<TVector3>() - trackPoint) /
389 trackDir * (trackDir.Dot(track->
End<TVector3>() - trackPoint) / trackDir.Mag2());
392 trackDir * (trackDir.Dot(track->
Vertex<TVector3>() - trackPoint) / trackDir.Mag2());
396 if (trackDir.Z() <= 0) {
397 trackDir.SetX(-trackDir.X());
398 trackDir.SetY(-trackDir.Y());
399 trackDir.SetZ(-trackDir.Z());
403 sortedTrack.
start = nearestPointStart;
404 sortedTrack.
end = nearestPointEnd;
405 sortedTrack.
dir = trackDir;
406 sortedTrack.
length = (nearestPointEnd - nearestPointStart).Mag();
410 for (
auto hitIter = hits.begin(); hitIter != hits.end(); ++hitIter) {
415 TVector3 nearestPoint =
416 trackPoint + trackDir * (trackDir.Dot(TVector3(sp->
XYZ()) - trackPoint) / trackDir.Mag2());
417 double lengthAlongTrack = (nearestPointStart - nearestPoint).Mag();
428 sortedShower.
hitMap.clear();
432 TVector3 showerEnd(0, 0, 0);
433 double furthestHitFromStart = -999.9;
434 for (
auto hitIter = hits.begin(); hitIter != hits.end(); ++hitIter) {
438 if ((TVector3(sp->
XYZ()) - shower->
ShowerStart()).Mag() > furthestHitFromStart) {
439 showerEnd = TVector3(sp->
XYZ());
440 furthestHitFromStart = (TVector3(sp->
XYZ()) - shower->
ShowerStart()).Mag();
444 TVector3 showerPoint, showerDir;
447 TVector3 nearestPointStart, nearestPointEnd;
456 showerDir * (showerDir.Dot(shower->
ShowerStart() - showerPoint) / showerDir.Mag2());
458 showerPoint + showerDir * (showerDir.Dot(showerEnd - showerPoint) / showerDir.Mag2());
463 showerPoint + showerDir * (showerDir.Dot(showerEnd - showerPoint) / showerDir.Mag2());
466 showerDir * (showerDir.Dot(shower->
ShowerStart() - showerPoint) / showerDir.Mag2());
475 showerDir * (showerDir.Dot(shower->
ShowerStart() - showerPoint) / showerDir.Mag2());
477 showerPoint + showerDir * (showerDir.Dot(showerEnd - showerPoint) / showerDir.Mag2());
482 showerPoint + showerDir * (showerDir.Dot(showerEnd - showerPoint) / showerDir.Mag2());
485 showerDir * (showerDir.Dot(shower->
ShowerStart() - showerPoint) / showerDir.Mag2());
489 if (showerDir.Z() <= 0) {
490 showerDir.SetX(-showerDir.X());
491 showerDir.SetY(-showerDir.Y());
492 showerDir.SetZ(-showerDir.Z());
496 sortedShower.
start = nearestPointStart;
497 sortedShower.
end = nearestPointEnd;
498 sortedShower.
dir = showerDir;
499 sortedShower.
length = (nearestPointEnd - nearestPointStart).Mag();
501 for (
auto hitIter = hits.begin(); hitIter != hits.end(); ++hitIter) {
506 TVector3 nearestPoint =
508 showerDir * (showerDir.Dot(TVector3(sp->
XYZ()) - showerPoint) / showerDir.Mag2());
509 double lengthAlongShower = (nearestPointStart - nearestPoint).Mag();
510 sortedShower.
hitMap.insert(
515 std::vector<double>& eVals,
516 std::vector<double>& eVecs)
518 TPrincipal* principal =
new TPrincipal(3,
"D");
520 for (
auto hitIter =
hits.begin(); hitIter !=
hits.end(); ++hitIter) {
528 principal->MakePrincipals();
530 for (
unsigned int i = 0; i < 3; ++i) {
531 eVals.push_back(principal->GetEigenValues()->GetMatrixArray()[i]);
534 for (
unsigned int i = 0; i < 9; ++i) {
535 eVecs.push_back(principal->GetEigenVectors()->GetMatrixArray()[i]);
539 double& coreHaloRatio,
540 double& concentration,
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;
550 double totalCharge = 0;
551 double totalChargeStart = 0;
552 double totalChargeEnd = 0;
554 double chargeCore = 0;
555 double chargeHalo = 0;
556 double chargeCon = 0;
557 unsigned int nHits = 0;
560 double chargeConStart = 0;
561 double chargeConEnd = 0;
562 unsigned int nHitsConStart = 0;
563 unsigned int nHitsConEnd = 0;
565 for (
auto hitIter = track.
hitMap.begin(); hitIter != track.
hitMap.end(); ++hitIter) {
569 double distFromTrackFit = ((TVector3(sp->
XYZ()) - track.
start).Cross(track.
dir)).Mag();
573 if (distFromTrackFit < MoliereRadiusFraction * MoliereRadius)
574 chargeCore += hitIter->second->Integral();
576 chargeHalo += hitIter->second->Integral();
578 totalCharge += hitIter->second->Integral();
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();
584 totalChargeStart += hitIter->second->Integral();
586 else if (1. - hitIter->first / track.
length < conFracRange) {
587 chargeConEnd += distFromTrackFit * distFromTrackFit * hitIter->second->Integral();
589 totalChargeEnd += hitIter->second->Integral();
594 coreHaloRatio = chargeHalo / TMath::Max(1.0
E-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);
613 double trackLength = (track.
end - track.
start).Mag();
614 return CalcSegmentdEdxDist(clock_data, det_prop, track, start * trackLength, end * trackLength);
623 double trackLength = (track.
end - track.
start).Mag();
624 return CalcSegmentdEdxDist(clock_data, det_prop, track, trackLength - distAtEnd, trackLength);
635 double totaldEdx = 0;
636 unsigned int nHits = 0;
639 for (
auto hitIter = track.
hitMap.begin(); hitIter != track.
hitMap.end(); ++hitIter) {
641 if (hitIter->first < start)
continue;
642 if (hitIter->first >= end)
break;
649 double xComponent, pitch3D;
662 xComponent = yzPitch * dir[0] / sqrt(dir[1] * dir[1] + dir[2] * dir[2]);
663 pitch3D = sqrt(xComponent * xComponent + yzPitch * yzPitch);
672 return nHits ? totaldEdx / nHits : 0;
676 TVector3& trackPoint,
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());
690 ROOT::Fit::Fitter fitter;
694 ROOT::Math::Functor fcn(sdist, 6);
697 TVector3 trackStart = track->
Vertex<TVector3>();
698 TVector3 trackEnd = track->
End<TVector3>();
699 trackDir = (trackEnd - trackStart).Unit();
701 TVector3 x0 = trackStart - trackDir;
702 TVector3 u = trackDir;
704 double pStart[6] = {x0.X(), u.X(), x0.Y(), u.Y(), x0.Z(), u.Z()};
706 fitter.SetFCN(fcn, pStart);
708 bool ok = fitter.FitFCN();
710 trackPoint.SetXYZ(x0.X(), x0.Y(), x0.Z());
711 trackDir.SetXYZ(u.X(), u.Y(), u.Z());
712 trackDir = trackDir.Unit();
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();
726 TVector3& showerPoint,
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());
740 ROOT::Fit::Fitter fitter;
744 ROOT::Math::Functor fcn(sdist, 6);
750 TVector3 x0 = showerStart - showerDir;
751 TVector3 u = showerDir;
753 double pStart[6] = {x0.X(), u.X(), x0.Y(), u.Y(), x0.Z(), u.Z()};
755 fitter.SetFCN(fcn, pStart);
757 bool ok = fitter.FitFCN();
759 showerPoint.SetXYZ(x0.X(), x0.Y(), x0.Z());
760 showerDir.SetXYZ(u.X(), u.Y(), u.Z());
761 showerDir = showerDir.Unit();
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();
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
std::vector< art::Ptr< recob::Shower > > fShowers
const TVector3 & ShowerStart() const
unsigned int NTPC(CryostatID const &cryoid=cryostat_zero) const
Returns the total number of TPCs in the specified cryostat.
anab::MVAPIDResult fResHolder
double CalcSegmentdEdxDistAtEnd(const detinfo::DetectorClocksData &clock_data, const detinfo::DetectorPropertiesData &det_prop, const mvapid::MVAAlg::SortedObj &track, double distAtEnd)
const calo::CalorimetryAlg fCaloAlg
int LinFit(const art::Ptr< recob::Track > track, TVector3 &trackPoint, TVector3 &trackDir)
void RunPCA(std::vector< art::Ptr< recob::Hit >> &hits, std::vector< double > &eVals, std::vector< double > &eVecs)
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::SpacePoint > > fHitsToSpacePoints
std::map< std::string, double > mvaOutput
std::vector< std::string > fMVAMethods
Declaration of signal hit object.
Geometry information for a single TPC.
CryostatID_t Cryostat
Index of cryostat.
int IsInActiveVol(const TVector3 &pos)
std::map< art::Ptr< recob::Track >, std::vector< art::Ptr< recob::Hit > > > fTracksToHits
geo::WireID const & WireID() const
Initial tdc tick for hit.
std::string fSpacePointLabel
void PrepareEvent(const art::Event &event, const detinfo::DetectorClocksData &clockData)
MVAAlg(fhicl::ParameterSet const &pset)
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)
std::vector< art::Ptr< recob::Hit > > fHits
int LinFitShower(const art::Ptr< recob::Shower > shower, TVector3 &showerPoint, TVector3 &showerDir)
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
void SortShower(art::Ptr< recob::Shower > shower, int &isStoppingReco, mvapid::MVAAlg::SortedObj &sortedShower)
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
std::map< double, const art::Ptr< recob::Hit > > hitMap
std::map< int, double > fNormToWiresY
T get(std::string const &key) const
std::vector< art::Ptr< recob::SpacePoint > > fSpacePoints
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Point_t const & Vertex() const
Access to track position at different points.
Provides recob::Track data product.
const Double32_t * XYZ() const
TLorentzVector fVertex4Vect
const TVector3 & Direction() const
constexpr PlaneID const & asPlaneID() const
Conversion to PlaneID (for convenience of notation).
std::map< art::Ptr< recob::Shower >, std::vector< art::Ptr< recob::Hit > > > fShowersToHits
PlaneID_t Plane
Index of the plane within its TPC.
float dEdx(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, const TCSlice &slc, TP3D &tp3d)
Detector simulation of raw signals on wires.
std::string fTrackingLabel
double CalcSegmentdEdxFrac(const detinfo::DetectorClocksData &clock_data, const detinfo::DetectorPropertiesData &det_prop, const SortedObj &track, double start, double end)
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.
decltype(auto) get(T &&obj)
ADL-aware version of std::to_string.
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.
std::map< art::Ptr< recob::Track >, std::vector< art::Ptr< recob::SpacePoint > > > fTracksToSpacePoints
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)
std::vector< art::Ptr< recob::Track > > fTracks
int trigger_offset(DetectorClocksData const &data)
Point_t const & End() const
Access to track position at different points.
unsigned int Nplanes(TPCID const &tpcid=tpc_zero) const
Returns the total number of planes in the specified TPC.
std::map< art::Ptr< recob::Shower >, std::vector< art::Ptr< recob::SpacePoint > > > fShowersToSpacePoints
TPCID_t TPC
Index of the TPC within its cryostat.
std::map< art::Ptr< recob::SpacePoint >, art::Ptr< recob::Hit > > fSpacePointsToHits
void FitAndSortTrack(art::Ptr< recob::Track > track, int &isStoppingReco, SortedObj &sortedObj)
std::map< int, double > fNormToWiresZ
Length_t WirePitch(PlaneID const &planeid=plane_zero) const
Returns the distance between two consecutive wires.
double CalcSegmentdEdxDist(const detinfo::DetectorClocksData &clock_data, const detinfo::DetectorPropertiesData &det_prop, const SortedObj &track, double start, double end)
art framework interface to geometry description
cet::coded_exception< error, detail::translate > exception
std::vector< std::string > fWeightFiles