LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
mvapid::MVAAlg Class Reference

#include "MVAAlg.h"

Classes

struct  SortedObj
 
struct  SumDistance2
 

Public Member Functions

 MVAAlg (fhicl::ParameterSet const &pset, const art::EDProducer *parentModule)
 
virtual ~MVAAlg ()
 
void GetDetectorEdges ()
 
void GetWireNormals ()
 
void reconfigure (fhicl::ParameterSet const &p)
 
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)
 

Private Member Functions

int IsInActiveVol (const TVector3 &pos)
 
void PrepareEvent (const art::Event &event)
 
void FitAndSortTrack (art::Ptr< recob::Track > track, int &isStoppingReco, SortedObj &sortedObj)
 
void SortShower (art::Ptr< recob::Shower > shower, int &isStoppingReco, mvapid::MVAAlg::SortedObj &sortedShower)
 
void RunPCA (std::vector< art::Ptr< recob::Hit > > &hits, std::vector< double > &eVals, std::vector< double > &eVecs)
 
void _Var_Shape (const SortedObj &track, double &coreHaloRatio, double &concentration, double &conicalness)
 
double CalcSegmentdEdxFrac (const SortedObj &track, double start, double end)
 
double CalcSegmentdEdxDist (const SortedObj &track, double start, double end)
 
double CalcSegmentdEdxDistAtEnd (const mvapid::MVAAlg::SortedObj &track, double distAtEnd)
 
int LinFit (const art::Ptr< recob::Track > track, TVector3 &trackPoint, TVector3 &trackDir)
 
int LinFitShower (const art::Ptr< recob::Shower > shower, TVector3 &showerPoint, TVector3 &showerDir)
 

Private Attributes

const calo::CalorimetryAlg fCaloAlg
 
double fEventT0
 
double fDetMinX
 
double fDetMaxX
 
double fDetMinY
 
double fDetMaxY
 
double fDetMinZ
 
double fDetMaxZ
 
std::map< int, double > fNormToWiresY
 
std::map< int, double > fNormToWiresZ
 
const art::EDProducerfParentModule
 
std::string fTrackLabel
 
std::string fShowerLabel
 
std::string fHitLabel
 
std::string fSpacePointLabel
 
std::string fTrackingLabel
 
std::vector< art::Ptr< recob::Track > > fTracks
 
std::vector< art::Ptr< recob::Shower > > fShowers
 
std::vector< art::Ptr< recob::SpacePoint > > fSpacePoints
 
std::vector< art::Ptr< recob::Hit > > fHits
 
std::map< art::Ptr< recob::Track >, std::vector< art::Ptr< recob::Hit > > > fTracksToHits
 
std::map< art::Ptr< recob::Track >, std::vector< art::Ptr< recob::SpacePoint > > > fTracksToSpacePoints
 
std::map< art::Ptr< recob::Shower >, std::vector< art::Ptr< recob::Hit > > > fShowersToHits
 
std::map< art::Ptr< recob::Shower >, std::vector< art::Ptr< recob::SpacePoint > > > fShowersToSpacePoints
 
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::SpacePoint > > fHitsToSpacePoints
 
std::map< art::Ptr< recob::SpacePoint >, art::Ptr< recob::Hit > > fSpacePointsToHits
 
anab::MVAPIDResult fResHolder
 
TMVA::Reader fReader
 
std::vector< std::string > fMVAMethods
 
std::vector< std::string > fWeightFiles
 
bool fCheatVertex
 
TLorentzVector fVertex4Vect
 

Detailed Description

Definition at line 47 of file MVAAlg.h.

Constructor & Destructor Documentation

mvapid::MVAAlg::MVAAlg ( fhicl::ParameterSet const &  pset,
const art::EDProducer parentModule 
)

Definition at line 28 of file MVAAlg.cxx.

References anab::MVAPIDResult::concentration, anab::MVAPIDResult::conicalness, anab::MVAPIDResult::coreHaloRatio, anab::MVAPIDResult::dEdxEnd, anab::MVAPIDResult::dEdxEndRatio, anab::MVAPIDResult::dEdxStart, anab::MVAPIDResult::evalRatio, fCheatVertex, fHitLabel, fMVAMethods, fReader, fResHolder, fShowerLabel, fSpacePointLabel, fTrackingLabel, fTrackLabel, fWeightFiles, and fhicl::ParameterSet::get().

Referenced by mvapid::MVAAlg::SumDistance2::operator()().

28  :
29  fCaloAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg")),
30  fParentModule(parentModule),fReader(""){
31  fHitLabel=pset.get<std::string>("HitLabel");
32  fTrackLabel=pset.get<std::string>("TrackLabel");
33  fShowerLabel=pset.get<std::string>("ShowerLabel");
34  fSpacePointLabel=pset.get<std::string>("SpacePointLabel");
35  fTrackingLabel=pset.get<std::string>("TrackingLabel","");
36 
37  fCheatVertex=pset.get<bool>("CheatVertex",false);
38 
39  fReader.AddVariable("evalRatio",&fResHolder.evalRatio);
40  fReader.AddVariable("coreHaloRatio",&fResHolder.coreHaloRatio);
41  fReader.AddVariable("concentration",&fResHolder.concentration);
42  fReader.AddVariable("conicalness",&fResHolder.conicalness);
43  fReader.AddVariable("dEdxStart",&fResHolder.dEdxStart);
44  fReader.AddVariable("dEdxEnd",&fResHolder.dEdxEnd);
45  fReader.AddVariable("dEdxEndRatio",&fResHolder.dEdxEndRatio);
46 
47  fMVAMethods=pset.get<std::vector<std::string> >("MVAMethods");
48  std::vector<std::string> weightFileBnames=pset.get<std::vector<std::string> >("WeightFiles");
49 
50  cet::search_path searchPath("FW_SEARCH_PATH");
51  for(auto fileIter=weightFileBnames.begin();fileIter!=weightFileBnames.end();++fileIter){
52  std::string fileWithPath;
53  if(!searchPath.find_file(*fileIter,fileWithPath)){
54  fWeightFiles.clear();
55  fMVAMethods.clear();
56  throw cet::exception("MVAPID")<<"Unable to find weight file "<<*fileIter<<" in search path "
57  <<searchPath.to_string()<<std::endl;
58  }
59  fWeightFiles.push_back(fileWithPath);
60  }
61 
62  if(fMVAMethods.size()!=fWeightFiles.size()){
63  std::cerr<<"Mismatch in number of MVA methods and weight files!"<<std::endl;
64  exit(1);
65  }
66 
67  for(unsigned int iMethod=0;iMethod!=fMVAMethods.size();++iMethod){
68  fReader.BookMVA(fMVAMethods[iMethod], fWeightFiles[iMethod]);
69  }
70 }
anab::MVAPIDResult fResHolder
Definition: MVAAlg.h:154
const calo::CalorimetryAlg fCaloAlg
Definition: MVAAlg.h:125
std::vector< std::string > fMVAMethods
Definition: MVAAlg.h:158
std::string fSpacePointLabel
Definition: MVAAlg.h:139
TMVA::Reader fReader
Definition: MVAAlg.h:156
bool fCheatVertex
Definition: MVAAlg.h:161
const art::EDProducer * fParentModule
Definition: MVAAlg.h:134
std::string fTrackingLabel
Definition: MVAAlg.h:140
std::string fTrackLabel
Definition: MVAAlg.h:136
std::string fHitLabel
Definition: MVAAlg.h:138
std::string fShowerLabel
Definition: MVAAlg.h:137
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
std::vector< std::string > fWeightFiles
Definition: MVAAlg.h:159
mvapid::MVAAlg::~MVAAlg ( )
virtual

Definition at line 72 of file MVAAlg.cxx.

Referenced by mvapid::MVAAlg::SumDistance2::operator()().

72 {}

Member Function Documentation

void mvapid::MVAAlg::_Var_Shape ( const SortedObj track,
double &  coreHaloRatio,
double &  concentration,
double &  conicalness 
)
private

Definition at line 539 of file MVAAlg.cxx.

References mvapid::MVAAlg::SortedObj::dir, E, fHitsToSpacePoints, mvapid::MVAAlg::SortedObj::hitMap, mvapid::MVAAlg::SortedObj::length, max, mvapid::MVAAlg::SortedObj::start, and recob::SpacePoint::XYZ().

Referenced by mvapid::MVAAlg::SumDistance2::operator()(), and RunPID().

541  {
542 
543  static const unsigned int conMinHits=3;
544  static const double minCharge = 0.1;
545  static const double conFracRange=0.2;
546  static const double MoliereRadius = 10.1;
547  static const double MoliereRadiusFraction = 0.2;
548 
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 
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&&sqrt(chargeConStart)>minCharge&&totalChargeStart>minCharge){
599  conicalness=(sqrt(chargeConEnd)/totalChargeEnd) / (sqrt(chargeConStart)/totalChargeStart);
600  }
601  else{
602  conicalness=1.;
603  }
604 }
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::SpacePoint > > fHitsToSpacePoints
Definition: MVAAlg.h:151
Float_t E
Definition: plot.C:23
Int_t max
Definition: plot.C:27
const Double32_t * XYZ() const
Definition: SpacePoint.h:65
Float_t track
Definition: plot.C:34
double mvapid::MVAAlg::CalcSegmentdEdxDist ( const SortedObj track,
double  start,
double  end 
)
private

Definition at line 618 of file MVAAlg.cxx.

References geo::CryostatID::Cryostat, calo::CalorimetryAlg::dEdx_AREA(), dir, mvapid::MVAAlg::SortedObj::dir, fCaloAlg, fEventT0, fNormToWiresY, fNormToWiresZ, mvapid::MVAAlg::SortedObj::hitMap, geo::GeometryCore::Nplanes(), geo::GeometryCore::NTPC(), geo::PlaneID::Plane, geo::TPCID::TPC, recob::Hit::WireID(), and geo::GeometryCore::WirePitch().

Referenced by CalcSegmentdEdxDistAtEnd(), CalcSegmentdEdxFrac(), and mvapid::MVAAlg::SumDistance2::operator()().

618  {
619 
621 
622  double totaldEdx=0;
623  unsigned int nHits=0;
624 
625  //Loop over hits again to calculate average dE/dx and shape variables
626  for ( auto hitIter = track.hitMap.begin(); hitIter!=track.hitMap.end(); ++hitIter ){
627 
628  if(hitIter->first<start) continue;
629  if(hitIter->first>=end) break;
630 
631  art::Ptr<recob::Hit> hit=hitIter->second;
632 
633  //Pitch to use in dEdx calculation
634  double yzPitch = geom->WirePitch(hit->WireID().Plane, hit->WireID().TPC); //pitch not taking into account angle of track or shower
635  double xComponent, pitch3D;
636 
637  TVector3 dir=track.dir;
638 
639  //This assumes equal numbers of TPCs in each cryostat and equal numbers of planes in each TPC
640  int planeKey = hit->WireID().Cryostat * geom->NTPC(0) * geom->Nplanes(0, 0) + hit->WireID().TPC * geom->Nplanes(0, 0) + hit->WireID().Plane;
641 
642  if(fNormToWiresY.count(planeKey) && fNormToWiresZ.count(planeKey))
643  {
644  TVector3 normToWires(0.0, fNormToWiresY.at(planeKey), fNormToWiresZ.at(planeKey));
645  yzPitch = geom->WirePitch(hit->WireID().Plane, hit->WireID().TPC) / fabs(dir.Dot(normToWires));
646  }
647 
648  xComponent = yzPitch * dir[0] / sqrt(dir[1] * dir[1] + dir[2] * dir[2]);
649  pitch3D = sqrt(xComponent * xComponent + yzPitch * yzPitch);
650 
651  double dEdx=fCaloAlg.dEdx_AREA(*hit, pitch3D, fEventT0);
652  if( dEdx < 50.){
653  ++nHits;
654  totaldEdx += dEdx;
655  }
656  }
657 
658  return nHits?totaldEdx/nHits:0;
659 }
const calo::CalorimetryAlg fCaloAlg
Definition: MVAAlg.h:125
double dEdx_AREA(art::Ptr< recob::Hit > hit, double pitch, double T0=0) const
geo::WireID WireID() const
Initial tdc tick for hit.
Definition: Hit.h:234
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:130
geo::Length_t WirePitch(geo::PlaneID const &planeid) const
Returns the distance between two consecutive wires.
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:258
Detector simulation of raw signals on wires.
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
TDirectory * dir
Definition: macro.C:5
std::map< int, double > fNormToWiresY
Definition: MVAAlg.h:131
std::vector< evd::details::RawDigitInfo_t >::const_iterator end(RawDigitCacheDataClass const &cache)
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:203
Float_t track
Definition: plot.C:34
double fEventT0
Definition: MVAAlg.h:127
std::map< int, double > fNormToWiresZ
Definition: MVAAlg.h:132
double mvapid::MVAAlg::CalcSegmentdEdxDistAtEnd ( const mvapid::MVAAlg::SortedObj track,
double  distAtEnd 
)
private

Definition at line 612 of file MVAAlg.cxx.

References CalcSegmentdEdxDist(), mvapid::MVAAlg::SortedObj::end, and mvapid::MVAAlg::SortedObj::start.

Referenced by mvapid::MVAAlg::SumDistance2::operator()().

612  {
613 
614  double trackLength=(track.end-track.start).Mag();
615  return CalcSegmentdEdxDist(track,trackLength-distAtEnd,trackLength);
616 }
double CalcSegmentdEdxDist(const SortedObj &track, double start, double end)
Definition: MVAAlg.cxx:618
double mvapid::MVAAlg::CalcSegmentdEdxFrac ( const SortedObj track,
double  start,
double  end 
)
private

Definition at line 606 of file MVAAlg.cxx.

References CalcSegmentdEdxDist(), mvapid::MVAAlg::SortedObj::end, and mvapid::MVAAlg::SortedObj::start.

Referenced by mvapid::MVAAlg::SumDistance2::operator()(), and RunPID().

606  {
607 
608  double trackLength=(track.end-track.start).Mag();
609  return CalcSegmentdEdxDist(track,start*trackLength,end*trackLength);
610 }
double CalcSegmentdEdxDist(const SortedObj &track, double start, double end)
Definition: MVAAlg.cxx:618
std::vector< evd::details::RawDigitInfo_t >::const_iterator end(RawDigitCacheDataClass const &cache)
Float_t track
Definition: plot.C:34
void mvapid::MVAAlg::FitAndSortTrack ( art::Ptr< recob::Track track,
int &  isStoppingReco,
SortedObj sortedObj 
)
private

Definition at line 374 of file MVAAlg.cxx.

References mvapid::MVAAlg::SortedObj::dir, mvapid::MVAAlg::SortedObj::end, recob::Track::End(), fCheatVertex, fHitsToSpacePoints, fTracksToHits, fVertex4Vect, mvapid::MVAAlg::SortedObj::hitMap, hits(), IsInActiveVol(), mvapid::MVAAlg::SortedObj::length, LinFit(), mvapid::MVAAlg::SortedObj::start, track, recob::Track::Vertex(), recob::SpacePoint::XYZ(), and Z.

Referenced by mvapid::MVAAlg::SumDistance2::operator()(), and RunPID().

375  {
376 
377  sortedTrack.hitMap.clear();
378  TVector3 trackPoint,trackDir;
379  this->LinFit(track,trackPoint,trackDir);
380 
381  TVector3 nearestPointStart, nearestPointEnd;
382 
383  //For single-particle events can opt to cheat vertex from start of primary trajectory.
384  //Ok since in real events it should be possible to identify the true vertex.
385  if(fCheatVertex){
386  if((track->End<TVector3>()-fVertex4Vect.Vect()).Mag()>(track->Vertex<TVector3>()-fVertex4Vect.Vect()).Mag()){
387  nearestPointStart = trackPoint+trackDir*(trackDir.Dot(track->Vertex<TVector3>()-trackPoint)/trackDir.Mag2());
388  nearestPointEnd = trackPoint+trackDir*(trackDir.Dot(track->End<TVector3>()-trackPoint)/trackDir.Mag2());
389  isStoppingReco = this->IsInActiveVol(track->End<TVector3>());
390  }
391  else{
392  nearestPointStart = trackPoint+trackDir*(trackDir.Dot(track->End<TVector3>()-trackPoint)/trackDir.Mag2());
393  nearestPointEnd = trackPoint+trackDir*(trackDir.Dot(track->Vertex<TVector3>()-trackPoint)/trackDir.Mag2());
394  isStoppingReco = this->IsInActiveVol(track->Vertex<TVector3>());
395  trackDir*=-1.;
396  }
397  }
398  else{
399  if(track->End<TVector3>().Z() >= track->Vertex<TVector3>().Z()){ //Otherwise assume particle is forward-going for now...
400  nearestPointStart = trackPoint+trackDir*(trackDir.Dot(track->Vertex<TVector3>()-trackPoint)/trackDir.Mag2());
401  nearestPointEnd = trackPoint+trackDir*(trackDir.Dot(track->End<TVector3>()-trackPoint)/trackDir.Mag2());
402  isStoppingReco = this->IsInActiveVol(track->End<TVector3>());
403  }
404  else{
405  nearestPointStart = trackPoint+trackDir*(trackDir.Dot(track->End<TVector3>()-trackPoint)/trackDir.Mag2());
406  nearestPointEnd = trackPoint+trackDir*(trackDir.Dot(track->Vertex<TVector3>()-trackPoint)/trackDir.Mag2());
407  isStoppingReco = this->IsInActiveVol(track->Vertex<TVector3>());
408  }
409 
410  if(trackDir.Z() <= 0){
411  trackDir.SetX(-trackDir.X());
412  trackDir.SetY(-trackDir.Y());
413  trackDir.SetZ(-trackDir.Z());
414  }
415  }
416 
417  sortedTrack.start=nearestPointStart;
418  sortedTrack.end=nearestPointEnd;
419  sortedTrack.dir=trackDir;
420  sortedTrack.length=(nearestPointEnd-nearestPointStart).Mag();
421 
422  std::vector<art::Ptr<recob::Hit>> hits=fTracksToHits[track];
423 
424  for(auto hitIter=hits.begin();hitIter!=hits.end();++hitIter){
425 
426  if(!fHitsToSpacePoints.count(*hitIter)) continue;
428 
429  TVector3 nearestPoint = trackPoint+trackDir*(trackDir.Dot(TVector3(sp->XYZ())-trackPoint)/trackDir.Mag2());
430  double lengthAlongTrack=(nearestPointStart-nearestPoint).Mag();
431  sortedTrack.hitMap.insert(std::pair<double,art::Ptr<recob::Hit> >(lengthAlongTrack,*hitIter));
432  }
433 }
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::SpacePoint > > fHitsToSpacePoints
Definition: MVAAlg.h:151
int LinFit(const art::Ptr< recob::Track > track, TVector3 &trackPoint, TVector3 &trackDir)
Definition: MVAAlg.cxx:661
int IsInActiveVol(const TVector3 &pos)
Definition: MVAAlg.cxx:76
bool fCheatVertex
Definition: MVAAlg.h:161
Float_t Z
Definition: plot.C:39
void hits()
Definition: readHits.C:15
std::map< art::Ptr< recob::Track >, std::vector< art::Ptr< recob::Hit > > > fTracksToHits
Definition: MVAAlg.h:147
Point_t const & Vertex() const
Access to track position at different points.
Definition: Track.h:127
const Double32_t * XYZ() const
Definition: SpacePoint.h:65
TLorentzVector fVertex4Vect
Definition: MVAAlg.h:163
Point_t const & End() const
Access to track position at different points.
Definition: Track.h:128
Float_t track
Definition: plot.C:34
void mvapid::MVAAlg::GetDetectorEdges ( )

Definition at line 89 of file MVAAlg.cxx.

References fDetMaxX, fDetMaxY, fDetMaxZ, fDetMinX, fDetMinY, fDetMinZ, geo::BoxBoundedGeo::MaxX(), geo::BoxBoundedGeo::MaxY(), geo::BoxBoundedGeo::MaxZ(), geo::BoxBoundedGeo::MinX(), geo::BoxBoundedGeo::MinY(), geo::BoxBoundedGeo::MinZ(), geo::GeometryCore::TotalNTPC(), and geo::GeometryCore::TPC().

Referenced by mvapid::MVAPID::beginJob(), and mvapid::MVAAlg::SumDistance2::operator()().

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(unsigned int t=0; t<geom->TotalNTPC(); t++)
102  {
103  if(geom->TPC(t).MinX() < fDetMinX)
104  fDetMinX = geom->TPC(t).MinX();
105  if(geom->TPC(t).MaxX() > fDetMaxX)
106  fDetMaxX = geom->TPC(t).MaxX();
107  if(geom->TPC(t).MinY() < fDetMinY)
108  fDetMinY = geom->TPC(t).MinY();
109  if(geom->TPC(t).MaxY() > fDetMaxY)
110  fDetMaxY = geom->TPC(t).MaxY();
111  if(geom->TPC(t).MinZ() < fDetMinZ)
112  fDetMinZ = geom->TPC(t).MinZ();
113  if(geom->TPC(t).MaxZ() > fDetMaxZ)
114  fDetMaxZ = geom->TPC(t).MaxZ();
115  }
116 }
unsigned int TotalNTPC() const
Returns the total number of TPCs in the detector.
double fDetMaxX
Definition: MVAAlg.h:129
double MinX() const
Returns the world x coordinate of the start of the box.
Definition: BoxBoundedGeo.h:90
double MaxX() const
Returns the world x coordinate of the end of the box.
Definition: BoxBoundedGeo.h:93
double fDetMinY
Definition: MVAAlg.h:129
double fDetMinZ
Definition: MVAAlg.h:129
double MinZ() const
Returns the world z coordinate of the start of the box.
double fDetMinX
Definition: MVAAlg.h:129
double MaxY() const
Returns the world y coordinate of the end of the box.
double MaxZ() const
Returns the world z coordinate of the end of the box.
double fDetMaxZ
Definition: MVAAlg.h:129
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
double fDetMaxY
Definition: MVAAlg.h:129
double MinY() const
Returns the world y coordinate of the start of the box.
void mvapid::MVAAlg::GetWireNormals ( )

Definition at line 118 of file MVAAlg.cxx.

References fNormToWiresY, fNormToWiresZ, geo::GeometryCore::IteratePlanes(), geo::GeometryCore::Nplanes(), and geo::GeometryCore::NTPC().

Referenced by mvapid::MVAPID::beginJob(), and mvapid::MVAAlg::SumDistance2::operator()().

119 {
120 
122 
123  fNormToWiresY.clear();
124  fNormToWiresZ.clear();
125 
126  int planeKey;
127 
128  //Get normals to wires for each plane in the detector
129  //This assumes equal numbers of TPCs in each cryostat and equal numbers of planes in each TPC
130  for (geo::PlaneGeo const& plane: geom->IteratePlanes()) {
131  std::string id = std::string(plane.ID());
132  int pcryo = id.find("C");
133  int ptpc = id.find("T");
134  int pplane = id.find("P");
135  std::string scryo = id.substr(pcryo+2,2);
136  std::string stpc = id.substr(ptpc+2,2);
137  std::string splane = id.substr(pplane+2,2);
138  int icryo = std::stoi(scryo);
139  int itpc = std::stoi(stpc);
140  int iplane = std::stoi(splane);
141  planeKey = icryo * geom->NTPC(0) * geom->Nplanes(0, 0) + itpc * geom->Nplanes(0, 0) + iplane; //single index for all planes in detector
142  fNormToWiresY.insert(std::make_pair(planeKey, -plane.Wire(0).Direction().Z())); //y component of normal
143  fNormToWiresZ.insert(std::make_pair(planeKey, plane.Wire(0).Direction().Y())); //z component of normal
144 
145  }
146 }
IteratorBox< plane_iterator,&GeometryCore::begin_plane,&GeometryCore::end_plane > IteratePlanes() const
Enables ranged-for loops on all planes of the detector.
unsigned int Nplanes(unsigned int tpc=0, unsigned int cstat=0) const
Returns the total number of wire planes in the specified TPC.
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:78
unsigned int NTPC(unsigned int cstat=0) const
Returns the total number of TPCs in the specified cryostat.
std::map< int, double > fNormToWiresY
Definition: MVAAlg.h:131
std::map< int, double > fNormToWiresZ
Definition: MVAAlg.h:132
int mvapid::MVAAlg::IsInActiveVol ( const TVector3 &  pos)
private

Definition at line 76 of file MVAAlg.cxx.

References fDetMaxX, fDetMaxY, fDetMaxZ, fDetMinX, fDetMinY, and fDetMinZ.

Referenced by FitAndSortTrack(), mvapid::MVAAlg::SumDistance2::operator()(), and SortShower().

77 {
78  const double fiducialDist = 5.0;
79 
80  if(pos.X() > (fDetMinX + fiducialDist) && pos.X() < (fDetMaxX - fiducialDist)
81  && pos.Y() > (fDetMinY + fiducialDist) && pos.Y() < (fDetMaxY - fiducialDist)
82  && pos.Z() > (fDetMinZ + fiducialDist) && pos.Z() < (fDetMaxZ - fiducialDist))
83  return 1;
84  else
85  return 0;
86 
87 }
double fDetMaxX
Definition: MVAAlg.h:129
double fDetMinY
Definition: MVAAlg.h:129
double fDetMinZ
Definition: MVAAlg.h:129
double fDetMinX
Definition: MVAAlg.h:129
double fDetMaxZ
Definition: MVAAlg.h:129
double fDetMaxY
Definition: MVAAlg.h:129
int mvapid::MVAAlg::LinFit ( const art::Ptr< recob::Track track,
TVector3 &  trackPoint,
TVector3 &  trackDir 
)
private

Definition at line 661 of file MVAAlg.cxx.

References recob::Track::End(), fTracksToSpacePoints, and recob::Track::Vertex().

Referenced by FitAndSortTrack(), and mvapid::MVAAlg::SumDistance2::operator()().

661  {
662 
663  const std::vector<art::Ptr<recob::SpacePoint> >& sp = fTracksToSpacePoints.at(track);
664 
665  TGraph2D grFit(1);
666  unsigned int iPt=0;
667  for(auto spIter=sp.begin();spIter!=sp.end();++spIter){
668  TVector3 point=(*spIter)->XYZ();
669  grFit.SetPoint(iPt++,point.X(),point.Y(),point.Z());
670  }
671 
672  //Lift from the ROOT line3Dfit.C tutorial
673  ROOT::Fit::Fitter fitter;
674  // make the functor object
675  mvapid::MVAAlg::SumDistance2 sdist(&grFit);
676 
677  ROOT::Math::Functor fcn(sdist,6);
678 
679  //Initial fit parameters from track start and end...
680  TVector3 trackStart=track->Vertex<TVector3>();
681  TVector3 trackEnd=track->End<TVector3>();
682  trackDir=(trackEnd-trackStart).Unit();
683 
684  TVector3 x0=trackStart-trackDir;
685  TVector3 u=trackDir;
686 
687  double pStart[6] = {x0.X(),u.X(),x0.Y(),u.Y(),x0.Z(),u.Z()};
688 
689  fitter.SetFCN(fcn,pStart);
690 
691  bool ok = fitter.FitFCN();
692  if (!ok) {
693  trackPoint.SetXYZ(x0.X(),x0.Y(),x0.Z());
694  trackDir.SetXYZ(u.X(),u.Y(),u.Z());
695  trackDir=trackDir.Unit();
696  return 1;
697  }
698  else{
699  const ROOT::Fit::FitResult & result = fitter.Result();
700  const double * parFit = result.GetParams();
701  trackPoint.SetXYZ(parFit[0],parFit[2],parFit[4]);
702  trackDir.SetXYZ(parFit[1],parFit[3],parFit[5]);
703  trackDir=trackDir.Unit();
704  return 0;
705  }
706 }
std::map< art::Ptr< recob::Track >, std::vector< art::Ptr< recob::SpacePoint > > > fTracksToSpacePoints
Definition: MVAAlg.h:148
Point_t const & Vertex() const
Access to track position at different points.
Definition: Track.h:127
Point_t const & End() const
Access to track position at different points.
Definition: Track.h:128
int mvapid::MVAAlg::LinFitShower ( const art::Ptr< recob::Shower shower,
TVector3 &  showerPoint,
TVector3 &  showerDir 
)
private

Definition at line 708 of file MVAAlg.cxx.

References recob::Shower::Direction(), fShowersToSpacePoints, and recob::Shower::ShowerStart().

Referenced by mvapid::MVAAlg::SumDistance2::operator()(), and SortShower().

708  {
709 
710  const std::vector<art::Ptr<recob::SpacePoint> >& sp = fShowersToSpacePoints.at(shower);
711 
712  TGraph2D grFit(1);
713  unsigned int iPt=0;
714  for(auto spIter=sp.begin();spIter!=sp.end();++spIter){
715  TVector3 point=(*spIter)->XYZ();
716  grFit.SetPoint(iPt++,point.X(),point.Y(),point.Z());
717  }
718 
719  //Lift from the ROOT line3Dfit.C tutorial
720  ROOT::Fit::Fitter fitter;
721  // make the functor object
722  mvapid::MVAAlg::SumDistance2 sdist(&grFit);
723 
724  ROOT::Math::Functor fcn(sdist,6);
725 
726  //Initial fit parameters from shower start and end...
727  TVector3 showerStart=shower->ShowerStart();
728  showerDir = shower->Direction().Unit();
729 
730  TVector3 x0=showerStart-showerDir;
731  TVector3 u=showerDir;
732 
733  double pStart[6] = {x0.X(),u.X(),x0.Y(),u.Y(),x0.Z(),u.Z()};
734 
735  fitter.SetFCN(fcn,pStart);
736 
737  bool ok = fitter.FitFCN();
738  if (!ok) {
739  showerPoint.SetXYZ(x0.X(),x0.Y(),x0.Z());
740  showerDir.SetXYZ(u.X(),u.Y(),u.Z());
741  showerDir=showerDir.Unit();
742  return 1;
743  }
744  else{
745  const ROOT::Fit::FitResult & result = fitter.Result();
746  const double * parFit = result.GetParams();
747  showerPoint.SetXYZ(parFit[0],parFit[2],parFit[4]);
748  showerDir.SetXYZ(parFit[1],parFit[3],parFit[5]);
749  showerDir=showerDir.Unit();
750  return 0;
751  }
752 }
const TVector3 & ShowerStart() const
Definition: Shower.h:192
const TVector3 & Direction() const
Definition: Shower.h:189
std::map< art::Ptr< recob::Shower >, std::vector< art::Ptr< recob::SpacePoint > > > fShowersToSpacePoints
Definition: MVAAlg.h:150
void mvapid::MVAAlg::PrepareEvent ( const art::Event event)
private

Definition at line 272 of file MVAAlg.cxx.

References fCheatVertex, fEventT0, fHitLabel, fHits, fHitsToSpacePoints, fShowerLabel, fShowers, fShowersToHits, fShowersToSpacePoints, fSpacePointLabel, fSpacePoints, fSpacePointsToHits, fTrackingLabel, fTrackLabel, fTracks, fTracksToHits, fTracksToSpacePoints, fVertex4Vect, art::DataViewImpl::getByLabel(), and track.

Referenced by mvapid::MVAAlg::SumDistance2::operator()(), and RunPID().

272  {
273 
274  fHits.clear();
275  fSpacePoints.clear();
276  fTracks.clear();
277  fShowers.clear();
278  fSpacePointsToHits.clear();
279  fHitsToSpacePoints.clear();
280  fTracksToHits.clear();
281  fTracksToSpacePoints.clear();
282  fShowersToHits.clear();
283  fShowersToSpacePoints.clear();
284 
285  auto const* detProp = lar::providerFrom<detinfo::DetectorPropertiesService>();
286  fEventT0=detProp->TriggerOffset();
287 
289  evt.getByLabel(fHitLabel, hitsHandle);
290 
291  for (unsigned int iHit = 0; iHit < hitsHandle->size(); ++iHit){
292  const art::Ptr<recob::Hit> hit(hitsHandle, iHit);
293  fHits.push_back(hit);
294  }
295 
297  evt.getByLabel(fTrackLabel, tracksHandle);
298 
299  for (unsigned int iTrack = 0; iTrack < tracksHandle->size(); ++iTrack){
300  const art::Ptr<recob::Track> track(tracksHandle, iTrack);
301  fTracks.push_back(track);
302  }
303 
305  evt.getByLabel(fShowerLabel, showersHandle);
306 
307  for (unsigned int iShower = 0; iShower < showersHandle->size(); ++iShower){
308  const art::Ptr<recob::Shower> shower(showersHandle, iShower);
309  fShowers.push_back(shower);
310  }
311 
313  evt.getByLabel(fSpacePointLabel, spHandle);
314 
315  for (unsigned int iSP = 0; iSP < spHandle->size(); ++iSP){
316  const art::Ptr<recob::SpacePoint> spacePoint(spHandle, iSP);
317  fSpacePoints.push_back(spacePoint);
318  }
319 
323 
324  for(unsigned int iSP = 0; iSP < fSpacePoints.size(); ++iSP)
325  {
326  const art::Ptr<recob::SpacePoint> spacePoint=fSpacePoints.at(iSP);
327 
328  const art::Ptr<recob::Hit> hit = findSPToHits.at(iSP);
329  fSpacePointsToHits[spacePoint]=hit;
330  fHitsToSpacePoints[hit]=spacePoint;
331  }
332 
333  for(unsigned int iTrack = 0; iTrack < fTracks.size(); ++iTrack){
334  const art::Ptr<recob::Track> track=fTracks.at(iTrack);
335 
336  const std::vector< art::Ptr<recob::Hit> > trackHits = findTracksToHits.at(iTrack);
337 
338  for (unsigned int iHit=0; iHit<trackHits.size(); ++iHit)
339  {
340  const art::Ptr<recob::Hit> hit = trackHits.at(iHit);
341  fTracksToHits[track].push_back(hit);
342  if(fHitsToSpacePoints.count(hit)){
343  fTracksToSpacePoints[track].push_back(fHitsToSpacePoints.at(hit));
344  }
345  }
346  }
347 
348  for(unsigned int iShower = 0; iShower < fShowers.size(); ++iShower){
349  const art::Ptr<recob::Shower> shower=fShowers.at(iShower);
350  const std::vector< art::Ptr<recob::Hit> > showerHits = findShowersToHits.at(iShower);
351 
352  for (unsigned int iHit=0; iHit<showerHits.size(); ++iHit)
353  {
354  const art::Ptr<recob::Hit> hit = showerHits.at(iHit);
355  fShowersToHits[shower].push_back(hit);
356  if(fHitsToSpacePoints.count(hit)){
357  fShowersToSpacePoints[shower].push_back(fHitsToSpacePoints.at(hit));
358  }
359  }
360  }
361 
362  if(fCheatVertex){
364  evt.getByLabel(fTrackingLabel, partHandle);
365 
366  if(partHandle->size()==0||partHandle->at(0).TrackId()!=1){
367  std::cout<<"Error, ID of first track in largeant list is not 0"<<std::endl;
368  exit(1);
369  }
370  fVertex4Vect=partHandle->at(0).Position();
371  }
372 }
std::vector< art::Ptr< recob::Shower > > fShowers
Definition: MVAAlg.h:143
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::SpacePoint > > fHitsToSpacePoints
Definition: MVAAlg.h:151
std::map< art::Ptr< recob::Track >, std::vector< art::Ptr< recob::SpacePoint > > > fTracksToSpacePoints
Definition: MVAAlg.h:148
std::string fSpacePointLabel
Definition: MVAAlg.h:139
bool fCheatVertex
Definition: MVAAlg.h:161
std::vector< art::Ptr< recob::Hit > > fHits
Definition: MVAAlg.h:145
std::map< art::Ptr< recob::Shower >, std::vector< art::Ptr< recob::Hit > > > fShowersToHits
Definition: MVAAlg.h:149
std::map< art::Ptr< recob::Track >, std::vector< art::Ptr< recob::Hit > > > fTracksToHits
Definition: MVAAlg.h:147
std::vector< art::Ptr< recob::SpacePoint > > fSpacePoints
Definition: MVAAlg.h:144
TLorentzVector fVertex4Vect
Definition: MVAAlg.h:163
Detector simulation of raw signals on wires.
std::string fTrackingLabel
Definition: MVAAlg.h:140
std::string fTrackLabel
Definition: MVAAlg.h:136
std::string fHitLabel
Definition: MVAAlg.h:138
std::map< art::Ptr< recob::Shower >, std::vector< art::Ptr< recob::SpacePoint > > > fShowersToSpacePoints
Definition: MVAAlg.h:150
std::vector< art::Ptr< recob::Track > > fTracks
Definition: MVAAlg.h:142
std::string fShowerLabel
Definition: MVAAlg.h:137
std::map< art::Ptr< recob::SpacePoint >, art::Ptr< recob::Hit > > fSpacePointsToHits
Definition: MVAAlg.h:152
TCEvent evt
Definition: DataStructs.cxx:5
Float_t track
Definition: plot.C:34
double fEventT0
Definition: MVAAlg.h:127
void mvapid::MVAAlg::reconfigure ( fhicl::ParameterSet const &  p)

Definition at line 74 of file MVAAlg.cxx.

Referenced by mvapid::MVAAlg::SumDistance2::operator()().

74 {}
void mvapid::MVAAlg::RunPCA ( std::vector< art::Ptr< recob::Hit > > &  hits,
std::vector< double > &  eVals,
std::vector< double > &  eVecs 
)
private

Definition at line 513 of file MVAAlg.cxx.

References fHitsToSpacePoints, and hits().

Referenced by mvapid::MVAAlg::SumDistance2::operator()(), and RunPID().

513  {
514 
515  // Define the TPrincipal
516  TPrincipal* principal = new TPrincipal(3,"D");
517  // Define variables to hold the eigenvalues and eigenvectors
518  //const TMatrixD* covar = new TMatrixD();
519  //const TVectorD* meanval = new TVectorD();
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 }
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::SpacePoint > > fHitsToSpacePoints
Definition: MVAAlg.h:151
void mvapid::MVAAlg::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 at line 148 of file MVAAlg.cxx.

References _Var_Shape(), CalcSegmentdEdxFrac(), anab::MVAPIDResult::concentration, anab::MVAPIDResult::conicalness, anab::MVAPIDResult::coreHaloRatio, util::CreateAssn(), anab::MVAPIDResult::dEdxEnd, anab::MVAPIDResult::dEdxEndRatio, anab::MVAPIDResult::dEdxStart, anab::MVAPIDResult::evalRatio, FitAndSortTrack(), fMVAMethods, fParentModule, fReader, fResHolder, fShowers, fShowersToHits, fTracks, fTracksToHits, mvapid::MVAAlg::SortedObj::hitMap, anab::MVAPIDResult::isStoppingReco, anab::MVAPIDResult::isTrack, anab::MVAPIDResult::length, mvapid::MVAAlg::SortedObj::length, anab::MVAPIDResult::mvaOutput, anab::MVAPIDResult::nSpacePoints, PrepareEvent(), RunPCA(), SortShower(), and anab::MVAPIDResult::trackID.

Referenced by mvapid::MVAAlg::SumDistance2::operator()(), and mvapid::MVAPID::produce().

150  {
151 
152  this->PrepareEvent(evt);
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(sortedObj,0.,0.05);
169  double dEdxEnd = CalcSegmentdEdxFrac(sortedObj,0.9,1.0);
170  double dEdxPenultimate = CalcSegmentdEdxFrac(sortedObj,0.8,0.9);
171 
172  /*
173  std::cout<<"coreHaloRatio: "<<coreHaloRatio<<std::endl;
174  std::cout<<"concentration: "<<concentration<<std::endl;
175  std::cout<<"conicalness: "<<conicalness<<std::endl;
176  std::cout<<"dEdxStart: "<<dEdxStart<<std::endl;
177  std::cout<<"dEdxEnd: "<<dEdxEnd<<std::endl;
178  std::cout<<"dEdxEndRatio: ";
179  if(dEdxPenultimate < 0.1)
180  std::cout<<"1.0";
181  else
182  std::cout<<dEdxEnd/dEdxPenultimate;
183  std::cout<<std::endl;
184  */
185 
187  fResHolder.isStoppingReco=isStoppingReco;
188  fResHolder.nSpacePoints=sortedObj.hitMap.size();
189  fResHolder.trackID=(*trackIter)->ID();
190  fResHolder.evalRatio=evalRatio;
191  fResHolder.concentration=concentration;
192  fResHolder.coreHaloRatio=coreHaloRatio;
193  fResHolder.conicalness=conicalness;
194  fResHolder.dEdxStart=dEdxStart;
195  fResHolder.dEdxEnd=dEdxEnd;
196  if(dEdxPenultimate < 0.1)
198  else
199  fResHolder.dEdxEndRatio=dEdxEnd/dEdxPenultimate;
200  fResHolder.length=sortedObj.length;
201 
202  for(auto methodIter=fMVAMethods.begin();methodIter!=fMVAMethods.end();++methodIter){
203  fResHolder.mvaOutput[*methodIter]=fReader.EvaluateMVA(*methodIter);
204  }
205  result.push_back(fResHolder);
206  util::CreateAssn(*fParentModule, evt, result, *trackIter, trackAssns);
207  }
208 
209  for(auto showerIter=fShowers.begin();showerIter!=fShowers.end();++showerIter){
210  mvapid::MVAAlg::SortedObj sortedObj;
211 
212  std::vector<double> eVals,eVecs;
213  int isStoppingReco;
214 
215  this->RunPCA(fShowersToHits[*showerIter],eVals,eVecs);
216 
217  double evalRatio;
218  if(eVals[0]< 0.0001)
219  evalRatio=0.0;
220  else
221  evalRatio=std::sqrt(eVals[1]*eVals[1]+eVals[2]*eVals[2])/eVals[0];
222 
223  //this->SortShower(*showerIter,TVector3(eVecs[0],eVecs[3],eVecs[6]),isStoppingReco,sortedObj);
224  this->SortShower(*showerIter,isStoppingReco,sortedObj);
225 
226  double coreHaloRatio,concentration,conicalness;
227  this->_Var_Shape(sortedObj,coreHaloRatio,concentration,conicalness);
228  double dEdxStart = CalcSegmentdEdxFrac(sortedObj,0.,0.05);
229  double dEdxEnd = CalcSegmentdEdxFrac(sortedObj,0.9,1.0);
230  double dEdxPenultimate = CalcSegmentdEdxFrac(sortedObj,0.8,0.9);
231 
232  /*
233  std::cout<<"coreHaloRatio: "<<coreHaloRatio<<std::endl;
234  std::cout<<"concentration: "<<concentration<<std::endl;
235  std::cout<<"conicalness: "<<conicalness<<std::endl;
236  std::cout<<"dEdxStart: "<<dEdxStart<<std::endl;
237  std::cout<<"dEdxEnd: "<<dEdxEnd<<std::endl;
238  std::cout<<"dEdxEndRatio: ";
239  if(dEdxPenultimate < 0.1)
240  std::cout<<"1.0";
241  else
242  std::cout<<dEdxEnd/dEdxPenultimate;
243  std::cout<<std::endl;
244  */
245 
247  fResHolder.isStoppingReco=isStoppingReco;
248  fResHolder.nSpacePoints=sortedObj.hitMap.size();
249  fResHolder.trackID=(*showerIter)->ID()+1000; //For the moment label showers by adding 1000 to ID
250 
251  fResHolder.evalRatio=evalRatio;
252  fResHolder.concentration=concentration;
253  fResHolder.coreHaloRatio=coreHaloRatio;
254  fResHolder.conicalness=conicalness;
255  fResHolder.dEdxStart=dEdxStart;
256  fResHolder.dEdxEnd=dEdxEnd;
257  if(dEdxPenultimate < 0.1)
259  else
260  fResHolder.dEdxEndRatio=dEdxEnd/dEdxPenultimate;
261  fResHolder.length=sortedObj.length;
262 
263  for(auto methodIter=fMVAMethods.begin();methodIter!=fMVAMethods.end();++methodIter){
264  fResHolder.mvaOutput[*methodIter]=fReader.EvaluateMVA(*methodIter);
265  }
266  result.push_back(fResHolder);
267  util::CreateAssn(*fParentModule, evt, result, *showerIter, showerAssns);
268  }
269 
270 }
std::vector< art::Ptr< recob::Shower > > fShowers
Definition: MVAAlg.h:143
anab::MVAPIDResult fResHolder
Definition: MVAAlg.h:154
std::vector< std::string > fMVAMethods
Definition: MVAAlg.h:158
double CalcSegmentdEdxFrac(const SortedObj &track, double start, double end)
Definition: MVAAlg.cxx:606
TMVA::Reader fReader
Definition: MVAAlg.h:156
void RunPCA(std::vector< art::Ptr< recob::Hit > > &hits, std::vector< double > &eVals, std::vector< double > &eVecs)
Definition: MVAAlg.cxx:513
void SortShower(art::Ptr< recob::Shower > shower, int &isStoppingReco, mvapid::MVAAlg::SortedObj &sortedShower)
Definition: MVAAlg.cxx:437
std::map< art::Ptr< recob::Shower >, std::vector< art::Ptr< recob::Hit > > > fShowersToHits
Definition: MVAAlg.h:149
void PrepareEvent(const art::Event &event)
Definition: MVAAlg.cxx:272
std::map< art::Ptr< recob::Track >, std::vector< art::Ptr< recob::Hit > > > fTracksToHits
Definition: MVAAlg.h:147
const art::EDProducer * fParentModule
Definition: MVAAlg.h:134
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
std::map< double, const art::Ptr< recob::Hit > > hitMap
Definition: MVAAlg.h:53
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:142
std::map< std::string, double > mvaOutput
Definition: MVAPIDResult.h:27
void FitAndSortTrack(art::Ptr< recob::Track > track, int &isStoppingReco, SortedObj &sortedObj)
Definition: MVAAlg.cxx:374
unsigned int trackID
Definition: MVAPIDResult.h:21
void mvapid::MVAAlg::SortShower ( art::Ptr< recob::Shower shower,
int &  isStoppingReco,
mvapid::MVAAlg::SortedObj sortedShower 
)
private

Definition at line 437 of file MVAAlg.cxx.

References mvapid::MVAAlg::SortedObj::dir, mvapid::MVAAlg::SortedObj::end, fCheatVertex, fHitsToSpacePoints, fShowersToHits, fVertex4Vect, mvapid::MVAAlg::SortedObj::hitMap, hits(), IsInActiveVol(), mvapid::MVAAlg::SortedObj::length, LinFitShower(), recob::Shower::ShowerStart(), mvapid::MVAAlg::SortedObj::start, and recob::SpacePoint::XYZ().

Referenced by mvapid::MVAAlg::SumDistance2::operator()(), and RunPID().

438  {
439  sortedShower.hitMap.clear();
440 
441  std::vector<art::Ptr<recob::Hit>> hits=fShowersToHits[shower];
442 
443  TVector3 showerEnd(0, 0, 0);
444  double furthestHitFromStart = -999.9;
445  for(auto hitIter=hits.begin();hitIter!=hits.end();++hitIter){
446 
447  if(!fHitsToSpacePoints.count(*hitIter)) continue;
449  if((TVector3(sp->XYZ()) - shower->ShowerStart()).Mag() > furthestHitFromStart)
450  {
451  showerEnd = TVector3(sp->XYZ());
452  furthestHitFromStart = (TVector3(sp->XYZ()) - shower->ShowerStart()).Mag();
453  }
454  }
455 
456  TVector3 showerPoint,showerDir;
457  this->LinFitShower(shower,showerPoint,showerDir);
458 
459  TVector3 nearestPointStart, nearestPointEnd;
460 
461  //Ensure that shower is fitted in correct direction (assuming for now that particle moves in +z direction)
462 
463  if(fCheatVertex){
464  if((showerEnd-fVertex4Vect.Vect()).Mag()>(shower->ShowerStart()-fVertex4Vect.Vect()).Mag()){
465  nearestPointStart = showerPoint+showerDir*(showerDir.Dot(shower->ShowerStart()-showerPoint)/showerDir.Mag2());
466  nearestPointEnd = showerPoint+showerDir*(showerDir.Dot(showerEnd-showerPoint)/showerDir.Mag2());
467  isStoppingReco = this->IsInActiveVol(showerEnd);
468  }
469  else
470  {
471  nearestPointStart = showerPoint+showerDir*(showerDir.Dot(showerEnd-showerPoint)/showerDir.Mag2());
472  nearestPointEnd = showerPoint+showerDir*(showerDir.Dot(shower->ShowerStart()-showerPoint)/showerDir.Mag2());
473  isStoppingReco = this->IsInActiveVol(shower->ShowerStart());
474  showerDir*=-1.;
475  }
476  }
477  else{
478  if(showerEnd.Z() >= shower->ShowerStart().Z()){
479  nearestPointStart = showerPoint+showerDir*(showerDir.Dot(shower->ShowerStart()-showerPoint)/showerDir.Mag2());
480  nearestPointEnd = showerPoint+showerDir*(showerDir.Dot(showerEnd-showerPoint)/showerDir.Mag2());
481  isStoppingReco = this->IsInActiveVol(showerEnd);
482  }
483  else{
484  nearestPointStart = showerPoint+showerDir*(showerDir.Dot(showerEnd-showerPoint)/showerDir.Mag2());
485  nearestPointEnd = showerPoint+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=dir;
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 = showerPoint+showerDir*(showerDir.Dot(TVector3(sp->XYZ())-showerPoint)/showerDir.Mag2());
508  double lengthAlongShower=(nearestPointStart-nearestPoint).Mag();
509  sortedShower.hitMap.insert(std::pair<double,art::Ptr<recob::Hit> >(lengthAlongShower,*hitIter));
510  }
511 
512 }
const TVector3 & ShowerStart() const
Definition: Shower.h:192
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::SpacePoint > > fHitsToSpacePoints
Definition: MVAAlg.h:151
int IsInActiveVol(const TVector3 &pos)
Definition: MVAAlg.cxx:76
bool fCheatVertex
Definition: MVAAlg.h:161
int LinFitShower(const art::Ptr< recob::Shower > shower, TVector3 &showerPoint, TVector3 &showerDir)
Definition: MVAAlg.cxx:708
std::map< art::Ptr< recob::Shower >, std::vector< art::Ptr< recob::Hit > > > fShowersToHits
Definition: MVAAlg.h:149
void hits()
Definition: readHits.C:15
const Double32_t * XYZ() const
Definition: SpacePoint.h:65
TLorentzVector fVertex4Vect
Definition: MVAAlg.h:163
std::map< double, const art::Ptr< recob::Hit > > hitMap
Definition: MVAAlg.h:53

Member Data Documentation

const calo::CalorimetryAlg mvapid::MVAAlg::fCaloAlg
private

Definition at line 125 of file MVAAlg.h.

Referenced by CalcSegmentdEdxDist().

bool mvapid::MVAAlg::fCheatVertex
private

Definition at line 161 of file MVAAlg.h.

Referenced by FitAndSortTrack(), MVAAlg(), PrepareEvent(), and SortShower().

double mvapid::MVAAlg::fDetMaxX
private

Definition at line 129 of file MVAAlg.h.

Referenced by GetDetectorEdges(), and IsInActiveVol().

double mvapid::MVAAlg::fDetMaxY
private

Definition at line 129 of file MVAAlg.h.

Referenced by GetDetectorEdges(), and IsInActiveVol().

double mvapid::MVAAlg::fDetMaxZ
private

Definition at line 129 of file MVAAlg.h.

Referenced by GetDetectorEdges(), and IsInActiveVol().

double mvapid::MVAAlg::fDetMinX
private

Definition at line 129 of file MVAAlg.h.

Referenced by GetDetectorEdges(), and IsInActiveVol().

double mvapid::MVAAlg::fDetMinY
private

Definition at line 129 of file MVAAlg.h.

Referenced by GetDetectorEdges(), and IsInActiveVol().

double mvapid::MVAAlg::fDetMinZ
private

Definition at line 129 of file MVAAlg.h.

Referenced by GetDetectorEdges(), and IsInActiveVol().

double mvapid::MVAAlg::fEventT0
private

Definition at line 127 of file MVAAlg.h.

Referenced by CalcSegmentdEdxDist(), and PrepareEvent().

std::string mvapid::MVAAlg::fHitLabel
private

Definition at line 138 of file MVAAlg.h.

Referenced by MVAAlg(), and PrepareEvent().

std::vector<art::Ptr<recob::Hit> > mvapid::MVAAlg::fHits
private

Definition at line 145 of file MVAAlg.h.

Referenced by PrepareEvent().

std::map<art::Ptr<recob::Hit>,art::Ptr<recob::SpacePoint> > mvapid::MVAAlg::fHitsToSpacePoints
private

Definition at line 151 of file MVAAlg.h.

Referenced by _Var_Shape(), FitAndSortTrack(), PrepareEvent(), RunPCA(), and SortShower().

std::vector<std::string> mvapid::MVAAlg::fMVAMethods
private

Definition at line 158 of file MVAAlg.h.

Referenced by MVAAlg(), and RunPID().

std::map<int,double> mvapid::MVAAlg::fNormToWiresY
private

Definition at line 131 of file MVAAlg.h.

Referenced by CalcSegmentdEdxDist(), and GetWireNormals().

std::map<int,double> mvapid::MVAAlg::fNormToWiresZ
private

Definition at line 132 of file MVAAlg.h.

Referenced by CalcSegmentdEdxDist(), and GetWireNormals().

const art::EDProducer* mvapid::MVAAlg::fParentModule
private

Definition at line 134 of file MVAAlg.h.

Referenced by RunPID().

TMVA::Reader mvapid::MVAAlg::fReader
private

Definition at line 156 of file MVAAlg.h.

Referenced by MVAAlg(), and RunPID().

anab::MVAPIDResult mvapid::MVAAlg::fResHolder
private

Definition at line 154 of file MVAAlg.h.

Referenced by MVAAlg(), and RunPID().

std::string mvapid::MVAAlg::fShowerLabel
private

Definition at line 137 of file MVAAlg.h.

Referenced by MVAAlg(), and PrepareEvent().

std::vector<art::Ptr<recob::Shower> > mvapid::MVAAlg::fShowers
private

Definition at line 143 of file MVAAlg.h.

Referenced by PrepareEvent(), and RunPID().

std::map<art::Ptr<recob::Shower>,std::vector<art::Ptr<recob::Hit> > > mvapid::MVAAlg::fShowersToHits
private

Definition at line 149 of file MVAAlg.h.

Referenced by PrepareEvent(), RunPID(), and SortShower().

std::map<art::Ptr<recob::Shower>,std::vector<art::Ptr<recob::SpacePoint> > > mvapid::MVAAlg::fShowersToSpacePoints
private

Definition at line 150 of file MVAAlg.h.

Referenced by LinFitShower(), and PrepareEvent().

std::string mvapid::MVAAlg::fSpacePointLabel
private

Definition at line 139 of file MVAAlg.h.

Referenced by MVAAlg(), and PrepareEvent().

std::vector<art::Ptr<recob::SpacePoint> > mvapid::MVAAlg::fSpacePoints
private

Definition at line 144 of file MVAAlg.h.

Referenced by PrepareEvent().

std::map<art::Ptr<recob::SpacePoint>,art::Ptr<recob::Hit > > mvapid::MVAAlg::fSpacePointsToHits
private

Definition at line 152 of file MVAAlg.h.

Referenced by PrepareEvent().

std::string mvapid::MVAAlg::fTrackingLabel
private

Definition at line 140 of file MVAAlg.h.

Referenced by MVAAlg(), and PrepareEvent().

std::string mvapid::MVAAlg::fTrackLabel
private

Definition at line 136 of file MVAAlg.h.

Referenced by MVAAlg(), and PrepareEvent().

std::vector<art::Ptr<recob::Track> > mvapid::MVAAlg::fTracks
private

Definition at line 142 of file MVAAlg.h.

Referenced by PrepareEvent(), and RunPID().

std::map<art::Ptr<recob::Track>,std::vector<art::Ptr<recob::Hit> > > mvapid::MVAAlg::fTracksToHits
private

Definition at line 147 of file MVAAlg.h.

Referenced by FitAndSortTrack(), PrepareEvent(), and RunPID().

std::map<art::Ptr<recob::Track>,std::vector<art::Ptr<recob::SpacePoint> > > mvapid::MVAAlg::fTracksToSpacePoints
private

Definition at line 148 of file MVAAlg.h.

Referenced by LinFit(), and PrepareEvent().

TLorentzVector mvapid::MVAAlg::fVertex4Vect
private

Definition at line 163 of file MVAAlg.h.

Referenced by FitAndSortTrack(), PrepareEvent(), and SortShower().

std::vector<std::string> mvapid::MVAAlg::fWeightFiles
private

Definition at line 159 of file MVAAlg.h.

Referenced by MVAAlg().


The documentation for this class was generated from the following files: